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