8265fa96 |
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 | |
8265fa96 |
16 | |
17 | #include "AliRICHClusterFinder.h" |
543d5224 |
18 | #include "AliRICHMap.h" |
9d6f9427 |
19 | #include <TMinuit.h> |
ed3ceb24 |
20 | #include <TParticle.h> |
9d6f9427 |
21 | #include <TVector3.h> |
cbaf35fb |
22 | #include <AliLoader.h> |
ed3ceb24 |
23 | #include <AliStack.h> |
cbaf35fb |
24 | #include <AliRun.h> |
8265fa96 |
25 | |
9d6f9427 |
26 | void RICHMinMathieson(Int_t &npar, Double_t *gin, Double_t &chi2, Double_t *par, Int_t iflag); |
8265fa96 |
27 | |
28 | ClassImp(AliRICHClusterFinder) |
543d5224 |
29 | //__________________________________________________________________________________________________ |
30 | AliRICHClusterFinder::AliRICHClusterFinder(AliRICH *pRICH) |
31 | {//main ctor |
e42a7b46 |
32 | fRICH = pRICH;//this goes first as GetDebug depends on it. |
33 | if(GetDebug())Info("main ctor","Start."); |
543d5224 |
34 | |
e42a7b46 |
35 | fDigitMap = 0; |
6865fb8d |
36 | fRawCluster.Reset(); |
37 | fResolvedCluster.Reset(); |
e42a7b46 |
38 | if(GetDebug())Info("main ctor","Stop."); |
543d5224 |
39 | }//main ctor |
40 | //__________________________________________________________________________________________________ |
543d5224 |
41 | void AliRICHClusterFinder::Exec() |
237c933d |
42 | { |
e42a7b46 |
43 | //Main method of cluster finder. Loops on events and chambers, everything else is done in FindClusters() |
89924db8 |
44 | if(GetDebug()) Info("Exec","Start."); |
e42a7b46 |
45 | |
46 | R()->GetLoader() ->LoadDigits(); |
47 | R()->GetLoader()->GetRunLoader()->LoadHeader(); |
48 | R()->GetLoader()->GetRunLoader()->LoadKinematics(); |
ed3ceb24 |
49 | |
543d5224 |
50 | for(Int_t iEventN=0;iEventN<gAlice->GetEventsPerRun();iEventN++){//events loop |
e42a7b46 |
51 | if(GetDebug()) Info("Exec","Processing event %i...",iEventN); |
52 | R()->GetLoader()->GetRunLoader()->GetEvent(iEventN); |
237c933d |
53 | |
e42a7b46 |
54 | R()->GetLoader()->MakeTree("R"); R()->MakeBranch("R"); |
55 | R()->ResetDigits(); R()->ResetClusters(); |
543d5224 |
56 | |
e42a7b46 |
57 | R()->GetLoader()->TreeD()->GetEntry(0); |
cbaf35fb |
58 | for(Int_t iChamber=1;iChamber<=kNCH;iChamber++){//chambers loop |
ed3ceb24 |
59 | FindClusters(iChamber); |
543d5224 |
60 | }//chambers loop |
e42a7b46 |
61 | R()->GetLoader()->TreeR()->Fill(); R()->GetLoader()->WriteRecPoints("OVERWRITE");//write out clusters for current event |
543d5224 |
62 | }//events loop |
ed3ceb24 |
63 | |
e42a7b46 |
64 | R()->ResetDigits();//reset and unload everything |
65 | R()->ResetClusters(); |
66 | R()->GetLoader() ->UnloadDigits(); |
67 | R()->GetLoader() ->UnloadRecPoints(); |
68 | R()->GetLoader()->GetRunLoader()->UnloadHeader(); |
69 | R()->GetLoader()->GetRunLoader()->UnloadKinematics(); |
ed3ceb24 |
70 | |
89924db8 |
71 | if(GetDebug()) Info("Exec","Stop."); |
543d5224 |
72 | }//Exec() |
73 | //__________________________________________________________________________________________________ |
ed3ceb24 |
74 | void AliRICHClusterFinder::FindClusters(Int_t iChamber) |
75 | { |
e42a7b46 |
76 | //Loops on digits for a given chamber, forms raw clusters, then tries to resolve them if requested |
77 | Int_t iNdigits=R()->Digits(iChamber)->GetEntriesFast(); |
78 | if(GetDebug())Info("FindClusters","Start for chamber %i with %i digits.",iChamber,iNdigits); |
79 | |
80 | if(iNdigits==0)return;//no digits for a given chamber, nothing to do |
cbaf35fb |
81 | |
e42a7b46 |
82 | fDigitMap=new AliRICHMap(R()->Digits(iChamber));//create digit map for the given chamber |
cbaf35fb |
83 | |
e42a7b46 |
84 | for(Int_t iDigN=0;iDigN<iNdigits;iDigN++){//digits loop for a given chamber |
85 | AliRICHdigit *dig=(AliRICHdigit*)R()->Digits(iChamber)->At(iDigN); |
cbaf35fb |
86 | Int_t i=dig->X(); Int_t j=dig->Y(); |
e42a7b46 |
87 | if(fDigitMap->TestHit(i,j)==kUsed) continue;//this digit is already taken, go after next digit |
cbaf35fb |
88 | |
e42a7b46 |
89 | FormRawCluster(i,j);//form raw cluster starting from (i,j) pad |
90 | if(GetDebug()) {Info("FindClusters","After FormRawCluster:");fRawCluster.Print();} |
91 | FindLocalMaxima(); //find number of local maxima and initial center of gravity |
92 | if(GetDebug()) {Info("FindClusters","After FindLocalMaxima:");fRawCluster.Print();} |
6865fb8d |
93 | |
e42a7b46 |
94 | if(AliRICHParam::IsResolveClusters()&&fRawCluster.Size()>1&&fRawCluster.Size()<6){ |
95 | FitCoG(); //serialization of resolved clusters will happen inside |
96 | }else{//cluster size=1 or resolving is switched off |
97 | WriteRawCluster();//simply output the formed raw cluster without deconvolution |
ed3ceb24 |
98 | } |
e42a7b46 |
99 | fRawCluster.Reset(); fResolvedCluster.Reset(); |
100 | }//digits loop for a given chamber |
cbaf35fb |
101 | |
e42a7b46 |
102 | delete fDigitMap; |
cbaf35fb |
103 | |
e42a7b46 |
104 | if(GetDebug())Info("FindClusters","Stop."); |
ed3ceb24 |
105 | }//FindClusters() |
106 | //__________________________________________________________________________________________________ |
6865fb8d |
107 | void AliRICHClusterFinder::FindClusterContribs(AliRICHcluster *pCluster) |
ed3ceb24 |
108 | { |
e42a7b46 |
109 | //Finds cerenkov-feedback-mip mixture for a given cluster |
110 | if(GetDebug()) {Info("FindClusterContribs","Start.");pCluster->Print();} |
ed3ceb24 |
111 | |
6865fb8d |
112 | TObjArray *pDigits = pCluster->Digits(); |
e42a7b46 |
113 | if(!pDigits) return; //?????????? |
ed3ceb24 |
114 | Int_t iNmips=0,iNckovs=0,iNfeeds=0; |
6865fb8d |
115 | TArrayI contribs(3*pCluster->Size()); |
116 | Int_t *pindex = new Int_t[3*pCluster->Size()]; |
117 | for(Int_t iDigN=0;iDigN<pCluster->Size();iDigN++) {//loop on digits of a given cluster |
e42a7b46 |
118 | contribs[3*iDigN] =((AliRICHdigit*)pDigits->At(iDigN))->GetTrack(0); |
89924db8 |
119 | if (contribs[3*iDigN] >= 10000000) contribs[3*iDigN] = 0; |
e42a7b46 |
120 | contribs[3*iDigN+1]=((AliRICHdigit*)pDigits->At(iDigN))->GetTrack(1); |
89924db8 |
121 | if (contribs[3*iDigN+1] >= 10000000) contribs[3*iDigN+1] = 0; |
e42a7b46 |
122 | contribs[3*iDigN+2]=((AliRICHdigit*)pDigits->At(iDigN))->GetTrack(2); |
89924db8 |
123 | if (contribs[3*iDigN+2] >= 10000000) contribs[3*iDigN+2] = 0; |
ed3ceb24 |
124 | }//loop on digits of a given cluster |
125 | TMath::Sort(contribs.GetSize(),contribs.GetArray(),pindex); |
e42a7b46 |
126 | for(Int_t iDigN=0;iDigN<3*pCluster->Size()-1;iDigN++) {//loop on digits to sort tids |
127 | if(GetDebug()) Info("FindClusterContribs","%4i for digit n. %4i",contribs[pindex[iDigN]],iDigN); |
ed3ceb24 |
128 | if(contribs[pindex[iDigN]]!=contribs[pindex[iDigN+1]]) { |
e42a7b46 |
129 | TParticle* particle = R()->GetLoader()->GetRunLoader()->Stack()->Particle(contribs[pindex[iDigN]]); |
89924db8 |
130 | Int_t code = particle->GetPdgCode(); |
131 | Double_t charge = 0; |
e42a7b46 |
132 | if(particle->GetPDG()) charge=particle->GetPDG()->Charge(); |
133 | if(GetDebug()) Info("FindClusterContribs"," charge of particle %f",charge); |
ed3ceb24 |
134 | |
135 | if(code==50000050) iNckovs++; |
e42a7b46 |
136 | if(code==50000051) iNfeeds++; |
137 | if(charge!=0) iNmips++; |
ed3ceb24 |
138 | } |
139 | }//loop on digits to sort Tid |
e42a7b46 |
140 | |
141 | if (contribs[pindex[3*pCluster->Size()-1]]!=kBad) { |
142 | |
143 | TParticle* particle = R()->GetLoader()->GetRunLoader()->Stack()->Particle(contribs[pindex[3*pCluster->Size()-1]]); |
144 | Int_t code = particle->GetPdgCode(); |
145 | Double_t charge = 0; |
146 | if(particle->GetPDG()) charge=particle->GetPDG()->Charge(); |
147 | if(GetDebug()) Info("FindClusterContribs"," charge of particle %f",charge); |
148 | if(code==50000050) iNckovs++; |
149 | if(code==50000051) iNfeeds++; |
150 | if(charge!=0) iNmips++; |
151 | } |
152 | |
153 | pCluster->CFM(iNckovs,iNfeeds,iNmips); |
154 | // |
ed3ceb24 |
155 | delete [] pindex; |
e42a7b46 |
156 | if(GetDebug()) {pCluster->Print();Info("FindClusterContribs","Stop.");} |
3582c1f9 |
157 | }//FindClusterContribs() |
cbaf35fb |
158 | //__________________________________________________________________________________________________ |
9d6f9427 |
159 | void AliRICHClusterFinder::FormRawCluster(Int_t i, Int_t j) |
ed3ceb24 |
160 | { |
e42a7b46 |
161 | //Builds the raw cluster (before deconvolution). Starts from the first pad (i,j) then calls itself recursevly for all neighbours. |
162 | if(GetDebug()) Info("FormRawCluster","Start with digit(%i,%i) Q=%f",i,j,((AliRICHdigit*)fDigitMap->GetHit(i,j))->Q()); |
cbaf35fb |
163 | |
e42a7b46 |
164 | fRawCluster.AddDigit((AliRICHdigit*) fDigitMap->GetHit(i,j));//take this pad in cluster |
165 | fDigitMap->FlagHit(i,j);//flag this pad as taken |
cbaf35fb |
166 | |
167 | Int_t listX[4], listY[4]; // Now look recursively for all neighbours |
e42a7b46 |
168 | for (Int_t iNei=0;iNei<R()->P()->PadNeighbours(i,j,listX,listY);iNei++) |
169 | if(fDigitMap->TestHit(listX[iNei],listY[iNei])==kUnused) FormRawCluster(listX[iNei],listY[iNei]); |
3582c1f9 |
170 | }//FormRawCluster() |
cbaf35fb |
171 | //__________________________________________________________________________________________________ |
e42a7b46 |
172 | void AliRICHClusterFinder::FindLocalMaxima() |
173 | { |
174 | //find number of local maxima in the current raw cluster and then calculates initial center of gravity |
175 | fNlocals=0; |
176 | for(Int_t iDig1=0;iDig1<fRawCluster.Size();iDig1++) { |
177 | Int_t iNotMax = 0; |
178 | AliRICHdigit *pDig1 = (AliRICHdigit *)fRawCluster.Digits()->At(iDig1); |
179 | TVector pad1 = pDig1->Pad(); |
180 | Int_t padQ1 = (Int_t)(pDig1->Q()+0.1); |
181 | Int_t padC1 = pDig1->ChFbMi(); |
182 | for(Int_t iDig2=0;iDig2<fRawCluster.Size();iDig2++) { |
183 | AliRICHdigit *pDig2 = (AliRICHdigit *)fRawCluster.Digits()->At(iDig2); |
184 | TVector pad2 = pDig2->Pad(); |
185 | Int_t padQ2 = (Int_t)(pDig2->Q()+0.1); |
186 | if(iDig1==iDig2) continue; |
187 | Int_t diffx = TMath::Sign(Int_t(pad1[0]-pad2[0]),1); |
188 | Int_t diffy = TMath::Sign(Int_t(pad1[1]-pad2[1]),1); |
189 | if((diffx+diffy)<=1) { |
190 | if(padQ2>=padQ1) iNotMax++; |
191 | } |
192 | } |
193 | if(iNotMax==0) { |
194 | TVector2 x2=AliRICHParam::Pad2Loc(pad1); |
195 | fLocalX[fNlocals]=x2.X();fLocalY[fNlocals]=x2.Y(); |
196 | fLocalQ[fNlocals] = (Double_t)padQ1; |
197 | fLocalC[fNlocals] = padC1; |
198 | fNlocals++; |
199 | } |
200 | } |
201 | fRawCluster.CoG(fNlocals); //first initial approximation of the CoG...to start minimization. |
202 | }//FindLocalMaxima() |
cbaf35fb |
203 | //__________________________________________________________________________________________________ |
9d6f9427 |
204 | void AliRICHClusterFinder::WriteRawCluster() |
6865fb8d |
205 | { |
e42a7b46 |
206 | //Add the current raw cluster to the list of clusters |
89924db8 |
207 | if(GetDebug()) Info("WriteRawCluster","Start."); |
cbaf35fb |
208 | |
e42a7b46 |
209 | FindClusterContribs(&fRawCluster); |
210 | R()->AddCluster(fRawCluster); |
211 | |
212 | if(GetDebug()){fRawCluster.Print(); Info("WriteRawCluster","Stop.");} |
cbaf35fb |
213 | }//WriteRawCluster() |
214 | //__________________________________________________________________________________________________ |
6865fb8d |
215 | void AliRICHClusterFinder::WriteResolvedCluster() |
216 | { |
e42a7b46 |
217 | //Add the current resolved cluster to the list of clusters |
89924db8 |
218 | if(GetDebug()) Info("WriteResolvedCluster","Start."); |
6865fb8d |
219 | |
e42a7b46 |
220 | FindClusterContribs(&fResolvedCluster); |
221 | R()->AddCluster(fResolvedCluster); |
6865fb8d |
222 | |
e42a7b46 |
223 | if(GetDebug()){fResolvedCluster.Print(); Info("WriteResolvedCluster","Stop.");} |
6865fb8d |
224 | }//WriteResolvedCluster() |
225 | //__________________________________________________________________________________________________ |
9d6f9427 |
226 | void AliRICHClusterFinder::FitCoG() |
e42a7b46 |
227 | { |
228 | //Fits cluster of size by the corresponding number of Mathieson shapes. |
229 | //This methode is only invoked in case everything is ok to start deconvolution |
230 | if(GetDebug()) {Info("FitCoG","Start with:"); fRawCluster.Print();} |
6865fb8d |
231 | |
0f72f306 |
232 | Double_t arglist; |
233 | Int_t ierflag = 0; |
6865fb8d |
234 | |
6865fb8d |
235 | TMinuit *pMinuit = new TMinuit(3*fNlocals-1); |
0f72f306 |
236 | pMinuit->mninit(5,10,7); |
237 | |
238 | arglist = -1; |
239 | pMinuit->mnexcm("SET PRI",&arglist, 1, ierflag); |
19dd5b2f |
240 | pMinuit->mnexcm("SET NOW",&arglist, 0, ierflag); |
9d6f9427 |
241 | |
9d6f9427 |
242 | TString chname; |
243 | Int_t ierflg; |
244 | |
245 | pMinuit->SetObjectFit((TObject*)this); |
246 | pMinuit->SetFCN(RICHMinMathieson); |
6865fb8d |
247 | |
248 | Double_t vstart,lower, upper; |
249 | Double_t stepX= 0.001; |
250 | Double_t stepY= 0.001; |
251 | Double_t stepQ= 0.0001; |
252 | |
253 | for(Int_t i=0;i<fNlocals;i++) { |
254 | vstart = fLocalX[i]; |
255 | lower = vstart - 2*AliRICHParam::PadSizeX(); |
256 | upper = vstart + 2*AliRICHParam::PadSizeX(); |
257 | pMinuit->mnparm(3*i ,Form("xCoG %i",i),vstart,stepX,lower,upper,ierflag); |
6865fb8d |
258 | vstart = fLocalY[i]; |
259 | lower = vstart - 2*AliRICHParam::PadSizeY(); |
260 | upper = vstart + 2*AliRICHParam::PadSizeY(); |
261 | pMinuit->mnparm(3*i+1,Form("yCoG %i",i),vstart,stepY,lower,upper,ierflag); |
6865fb8d |
262 | if(i==fNlocals-1) break; // last parameter is constrained |
263 | vstart = fLocalQ[i]/fRawCluster.Q(); |
264 | lower = 0; |
265 | upper = 1; |
266 | pMinuit->mnparm(3*i+2,Form("qfrac %i",i),vstart,stepQ,lower,upper,ierflag); |
6865fb8d |
267 | } |
268 | |
9d6f9427 |
269 | arglist = -1; |
9d6f9427 |
270 | pMinuit->mnexcm("SET NOGR",&arglist, 1, ierflag); |
9d6f9427 |
271 | arglist = 1; |
272 | pMinuit->mnexcm("SET ERR", &arglist, 1,ierflg); |
273 | arglist = -1; |
274 | pMinuit->mnexcm("SIMPLEX",&arglist, 0, ierflag); |
275 | pMinuit->mnexcm("MIGRAD",&arglist, 0, ierflag); |
6865fb8d |
276 | pMinuit->mnexcm("EXIT" ,&arglist, 0, ierflag); |
277 | |
278 | Double_t xCoG[50],yCoG[50],qfracCoG[50]; |
9d6f9427 |
279 | Double_t eps, b1, b2; |
6865fb8d |
280 | |
281 | Double_t qfraclast=0; |
282 | for(Int_t i=0;i<fNlocals;i++) { |
283 | pMinuit->mnpout(3*i ,chname, xCoG[i], eps , b1, b2, ierflg); |
284 | pMinuit->mnpout(3*i+1,chname, yCoG[i], eps , b1, b2, ierflg); |
285 | if(i==fNlocals-1) break; |
286 | pMinuit->mnpout(3*i+2,chname, qfracCoG[i], eps , b1, b2, ierflg); |
287 | qfraclast+=qfracCoG[i]; |
288 | } |
289 | qfracCoG[fNlocals-1] = 1 - qfraclast; |
9d6f9427 |
290 | delete pMinuit; |
6865fb8d |
291 | |
e42a7b46 |
292 | for(Int_t i=0;i<fNlocals;i++){//resolved positions loop |
6865fb8d |
293 | fResolvedCluster.Fill(&fRawCluster,xCoG[i],yCoG[i],qfracCoG[i],fLocalC[i]); |
6865fb8d |
294 | WriteResolvedCluster(); |
295 | } |
e42a7b46 |
296 | if(GetDebug()) Info("CoG","Stop."); |
09c52ebc |
297 | }//FitCoG() |
9d6f9427 |
298 | //__________________________________________________________________________________________________ |
6865fb8d |
299 | void RICHMinMathieson(Int_t &npar, Double_t *, Double_t &chi2, Double_t *par, Int_t ) |
09c52ebc |
300 | { |
e42a7b46 |
301 | //Mathieson minimization function |
9d6f9427 |
302 | |
6865fb8d |
303 | AliRICHcluster *pRawCluster = ((AliRICHClusterFinder*)gMinuit->GetObjectFit())->GetRawCluster(); |
9d6f9427 |
304 | |
3582c1f9 |
305 | TVector2 centroid[50]; |
6865fb8d |
306 | Double_t q[50]; |
307 | Int_t nFunctions = (npar+1)/3; |
308 | Double_t qfract = 0; |
309 | for(Int_t i=0;i<nFunctions;i++) { |
3582c1f9 |
310 | centroid[i].Set(par[3*i],par[3*i+1]); |
6865fb8d |
311 | if(i==nFunctions-1) break; |
312 | q[i]=par[3*i+2]; |
313 | qfract+=q[i]; |
314 | } |
315 | q[nFunctions-1] = 1 - qfract; |
316 | |
9d6f9427 |
317 | chi2 = 0; |
318 | Int_t qtot = pRawCluster->Q(); |
319 | for(Int_t i=0;i<pRawCluster->Size();i++) { |
e42a7b46 |
320 | TVector pad=((AliRICHdigit *)pRawCluster->Digits()->At(i))->Pad(); |
9d6f9427 |
321 | Double_t padQ = ((AliRICHdigit *)pRawCluster->Digits()->At(i))->Q(); |
6865fb8d |
322 | Double_t qfracpar=0; |
323 | for(Int_t j=0;j<nFunctions;j++) { |
e42a7b46 |
324 | qfracpar += q[j]*AliRICHParam::FracQdc(centroid[j],pad); |
9d6f9427 |
325 | } |
6865fb8d |
326 | chi2 += TMath::Power((qtot*qfracpar-padQ),2)/padQ; |
6865fb8d |
327 | } |
9d6f9427 |
328 | }//RICHMinMathieson() |
e42a7b46 |
329 | //__________________________________________________________________________________________________ |