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