Introduction of the Copyright and cvs Log
[u/mrichter/AliRoot.git] / MUON / AliMUONClusterFinderv0.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /*
17 $Log$
18 */
19
20 #include "AliMUONClusterFinderv0.h"
21 #include "AliMUONSegResV1.h"
22 //#include "TTree.h"
23 //#include "AliRun.h"
24 //#include <TCanvas.h>
25 //#include <TH1.h>
26 //#include <TPad.h>
27 //#include <TGraph.h> 
28
29 //----------------------------------------------------------
30 ClassImp(AliMUONClusterFinderv0)
31
32     AliMUONClusterFinderv0::AliMUONClusterFinderv0
33 (AliMUONsegmentation *segmentation, AliMUONresponse *response, 
34  TClonesArray *digits, Int_t chamber) : AliMUONClusterFinder(segmentation,response,digits,chamber)
35 {;}
36
37     AliMUONClusterFinderv0::AliMUONClusterFinderv0():AliMUONClusterFinder()
38 {;}
39
40 /*
41 void AliMUONClusterFinder::AddRawCluster(const AliMUONRawCluster c)
42 {
43   //
44   // Add a raw cluster copy to the list
45   //
46     AliMUON *MUON=(AliMUON*)gAlice->GetModule("MUON");
47     MUON->AddRawCluster(fChamber,c); 
48     fNRawClusters++;
49 }
50 */
51
52
53
54 void AliMUONClusterFinderv0::Decluster(AliMUONRawCluster *cluster)
55 {
56 //    AliMUONdigit *dig;
57 //    Int_t q;
58     static int done=0;
59     if (!done) {
60         printf("Calling decluster\n");
61         done=1;
62     }
63     
64
65     
66     Int_t mul = cluster->fMultiplicity;
67 //    printf("Decluster - multiplicity   %d \n",mul);
68
69     if (mul == 1) {
70 //      printf("\n Nothing special for 1-clusters \n");
71 //
72 // Nothing special for 1-clusters
73 //
74         AddRawCluster(*cluster); 
75     } else if (mul ==2) {
76 //
77 // 2-cluster, compute offset
78 //
79         SetOffset(cluster);
80         FillCluster(cluster);
81         AddRawCluster(*cluster); 
82     } else if (mul ==3) {
83 //
84 // 3-cluster, check topology
85 //      printf("\n 3-cluster, check topology \n");
86 //
87         if (Centered(cluster)) {
88 //
89 // ok, cluster is centered 
90 //          printf("\n ok, cluster is centered \n");
91         } else {
92 //
93 // cluster is not centered, split into 2+1
94 //          printf("\n cluster is not centered, split into 2+1 \n");
95         }
96             
97     } else {
98         if (mul >(50-5)) printf("Decluster - multiplicity %d approaching 50\n",mul);
99 // 
100 // 4-and more-pad clusters
101 //
102         SplitByLocalMaxima(cluster);
103     } // multiplicity 
104 }
105
106 Int_t AliMUONClusterFinderv0::PeakOffsetAndCoordinates(Int_t DigitIndex, Float_t *X, Float_t *Y)
107 //
108 // Computes for which allowed offsets the digit has the highest neighbouring charge
109 // Returns the value of the offset, and sets the pyisical coordinates of that pad
110 // Loop on physical neighbours is specific to AliMUONsegmentationV1
111 {
112 Int_t nPara, offset, returnOffset=0 ;
113 AliMUONdigit* dig= (AliMUONdigit*)fDigits->UncheckedAt(DigitIndex);
114 AliMUONsegmentationV1* seg = (AliMUONsegmentationV1*) fSegmentation;
115 seg->GetNParallelAndOffset(dig->fPadX,dig->fPadY,&nPara,&offset);
116 if (nPara>1)
117   {
118   Float_t qMax=0;
119   for (Int_t i=0;i<nPara; i++)
120     {
121     // Compute the charge on the 9 neighbouring pads
122     // We assume that there are no pads connected in parallel in the neighbourhood
123     //
124     Float_t q=0;
125     for (Int_t dx=-1;dx<2;dx++)
126       for (Int_t dy=-1;dy<2;dy++)
127         {
128         if (dx==dy && dy==0)
129           continue; 
130         Int_t padY=dig->fPadY+dy;
131         Int_t padX=seg->Ix((Int_t) (dig->fPadX+dx+i*offset) , padY);
132         if (fHitMap->TestHit(padX, padY)==empty) 
133            continue;
134         AliMUONdigit* digt = (AliMUONdigit*) fHitMap->GetHit(padX,padY);
135         q += digt->fSignal;
136         }
137     if (q>qMax)
138       {
139       returnOffset=i*offset;
140       qMax=q;
141       }
142     }
143   }
144 fSegmentation->GetPadCxy(dig->fPadX+returnOffset,dig->fPadY,*X,*Y);
145 return returnOffset;
146 }
147
148
149 void AliMUONClusterFinderv0::SetOffset(AliMUONRawCluster *cluster)
150 // compute the offsets assuming that there is only one peak !
151 {
152 //DumpCluster(cluster);
153 Float_t X,Y;
154 cluster->fOffsetMap[0]=PeakOffsetAndCoordinates(cluster->fIndexMap[0],&X,&Y);
155 for (Int_t i=1;i<cluster->fMultiplicity;i++) {
156   AliMUONdigit* dig= (AliMUONdigit*)fDigits->UncheckedAt(cluster->fIndexMap[i]);
157   fSegmentation->Distance2AndOffset(dig->fPadX,dig->fPadY,X,Y,&(cluster->fOffsetMap[i]));
158   }
159 }
160
161 void AliMUONClusterFinderv0::DumpCluster(AliMUONRawCluster *cluster)
162 {    
163 printf ("other cluster\n");
164 for (Int_t i=0; i<cluster->fMultiplicity; i++)
165     {
166         AliMUONdigit* dig= (AliMUONdigit*)fDigits->UncheckedAt(cluster->fIndexMap[i]);
167         Int_t nPara, offset;
168         fSegmentation->GetNParallelAndOffset(dig->fPadX,dig->fPadY,&nPara,&offset);
169
170         printf("X %d Y %d Q %d NPara %d \n",dig->fPadX, dig->fPadY,dig->fSignal, nPara);
171     }
172 }
173
174 Bool_t AliMUONClusterFinderv0::Centered(AliMUONRawCluster *cluster)
175 {
176     AliMUONdigit* dig;
177     dig= (AliMUONdigit*)fDigits->UncheckedAt(cluster->fIndexMap[0]);
178     Int_t ix=dig->fPadX;
179     Int_t iy=dig->fPadY;
180     Int_t nn;
181     Int_t X[kMaxNeighbours], Y[kMaxNeighbours], XN[kMaxNeighbours], YN[kMaxNeighbours];
182     fSegmentation->Neighbours(ix,iy,&nn,X,Y);
183
184     Int_t nd=0;
185     for (Int_t i=0; i<nn; i++) {
186         if (fHitMap->TestHit(X[i],Y[i]) == used) {
187             XN[nd]=X[i];
188             YN[nd]=Y[i];
189             nd++;
190         }
191     }
192     if (nd==2) {
193 //
194 // cluster is centered !
195         SetOffset(cluster);
196         FillCluster(cluster);
197         AddRawCluster(*cluster);
198         return kTRUE;
199     } else if (nd ==1) {
200 //
201 // Highest signal on an edge, split cluster into 2+1
202 //
203 // who is the neighbour ?
204         Int_t nind=fHitMap->GetHitIndex(XN[0], YN[0]);
205         Int_t i1= (nind==cluster->fIndexMap[1]) ? 1:2;
206         Int_t i2= (nind==cluster->fIndexMap[1]) ? 2:1;    
207 //
208 // 2-cluster
209         AliMUONRawCluster cnew;
210         cnew.fMultiplicity=2;
211         cnew.fIndexMap[0]=cluster->fIndexMap[0];
212         cnew.fIndexMap[1]=cluster->fIndexMap[i1];
213         SetOffset(&cnew);
214         FillCluster(&cnew);
215         AddRawCluster(cnew);
216 //
217 // 1-cluster
218         cluster->fMultiplicity=1;
219         cluster->fIndexMap[0]=cluster->fIndexMap[i2];
220         cluster->fIndexMap[1]=0;
221         cluster->fIndexMap[2]=0;        
222         FillCluster(cluster);
223         AddRawCluster(*cluster);
224         return kFALSE;
225     } else {
226         printf("\n Completely screwed up %d !! \n",nd);
227         
228     }
229     
230         return kFALSE;
231 }
232
233
234 void AliMUONClusterFinderv0::SplitByLocalMaxima(AliMUONRawCluster *c)
235 {
236     AliMUONdigit* dig[50], *digt;
237     Int_t ix[50], iy[50], q[50];
238     Float_t x[50], y[50];
239     Int_t i; // loops over digits
240     Int_t j; // loops over local maxima
241     
242     Int_t mul=c->fMultiplicity;
243 //
244 //  dump digit information into arrays
245 //
246     for (i=0; i<mul; i++)
247     {
248         dig[i]= (AliMUONdigit*)fDigits->UncheckedAt(c->fIndexMap[i]);
249         ix[i]= dig[i]->fPadX;
250         iy[i]= dig[i]->fPadY;
251         q[i] = dig[i]->fSignal;
252         fSegmentation->GetPadCxy(ix[i], iy[i], x[i], y[i]);
253     }
254 //
255 //  Find local maxima
256 //
257     Bool_t IsLocal[50];
258     Int_t NLocal=0;
259     Int_t AssocPeak[50];
260     Int_t IndLocal[50];
261     Int_t nn;
262     Int_t X[kMaxNeighbours], Y[kMaxNeighbours];
263     for (i=0; i<mul; i++) {
264         fSegmentation->Neighbours(ix[i], iy[i], &nn, X, Y);
265         IsLocal[i]=kTRUE;
266         for (j=0; j<nn; j++) {
267             if (fHitMap->TestHit(X[j], Y[j])==empty) continue;
268             digt=(AliMUONdigit*) fHitMap->GetHit(X[j], Y[j]);
269             if (digt->fSignal > q[i]) {
270                 IsLocal[i]=kFALSE;
271                 break;
272 //
273 // handle special case of neighbouring pads with equal signal
274             } else if (digt->fSignal == q[i]) {
275                 if (NLocal >0) {
276                     for (Int_t k=0; k<NLocal; k++) {
277                         if (X[j]==ix[IndLocal[k]] && Y[j]==iy[IndLocal[k]]){
278                             IsLocal[i]=kFALSE;
279                         }
280                     }
281                 }
282             } 
283         } // loop over next neighbours
284         if (IsLocal[i]) {
285             IndLocal[NLocal]=i;
286             // New for LYON : we guess which is the actual position of the pad hit
287             // But this would run like that for normal chamber !
288             c->fOffsetMap[i]=PeakOffsetAndCoordinates(c->fIndexMap[i], &(x[i]), &(y[i]));
289             NLocal++;
290         } 
291     } // loop over all digits
292 //    printf("Found %d local Maxima",NLocal);
293 // 
294 // Check if enough local clusters have been found,
295 // if not add global maxima to the list 
296 //
297 // But what the hell is that ? (Manu)
298 //
299 //    Int_t nPerMax=mul/NLocal;
300 //    if (nPerMax > 5) {
301 //      Int_t nGlob=mul/5-NLocal+1;
302 //      if (nGlob > 0) {
303 //          Int_t nnew=0;
304 //          for (i=0; i<mul; i++) {
305 //              if (!IsLocal[i]) {
306 //                  IndLocal[NLocal]=i;
307 //                  IsLocal[i]=kTRUE;
308 //                  NLocal++;
309 //                  nnew++;
310 //              }
311 //              if (nnew==nGlob) break;
312 //          }
313 //      }
314 //   }
315 //
316 //
317 // Associate hits to peaks
318 //
319     for (i=0; i<mul; i++) {
320         //
321         // loop on digits
322         //
323         Float_t dmin=1.E10;
324         Float_t qmax=0;
325         Int_t offset;
326         if (IsLocal[i]) continue;
327         for (j=0; j<NLocal; j++) {
328             //
329             // Loop on peaks
330             //
331             Int_t il=IndLocal[j];
332 //          Float_t d=TMath::Sqrt((x[i]-x[il])*(x[i]-x[il])
333 //                                +(y[i]-y[il])*(y[i]-y[il]));
334             // Can run like that for non-Lyon chambers
335             Float_t d = fSegmentation->Distance2AndOffset(ix[i],iy[i],x[il],y[il], &offset);
336             Float_t ql=q[il];
337 //
338 // Select nearest peak
339 //
340             if (d<dmin) {
341                 dmin=d;
342                 qmax=ql;
343                 AssocPeak[i]=j;
344                 c->fOffsetMap[i]=offset;
345             } else if (d==dmin) {
346 //
347 // If more than one take highest peak
348 //
349                 if (ql>qmax) {
350                     dmin=d;
351                     qmax=ql;
352                     AssocPeak[i]=j;
353                     c->fOffsetMap[i]=offset;
354                 }
355             } // end if
356         } // End loop on peaks
357     } // end loop on digits
358 //
359 // One cluster for each maximum
360 //
361     for (j=0; j<NLocal; j++) {
362         AliMUONRawCluster cnew;
363         cnew.fIndexMap[0]=c->fIndexMap[IndLocal[j]];
364         cnew.fOffsetMap[0]=c->fOffsetMap[IndLocal[j]];
365         cnew.fMultiplicity=1;
366         for (i=0; i<mul; i++) {
367             if (IsLocal[i]) continue;
368             if (AssocPeak[i]==j) {
369                 cnew.fIndexMap[cnew.fMultiplicity]=c->fIndexMap[i];
370                 cnew.fOffsetMap[cnew.fMultiplicity]=c->fOffsetMap[i];
371                 cnew.fMultiplicity++;
372             }
373         }
374         FillCluster(&cnew);
375         AddRawCluster(cnew);
376     }
377 }
378
379 /*
380 void  AliMUONClusterFinderv0::FillCluster(AliMUONRawCluster* c)
381 {
382 //
383 //  Completes cluster information starting from list of digits
384 //
385     AliMUONdigit* dig;
386     Float_t x, y;
387     Int_t  ix, iy;
388     
389     c->fPeakSignal=0;
390     c->fX=0;
391     c->fY=0;
392     c->fQ=0;
393     for (Int_t i=0; i<c->fMultiplicity; i++)
394     {
395         dig= (AliMUONdigit*)fDigits->UncheckedAt(c->fIndexMap[i]);
396         ix=dig->fPadX + c.fOffsetMap[i]; // should be 0 for non-LYON
397         iy=dig->fPadY;
398         Int_t q=dig->fSignal;
399 //
400 //
401 // peak signal and track list
402         if (q>c->fPeakSignal) {
403             c->fPeakSignal=0;
404             c->fTracks[0]=dig->fTracks[0];
405             c->fTracks[1]=dig->fTracks[1];
406             c->fTracks[2]=dig->fTracks[2];
407         }
408 //
409 // centre of gravity
410         fSegmentation->GetPadCxy(ix, iy, x, y);
411         c->fX += q*x;
412         c->fY += q*y;
413         c->fQ += q;
414     }
415     c->fX/=c->fQ;
416 // Not valid for inclined tracks in X !!! (Manu)
417 //    c->fX=fSegmentation->GetAnod(c->fX);
418     c->fY/=c->fQ; 
419 //
420 //  apply correction to the coordinate along the anode wire
421 //
422     if (fCogCorr) {
423         x=c->fX;   
424         y=c->fY;
425         fSegmentation->GetPadIxy(x, y, ix, iy);
426         fSegmentation->GetPadCxy(ix, iy, x, y);
427         Float_t YonPad=(c->fY-y)/fSegmentation->Dpy();
428         c->fY=y-fCogCorr->Eval(YonPad, 0, 0);
429     }
430
431 }
432 */
433 /*
434 void  AliMUONClusterFinder::FindCluster(Int_t i, Int_t j, AliMUONRawCluster &c){
435 //
436 //  Find clusters
437 //
438 //
439 //  Add i,j as element of the cluster
440 //    
441     
442     Int_t idx = fHitMap->GetHitIndex(i,j);
443     AliMUONdigit* dig = (AliMUONdigit*) fHitMap->GetHit(i,j);
444     Int_t q=dig->fSignal;
445     if (q > TMath::Abs(c.fPeakSignal)) {
446         c.fPeakSignal=q;
447         c.fTracks[0]=dig->fTracks[0];
448         c.fTracks[1]=dig->fTracks[1];
449         c.fTracks[2]=dig->fTracks[2];
450     }
451 //
452 //  Make sure that list of digits is ordered 
453 // 
454     Int_t mu=c.fMultiplicity;
455     c.fIndexMap[mu]=idx;
456
457     if (mu > 0) {
458         for (Int_t ind=mu-1; ind>=0; ind--) {
459             Int_t ist=(c.fIndexMap)[ind];
460             Int_t ql=((AliMUONdigit*)fDigits
461                       ->UncheckedAt(ist))->fSignal;
462             if (q>ql) {
463                 c.fIndexMap[ind]=idx;
464                 c.fIndexMap[ind+1]=ist;
465             } else {
466                 break;
467             }
468         }
469     }
470     
471     c.fMultiplicity++;
472     
473     if (c.fMultiplicity >= 50 ) {
474         printf("FindCluster - multiplicity >50  %d \n",c.fMultiplicity);
475         c.fMultiplicity=50-1;
476     }
477
478 // Prepare center of gravity calculation
479     Float_t x, y;
480     fSegmentation->GetPadCxy(i, j, x, y);
481     c.fX += q*x;
482     c.fY += q*y;
483     c.fQ += q;
484 // Flag hit as taken  
485     fHitMap->FlagHit(i,j);
486 //
487 //  Now look recursively for all neighbours
488 //  
489     Int_t nn;
490     Int_t Xlist[kMaxNeighbours], Ylist[kMaxNeighbours];
491     fSegmentation->Neighbours(i,j,&nn,Xlist,Ylist);
492     for (Int_t in=0; in<nn; in++) {
493         Int_t ix=Xlist[in];
494         Int_t iy=Ylist[in];
495         if (fHitMap->TestHit(ix,iy)==unused) FindCluster(ix, iy, c);
496     }
497 }
498 */
499
500 //_____________________________________________________________________________
501
502 void AliMUONClusterFinderv0::FindRawClusters()
503 {
504   //
505   // simple MUON cluster finder from digits -- finds neighbours and 
506   // fill the tree with raw clusters
507   //
508     if (!fNdigits) return;
509
510     fHitMap = new AliMUONHitMapA1(fSegmentation, fDigits);
511
512     AliMUONdigit *dig;
513
514     int ndig;
515     int nskip=0;
516
517     fHitMap->FillHits();
518     for (ndig=0; ndig<fNdigits; ndig++) {
519         dig = (AliMUONdigit*)fDigits->UncheckedAt(ndig);
520         Int_t i=dig->fPadX;
521         Int_t j=dig->fPadY;
522         if (fHitMap->TestHit(i,j)==used ||fHitMap->TestHit(i,j)==empty) {
523             nskip++;
524             continue;
525         }
526         AliMUONRawCluster c;
527         c.fMultiplicity=0;
528 //      c.fPeakSignal=dig->fSignal;
529 //      c.fTracks[0]=dig->fTracks[0];
530 //      c.fTracks[1]=dig->fTracks[1];
531 //      c.fTracks[2]=dig->fTracks[2];
532         c.fPeakSignal=0;
533         FindCluster(i,j, c);
534         // center of gravity
535         c.fX /= c.fQ;
536         c.fX=fSegmentation->GetAnod(c.fX);
537         c.fY /= c.fQ;
538 //
539 //  apply correction to the coordinate along the anode wire
540 //
541
542         
543     if (fCogCorr) {
544         Int_t ix,iy;
545         Float_t x=c.fX;   
546         Float_t y=c.fY;
547         fSegmentation->GetPadIxy(x, y, ix, iy);
548         fSegmentation->GetPadCxy(ix, iy, x, y);
549         Float_t YonPad=(c.fY-y)/fSegmentation->Dpy();
550         c.fY=y-fCogCorr->Eval(YonPad,0,0);
551     }
552 //
553 //      Analyse cluster and decluster if necessary
554 //      
555         Decluster(&c);
556 //
557 //
558 //
559 //      reset Cluster object
560         for (int k=0;k<c.fMultiplicity;k++) {
561             c.fIndexMap[k]=0;
562             c.fOffsetMap[k]=0;
563         }
564         c.fMultiplicity=0;
565     } // end loop ndig    
566     delete fHitMap;
567 }
568
569 /*
570
571 void AliMUONClusterFinder::
572 CalibrateCOG()
573 {
574     Float_t x[5];
575     Float_t y[5];
576     Int_t n, i;
577     TF1 func;
578     
579     if (fSegmentation) {
580         fSegmentation->GiveTestPoints(n, x, y);
581         for (i=0; i<n; i++) {
582             Float_t xtest=x[i];
583             Float_t ytest=y[i];     
584             SinoidalFit(xtest, ytest, func);
585         }
586         fCogCorr = new TF1(func);
587     }
588 }
589 */
590 /*
591
592 void AliMUONClusterFinder::
593 SinoidalFit(Float_t x, Float_t y, TF1 &func)
594 {
595 //
596     static Int_t count=0;
597     char canvasname[3];
598     count++;
599     sprintf(canvasname,"c%d",count);
600
601 // MANU : without const, error on HP
602     const Int_t ns=101;
603     Float_t xg[ns], yg[ns], xrg[ns], yrg[ns];
604     Float_t xsig[ns], ysig[ns];
605    
606     AliMUONsegmentation *segmentation=fSegmentation;
607
608     Int_t ix,iy;
609     segmentation->GetPadIxy(x,y,ix,iy);   
610     segmentation->GetPadCxy(ix,iy,x,y);   
611     Int_t isec=segmentation->Sector(ix,iy);
612 // Pad Limits    
613     Float_t xmin = x-segmentation->GetRealDpx(isec)/2;
614     Float_t ymin = y-segmentation->Dpy()/2;
615 //              
616 //      Integration Limits
617     Float_t dxI=fResponse->Nsigma()*fResponse->ChwX();
618     Float_t dyI=fResponse->Nsigma()*fResponse->ChwY();
619
620 //
621 //  Scanning
622 //
623     Int_t i;
624     Float_t qp;
625 //
626 //  y-position
627     Float_t yscan=ymin;
628     Float_t dy=segmentation->Dpy()/(ns-1);
629
630     for (i=0; i<ns; i++) {
631 //
632 //      Pad Loop
633 //      
634         Float_t sum=0;
635         Float_t qcheck=0;
636         segmentation->SigGenInit(x, yscan, 0);
637         
638         for (segmentation->FirstPad(x, yscan, dxI, dyI); 
639              segmentation->MorePads(); 
640              segmentation->NextPad()) 
641         {
642             qp=fResponse->IntXY(segmentation);
643             qp=TMath::Abs(qp);
644 //
645 //
646             if (qp > 1.e-4) {
647                 qcheck+=qp;
648                 Int_t ixs=segmentation->Ix();
649                 Int_t iys=segmentation->Iy();
650                 Float_t xs,ys;
651                 segmentation->GetPadCxy(ixs,iys,xs,ys);
652                 sum+=qp*ys;
653             }
654         } // Pad loop
655         Float_t ycog=sum/qcheck;
656         yg[i]=(yscan-y)/segmentation->Dpy();
657         yrg[i]=(ycog-y)/segmentation->Dpy();
658         ysig[i]=ycog-yscan;
659         yscan+=dy;
660     } // scan loop
661 //
662 //  x-position
663     Float_t xscan=xmin;
664     Float_t dx=segmentation->GetRealDpx(isec)/(ns-1);
665
666     for (i=0; i<ns; i++) {
667 //
668 //      Pad Loop
669 //      
670         Float_t sum=0;
671         Float_t qcheck=0;
672         segmentation->SigGenInit(xscan, y, 0);
673         
674         for (segmentation->FirstPad(xscan, y, dxI, dyI); 
675              segmentation->MorePads(); 
676              segmentation->NextPad()) 
677         {
678             qp=fResponse->IntXY(segmentation);
679             qp=TMath::Abs(qp);
680 //
681 //
682             if (qp > 1.e-2) {
683                 qcheck+=qp;
684                 Int_t ixs=segmentation->Ix();
685                 Int_t iys=segmentation->Iy();
686                 Float_t xs,ys;
687                 segmentation->GetPadCxy(ixs,iys,xs,ys);
688                 sum+=qp*xs;
689             }
690         } // Pad loop
691         Float_t xcog=sum/qcheck;
692         xcog=segmentation->GetAnod(xcog);
693         
694         xg[i]=(xscan-x)/segmentation->GetRealDpx(isec);
695         xrg[i]=(xcog-x)/segmentation->GetRealDpx(isec);
696         xsig[i]=xcog-xscan;
697         xscan+=dx;
698     }
699
700    TCanvas *c1=new TCanvas(canvasname,canvasname,400,10,600,700);
701    TPad* pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99);
702    TPad* pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99);
703    TPad* pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49);
704    TPad* pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49);
705    pad11->SetFillColor(11);
706    pad12->SetFillColor(11);
707    pad13->SetFillColor(11);
708    pad14->SetFillColor(11);
709    pad11->Draw();
710    pad12->Draw();
711    pad13->Draw();
712    pad14->Draw();
713    TGraph *graphx = new TGraph(ns,xg ,xsig);
714    TGraph *graphxr= new TGraph(ns,xrg,xsig);   
715    TGraph *graphy = new TGraph(ns,yg ,ysig);
716    TGraph *graphyr= new TGraph(ns,yrg,ysig);
717 //
718 // Creates a Root function based on function sinoid above
719 // and perform the fit
720 //
721    Double_t sinoid(Double_t *x, Double_t *par);
722    TF1 *sinoidf = new TF1("sinoidf",sinoid,0.5,0.5,5);
723    graphyr->Fit("sinoidf","V");
724    sinoidf->Copy(func);
725    func.Eval(0,0,0);
726 //
727    pad11->cd();
728    graphx->SetFillColor(42);
729    graphx->SetMarkerColor(4);
730    graphx->SetMarkerStyle(21);
731    graphx->Draw("AC");
732    graphx->GetHistogram()->SetXTitle("x on pad");
733    graphx->GetHistogram()->SetYTitle("xcog-x");   
734
735
736    pad12->cd();
737    graphxr->SetFillColor(42);
738    graphxr->SetMarkerColor(4);
739    graphxr->SetMarkerStyle(21);
740    graphxr->Draw("AP");
741    graphxr->GetHistogram()->SetXTitle("xcog on pad");
742    graphxr->GetHistogram()->SetYTitle("xcog-x");   
743
744
745    pad13->cd();
746    graphy->SetFillColor(42);
747    graphy->SetMarkerColor(4);
748    graphy->SetMarkerStyle(21);
749    graphy->Draw("AF");
750    graphy->GetHistogram()->SetXTitle("y on pad");
751    graphy->GetHistogram()->SetYTitle("ycog-y");   
752
753
754
755    pad14->cd();
756    graphyr->SetFillColor(42);
757    graphyr->SetMarkerColor(4);
758    graphyr->SetMarkerStyle(21);
759    graphyr->Draw("AF");
760    graphyr->GetHistogram()->SetXTitle("ycog on pad");
761    graphyr->GetHistogram()->SetYTitle("ycog-y");   
762
763    c1->Update();
764    
765 }
766 */
767 /*
768 Double_t sinoid(Double_t *x, Double_t *par)
769 {
770     Double_t arg = -2*TMath::Pi()*x[0];
771     Double_t fitval= par[0]*TMath::Sin(arg)+
772         par[1]*TMath::Sin(2*arg)+
773         par[2]*TMath::Sin(3*arg)+
774         par[3]*TMath::Sin(4*arg)+
775         par[4]*TMath::Sin(5*arg);
776     return fitval;
777  }
778 */
779
780
781
782
783
784
785