472c8a13cb4e0dbe88a4b1c2999a39042466b9d3
[u/mrichter/AliRoot.git] / MUON / AliMUONClusterFinder.cxx
1 #include "AliMUONClusterFinder.h"
2 #include "TTree.h"
3 #include "AliRun.h"
4 #include <TCanvas.h>
5 #include <TH1.h>
6 #include <TPad.h>
7 #include <TGraph.h> 
8 #include <TPostScript.h> 
9 #include <TMinuit.h> 
10
11 ClassImp(AliMUONRecCluster)
12 //_____________________________________________________________________
13 static AliMUONsegmentation* gSegmentation;
14 static AliMUONresponse*     gResponse;
15 static Int_t                gix[500];
16 static Int_t                giy[500];
17 static Float_t              gCharge[500];
18 static Int_t                gNbins;
19 static Int_t                gFirst=kTRUE;
20 static TMinuit *gMyMinuit ;
21 void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
22 static Int_t                gChargeTot;
23
24
25 AliMUONRecCluster::AliMUONRecCluster()
26
27 {
28     fDigits=0;
29     fNdigit=-1;
30 }
31
32 AliMUONRecCluster::
33 AliMUONRecCluster(Int_t FirstDigit,Int_t Ichamber, Int_t Icathod)
34 {
35     fX = 0.;
36     fY = 0.;
37     fDigits = new TArrayI(10);
38     fNdigit=0;
39     AddDigit(FirstDigit);
40     fChamber=Ichamber;
41     fCathod=Icathod;
42 }
43
44 void AliMUONRecCluster::AddDigit(Int_t Digit)
45 {
46     if (fNdigit==fDigits->GetSize()) {
47         //enlarge the list by hand!
48         Int_t *array= new Int_t[fNdigit*2];
49         for(Int_t i=0;i<fNdigit;i++)
50             array[i] = fDigits->At(i);
51         fDigits->Adopt(fNdigit*2,array);
52     }
53     fDigits->AddAt(Digit,fNdigit);
54     fNdigit++;
55 }
56
57
58 AliMUONRecCluster::~AliMUONRecCluster()
59 {
60     if (fDigits)
61         delete fDigits;
62 }
63
64 Int_t AliMUONRecCluster::FirstDigitIndex()
65 {
66     fCurrentDigit=0;
67     return fDigits->At(fCurrentDigit);
68 }
69
70 Int_t AliMUONRecCluster::NextDigitIndex()
71 {
72     fCurrentDigit++;
73     if (fCurrentDigit<fNdigit)
74         return fDigits->At(fCurrentDigit);
75     else 
76         return InvalidDigitIndex();
77 }
78
79 Int_t AliMUONRecCluster::NDigits()
80 {
81     return fNdigit;
82 }
83 void AliMUONRecCluster::Finish()
84 {
85     // In order to reconstruct coordinates, one has to
86     // get back to the digits which is not trivial here,
87     // because we don't know where digits are stored!
88     // Center Of Gravity, or other method should be
89     // a property of AliMUON class!
90 }
91
92
93 //----------------------------------------------------------
94 ClassImp(AliMUONClusterFinder)
95
96     AliMUONClusterFinder::AliMUONClusterFinder
97 (AliMUONsegmentation *segmentation, AliMUONresponse *response, 
98  TClonesArray *digits, Int_t chamber)   
99 {
100     fSegmentation=segmentation;
101     fResponse=response;
102     
103     fDigits=digits;
104     fNdigits = fDigits->GetEntriesFast();
105     fChamber=chamber;
106     fRawClusters=new TClonesArray("AliMUONRawCluster",10000);
107     fNRawClusters=0;
108     fCogCorr = 0;
109     SetNperMax();
110     SetClusterSize();
111     SetDeclusterFlag();
112     fNPeaks=-1;
113 }
114
115     AliMUONClusterFinder::AliMUONClusterFinder()
116 {
117     fSegmentation=0;
118     fResponse=0;
119     
120     fDigits=0;
121     fNdigits = 0;
122     fChamber=-1;
123     fRawClusters=new TClonesArray("AliMUONRawCluster",10000);
124     fNRawClusters=0;
125     fHitMap = 0;
126     fCogCorr = 0;
127     SetNperMax();
128     SetClusterSize();
129     SetDeclusterFlag();
130     fNPeaks=-1;
131 }
132
133 void AliMUONClusterFinder::AddRawCluster(const AliMUONRawCluster c)
134 {
135   //
136   // Add a raw cluster copy to the list
137   //
138     AliMUON *MUON=(AliMUON*)gAlice->GetModule("MUON");
139     MUON->AddRawCluster(fChamber,c); 
140     fNRawClusters++;
141 }
142
143
144
145 void AliMUONClusterFinder::Decluster(AliMUONRawCluster *cluster)
146 {
147 //    AliMUONdigit *dig;
148 //    Int_t q;
149
150     
151     Int_t mul = cluster->fMultiplicity;
152 //    printf("Decluster - multiplicity   %d \n",mul);
153
154     if (mul == 1 || mul ==2) {
155 //      printf("\n Nothing special for 1- and 2-clusters \n");
156 //
157 // Nothing special for 1- and 2-clusters
158         if (fNPeaks != 0) {
159             cluster->fNcluster[0]=fNPeaks;
160             cluster->fNcluster[1]=0;
161         } 
162         AddRawCluster(*cluster); 
163         fNPeaks++;
164     } else if (mul ==3) {
165 //
166 // 3-cluster, check topology
167 //      printf("\n 3-cluster, check topology \n");
168         if (fDeclusterFlag) {
169           if (Centered(cluster)) {
170             // ok, cluster is centered 
171             //      printf("\n ok, cluster is centered \n");
172           } else {
173             // cluster is not centered, split into 2+1
174             //      printf("\n cluster is not centered, split into 2+1 \n");
175           }
176         } else {
177           if (fNPeaks != 0) {
178             cluster->fNcluster[0]=fNPeaks;
179             cluster->fNcluster[1]=0;
180           } 
181           AddRawCluster(*cluster); 
182           fNPeaks++;
183         }           
184     } else {
185       //if (mul < 12) {
186           //      printf("Decluster - multiplicity > 45   %d \n",mul);
187           //printf("Decluster - multiplicity < 25   %d \n",mul);
188 // 
189 // 4-and more-pad clusters
190 //
191       if (mul <= fClusterSize) {
192         if (fDeclusterFlag) {
193             SplitByLocalMaxima(cluster);
194         } else {
195             if (fNPeaks != 0) {
196                 cluster->fNcluster[0]=fNPeaks;
197                 cluster->fNcluster[1]=0;
198             } 
199             AddRawCluster(*cluster);
200             fNPeaks++;
201         }       
202       }
203         //}
204     } // multiplicity 
205 }
206
207
208 Bool_t AliMUONClusterFinder::Centered(AliMUONRawCluster *cluster)
209 {
210     AliMUONdigit* dig;
211     dig= (AliMUONdigit*)fDigits->UncheckedAt(cluster->fIndexMap[0]);
212     Int_t ix=dig->fPadX;
213     Int_t iy=dig->fPadY;
214     Int_t nn;
215     Int_t X[kMaxNeighbours], Y[kMaxNeighbours], XN[kMaxNeighbours], YN[kMaxNeighbours];
216     
217     fSegmentation->Neighbours(ix,iy,&nn,X,Y);
218     Int_t nd=0;
219     for (Int_t i=0; i<nn; i++) {
220         if (fHitMap->TestHit(X[i],Y[i]) == used) {
221             XN[nd]=X[i];
222             YN[nd]=Y[i];
223             nd++;
224         }
225     }
226     if (nd==2) {
227 //
228 // cluster is centered !
229         if (fNPeaks != 0) {
230             cluster->fNcluster[0]=fNPeaks;
231             cluster->fNcluster[1]=0;
232         }  
233         AddRawCluster(*cluster);
234         fNPeaks++;
235         return kTRUE;
236     } else if (nd ==1) {
237 //
238 // Highest signal on an edge, split cluster into 2+1
239 //
240 // who is the neighbour ?
241         Int_t nind=fHitMap->GetHitIndex(XN[0], YN[0]);
242         Int_t i1= (nind==cluster->fIndexMap[1]) ? 1:2;
243         Int_t i2= (nind==cluster->fIndexMap[1]) ? 2:1;    
244 //
245 // 2-cluster
246         AliMUONRawCluster cnew;
247         if (fNPeaks == 0) {
248             cnew.fNcluster[0]=-1;
249             cnew.fNcluster[1]=fNRawClusters;
250         } else {
251             cnew.fNcluster[0]=fNPeaks;
252             cnew.fNcluster[1]=0;
253         }
254         cnew.fMultiplicity=2;
255         cnew.fIndexMap[0]=cluster->fIndexMap[0];
256         cnew.fIndexMap[1]=cluster->fIndexMap[i1];
257         FillCluster(&cnew);
258         cnew.fClusterType=cnew.PhysicsContribution();
259         AddRawCluster(cnew);
260         fNPeaks++;
261 //
262 // 1-cluster
263         cluster->fMultiplicity=1;
264         cluster->fIndexMap[0]=cluster->fIndexMap[i2];
265         cluster->fIndexMap[1]=0;
266         cluster->fIndexMap[2]=0;        
267         FillCluster(cluster);
268         if (fNPeaks != 0) {
269             cluster->fNcluster[0]=fNPeaks;
270             cluster->fNcluster[1]=0;
271         }  
272         cluster->fClusterType=cluster->PhysicsContribution();
273         AddRawCluster(*cluster);
274         fNPeaks++;
275         return kFALSE;
276     } else {
277         printf("\n Completely screwed up %d !! \n",nd);
278         
279     }
280     
281         return kFALSE;
282 }
283 void AliMUONClusterFinder::SplitByLocalMaxima(AliMUONRawCluster *c)
284 {
285     AliMUONdigit* dig[100], *digt;
286     Int_t ix[100], iy[100], q[100];
287     Float_t x[100], y[100];
288     Int_t i; // loops over digits
289     Int_t j; // loops over local maxima
290     //    Float_t xPeak[2];
291     //    Float_t yPeak[2];
292     //    Int_t threshold=500;
293     Int_t mul=c->fMultiplicity;
294 //
295 //  dump digit information into arrays
296 //
297     for (i=0; i<mul; i++)
298     {
299         dig[i]= (AliMUONdigit*)fDigits->UncheckedAt(c->fIndexMap[i]);
300         ix[i]= dig[i]->fPadX;
301         iy[i]= dig[i]->fPadY;
302         q[i] = dig[i]->fSignal;
303         fSegmentation->GetPadCxy(ix[i], iy[i], x[i], y[i]);
304     }
305 //
306 //  Find local maxima
307 //
308     Bool_t IsLocal[100];
309     Int_t NLocal=0;
310     Int_t AssocPeak[100];
311     Int_t IndLocal[100];
312     Int_t nn;
313     Int_t X[kMaxNeighbours], Y[kMaxNeighbours];
314     for (i=0; i<mul; i++) {
315         fSegmentation->Neighbours(ix[i], iy[i], &nn, X, Y);
316         IsLocal[i]=kTRUE;
317         for (j=0; j<nn; j++) {
318             if (fHitMap->TestHit(X[j], Y[j])==empty) continue;
319             digt=(AliMUONdigit*) fHitMap->GetHit(X[j], Y[j]);
320             if (digt->fSignal > q[i]) {
321                 IsLocal[i]=kFALSE;
322                 break;
323 //
324 // handle special case of neighbouring pads with equal signal
325             } else if (digt->fSignal == q[i]) {
326                 if (NLocal >0) {
327                     for (Int_t k=0; k<NLocal; k++) {
328                         if (X[j]==ix[IndLocal[k]] && Y[j]==iy[IndLocal[k]]){
329                             IsLocal[i]=kFALSE;
330                         }
331                     }
332                 }
333             } 
334         } // loop over next neighbours
335         // Maxima should not be on the edge
336         if (IsLocal[i]) {
337             IndLocal[NLocal]=i;
338             NLocal++;
339         } 
340     } // loop over all digits
341 //    printf("Found %d local Maxima",NLocal);
342 //
343 // If only one local maximum found but multiplicity is high 
344 // take global maximum from the list of digits.    
345     if (NLocal==1 && mul>12) {
346         Int_t nnew=0;
347         for (i=0; i<mul; i++) {
348             if (!IsLocal[i]) {
349                 IndLocal[NLocal]=i;
350                 IsLocal[i]=kTRUE;
351                 NLocal++;
352                 nnew++;
353             }
354             if (nnew==1) break;
355         }
356     }
357     
358 // If number of local maxima is 2 try to fit a double gaussian
359     if (NLocal==2) {
360 //
361 //  Initialise global variables for fit
362         gFirst=1;
363         gSegmentation=fSegmentation;
364         gResponse    =fResponse;
365         gNbins=mul;
366         
367         for (i=0; i<mul; i++) {
368             gix[i]=ix[i];
369             giy[i]=iy[i];
370             gCharge[i]=Float_t(q[i]);
371         }
372 //
373         if (gFirst) {
374           gFirst=kFALSE;
375           gMyMinuit = new TMinuit(5);
376         }
377         gMyMinuit->SetFCN(fcn);
378         gMyMinuit->mninit(5,10,7);
379         Double_t arglist[20];
380         Int_t ierflag=0;
381         arglist[0]=1;
382 //      gMyMinuit->mnexcm("SET ERR",arglist,1,ierflag);
383 // Set starting values 
384         static Double_t vstart[5];
385         vstart[0]=x[IndLocal[0]];
386         vstart[1]=y[IndLocal[0]];       
387         vstart[2]=x[IndLocal[1]];
388         vstart[3]=y[IndLocal[1]];       
389         vstart[4]=Float_t(q[IndLocal[0]])/
390             Float_t(q[IndLocal[0]]+q[IndLocal[1]]);
391 // lower and upper limits
392         static Double_t lower[5], upper[5];
393         Int_t isec=fSegmentation->Sector(ix[IndLocal[0]], iy[IndLocal[0]]);
394         lower[0]=vstart[0]-fSegmentation->Dpx(isec)/2;
395         lower[1]=vstart[1]-fSegmentation->Dpy(isec)/2;
396 //      lower[1]=vstart[1];
397         
398         upper[0]=lower[0]+fSegmentation->Dpx(isec);
399         upper[1]=lower[1]+fSegmentation->Dpy(isec);
400 //      upper[1]=vstart[1];
401         
402         isec=fSegmentation->Sector(ix[IndLocal[1]], iy[IndLocal[1]]);
403         lower[2]=vstart[2]-fSegmentation->Dpx(isec)/2;
404         lower[3]=vstart[3]-fSegmentation->Dpy(isec)/2;
405 //      lower[3]=vstart[3];
406         
407         upper[2]=lower[2]+fSegmentation->Dpx(isec);
408         upper[3]=lower[3]+fSegmentation->Dpy(isec);
409 //      upper[3]=vstart[3];
410         
411         lower[4]=0.;
412         upper[4]=1.;
413 // step sizes
414         static Double_t step[5]={0.005, 0.03, 0.005, 0.03, 0.01};
415         
416         gMyMinuit->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
417         gMyMinuit->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
418         gMyMinuit->mnparm(2,"x2",vstart[2],step[2],lower[2],upper[2],ierflag);
419         gMyMinuit->mnparm(3,"y2",vstart[3],step[3],lower[3],upper[3],ierflag);
420         gMyMinuit->mnparm(4,"a0",vstart[4],step[4],lower[4],upper[4],ierflag);
421 // ready for minimisation       
422         gMyMinuit->SetPrintLevel(-1);
423         gMyMinuit->mnexcm("SET OUT", arglist, 0, ierflag);
424         arglist[0]= -1;
425         arglist[1]= 0;
426         
427         gMyMinuit->mnexcm("SET NOGR", arglist, 0, ierflag);
428         gMyMinuit->mnexcm("SCAN", arglist, 0, ierflag);
429         gMyMinuit->mnexcm("EXIT" , arglist, 0, ierflag);
430 // Print results
431 //      Double_t amin,edm,errdef;
432 //      Int_t nvpar,nparx,icstat;
433 //      gMyMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
434 //      gMyMinuit->mnprin(3,amin);
435 // Get fitted parameters
436
437         Double_t xrec[2], yrec[2], qfrac;
438         TString chname;
439         Double_t epxz, b1, b2;
440         Int_t ierflg;
441         gMyMinuit->mnpout(0, chname, xrec[0], epxz, b1, b2, ierflg);    
442         gMyMinuit->mnpout(1, chname, yrec[0], epxz, b1, b2, ierflg);    
443         gMyMinuit->mnpout(2, chname, xrec[1], epxz, b1, b2, ierflg);    
444         gMyMinuit->mnpout(3, chname, yrec[1], epxz, b1, b2, ierflg);    
445         gMyMinuit->mnpout(4, chname, qfrac,   epxz, b1, b2, ierflg);    
446         printf("\n %f %f %f %f %f\n", xrec[0], yrec[0], xrec[1], yrec[1],qfrac);
447 //      delete gMyMinuit;
448         
449         
450  //
451  // One cluster for each maximum
452  //
453          for (j=0; j<2; j++) {
454              AliMUONRawCluster cnew;
455              if (fNPeaks == 0) {
456                  cnew.fNcluster[0]=-1;
457                  cnew.fNcluster[1]=fNRawClusters;
458              } else {
459                  cnew.fNcluster[0]=fNPeaks;
460                  cnew.fNcluster[1]=0;
461              }
462              cnew.fMultiplicity=0;
463              cnew.fX=Float_t(xrec[j]);
464              cnew.fY=Float_t(yrec[j]);
465              if (j==0) {
466                  cnew.fQ=Int_t(gChargeTot*qfrac);
467              } else {
468                  cnew.fQ=Int_t(gChargeTot*(1-qfrac));
469              }
470              gSegmentation->SetHit(xrec[j],yrec[j]);
471              for (i=0; i<mul; i++) {
472                  cnew.fIndexMap[cnew.fMultiplicity]=c->fIndexMap[i];
473                  gSegmentation->SetPad(gix[i], giy[i]);
474                  Float_t q1=gResponse->IntXY(gSegmentation);
475                  cnew.fContMap[cnew.fMultiplicity]=Float_t(q[i])/(q1*cnew.fQ);
476                  cnew.fMultiplicity++;
477              }
478              FillCluster(&cnew,0);
479              //printf("\n x,y %f %f ", cnew.fX, cnew.fY);
480              cnew.fClusterType=cnew.PhysicsContribution();
481              AddRawCluster(cnew);
482              fNPeaks++;
483          }
484      }
485
486      Bool_t fitted=kTRUE;
487
488      if (NLocal !=2 || !fitted) {
489  // Check if enough local clusters have been found,
490  // if not add global maxima to the list 
491  //
492          Int_t nPerMax;
493          if (NLocal!=0) {
494              nPerMax=mul/NLocal;
495          } else {
496              printf("\n Warning, no local maximum found \n");
497              nPerMax=fNperMax+1;
498          }
499
500          if (nPerMax > fNperMax) {
501              Int_t nGlob=mul/fNperMax-NLocal+1;
502              if (nGlob > 0) {
503                  Int_t nnew=0;
504                  for (i=0; i<mul; i++) {
505                      if (!IsLocal[i]) {
506                          IndLocal[NLocal]=i;
507                          IsLocal[i]=kTRUE;
508                          NLocal++;
509                          nnew++;
510                      }
511                      if (nnew==nGlob) break;
512                  }
513              }
514          }
515  //
516  // Associate hits to peaks
517  //
518          for (i=0; i<mul; i++) {
519              Float_t dmin=1.E10;
520              Float_t qmax=0;
521              if (IsLocal[i]) continue;
522              for (j=0; j<NLocal; j++) {
523                  Int_t il=IndLocal[j];
524                  Float_t d=TMath::Sqrt((x[i]-x[il])*(x[i]-x[il])
525                                    +(y[i]-y[il])*(y[i]-y[il]));
526                  Float_t ql=q[il];
527  //
528  // Select nearest peak
529  //
530                  if (d<dmin) {
531                      dmin=d;
532                      qmax=ql;
533                      AssocPeak[i]=j;
534                  } else if (d==dmin) {
535  //
536  // If more than one take highest peak
537  //
538                      if (ql>qmax) {
539                          dmin=d;
540                          qmax=ql;
541                          AssocPeak[i]=j;
542                      }
543                  }
544              }
545          }
546
547
548  //
549  // One cluster for each maximum
550  //
551          for (j=0; j<NLocal; j++) {
552              AliMUONRawCluster cnew;
553              if (fNPeaks == 0) {
554                  cnew.fNcluster[0]=-1;
555                  cnew.fNcluster[1]=fNRawClusters;
556              } else {
557                  cnew.fNcluster[0]=fNPeaks;
558                  cnew.fNcluster[1]=0;
559              }
560              cnew.fIndexMap[0]=c->fIndexMap[IndLocal[j]];
561              cnew.fMultiplicity=1;
562              for (i=0; i<mul; i++) {
563                  if (IsLocal[i]) continue;
564                  if (AssocPeak[i]==j) {
565                      cnew.fIndexMap[cnew.fMultiplicity]=c->fIndexMap[i];
566                      cnew.fMultiplicity++;
567                  }
568              }
569              FillCluster(&cnew);
570              cnew.fClusterType=cnew.PhysicsContribution();
571              AddRawCluster(cnew);
572              fNPeaks++;
573          }
574      }
575      
576 }
577
578
579 void  AliMUONClusterFinder::FillCluster(AliMUONRawCluster* c, Int_t flag) 
580 {
581 //
582 //  Completes cluster information starting from list of digits
583 //
584     AliMUONdigit* dig;
585     Float_t x, y;
586     Int_t  ix, iy;
587     Float_t frac=0;
588     
589     c->fPeakSignal=0;
590     if (flag) {
591         c->fX=0;
592         c->fY=0;
593         c->fQ=0;
594     }
595     //c->fQ=0;
596  
597
598     for (Int_t i=0; i<c->fMultiplicity; i++)
599     {
600         dig= (AliMUONdigit*)fDigits->UncheckedAt(c->fIndexMap[i]);
601         ix=dig->fPadX+c->fOffsetMap[i];
602         iy=dig->fPadY;
603         Int_t q=dig->fSignal;
604         if (dig->fPhysics >= dig->fSignal) {
605           c->fPhysicsMap[i]=2;
606         } else if (dig->fPhysics == 0) {
607           c->fPhysicsMap[i]=0;
608         } else  c->fPhysicsMap[i]=1;
609 //
610 //
611 // peak signal and track list
612         if (flag) {
613            if (q>c->fPeakSignal) {
614               c->fPeakSignal=q;
615 /*
616             c->fTracks[0]=dig->fTracks[0];
617             c->fTracks[1]=dig->fTracks[1];
618             c->fTracks[2]=dig->fTracks[2];
619 */
620               //c->fTracks[0]=dig->fTrack;
621             c->fTracks[0]=dig->fHit;
622             c->fTracks[1]=dig->fTracks[0];
623             c->fTracks[2]=dig->fTracks[1];
624            }
625         } else {
626            if (c->fContMap[i] > frac) {
627               frac=c->fContMap[i];
628               c->fPeakSignal=q;
629 /*
630             c->fTracks[0]=dig->fTracks[0];
631             c->fTracks[1]=dig->fTracks[1];
632             c->fTracks[2]=dig->fTracks[2];
633 */
634               //c->fTracks[0]=dig->fTrack;
635             c->fTracks[0]=dig->fHit;
636             c->fTracks[1]=dig->fTracks[0];
637             c->fTracks[2]=dig->fTracks[1];
638            }
639         }
640 //
641         if (flag) {
642             fSegmentation->GetPadCxy(ix, iy, x, y);
643             c->fX += q*x;
644             c->fY += q*y;
645             c->fQ += q;
646         }
647
648     } // loop over digits
649
650  if (flag) {
651      
652      c->fX/=c->fQ;
653      c->fX=fSegmentation->GetAnod(c->fX);
654      c->fY/=c->fQ; 
655 //
656 //  apply correction to the coordinate along the anode wire
657 //
658      x=c->fX;   
659      y=c->fY;
660      fSegmentation->GetPadIxy(x, y, ix, iy);
661      fSegmentation->GetPadCxy(ix, iy, x, y);
662      Int_t isec=fSegmentation->Sector(ix,iy);
663      TF1* CogCorr = fSegmentation->CorrFunc(isec-1);
664      
665      if (CogCorr) {
666          Float_t YonPad=(c->fY-y)/fSegmentation->Dpy(isec);
667          c->fY=c->fY-CogCorr->Eval(YonPad, 0, 0);
668      }
669  }
670 }
671
672
673 void  AliMUONClusterFinder::FindCluster(Int_t i, Int_t j, AliMUONRawCluster &c){
674 //
675 //  Find clusters
676 //
677 //
678 //  Add i,j as element of the cluster
679 //    
680     
681     Int_t idx = fHitMap->GetHitIndex(i,j);
682     AliMUONdigit* dig = (AliMUONdigit*) fHitMap->GetHit(i,j);
683     Int_t q=dig->fSignal;
684     if (q > TMath::Abs(c.fPeakSignal)) {
685         c.fPeakSignal=q;
686 /*
687         c.fTracks[0]=dig->fTracks[0];
688         c.fTracks[1]=dig->fTracks[1];
689         c.fTracks[2]=dig->fTracks[2];
690 */
691         //c.fTracks[0]=dig->fTrack;
692         c.fTracks[0]=dig->fHit;
693         c.fTracks[1]=dig->fTracks[0];
694         c.fTracks[2]=dig->fTracks[1];
695     }
696 //
697 //  Make sure that list of digits is ordered 
698 // 
699     Int_t mu=c.fMultiplicity;
700     c.fIndexMap[mu]=idx;
701
702     if (dig->fPhysics >= dig->fSignal) {
703         c.fPhysicsMap[mu]=2;
704     } else if (dig->fPhysics == 0) {
705         c.fPhysicsMap[mu]=0;
706     } else  c.fPhysicsMap[mu]=1;
707
708     if (mu > 0) {
709         for (Int_t ind=mu-1; ind>=0; ind--) {
710             Int_t ist=(c.fIndexMap)[ind];
711             Int_t ql=((AliMUONdigit*)fDigits
712                       ->UncheckedAt(ist))->fSignal;
713             if (q>ql) {
714                 c.fIndexMap[ind]=idx;
715                 c.fIndexMap[ind+1]=ist;
716             } else {
717                 break;
718             }
719         }
720     }
721     
722     c.fMultiplicity++;
723     
724     if (c.fMultiplicity >= 50 ) {
725         printf("FindCluster - multiplicity >50  %d \n",c.fMultiplicity);
726         c.fMultiplicity=49;
727     }
728
729 // Prepare center of gravity calculation
730     Float_t x, y;
731     fSegmentation->GetPadCxy(i, j, x, y);
732     c.fX += q*x;
733     c.fY += q*y;
734     c.fQ += q;
735 // Flag hit as taken  
736     fHitMap->FlagHit(i,j);
737 //
738 //  Now look recursively for all neighbours
739 //  
740     Int_t nn;
741     Int_t Xlist[kMaxNeighbours], Ylist[kMaxNeighbours];
742     fSegmentation->Neighbours(i,j,&nn,Xlist,Ylist);
743     for (Int_t in=0; in<nn; in++) {
744         Int_t ix=Xlist[in];
745         Int_t iy=Ylist[in];
746         if (fHitMap->TestHit(ix,iy)==unused) FindCluster(ix, iy, c);
747     }
748 }
749
750 //_____________________________________________________________________________
751
752 void AliMUONClusterFinder::FindRawClusters()
753 {
754   //
755   // simple MUON cluster finder from digits -- finds neighbours and 
756   // fill the tree with raw clusters
757   //
758     if (!fNdigits) return;
759
760     fHitMap = new AliMUONHitMapA1(fSegmentation, fDigits);
761
762     AliMUONdigit *dig;
763
764     Int_t ndig;
765     Int_t nskip=0;
766     Int_t ncls=0;
767     fHitMap->FillHits();
768     for (ndig=0; ndig<fNdigits; ndig++) {
769         dig = (AliMUONdigit*)fDigits->UncheckedAt(ndig);
770         Int_t i=dig->fPadX;
771         Int_t j=dig->fPadY;
772         if (fHitMap->TestHit(i,j)==used ||fHitMap->TestHit(i,j)==empty) {
773             nskip++;
774             continue;
775         }
776         AliMUONRawCluster c;
777         c.fMultiplicity=0;
778         c.fPeakSignal=dig->fSignal;
779 /*
780         c.fTracks[0]=dig->fTracks[0];
781         c.fTracks[1]=dig->fTracks[1];
782         c.fTracks[2]=dig->fTracks[2];
783 */
784         //c.fTracks[0]=dig->fTrack;
785         c.fTracks[0]=dig->fHit;
786         c.fTracks[1]=dig->fTracks[0];
787         c.fTracks[2]=dig->fTracks[1];
788         // tag the beginning of cluster list in a raw cluster
789         c.fNcluster[0]=-1;
790         FindCluster(i,j, c);
791         // center of gravity
792         c.fX /= c.fQ;
793         c.fX=fSegmentation->GetAnod(c.fX);
794         c.fY /= c.fQ;
795 //
796 //  apply correction to the coordinate along the anode wire
797 //
798         Int_t ix,iy;
799         Float_t x=c.fX;   
800         Float_t y=c.fY;
801         fSegmentation->GetPadIxy(x, y, ix, iy);
802         fSegmentation->GetPadCxy(ix, iy, x, y);
803         Int_t isec=fSegmentation->Sector(ix,iy);
804         TF1* CogCorr=fSegmentation->CorrFunc(isec-1);
805         if (CogCorr) {
806             Float_t YonPad=(c.fY-y)/fSegmentation->Dpy(isec);
807             c.fY=c.fY-CogCorr->Eval(YonPad,0,0);
808         }
809
810 //
811 //      Analyse cluster and decluster if necessary
812 //      
813     ncls++;
814     c.fNcluster[1]=fNRawClusters;
815     c.fClusterType=c.PhysicsContribution();
816     Decluster(&c);
817     fNPeaks=0;
818 //
819 //
820 //
821 //      reset Cluster object
822         for (int k=0;k<c.fMultiplicity;k++) {
823             c.fIndexMap[k]=0;
824         }
825         c.fMultiplicity=0;
826     } // end loop ndig    
827     delete fHitMap;
828 }
829
830 void AliMUONClusterFinder::
831 CalibrateCOG()
832 {
833     Float_t x[5];
834     Float_t y[5];
835     Int_t n, i;
836     TF1 func;
837     if (fSegmentation) {
838         fSegmentation->GiveTestPoints(n, x, y);
839         for (i=0; i<n; i++) {
840             Float_t xtest=x[i];
841             Float_t ytest=y[i];     
842             SinoidalFit(xtest, ytest, func);
843             fSegmentation->SetCorrFunc(i, new TF1(func));
844         }
845     }
846 }
847
848
849 void AliMUONClusterFinder::
850 SinoidalFit(Float_t x, Float_t y, TF1 &func)
851 {
852 //
853     static Int_t count=0;
854     char canvasname[3];
855     count++;
856     sprintf(canvasname,"c%d",count);
857
858     const Int_t ns=101;
859     Float_t xg[ns], yg[ns], xrg[ns], yrg[ns];
860     Float_t xsig[ns], ysig[ns];
861    
862     AliMUONsegmentation *segmentation=fSegmentation;
863
864     Int_t ix,iy;
865     segmentation->GetPadIxy(x,y,ix,iy);   
866     segmentation->GetPadCxy(ix,iy,x,y);   
867     Int_t isec=segmentation->Sector(ix,iy);
868 // Pad Limits    
869     Float_t xmin = x-segmentation->Dpx(isec)/2;
870     Float_t ymin = y-segmentation->Dpy(isec)/2;
871 //              
872 //      Integration Limits
873     Float_t dxI=fResponse->SigmaIntegration()*fResponse->ChargeSpreadX();
874     Float_t dyI=fResponse->SigmaIntegration()*fResponse->ChargeSpreadY();
875
876 //
877 //  Scanning
878 //
879     Int_t i;
880     Float_t qp;
881 //
882 //  y-position
883     Float_t yscan=ymin;
884     Float_t dy=segmentation->Dpy(isec)/(ns-1);
885
886     for (i=0; i<ns; i++) {
887 //
888 //      Pad Loop
889 //      
890         Float_t sum=0;
891         Float_t qcheck=0;
892         segmentation->SigGenInit(x, yscan, 0);
893         
894         for (segmentation->FirstPad(x, yscan, dxI, dyI); 
895              segmentation->MorePads(); 
896              segmentation->NextPad()) 
897         {
898             qp=fResponse->IntXY(segmentation);
899             qp=TMath::Abs(qp);
900 //
901 //
902             if (qp > 1.e-4) {
903                 qcheck+=qp;
904                 Int_t ixs=segmentation->Ix();
905                 Int_t iys=segmentation->Iy();
906                 Float_t xs,ys;
907                 segmentation->GetPadCxy(ixs,iys,xs,ys);
908                 sum+=qp*ys;
909             }
910         } // Pad loop
911         Float_t ycog=sum/qcheck;
912         yg[i]=(yscan-y)/segmentation->Dpy(isec);
913         yrg[i]=(ycog-y)/segmentation->Dpy(isec);
914         ysig[i]=ycog-yscan;
915         yscan+=dy;
916     } // scan loop
917 //
918 //  x-position
919     Float_t xscan=xmin;
920     Float_t dx=segmentation->Dpx(isec)/(ns-1);
921
922     for (i=0; i<ns; i++) {
923 //
924 //      Pad Loop
925 //      
926         Float_t sum=0;
927         Float_t qcheck=0;
928         segmentation->SigGenInit(xscan, y, 0);
929         
930         for (segmentation->FirstPad(xscan, y, dxI, dyI); 
931              segmentation->MorePads(); 
932              segmentation->NextPad()) 
933         {
934             qp=fResponse->IntXY(segmentation);
935             qp=TMath::Abs(qp);
936 //
937 //
938             if (qp > 1.e-2) {
939                 qcheck+=qp;
940                 Int_t ixs=segmentation->Ix();
941                 Int_t iys=segmentation->Iy();
942                 Float_t xs,ys;
943                 segmentation->GetPadCxy(ixs,iys,xs,ys);
944                 sum+=qp*xs;
945             }
946         } // Pad loop
947         Float_t xcog=sum/qcheck;
948         xcog=segmentation->GetAnod(xcog);
949         
950         xg[i]=(xscan-x)/segmentation->Dpx(isec);
951         xrg[i]=(xcog-x)/segmentation->Dpx(isec);
952         xsig[i]=xcog-xscan;
953         xscan+=dx;
954     }
955 //
956 // Creates a Root function based on function sinoid above
957 // and perform the fit
958 //
959     //    TGraph *graphx = new TGraph(ns,xg ,xsig);
960     //    TGraph *graphxr= new TGraph(ns,xrg,xsig);   
961     //    TGraph *graphy = new TGraph(ns,yg ,ysig);
962     TGraph *graphyr= new TGraph(ns,yrg,ysig);
963
964     Double_t sinoid(Double_t *x, Double_t *par);
965     new TF1("sinoidf",sinoid,0.5,0.5,5);
966     graphyr->Fit("sinoidf","Q");
967     func = *((TF1*)((graphyr->GetListOfFunctions())->At(0)));
968 /*
969     
970     TCanvas *c1=new TCanvas(canvasname,canvasname,400,10,600,700);
971     TPad* pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99);
972     TPad* pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99);
973     TPad* pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49);
974     TPad* pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49);
975     pad11->SetFillColor(11);
976     pad12->SetFillColor(11);
977     pad13->SetFillColor(11);
978     pad14->SetFillColor(11);
979     pad11->Draw();
980     pad12->Draw();
981     pad13->Draw();
982     pad14->Draw();
983     
984 //
985     pad11->cd();
986     graphx->SetFillColor(42);
987     graphx->SetMarkerColor(4);
988     graphx->SetMarkerStyle(21);
989     graphx->Draw("AC");
990     graphx->GetHistogram()->SetXTitle("x on pad");
991     graphx->GetHistogram()->SetYTitle("xcog-x");   
992
993
994     pad12->cd();
995     graphxr->SetFillColor(42);
996     graphxr->SetMarkerColor(4);
997     graphxr->SetMarkerStyle(21);
998     graphxr->Draw("AP");
999     graphxr->GetHistogram()->SetXTitle("xcog on pad");
1000     graphxr->GetHistogram()->SetYTitle("xcog-x");   
1001     
1002
1003     pad13->cd();
1004     graphy->SetFillColor(42);
1005     graphy->SetMarkerColor(4);
1006     graphy->SetMarkerStyle(21);
1007     graphy->Draw("AF");
1008     graphy->GetHistogram()->SetXTitle("y on pad");
1009     graphy->GetHistogram()->SetYTitle("ycog-y");   
1010     
1011
1012
1013     pad14->cd();
1014     graphyr->SetFillColor(42);
1015     graphyr->SetMarkerColor(4);
1016     graphyr->SetMarkerStyle(21);
1017     graphyr->Draw("AF");
1018     graphyr->GetHistogram()->SetXTitle("ycog on pad");
1019     graphyr->GetHistogram()->SetYTitle("ycog-y");   
1020     
1021     c1->Update();
1022 */
1023 }
1024
1025 Double_t sinoid(Double_t *x, Double_t *par)
1026 {
1027     Double_t arg = -2*TMath::Pi()*x[0];
1028     Double_t fitval= par[0]*TMath::Sin(arg)+
1029         par[1]*TMath::Sin(2*arg)+
1030         par[2]*TMath::Sin(3*arg)+
1031         par[3]*TMath::Sin(4*arg)+
1032         par[4]*TMath::Sin(5*arg);
1033     return fitval;
1034  }
1035
1036
1037 Double_t DoubleGauss(Double_t *x, Double_t *par)
1038 {
1039     Double_t arg1 = (x[0]-par[1])/0.18;
1040     Double_t arg2 = (x[0]-par[3])/0.18;
1041     Double_t fitval= par[0]*TMath::Exp(-arg1*arg1/2)
1042         +par[2]*TMath::Exp(-arg2*arg2/2);
1043     return fitval;
1044  }
1045
1046 Float_t DiscrCharge(Int_t i,Double_t *par) 
1047 {
1048 // par[0]    x-position of first  cluster
1049 // par[1]    y-position of first  cluster
1050 // par[2]    x-position of second cluster
1051 // par[3]    y-position of second cluster
1052 // par[4]    charge fraction of first  cluster
1053 // 1-par[4]  charge fraction of second cluster
1054
1055     static Float_t qtot;
1056     if (gFirst) {
1057         qtot=0;
1058         for (Int_t jbin=0; jbin<gNbins; jbin++) {
1059             qtot+=gCharge[jbin];
1060         }
1061         gFirst=0;
1062         //printf("\n sum of charge from DiscrCharge %f\n", qtot);
1063         gChargeTot=Int_t(qtot);
1064         
1065     }
1066     gSegmentation->SetPad(gix[i], giy[i]);
1067 //  First Cluster
1068     gSegmentation->SetHit(par[0],par[1]);
1069     Float_t q1=gResponse->IntXY(gSegmentation);
1070     
1071 //  Second Cluster
1072     gSegmentation->SetHit(par[2],par[3]);
1073     Float_t q2=gResponse->IntXY(gSegmentation);
1074     
1075     Float_t value = qtot*(par[4]*q1+(1.-par[4])*q2);
1076     return value;
1077 }
1078
1079 //
1080 // Minimisation function
1081 void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
1082 {
1083     Int_t i;
1084     Float_t delta;
1085     Float_t chisq=0;
1086     Float_t qcont=0;
1087     Float_t qtot=0;
1088     
1089     for (i=0; i<gNbins; i++) {
1090         Float_t q0=gCharge[i];
1091         Float_t q1=DiscrCharge(i,par);
1092         delta=(q0-q1)/TMath::Sqrt(q0);
1093         chisq+=delta*delta;
1094         qcont+=q1;
1095         qtot+=q0;
1096     }
1097     chisq=chisq+=(qtot-qcont)*(qtot-qcont)*0.5;
1098     f=chisq;
1099 }
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110