]> git.uio.no Git - u/mrichter/AliRoot.git/blame - MUON/AliMUONClusterFinderv0.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Add all TNodes to the fNodes list
[u/mrichter/AliRoot.git] / MUON / AliMUONClusterFinderv0.cxx
CommitLineData
a897a37a 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//----------------------------------------------------------
11ClassImp(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/*
22void 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
35void 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
87Int_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{
93Int_t nPara, offset, returnOffset=0 ;
94AliMUONdigit* dig= (AliMUONdigit*)fDigits->UncheckedAt(DigitIndex);
95AliMUONsegmentationV1* seg = (AliMUONsegmentationV1*) fSegmentation;
96seg->GetNParallelAndOffset(dig->fPadX,dig->fPadY,&nPara,&offset);
97if (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 }
125fSegmentation->GetPadCxy(dig->fPadX+returnOffset,dig->fPadY,*X,*Y);
126return returnOffset;
127}
128
129
130void AliMUONClusterFinderv0::SetOffset(AliMUONRawCluster *cluster)
131// compute the offsets assuming that there is only one peak !
132{
133//DumpCluster(cluster);
134Float_t X,Y;
135cluster->fOffsetMap[0]=PeakOffsetAndCoordinates(cluster->fIndexMap[0],&X,&Y);
136for (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
142void AliMUONClusterFinderv0::DumpCluster(AliMUONRawCluster *cluster)
143{
144printf ("other cluster\n");
145for (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
155Bool_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
215void 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/*
361void 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/*
415void 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
483void 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
552void AliMUONClusterFinder::
553CalibrateCOG()
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
573void AliMUONClusterFinder::
574SinoidalFit(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/*
749Double_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