]>
Commit | Line | Data |
---|---|---|
0fe8fa07 | 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 "AliRICHCluster.h" | |
0422a446 | 17 | #include <TMinuit.h> //Solve() |
0fe8fa07 | 18 | |
19 | ClassImp(AliRICHCluster) | |
20 | //__________________________________________________________________________________________________ | |
21 | void AliRICHCluster::Print(Option_t*)const | |
22 | { | |
23 | //Print current cluster | |
24 | const char *status=0; | |
25 | switch(fStatus){ | |
0422a446 | 26 | case kFormed: status="formed" ;break; |
27 | case kUnfolded: status="unfolded" ;break; | |
28 | case kCoG: status="CoGed" ;break; | |
29 | case kEmpty: status="empty" ;break; | |
0fe8fa07 | 30 | } |
0422a446 | 31 | Int_t iNdigs=0; if(fDigits) iNdigs=fDigits->GetEntriesFast(); |
0fe8fa07 | 32 | |
0422a446 | 33 | Printf("cfm=%10i, cs=%2i, Size=%2i Maxima=%2i, Shape=%5i, pos=(%7.3f,%7.3f) Q=%6i, %s", |
34 | fCFM,fChamber,Size(),Nlocmax(),fShape,fX,fY,fQdc,status,iNdigs); | |
35 | for(Int_t i=0;i<iNdigs;i++) Digit(i)->Print(); | |
0fe8fa07 | 36 | }//Print() |
0422a446 | 37 | |
38 | ||
39 | //__________________________________________________________________________________________________ | |
40 | TMinuit* AliRICHCluster::Solve() | |
41 | { | |
42 | //At this point, cluster contains a list of digits, cluster charge is precalculated as a sum of digits charges (in AddDigit()), | |
43 | //position is preset to (-1,-1) (in ctor), status is preset to kFormed in (AddDigit()), chamber-sector info is preseted to actual value (in AddDigit()) | |
44 | //Here we decide what to do with this cluster: unfold or just calculate center of gravity | |
45 | //Arguments: none | |
46 | // Returns: pointer to fitter or 0 if no unfolding decided | |
47 | TMinuit *pMinuit=0; | |
48 | if(Size()>=2 && AliRICHParam::IsResolveClusters()) | |
49 | pMinuit=Unfold(); | |
50 | else | |
51 | CoG(0); | |
52 | return pMinuit; | |
53 | }//Solve() | |
54 | //__________________________________________________________________________________________________ | |
55 | void AliRICHCluster::FitFunc(Int_t &iNpars, Double_t *, Double_t &chi2, Double_t *aPar, Int_t ) | |
56 | { | |
57 | //Cluster fit function | |
58 | //par[0]=x par[1]=y par[2]=q for the first Mathieson shape | |
59 | //par[3]=x par[4]=y par[5]=q for the second Mathieson shape and so on up to iNpars/3 Mathieson shapes | |
60 | //We need to calculate Qpad - Qpadmath summup over all pads of the cluster | |
61 | //Here Qpad is a actual charge of the pad, Qpadmath is calculated charge of the pad induced by all Mathiesons | |
62 | //Arguments: iNpars - number of parameters which is number of local maxima of cluster * 3 | |
63 | // chi2 - function result to be minimised | |
64 | // aPar - parametrs array of size iNpars | |
65 | // Returns: none | |
66 | AliRICHCluster *pClu=(AliRICHCluster*)gMinuit->GetObjectFit(); | |
67 | Int_t iNmathiesons = iNpars/3; | |
68 | ||
69 | TVector2 curMathiesonPos; | |
70 | chi2 = 0; | |
71 | for(Int_t i=0;i<pClu->Size();i++){//digits loop | |
72 | TVector pad = pClu->Digit(i)->Pad(); | |
73 | Double_t dQpad = pClu->Digit(i)->Qdc(); | |
74 | Double_t dQpadmath = 0; | |
75 | for(Int_t j=0;j<iNmathiesons;j++){//all Mathiesons may contribute to this pad | |
76 | curMathiesonPos.Set(aPar[3*j],aPar[3*j+1]);//get current position for current Mathieson | |
77 | dQpadmath += aPar[3*j+2]*AliRICHParam::FracQdc(curMathiesonPos,pad);//sums up contributions to the current pad from all Mathiesons | |
78 | } | |
79 | chi2 += TMath::Power((dQpadmath-dQpad),2); | |
80 | }//digits loop | |
81 | }//RichClusterFitFunction() | |
82 | //__________________________________________________________________________________________________ | |
83 | TMinuit* AliRICHCluster::Unfold() | |
84 | { | |
85 | //This methode is invoked from Solve() when decided to unfold this cluster | |
86 | //Method first finds number of local maxima and if it's more then one tries to unfold this cluster into local maxima number of clusters | |
87 | //Arguments: none | |
88 | // Returns: pointer to fitter for retriving parameters | |
89 | ||
90 | TMinuit *pMinuit = new TMinuit(15); //init MINUIT with max 15 parameters (maxim 5 mathiesons, 3 params per matheson ) | |
91 | pMinuit->SetObjectFit((TObject*)this); | |
92 | pMinuit->SetFCN(AliRICHCluster::FitFunc);//set fit function | |
93 | Double_t aArg=-1,parStart,parStep,parLow,parHigh; Int_t iErrFlg; //tmp for MINUIT parameters definitions | |
94 | pMinuit->mnexcm("SET PRI" ,&aArg,1,iErrFlg); //suspend all printout from TMinuit | |
95 | ||
96 | Int_t iLocMaxCnt=0; | |
97 | //Strategy is to check if the current pad has QDC more then all neigbours | |
98 | for(Int_t iDig1=0;iDig1<Size();iDig1++) {//first digits loop | |
99 | AliRICHDigit *pDig1 = Digit(iDig1);//take the current digit | |
100 | Int_t iHowManyMoreCnt = 0;//counts how many neighbouring pads has QDC more then current one | |
101 | for(Int_t iDig2=0;iDig2<Size();iDig2++) {//loop on all digits again | |
102 | AliRICHDigit *pDig2 = Digit(iDig2); | |
103 | if(iDig1==iDig2) continue; //no need to compare | |
104 | Int_t dist = TMath::Sign(Int_t(pDig1->PadX()-pDig2->PadX()),1)+TMath::Sign(Int_t(pDig1->PadY()-pDig2->PadY()),1);//distance between pads | |
105 | if(dist==1)//means pads are neighbours | |
106 | if(pDig2->Qdc()>=pDig1->Qdc()) iHowManyMoreCnt++;//count number of pads with Q more then Q of current pad | |
107 | }//second digits loop | |
108 | if(iHowManyMoreCnt==0&&iLocMaxCnt<6){//this pad has Q more then any neighbour so it's local maximum | |
109 | TVector2 x2=AliRICHParam::Pad2Loc(pDig1->Pad());//take pad center position and use it as parameter for current Mathienson shape | |
110 | pMinuit->mnparm(3*iLocMaxCnt ,Form("x%i",iLocMaxCnt),parStart=x2.X() ,parStep=0.01,parLow=0,parHigh=0,iErrFlg); | |
111 | pMinuit->mnparm(3*iLocMaxCnt+1,Form("y%i",iLocMaxCnt),parStart=x2.Y() ,parStep=0.01,parLow=0,parHigh=0,iErrFlg); | |
112 | pMinuit->mnparm(3*iLocMaxCnt+2,Form("q%i",iLocMaxCnt),parStart=pDig1->Qdc(),parStep=0.01,parLow=0,parHigh=0,iErrFlg);// | |
113 | iLocMaxCnt++; | |
114 | }//if this pad is local maximum | |
115 | }//first digits loop | |
116 | ||
117 | fSize+=iLocMaxCnt; | |
118 | if(iLocMaxCnt>0&&iLocMaxCnt<6){ //resonable number of local maxima to fit | |
119 | Double_t aArg=0; | |
120 | pMinuit->mnexcm("MIGRAD",&aArg,0,iErrFlg);//start fitting | |
121 | fStatus=kUnfolded; | |
122 | }else{ | |
123 | delete pMinuit; | |
124 | pMinuit=0; | |
125 | CoG(0); | |
126 | } | |
127 | return pMinuit; | |
128 | }//Unfold() | |
0fe8fa07 | 129 | //__________________________________________________________________________________________________ |
0422a446 | 130 | void AliRICHCluster::CoG(Int_t nLocals) |
131 | { | |
132 | //Calculates naive cluster position as a center of gravity of its digits. | |
133 | //Also determines the box fully contaning this cluster | |
134 | //Arguments: | |
135 | Float_t xmin=999,ymin=999,xmax=0,ymax=0; | |
136 | fX=fY=0; | |
137 | for(Int_t iDig=0;iDig<Size();iDig++) { | |
138 | AliRICHDigit *pDig=Digit(iDig); | |
139 | TVector pad=pDig->Pad(); Double_t q=pDig->Qdc(); | |
140 | TVector2 x2=AliRICHParam::Pad2Loc(pad); | |
141 | fX += x2.X()*q;fY +=x2.Y()*q; | |
142 | if(pad[0]<xmin)xmin=pad[0];if(pad[0]>xmax)xmax=pad[0];if(pad[1]<ymin)ymin=pad[1];if(pad[1]>ymax)ymax=pad[1]; | |
143 | } | |
144 | fX/=fQdc;fY/=fQdc;//Center of Gravity | |
145 | ||
146 | TVector2 center = AliRICHParam::Pad2Loc(AliRICHParam::Loc2Pad(TVector2(fX,fY))); | |
147 | fX += AliRICHParam::CogCorr(fX-center.X());//correct cluster position for sinoid | |
148 | ||
149 | fShape=Int_t(100*(xmax-xmin+1)+ymax-ymin+1);//find box containing cluster | |
150 | fSize+=nLocals; | |
151 | fStatus=kCoG; | |
152 | }//CoG() | |
153 | //__________________________________________________________________________________________________ | |
154 | void AliRICHCluster::Test(const TVector2 &hitX2,Double_t dEloss) | |
155 | { | |
156 | //This is to test all cluster functionality | |
157 | //Method uses AddDigit() to add a predifined pad structure and then calls Solve | |
158 | Int_t iQtot=AliRICHParam::TotQdc(hitX2,dEloss); | |
159 | if(iQtot==0){ | |
160 | Printf("Provided hit position out of sensitive area"); | |
161 | return; | |
162 | } | |
163 | TVector area=AliRICHParam::Loc2Area(hitX2); | |
164 | TVector pad(2); | |
165 | for(pad[1]=area[1];pad[1]<=area[3];pad[1]++){//affected pads loop first y | |
166 | for(pad[0]=area[0];pad[0]<=area[2];pad[0]++){//then x | |
167 | Double_t dQpad=iQtot*AliRICHParam::FracQdc(hitX2,pad);//charge fraction from Mathieson centered at x to pad | |
168 | AddDigit(new AliRICHDigit(3,(Int_t)pad[0],(Int_t)pad[1],dQpad)); | |
169 | }//affected pads loop | |
170 | } | |
171 | TMinuit *pMinuit=Solve(); | |
172 | Print(); | |
173 | Printf("Initial hit (%.2f,%.2f) Qtot=%i Eloss=%.2f",hitX2.X(),hitX2.Y(),iQtot,dEloss); | |
174 | Double_t d1,d2,d3; Int_t iErrFlg;TString sName; //tmp vars for TMinuit | |
175 | Double_t x,y,q; | |
176 | for(Int_t i=0;i<Nlocmax();i++){//retrive fitting results | |
177 | pMinuit->mnpout(3*i ,sName, x, d1 , d2, d3, iErrFlg); | |
178 | pMinuit->mnpout(3*i+1 ,sName, y, d1 , d2, d3, iErrFlg); | |
179 | pMinuit->mnpout(3*i+2 ,sName, q, d1 , d2, d3, iErrFlg); | |
180 | } | |
181 | Printf(" Fitted hit (%.2f,%.2f) Qfit=%.0f",x,y,q); | |
182 | delete pMinuit; pMinuit=0; Reset(); | |
183 | }//Test() |