StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
digcluster.cxx
1 // //
3 // DIGCluster //
4 // //
5 // Class containing cluster information //
6 // (pixel list, digital charge) //
7 // //
8 // //
9 // //
10 // //
11 // //
13 #include <digcluster.h>
14 
15 
16 #include <TROOT.h> // for gROOT object
17 #include <TMath.h>
18 #include <TMatrixD.h>
19 #include <TCanvas.h>
20 #include <TGraph.h>
21 #include <TAxis.h>
22 #include <TRandom3.h>
23 #include <TFile.h>
24 #include <TTree.h>
25 #include <TBranch.h>
26 #include <TClonesArray.h>
27 
28 //include other classes.h:
29 #include "digplane.h"
30 #include "digadc.h"
31 
32 
33 using namespace std;
34 
35 //==============================================================================
36 ClassImp(DIGCluster)
37 //______________________________________________________________________________
38 //
39 DIGCluster::DIGCluster()
40 {
41  //
42  // default constructor
43  //
44  fNpixels=0;
45  fPixelMap.clear();
46  fDigitalChargeMap.clear();
47  Xposition_CoG =-2.0;
48  Yposition_CoG =-2.0;
49  fSeedPixelIndex = -1;
50 }
51 //______________________________________________________________________________
52 //
53 DIGCluster::~DIGCluster() {
54  //
55  // virtual destructor
56  //
57 }
58 //______________________________________________________________________________
59 //
60 DIGCluster::DIGCluster(DIGCluster & adigCluster) : TObject()
61 {
62 
63  //recopy constructor
64  fNpixels = adigCluster.GetNpixels();
65 
66  fPixelMap.resize((adigCluster.GetPixelMap()).size());
67  fDigitalChargeMap.resize((adigCluster.GetPixelMap()).size());
68  for (Int_t i=0 ; i<fNpixels ; i++){
69  fPixelMap[i] = adigCluster.GetPixelMap()[i];
70  fDigitalChargeMap[i] = adigCluster.GetDigitalCharge()[i];
71  }
72  Xposition_CoG = adigCluster.GetXposition_CoG();
73  Yposition_CoG = adigCluster.GetYposition_CoG();
74  fSeedPixelIndex = adigCluster.GetSeedPixelIndex();
75 
76 }
77 //______________________________________________________________________________
78 //
79 void DIGCluster::Clear(const Option_t *)
80 {
81  // delete pointers. fDIGParticleArray->Clear("C");
82 
83 }
84 //______________________________________________________________________________
85 //
86 Int_t DIGCluster::GetTotalCharge(){
87  Int_t TotalCharge = 0;
88  for (Int_t i=0 ; i<fNpixels ; i++){
89  TotalCharge+=fDigitalChargeMap[i];
90  }
91  return TotalCharge;
92 }
93 //______________________________________________________________________________
94 //
95 Int_t DIGCluster::Get1stCrownCharge(DIGPlane *myDIGPlane){
96  Int_t f1stCrownCharge = 0;
97  std::vector<Int_t> pixmapvector;
98  pixmapvector = Get1stCrownPixelsIndex(myDIGPlane);
99  Int_t imax = pixmapvector.size();
100  for ( Int_t i=0 ; i<imax ; i++ ) {
101  for ( Int_t j=0 ; j<fNpixels; j++ ) {
102  if(pixmapvector[i]==fPixelMap[j]){
103  f1stCrownCharge+=fDigitalChargeMap[j];
104  }
105  // cout<<" DIGCluster::Get1stCrownCharge i imax pixmapvector "<<i<<" "<<imax<<" "<<pixmapvector[i]<<endl;
106  // cout<<" j fNpixels fPixelMap "<< j<<" "<<fNpixels<<" "<<" "<<fPixelMap[j]<<endl;
107  // cout<<"f1stCrownCharge = "<<f1stCrownCharge<<endl;
108  }
109  }
110  //cout<<" DIGCluster::Get1stCrownCharge "<<f1stCrownCharge<<endl;
111  return f1stCrownCharge;
112 }
113 //______________________________________________________________________________
114 //
115 Int_t DIGCluster::Get2ndCrownCharge(DIGPlane *myDIGPlane){
116  Int_t f2ndCrownCharge = 0;
117  std::vector<Int_t> pixmapvector;
118  pixmapvector = Get2ndCrownPixelsIndex(myDIGPlane);
119  Int_t imax = pixmapvector.size();
120  for ( Int_t i=0 ; i<imax; i++ ) {
121  for ( Int_t j=0 ; j<fNpixels; j++ ) {
122  if(pixmapvector[i]==fPixelMap[j]){
123  f2ndCrownCharge+=fDigitalChargeMap[j];
124  }
125  }
126  }
127  return f2ndCrownCharge;
128 
129 }
130 //______________________________________________________________________________
131 //
132 Int_t DIGCluster::Get4NeigboursCharge(DIGPlane *myDIGPlane){
133  Int_t f4NeigboursCharge = 0;
134  std::vector<Int_t> pixmapvector;
135  pixmapvector = Get4NeigboursPixelsIndex(myDIGPlane);
136  Int_t imax = pixmapvector.size();
137  for ( Int_t i=0 ; i<imax ; i++ ) {
138  for ( Int_t j=0 ; j<fNpixels; j++ ) {
139  if(pixmapvector[i]==fPixelMap[j]){
140  f4NeigboursCharge+=fDigitalChargeMap[j];
141  }
142  }
143  }
144  return f4NeigboursCharge;
145 }
146 //______________________________________________________________________________
147 //
148 Int_t DIGCluster::GetMultiplicity(Int_t Threshold){
149  Int_t TotalMultiplicity = 0;
150  for (Int_t i=0 ; i<fNpixels ; i++){
151  if(fDigitalChargeMap[i]>=Threshold){
152  TotalMultiplicity++;
153  }
154  }
155  return TotalMultiplicity;
156 
157 
158 }
159 //______________________________________________________________________________
160 //
161 void DIGCluster::Compute_CoG(DIGPlane *myDIGPlane){
162 
163  Double_t XCoG = 0.0;
164  Double_t YCoG = 0.0;
165  Double_t TotalCharge = 0.0;
166  for (Int_t i = 0; i < fNpixels ; i++){
167  Int_t PixelNumber = fPixelMap[i];
168  Double_t Xpix = (myDIGPlane->GetPitchY()) * (0.5+PixelNumber%(myDIGPlane->GetNpixelsX()));
169  Double_t Ypix = (myDIGPlane->GetPitchX()) * (0.5+PixelNumber/(myDIGPlane->GetNpixelsX()));
170  Int_t Charge = fDigitalChargeMap[i];
171  XCoG += Charge * Xpix;
172  YCoG += Charge * Ypix;
173  TotalCharge += Charge;
174  }
175  if(TotalCharge>0){
176  XCoG = XCoG / TotalCharge;
177  YCoG = YCoG / TotalCharge;
178  }else{
179  XCoG = -1.0;
180  YCoG = -1.0;
181  }
182  SetXposition_CoG(XCoG);
183  SetYposition_CoG(YCoG);
184 }
185 //______________________________________________________________________________
186 //
187 
188 void DIGCluster::AddPixel(Int_t DigitalCharge, Int_t PixelNumber){
189  fNpixels++;
190  fPixelMap.push_back(PixelNumber);
191  fDigitalChargeMap.push_back(DigitalCharge);
192 }
193 //______________________________________________________________________________
194 //
195 void DIGCluster::PrintInfo() {
196  std::cout<<"---------DIGCluster properties------------- "<<endl;
197  std::cout<<" Total charge, mul1, mul2, mul3, mul4 "<<endl;
198  std::cout<<GetTotalCharge()
199  <<" "<<GetMultiplicity(1)
200  <<" "<<GetMultiplicity(2)
201  <<" "<<GetMultiplicity(3)
202  <<" "<<GetMultiplicity(4)<<endl;
203  for (Int_t i=0 ; i<fNpixels ; i++){
204  std::cout<<i<<" "<<fPixelMap[i]<<" "<<" "<<fDigitalChargeMap[i]<<endl;
205  }
206  std::cout<<" CoG "<<GetXposition_CoG() <<" "<< GetYposition_CoG()<<endl;
207  std::cout<<" SEED index, totalC seedcharge "<< GetSeedPixelIndex()<<" "<<GetTotalCharge()<<" "<<GetDigitalCharge()[GetSeedPixelIndex()]<<endl;
208 
209 }
210 //______________________________________________________________________________
211 //
212 void DIGCluster::Compute_SeedPixel(DIGPlane *myDIGPlane){
213  // compute the seed pixel.
214  // Find the highest charge and if there are several pixels with the maximum charge, chose the pixel which is the closest of the other pixels.
215 
216  Int_t CurrentDigSeedCharge = 0;
217  Int_t CurrentDigSeedIndex = -1;
218  Float_t CurrentDist = 0;
219  for (Int_t i = 0; i < fNpixels ; i++){
220  if(fDigitalChargeMap[i]>CurrentDigSeedCharge){
221  CurrentDigSeedCharge = fDigitalChargeMap[i];
222  CurrentDigSeedIndex = i;
223  for (Int_t j = 0; j < fNpixels ; j++){
224  if(fDigitalChargeMap[j]==CurrentDigSeedCharge){
225  Int_t SeedCandidatePixelNumber = fPixelMap[i];
226  Double_t SeedCandidateXpix = (myDIGPlane->GetPitchY()) * (0.5+SeedCandidatePixelNumber%(myDIGPlane->GetNpixelsX()));
227  Double_t SeedCandidateYpix = (myDIGPlane->GetPitchX()) * (0.5+SeedCandidatePixelNumber/(myDIGPlane->GetNpixelsX()));
228  Int_t PixelNumber = fPixelMap[j];
229  Double_t Xpix = (myDIGPlane->GetPitchY()) * (0.5+PixelNumber%(myDIGPlane->GetNpixelsX()));
230  Double_t Ypix = (myDIGPlane->GetPitchX()) * (0.5+PixelNumber/(myDIGPlane->GetNpixelsX()));
231  CurrentDist += TMath::Sqrt((Xpix-SeedCandidateXpix)*(Xpix-SeedCandidateXpix)+(Ypix-SeedCandidateYpix)*(Ypix-SeedCandidateYpix));
232  }
233  }
234  }else if(fDigitalChargeMap[i]==CurrentDigSeedCharge){
235  Float_t newdist = 0;
236  for (Int_t j = 0; j < fNpixels ; j++){
237  if(fDigitalChargeMap[j]==CurrentDigSeedCharge){
238  Int_t SeedCandidatePixelNumber = fPixelMap[i];
239  Double_t SeedCandidateXpix = (myDIGPlane->GetPitchY()) * (0.5+SeedCandidatePixelNumber%(myDIGPlane->GetNpixelsX()));
240  Double_t SeedCandidateYpix = (myDIGPlane->GetPitchX()) * (0.5+SeedCandidatePixelNumber/(myDIGPlane->GetNpixelsX()));
241  Int_t PixelNumber = fPixelMap[j];
242  Double_t Xpix = (myDIGPlane->GetPitchY()) * (0.5+PixelNumber%(myDIGPlane->GetNpixelsX()));
243  Double_t Ypix = (myDIGPlane->GetPitchX()) * (0.5+PixelNumber/(myDIGPlane->GetNpixelsX()));
244  newdist += TMath::Sqrt((Xpix-SeedCandidateXpix)*(Xpix-SeedCandidateXpix)+(Ypix-SeedCandidateYpix)*(Ypix-SeedCandidateYpix));
245  }
246  }
247  if(newdist<CurrentDist){
248  CurrentDigSeedCharge = fDigitalChargeMap[i];
249  CurrentDigSeedIndex = i;
250  CurrentDist = newdist;
251  }
252  }
253  }
254  SetSeedPixelIndex(CurrentDigSeedIndex);
255 
256 }
257 //______________________________________________________________________________
258 //
259 void DIGCluster::GetXYPixelNumber(Int_t &Xpix, Int_t &Ypix, DIGPlane *myDIGPlane, Int_t PixelNumber){
260 
261  Xpix = PixelNumber%(myDIGPlane->GetNpixelsX());
262  Ypix = PixelNumber/(myDIGPlane->GetNpixelsX());
263 
264 }
265 //______________________________________________________________________________
266 //
267 Int_t DIGCluster::GetSeedPixelIndex(){
268 //return the seed pixel index in the pixel map.
269  return fSeedPixelIndex;
270 }
271 //______________________________________________________________________________
272 //
273 void DIGCluster::SetSeedPixelIndex(Int_t Index)
274 {
275  fSeedPixelIndex = Index;
276 }
277 //______________________________________________________________________________
278 //
279 std::vector<Int_t> DIGCluster::Get4NeigboursPixelsIndex(DIGPlane *myDIGPlane){
280  vector< Int_t > PixelMap;
281  PixelMap.clear();
282  Int_t Xpix;
283  Int_t Ypix;
284  GetXYPixelNumber(Xpix,Ypix,myDIGPlane,fPixelMap[fSeedPixelIndex]); // <--------------------------------
285  if(IsPixelInThePlane(Xpix+1,Ypix,myDIGPlane)){
286  PixelMap.push_back(fPixelMap[fSeedPixelIndex]+1);
287  }
288  if(IsPixelInThePlane(Xpix-1,Ypix,myDIGPlane)){
289  PixelMap.push_back(fPixelMap[fSeedPixelIndex]-1);
290  }
291  if(IsPixelInThePlane(Xpix,Ypix+1,myDIGPlane)){
292  PixelMap.push_back(fPixelMap[fSeedPixelIndex]+myDIGPlane->GetNpixelsX());
293  }
294  if(IsPixelInThePlane(Xpix,Ypix-1,myDIGPlane)){
295  PixelMap.push_back(fPixelMap[fSeedPixelIndex]-myDIGPlane->GetNpixelsX());
296  }
297  return PixelMap;
298 }
299 //______________________________________________________________________________
300 //
301 std::vector<Int_t> DIGCluster::Get1stCrownPixelsIndex(DIGPlane *myDIGPlane){
302  vector< Int_t > PixelMap;
303  PixelMap.clear();
304  Int_t Xpix;
305  Int_t Ypix;
306  GetXYPixelNumber(Xpix,Ypix,myDIGPlane,fPixelMap[fSeedPixelIndex]);
307  if(IsPixelInThePlane(Xpix-1,Ypix-1,myDIGPlane)){
308  PixelMap.push_back(fPixelMap[fSeedPixelIndex]-1-myDIGPlane->GetNpixelsX());
309  }
310  if(IsPixelInThePlane(Xpix,Ypix-1,myDIGPlane)){
311  PixelMap.push_back(fPixelMap[fSeedPixelIndex]-myDIGPlane->GetNpixelsX());
312  }
313  if(IsPixelInThePlane(Xpix+1,Ypix-1,myDIGPlane)){
314  PixelMap.push_back(fPixelMap[fSeedPixelIndex]+1-myDIGPlane->GetNpixelsX());
315  }
316  if(IsPixelInThePlane(Xpix-1,Ypix,myDIGPlane)){
317  PixelMap.push_back(fPixelMap[fSeedPixelIndex]-1);
318  }
319  if(IsPixelInThePlane(Xpix+1,Ypix,myDIGPlane)){
320  PixelMap.push_back(fPixelMap[fSeedPixelIndex]+1);
321  }
322  if(IsPixelInThePlane(Xpix-1,Ypix+1,myDIGPlane)){
323  PixelMap.push_back(fPixelMap[fSeedPixelIndex]-1+myDIGPlane->GetNpixelsX());
324  }
325  if(IsPixelInThePlane(Xpix,Ypix+1,myDIGPlane)){
326  PixelMap.push_back(fPixelMap[fSeedPixelIndex]+myDIGPlane->GetNpixelsX());
327  }
328  if(IsPixelInThePlane(Xpix+1,Ypix+1,myDIGPlane)){
329  PixelMap.push_back(fPixelMap[fSeedPixelIndex]+1+myDIGPlane->GetNpixelsX());
330  }
331  return PixelMap;
332 }
333 //______________________________________________________________________________
334 //
335 std::vector<Int_t> DIGCluster::Get2ndCrownPixelsIndex(DIGPlane *myDIGPlane){
336  vector< Int_t > PixelMap;
337  PixelMap.clear();
338  Int_t Xpix;
339  Int_t Ypix;
340  GetXYPixelNumber(Xpix,Ypix,myDIGPlane,fPixelMap[fSeedPixelIndex]);
341  if(IsPixelInThePlane(Xpix-2,Ypix-2,myDIGPlane)){
342  PixelMap.push_back(fPixelMap[fSeedPixelIndex]-2-2*myDIGPlane->GetNpixelsX());
343  }
344  if(IsPixelInThePlane(Xpix-1,Ypix-2,myDIGPlane)){
345  PixelMap.push_back(fPixelMap[fSeedPixelIndex]-1-2*myDIGPlane->GetNpixelsX());
346  }
347  if(IsPixelInThePlane(Xpix,Ypix-2,myDIGPlane)){
348  PixelMap.push_back(fPixelMap[fSeedPixelIndex]-2*myDIGPlane->GetNpixelsX());
349  }
350  if(IsPixelInThePlane(Xpix+1,Ypix-2,myDIGPlane)){
351  PixelMap.push_back(fPixelMap[fSeedPixelIndex]+1-2*myDIGPlane->GetNpixelsX());
352  }
353  if(IsPixelInThePlane(Xpix+2,Ypix-2,myDIGPlane)){
354  PixelMap.push_back(fPixelMap[fSeedPixelIndex]+2-2*myDIGPlane->GetNpixelsX());
355  }
356 
357 
358  if(IsPixelInThePlane(Xpix-2,Ypix-1,myDIGPlane)){
359  PixelMap.push_back(fPixelMap[fSeedPixelIndex]-2-myDIGPlane->GetNpixelsX());
360  }
361  if(IsPixelInThePlane(Xpix+2,Ypix-1,myDIGPlane)){
362  PixelMap.push_back(fPixelMap[fSeedPixelIndex]+2-myDIGPlane->GetNpixelsX());
363  }
364 
365  if(IsPixelInThePlane(Xpix-2,Ypix,myDIGPlane)){
366  PixelMap.push_back(fPixelMap[fSeedPixelIndex]-2);
367  }
368  if(IsPixelInThePlane(Xpix+2,Ypix,myDIGPlane)){
369  PixelMap.push_back(fPixelMap[fSeedPixelIndex]+2);
370  }
371 
372 
373  if(IsPixelInThePlane(Xpix-2,Ypix+1,myDIGPlane)){
374  PixelMap.push_back(fPixelMap[fSeedPixelIndex]-2+myDIGPlane->GetNpixelsX());
375  }
376  if(IsPixelInThePlane(Xpix+2,Ypix+1,myDIGPlane)){
377  PixelMap.push_back(fPixelMap[fSeedPixelIndex]+2+myDIGPlane->GetNpixelsX());
378  }
379 
380  if(IsPixelInThePlane(Xpix-2,Ypix+2,myDIGPlane)){
381  PixelMap.push_back(fPixelMap[fSeedPixelIndex]-2+2*myDIGPlane->GetNpixelsX());
382  }
383  if(IsPixelInThePlane(Xpix-1,Ypix+2,myDIGPlane)){
384  PixelMap.push_back(fPixelMap[fSeedPixelIndex]-1+2*myDIGPlane->GetNpixelsX());
385  }
386  if(IsPixelInThePlane(Xpix,Ypix+2,myDIGPlane)){
387  PixelMap.push_back(fPixelMap[fSeedPixelIndex]+2*myDIGPlane->GetNpixelsX());
388  }
389  if(IsPixelInThePlane(Xpix+1,Ypix+2,myDIGPlane)){
390  PixelMap.push_back(fPixelMap[fSeedPixelIndex]+1+2*myDIGPlane->GetNpixelsX());
391  }
392  if(IsPixelInThePlane(Xpix+2,Ypix+2,myDIGPlane)){
393  PixelMap.push_back(fPixelMap[fSeedPixelIndex]+2+2*myDIGPlane->GetNpixelsX());
394  }
395 
396 
397  return PixelMap;
398 }
399 //______________________________________________________________________________
400 //
401 Bool_t DIGCluster::IsPixelInThePlane(Int_t Xpix, Int_t Ypix, DIGPlane *myDIGPlane){
402 
403  Bool_t IsIn = true;
404  Int_t NpixelsX = myDIGPlane->GetNpixelsX();
405  Int_t NpixelsY = myDIGPlane->GetNpixelsY();
406  if(Xpix<0){IsIn = false;}
407  if(Ypix<0){IsIn = false;}
408  if(Xpix>=NpixelsX){IsIn = false;}
409  if(Ypix>=NpixelsY){IsIn = false;}
410  return IsIn;
411 }
412 
413 
414 
StarPicoDst::Charge
Charge
Define charge.
Definition: StPicoCommon.h:20
Generated by   doxygen 1.8.5