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