1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
20 #include "AliMUONClusterFinder.h"
27 #include <TPostScript.h>
30 //----------------------------------------------------------
31 static AliMUONsegmentation* gSegmentation;
32 static AliMUONresponse* gResponse;
33 static Int_t gix[500];
34 static Int_t giy[500];
35 static Float_t gCharge[500];
37 static Int_t gFirst=kTRUE;
38 static TMinuit *gMyMinuit ;
39 void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
40 static Int_t gChargeTot;
42 //----------------------------------------------------------
44 ClassImp(AliMUONClusterFinder)
46 AliMUONClusterFinder::AliMUONClusterFinder
47 (AliMUONsegmentation *segmentation, AliMUONresponse *response,
48 TClonesArray *digits, Int_t chamber)
50 fSegmentation=segmentation;
54 fNdigits = fDigits->GetEntriesFast();
56 fRawClusters=new TClonesArray("AliMUONRawCluster",10000);
65 AliMUONClusterFinder::AliMUONClusterFinder()
73 fRawClusters=new TClonesArray("AliMUONRawCluster",10000);
83 void AliMUONClusterFinder::AddRawCluster(const AliMUONRawCluster c)
86 // Add a raw cluster copy to the list
88 AliMUON *MUON=(AliMUON*)gAlice->GetModule("MUON");
89 MUON->AddRawCluster(fChamber,c);
95 void AliMUONClusterFinder::Decluster(AliMUONRawCluster *cluster)
101 Int_t mul = cluster->fMultiplicity;
102 // printf("Decluster - multiplicity %d \n",mul);
104 if (mul == 1 || mul ==2) {
105 // printf("\n Nothing special for 1- and 2-clusters \n");
107 // Nothing special for 1- and 2-clusters
109 cluster->fNcluster[0]=fNPeaks;
110 cluster->fNcluster[1]=0;
112 AddRawCluster(*cluster);
114 } else if (mul ==3) {
116 // 3-cluster, check topology
117 // printf("\n 3-cluster, check topology \n");
118 if (fDeclusterFlag) {
119 if (Centered(cluster)) {
120 // ok, cluster is centered
121 // printf("\n ok, cluster is centered \n");
123 // cluster is not centered, split into 2+1
124 // printf("\n cluster is not centered, split into 2+1 \n");
128 cluster->fNcluster[0]=fNPeaks;
129 cluster->fNcluster[1]=0;
131 AddRawCluster(*cluster);
136 // printf("Decluster - multiplicity > 45 %d \n",mul);
137 //printf("Decluster - multiplicity < 25 %d \n",mul);
139 // 4-and more-pad clusters
141 if (mul <= fClusterSize) {
142 if (fDeclusterFlag) {
143 SplitByLocalMaxima(cluster);
146 cluster->fNcluster[0]=fNPeaks;
147 cluster->fNcluster[1]=0;
149 AddRawCluster(*cluster);
158 Bool_t AliMUONClusterFinder::Centered(AliMUONRawCluster *cluster)
161 dig= (AliMUONdigit*)fDigits->UncheckedAt(cluster->fIndexMap[0]);
165 Int_t X[kMaxNeighbours], Y[kMaxNeighbours], XN[kMaxNeighbours], YN[kMaxNeighbours];
167 fSegmentation->Neighbours(ix,iy,&nn,X,Y);
169 for (Int_t i=0; i<nn; i++) {
170 if (fHitMap->TestHit(X[i],Y[i]) == used) {
178 // cluster is centered !
180 cluster->fNcluster[0]=fNPeaks;
181 cluster->fNcluster[1]=0;
183 AddRawCluster(*cluster);
188 // Highest signal on an edge, split cluster into 2+1
190 // who is the neighbour ?
191 Int_t nind=fHitMap->GetHitIndex(XN[0], YN[0]);
192 Int_t i1= (nind==cluster->fIndexMap[1]) ? 1:2;
193 Int_t i2= (nind==cluster->fIndexMap[1]) ? 2:1;
196 AliMUONRawCluster cnew;
198 cnew.fNcluster[0]=-1;
199 cnew.fNcluster[1]=fNRawClusters;
201 cnew.fNcluster[0]=fNPeaks;
204 cnew.fMultiplicity=2;
205 cnew.fIndexMap[0]=cluster->fIndexMap[0];
206 cnew.fIndexMap[1]=cluster->fIndexMap[i1];
208 cnew.fClusterType=cnew.PhysicsContribution();
213 cluster->fMultiplicity=1;
214 cluster->fIndexMap[0]=cluster->fIndexMap[i2];
215 cluster->fIndexMap[1]=0;
216 cluster->fIndexMap[2]=0;
217 FillCluster(cluster);
219 cluster->fNcluster[0]=fNPeaks;
220 cluster->fNcluster[1]=0;
222 cluster->fClusterType=cluster->PhysicsContribution();
223 AddRawCluster(*cluster);
227 printf("\n Completely screwed up %d !! \n",nd);
233 void AliMUONClusterFinder::SplitByLocalMaxima(AliMUONRawCluster *c)
235 AliMUONdigit* dig[100], *digt;
236 Int_t ix[100], iy[100], q[100];
237 Float_t x[100], y[100];
238 Int_t i; // loops over digits
239 Int_t j; // loops over local maxima
242 // Int_t threshold=500;
243 Int_t mul=c->fMultiplicity;
245 // dump digit information into arrays
247 for (i=0; i<mul; i++)
249 dig[i]= (AliMUONdigit*)fDigits->UncheckedAt(c->fIndexMap[i]);
250 ix[i]= dig[i]->fPadX;
251 iy[i]= dig[i]->fPadY;
252 q[i] = dig[i]->fSignal;
253 fSegmentation->GetPadCxy(ix[i], iy[i], x[i], y[i]);
260 Int_t AssocPeak[100];
263 Int_t X[kMaxNeighbours], Y[kMaxNeighbours];
264 for (i=0; i<mul; i++) {
265 fSegmentation->Neighbours(ix[i], iy[i], &nn, X, Y);
267 for (j=0; j<nn; j++) {
268 if (fHitMap->TestHit(X[j],Y[j])==empty) continue;
269 digt=(AliMUONdigit*) fHitMap->GetHit(X[j], Y[j]);
270 if (digt->fSignal > q[i]) {
274 // handle special case of neighbouring pads with equal signal
275 } else if (digt->fSignal == q[i]) {
277 for (Int_t k=0; k<NLocal; k++) {
278 if (X[j]==ix[IndLocal[k]] && Y[j]==iy[IndLocal[k]]){
284 } // loop over next neighbours
285 // Maxima should not be on the edge
290 } // loop over all digits
291 // printf("Found %d local Maxima",NLocal);
293 // If only one local maximum found but multiplicity is high
294 // take global maximum from the list of digits.
295 if (NLocal==1 && mul>12) {
297 for (i=0; i<mul; i++) {
308 // If number of local maxima is 2 try to fit a double gaussian
311 // Initialise global variables for fit
313 gSegmentation=fSegmentation;
314 gResponse =fResponse;
317 for (i=0; i<mul; i++) {
320 gCharge[i]=Float_t(q[i]);
325 gMyMinuit = new TMinuit(5);
327 gMyMinuit->SetFCN(fcn);
328 gMyMinuit->mninit(5,10,7);
329 Double_t arglist[20];
332 // gMyMinuit->mnexcm("SET ERR",arglist,1,ierflag);
333 // Set starting values
334 static Double_t vstart[5];
335 vstart[0]=x[IndLocal[0]];
336 vstart[1]=y[IndLocal[0]];
337 vstart[2]=x[IndLocal[1]];
338 vstart[3]=y[IndLocal[1]];
339 vstart[4]=Float_t(q[IndLocal[0]])/
340 Float_t(q[IndLocal[0]]+q[IndLocal[1]]);
341 // lower and upper limits
342 static Double_t lower[5], upper[5];
343 Int_t isec=fSegmentation->Sector(ix[IndLocal[0]], iy[IndLocal[0]]);
344 lower[0]=vstart[0]-fSegmentation->Dpx(isec)/2;
345 lower[1]=vstart[1]-fSegmentation->Dpy(isec)/2;
346 // lower[1]=vstart[1];
348 upper[0]=lower[0]+fSegmentation->Dpx(isec);
349 upper[1]=lower[1]+fSegmentation->Dpy(isec);
350 // upper[1]=vstart[1];
352 isec=fSegmentation->Sector(ix[IndLocal[1]], iy[IndLocal[1]]);
353 lower[2]=vstart[2]-fSegmentation->Dpx(isec)/2;
354 lower[3]=vstart[3]-fSegmentation->Dpy(isec)/2;
355 // lower[3]=vstart[3];
357 upper[2]=lower[2]+fSegmentation->Dpx(isec);
358 upper[3]=lower[3]+fSegmentation->Dpy(isec);
359 // upper[3]=vstart[3];
364 static Double_t step[5]={0.005, 0.03, 0.005, 0.03, 0.01};
366 gMyMinuit->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
367 gMyMinuit->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
368 gMyMinuit->mnparm(2,"x2",vstart[2],step[2],lower[2],upper[2],ierflag);
369 gMyMinuit->mnparm(3,"y2",vstart[3],step[3],lower[3],upper[3],ierflag);
370 gMyMinuit->mnparm(4,"a0",vstart[4],step[4],lower[4],upper[4],ierflag);
371 // ready for minimisation
372 gMyMinuit->SetPrintLevel(-1);
373 gMyMinuit->mnexcm("SET OUT", arglist, 0, ierflag);
377 gMyMinuit->mnexcm("SET NOGR", arglist, 0, ierflag);
378 gMyMinuit->mnexcm("SCAN", arglist, 0, ierflag);
379 gMyMinuit->mnexcm("EXIT" , arglist, 0, ierflag);
381 // Double_t amin,edm,errdef;
382 // Int_t nvpar,nparx,icstat;
383 // gMyMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
384 // gMyMinuit->mnprin(3,amin);
385 // Get fitted parameters
387 Double_t xrec[2], yrec[2], qfrac;
389 Double_t epxz, b1, b2;
391 gMyMinuit->mnpout(0, chname, xrec[0], epxz, b1, b2, ierflg);
392 gMyMinuit->mnpout(1, chname, yrec[0], epxz, b1, b2, ierflg);
393 gMyMinuit->mnpout(2, chname, xrec[1], epxz, b1, b2, ierflg);
394 gMyMinuit->mnpout(3, chname, yrec[1], epxz, b1, b2, ierflg);
395 gMyMinuit->mnpout(4, chname, qfrac, epxz, b1, b2, ierflg);
396 printf("\n %f %f %f %f %f\n", xrec[0], yrec[0], xrec[1], yrec[1],qfrac);
401 // One cluster for each maximum
403 for (j=0; j<2; j++) {
404 AliMUONRawCluster cnew;
406 cnew.fNcluster[0]=-1;
407 cnew.fNcluster[1]=fNRawClusters;
409 cnew.fNcluster[0]=fNPeaks;
412 cnew.fMultiplicity=0;
413 cnew.fX=Float_t(xrec[j]);
414 cnew.fY=Float_t(yrec[j]);
416 cnew.fQ=Int_t(gChargeTot*qfrac);
418 cnew.fQ=Int_t(gChargeTot*(1-qfrac));
420 gSegmentation->SetHit(xrec[j],yrec[j]);
421 for (i=0; i<mul; i++) {
422 cnew.fIndexMap[cnew.fMultiplicity]=c->fIndexMap[i];
423 gSegmentation->SetPad(gix[i], giy[i]);
424 Float_t q1=gResponse->IntXY(gSegmentation);
425 cnew.fContMap[cnew.fMultiplicity]=Float_t(q[i])/(q1*cnew.fQ);
426 cnew.fMultiplicity++;
428 FillCluster(&cnew,0);
429 //printf("\n x,y %f %f ", cnew.fX, cnew.fY);
430 cnew.fClusterType=cnew.PhysicsContribution();
438 if (NLocal !=2 || !fitted) {
439 // Check if enough local clusters have been found,
440 // if not add global maxima to the list
446 printf("\n Warning, no local maximum found \n");
450 if (nPerMax > fNperMax) {
451 Int_t nGlob=mul/fNperMax-NLocal+1;
454 for (i=0; i<mul; i++) {
461 if (nnew==nGlob) break;
466 // Associate hits to peaks
468 for (i=0; i<mul; i++) {
471 if (IsLocal[i]) continue;
472 for (j=0; j<NLocal; j++) {
473 Int_t il=IndLocal[j];
474 Float_t d=TMath::Sqrt((x[i]-x[il])*(x[i]-x[il])
475 +(y[i]-y[il])*(y[i]-y[il]));
478 // Select nearest peak
484 } else if (d==dmin) {
486 // If more than one take highest peak
499 // One cluster for each maximum
501 for (j=0; j<NLocal; j++) {
502 AliMUONRawCluster cnew;
504 cnew.fNcluster[0]=-1;
505 cnew.fNcluster[1]=fNRawClusters;
507 cnew.fNcluster[0]=fNPeaks;
510 cnew.fIndexMap[0]=c->fIndexMap[IndLocal[j]];
511 cnew.fMultiplicity=1;
512 for (i=0; i<mul; i++) {
513 if (IsLocal[i]) continue;
514 if (AssocPeak[i]==j) {
515 cnew.fIndexMap[cnew.fMultiplicity]=c->fIndexMap[i];
516 cnew.fMultiplicity++;
520 cnew.fClusterType=cnew.PhysicsContribution();
529 void AliMUONClusterFinder::FillCluster(AliMUONRawCluster* c, Int_t flag)
532 // Completes cluster information starting from list of digits
548 for (Int_t i=0; i<c->fMultiplicity; i++)
550 dig= (AliMUONdigit*)fDigits->UncheckedAt(c->fIndexMap[i]);
551 ix=dig->fPadX+c->fOffsetMap[i];
553 Int_t q=dig->fSignal;
554 if (dig->fPhysics >= dig->fSignal) {
556 } else if (dig->fPhysics == 0) {
558 } else c->fPhysicsMap[i]=1;
561 // peak signal and track list
563 if (q>c->fPeakSignal) {
566 c->fTracks[0]=dig->fTracks[0];
567 c->fTracks[1]=dig->fTracks[1];
568 c->fTracks[2]=dig->fTracks[2];
570 //c->fTracks[0]=dig->fTrack;
571 c->fTracks[0]=dig->fHit;
572 c->fTracks[1]=dig->fTracks[0];
573 c->fTracks[2]=dig->fTracks[1];
576 if (c->fContMap[i] > frac) {
580 c->fTracks[0]=dig->fTracks[0];
581 c->fTracks[1]=dig->fTracks[1];
582 c->fTracks[2]=dig->fTracks[2];
584 //c->fTracks[0]=dig->fTrack;
585 c->fTracks[0]=dig->fHit;
586 c->fTracks[1]=dig->fTracks[0];
587 c->fTracks[2]=dig->fTracks[1];
592 fSegmentation->GetPadCxy(ix, iy, x, y);
598 } // loop over digits
603 c->fX=fSegmentation->GetAnod(c->fX);
606 // apply correction to the coordinate along the anode wire
610 fSegmentation->GetPadIxy(x, y, ix, iy);
611 fSegmentation->GetPadCxy(ix, iy, x, y);
612 Int_t isec=fSegmentation->Sector(ix,iy);
613 TF1* CogCorr = fSegmentation->CorrFunc(isec-1);
616 Float_t YonPad=(c->fY-y)/fSegmentation->Dpy(isec);
617 c->fY=c->fY-CogCorr->Eval(YonPad, 0, 0);
623 void AliMUONClusterFinder::FindCluster(Int_t i, Int_t j, AliMUONRawCluster &c){
628 // Add i,j as element of the cluster
631 Int_t idx = fHitMap->GetHitIndex(i,j);
632 AliMUONdigit* dig = (AliMUONdigit*) fHitMap->GetHit(i,j);
633 Int_t q=dig->fSignal;
634 if (q > TMath::Abs(c.fPeakSignal)) {
637 c.fTracks[0]=dig->fTracks[0];
638 c.fTracks[1]=dig->fTracks[1];
639 c.fTracks[2]=dig->fTracks[2];
641 //c.fTracks[0]=dig->fTrack;
642 c.fTracks[0]=dig->fHit;
643 c.fTracks[1]=dig->fTracks[0];
644 c.fTracks[2]=dig->fTracks[1];
647 // Make sure that list of digits is ordered
649 Int_t mu=c.fMultiplicity;
652 if (dig->fPhysics >= dig->fSignal) {
654 } else if (dig->fPhysics == 0) {
656 } else c.fPhysicsMap[mu]=1;
659 for (Int_t ind=mu-1; ind>=0; ind--) {
660 Int_t ist=(c.fIndexMap)[ind];
661 Int_t ql=((AliMUONdigit*)fDigits
662 ->UncheckedAt(ist))->fSignal;
664 c.fIndexMap[ind]=idx;
665 c.fIndexMap[ind+1]=ist;
674 if (c.fMultiplicity >= 50 ) {
675 printf("FindCluster - multiplicity >50 %d \n",c.fMultiplicity);
679 // Prepare center of gravity calculation
681 fSegmentation->GetPadCxy(i, j, x, y);
686 fHitMap->FlagHit(i,j);
688 // Now look recursively for all neighbours
691 Int_t Xlist[kMaxNeighbours], Ylist[kMaxNeighbours];
692 fSegmentation->Neighbours(i,j,&nn,Xlist,Ylist);
693 for (Int_t in=0; in<nn; in++) {
696 if (fHitMap->TestHit(ix,iy)==unused) FindCluster(ix, iy, c);
700 //_____________________________________________________________________________
702 void AliMUONClusterFinder::FindRawClusters()
705 // simple MUON cluster finder from digits -- finds neighbours and
706 // fill the tree with raw clusters
709 if (!fNdigits) return;
711 fHitMap = new AliMUONHitMapA1(fSegmentation, fDigits);
720 for (ndig=0; ndig<fNdigits; ndig++) {
721 dig = (AliMUONdigit*)fDigits->UncheckedAt(ndig);
724 if (fHitMap->TestHit(i,j)==used ||fHitMap->TestHit(i,j)==empty) {
730 c.fPeakSignal=dig->fSignal;
732 c.fTracks[0]=dig->fTracks[0];
733 c.fTracks[1]=dig->fTracks[1];
734 c.fTracks[2]=dig->fTracks[2];
736 //c.fTracks[0]=dig->fTrack;
737 c.fTracks[0]=dig->fHit;
738 c.fTracks[1]=dig->fTracks[0];
739 c.fTracks[2]=dig->fTracks[1];
740 // tag the beginning of cluster list in a raw cluster
745 c.fX=fSegmentation->GetAnod(c.fX);
748 // apply correction to the coordinate along the anode wire
753 fSegmentation->GetPadIxy(x, y, ix, iy);
754 fSegmentation->GetPadCxy(ix, iy, x, y);
755 Int_t isec=fSegmentation->Sector(ix,iy);
756 TF1* CogCorr=fSegmentation->CorrFunc(isec-1);
758 Float_t YonPad=(c.fY-y)/fSegmentation->Dpy(isec);
759 c.fY=c.fY-CogCorr->Eval(YonPad,0,0);
763 // Analyse cluster and decluster if necessary
766 c.fNcluster[1]=fNRawClusters;
767 c.fClusterType=c.PhysicsContribution();
773 // reset Cluster object
774 for (int k=0;k<c.fMultiplicity;k++) {
782 void AliMUONClusterFinder::
790 fSegmentation->GiveTestPoints(n, x, y);
791 for (i=0; i<n; i++) {
794 SinoidalFit(xtest, ytest, func);
795 fSegmentation->SetCorrFunc(i, new TF1(func));
801 void AliMUONClusterFinder::
802 SinoidalFit(Float_t x, Float_t y, TF1 &func)
805 static Int_t count=0;
808 sprintf(canvasname,"c%d",count);
811 Float_t xg[ns], yg[ns], xrg[ns], yrg[ns];
812 Float_t xsig[ns], ysig[ns];
814 AliMUONsegmentation *segmentation=fSegmentation;
817 segmentation->GetPadIxy(x,y,ix,iy);
818 segmentation->GetPadCxy(ix,iy,x,y);
819 Int_t isec=segmentation->Sector(ix,iy);
821 Float_t xmin = x-segmentation->Dpx(isec)/2;
822 Float_t ymin = y-segmentation->Dpy(isec)/2;
824 // Integration Limits
825 Float_t dxI=fResponse->SigmaIntegration()*fResponse->ChargeSpreadX();
826 Float_t dyI=fResponse->SigmaIntegration()*fResponse->ChargeSpreadY();
836 Float_t dy=segmentation->Dpy(isec)/(ns-1);
838 for (i=0; i<ns; i++) {
844 segmentation->SigGenInit(x, yscan, 0);
846 for (segmentation->FirstPad(x, yscan, dxI, dyI);
847 segmentation->MorePads();
848 segmentation->NextPad())
850 qp=fResponse->IntXY(segmentation);
856 Int_t ixs=segmentation->Ix();
857 Int_t iys=segmentation->Iy();
859 segmentation->GetPadCxy(ixs,iys,xs,ys);
863 Float_t ycog=sum/qcheck;
864 yg[i]=(yscan-y)/segmentation->Dpy(isec);
865 yrg[i]=(ycog-y)/segmentation->Dpy(isec);
872 Float_t dx=segmentation->Dpx(isec)/(ns-1);
874 for (i=0; i<ns; i++) {
880 segmentation->SigGenInit(xscan, y, 0);
882 for (segmentation->FirstPad(xscan, y, dxI, dyI);
883 segmentation->MorePads();
884 segmentation->NextPad())
886 qp=fResponse->IntXY(segmentation);
892 Int_t ixs=segmentation->Ix();
893 Int_t iys=segmentation->Iy();
895 segmentation->GetPadCxy(ixs,iys,xs,ys);
899 Float_t xcog=sum/qcheck;
900 xcog=segmentation->GetAnod(xcog);
902 xg[i]=(xscan-x)/segmentation->Dpx(isec);
903 xrg[i]=(xcog-x)/segmentation->Dpx(isec);
908 // Creates a Root function based on function sinoid above
909 // and perform the fit
911 // TGraph *graphx = new TGraph(ns,xg ,xsig);
912 // TGraph *graphxr= new TGraph(ns,xrg,xsig);
913 // TGraph *graphy = new TGraph(ns,yg ,ysig);
914 TGraph *graphyr= new TGraph(ns,yrg,ysig);
916 Double_t sinoid(Double_t *x, Double_t *par);
917 new TF1("sinoidf",sinoid,0.5,0.5,5);
918 graphyr->Fit("sinoidf","Q");
919 func = *((TF1*)((graphyr->GetListOfFunctions())->At(0)));
922 TCanvas *c1=new TCanvas(canvasname,canvasname,400,10,600,700);
923 TPad* pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99);
924 TPad* pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99);
925 TPad* pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49);
926 TPad* pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49);
927 pad11->SetFillColor(11);
928 pad12->SetFillColor(11);
929 pad13->SetFillColor(11);
930 pad14->SetFillColor(11);
938 graphx->SetFillColor(42);
939 graphx->SetMarkerColor(4);
940 graphx->SetMarkerStyle(21);
942 graphx->GetHistogram()->SetXTitle("x on pad");
943 graphx->GetHistogram()->SetYTitle("xcog-x");
947 graphxr->SetFillColor(42);
948 graphxr->SetMarkerColor(4);
949 graphxr->SetMarkerStyle(21);
951 graphxr->GetHistogram()->SetXTitle("xcog on pad");
952 graphxr->GetHistogram()->SetYTitle("xcog-x");
956 graphy->SetFillColor(42);
957 graphy->SetMarkerColor(4);
958 graphy->SetMarkerStyle(21);
960 graphy->GetHistogram()->SetXTitle("y on pad");
961 graphy->GetHistogram()->SetYTitle("ycog-y");
966 graphyr->SetFillColor(42);
967 graphyr->SetMarkerColor(4);
968 graphyr->SetMarkerStyle(21);
970 graphyr->GetHistogram()->SetXTitle("ycog on pad");
971 graphyr->GetHistogram()->SetYTitle("ycog-y");
977 Double_t sinoid(Double_t *x, Double_t *par)
979 Double_t arg = -2*TMath::Pi()*x[0];
980 Double_t fitval= par[0]*TMath::Sin(arg)+
981 par[1]*TMath::Sin(2*arg)+
982 par[2]*TMath::Sin(3*arg)+
983 par[3]*TMath::Sin(4*arg)+
984 par[4]*TMath::Sin(5*arg);
989 Double_t DoubleGauss(Double_t *x, Double_t *par)
991 Double_t arg1 = (x[0]-par[1])/0.18;
992 Double_t arg2 = (x[0]-par[3])/0.18;
993 Double_t fitval= par[0]*TMath::Exp(-arg1*arg1/2)
994 +par[2]*TMath::Exp(-arg2*arg2/2);
998 Float_t DiscrCharge(Int_t i,Double_t *par)
1000 // par[0] x-position of first cluster
1001 // par[1] y-position of first cluster
1002 // par[2] x-position of second cluster
1003 // par[3] y-position of second cluster
1004 // par[4] charge fraction of first cluster
1005 // 1-par[4] charge fraction of second cluster
1007 static Float_t qtot;
1010 for (Int_t jbin=0; jbin<gNbins; jbin++) {
1011 qtot+=gCharge[jbin];
1014 //printf("\n sum of charge from DiscrCharge %f\n", qtot);
1015 gChargeTot=Int_t(qtot);
1018 gSegmentation->SetPad(gix[i], giy[i]);
1020 gSegmentation->SetHit(par[0],par[1]);
1021 Float_t q1=gResponse->IntXY(gSegmentation);
1024 gSegmentation->SetHit(par[2],par[3]);
1025 Float_t q2=gResponse->IntXY(gSegmentation);
1027 Float_t value = qtot*(par[4]*q1+(1.-par[4])*q2);
1032 // Minimisation function
1033 void fcn(Int_t &, Double_t *, Double_t &f, Double_t *par, Int_t )
1041 for (i=0; i<gNbins; i++) {
1042 Float_t q0=gCharge[i];
1043 Float_t q1=DiscrCharge(i,par);
1044 delta=(q0-q1)/TMath::Sqrt(q0);
1049 chisq=chisq+=(qtot-qcont)*(qtot-qcont)*0.5;