]>
Commit | Line | Data |
---|---|---|
5a3482a0 | 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 | ////////////////////////////////////////////////////////////////////////// | |
17 | // // | |
18 | // AliHMPIDReconHTA // | |
19 | // // | |
20 | // HMPID class to perfom pattern recognition based on Hough transfrom // | |
21 | // for single chamber // | |
22 | ////////////////////////////////////////////////////////////////////////// | |
23 | ||
24 | #include "AliHMPIDReconHTA.h"//class header | |
25 | #include "AliHMPIDCluster.h" //CkovHiddenTrk() | |
26 | #include "AliHMPIDRecon.h" //FunMinPhot() | |
27 | #include <TMinuit.h> //FitFree() | |
28 | #include <TClonesArray.h> //CkovHiddenTrk() | |
29 | #include <AliESDtrack.h> //CkovHiddenTrk() | |
30 | #include <TH2I.h> //InitDatabase() | |
31 | #include <TGraph.h> //ShapeModel() | |
5a3482a0 | 32 | #include <TSpline.h> //ShapeModel() |
33 | #include "TStopwatch.h" // | |
34 | ||
aa00f952 | 35 | TH2I* AliHMPIDReconHTA::fgDatabase = 0x0; |
5a3482a0 | 36 | |
37 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
38 | AliHMPIDReconHTA::AliHMPIDReconHTA():TTask("RichRec","RichPat") | |
39 | { | |
40 | //.. | |
41 | //hidden algorithm | |
42 | //.. | |
43 | fMipX=fMipY=fThTrkFit=fPhTrkFit=fCkovFit=fMipQ=fRadX=fRadY=-999; | |
44 | fIdxMip=fNClu=0; | |
45 | fCkovSig2=0; | |
46 | fXClu = 0x0; | |
47 | fYClu = 0x0; | |
48 | fClCk = 0x0; | |
49 | ||
50 | fParam=AliHMPIDParam::Instance(); | |
51 | ||
52 | fParam->SetRefIdx(fParam->MeanIdxRad()); // initialization of ref index to a default one | |
aa00f952 | 53 | if(!fgDatabase) InitDatabase(); |
5a3482a0 | 54 | } |
55 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
56 | AliHMPIDReconHTA::~AliHMPIDReconHTA() | |
57 | { | |
58 | DeleteVars(); | |
59 | } | |
60 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
61 | void AliHMPIDReconHTA::InitVars(Int_t n) | |
62 | { | |
63 | //.. | |
64 | //Init some variables | |
65 | //.. | |
66 | fXClu = new Double_t[n]; | |
67 | fYClu = new Double_t[n]; | |
68 | fClCk = new Bool_t[n]; | |
69 | for(Int_t i=0;i<n;i++) fClCk[i] = kTRUE; | |
70 | // | |
71 | } | |
72 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
aa00f952 | 73 | void AliHMPIDReconHTA::DeleteVars()const |
5a3482a0 | 74 | { |
75 | //.. | |
76 | //Delete variables | |
77 | //.. | |
78 | if(fXClu) delete fXClu; | |
79 | if(fYClu) delete fYClu; | |
80 | if(fClCk) delete fClCk; | |
81 | } | |
82 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
83 | Bool_t AliHMPIDReconHTA::CkovHiddenTrk(AliESDtrack *pTrk,TClonesArray *pCluLst,Double_t nmean, Double_t qthre) | |
84 | { | |
85 | // Pattern recognition method without any infos from tracking:HTA (Hidden Track Algorithm)... | |
86 | // The method finds in the chmuber the cluster with the highest charge | |
87 | // compatibile with a MIP, then the strategy is applied | |
88 | // Arguments: pTrk - pointer to ESD track | |
89 | // pCluLs - list of clusters for a given chamber | |
90 | // nmean - mean freon ref. index | |
91 | // Returns: - 0=ok,1=not fitted | |
92 | ||
93 | AliHMPIDParam *pParam = AliHMPIDParam::Instance(); | |
94 | pParam->SetRefIdx(nmean); | |
95 | ||
96 | if(!CluPreFilter(pCluLst)) {return kFALSE;} | |
97 | ||
98 | Float_t mipX=-1,mipY=-1;Int_t mipId=-1,mipQ=-1; | |
99 | Double_t qRef = 0; | |
100 | Int_t nCh=0; | |
101 | for (Int_t iClu=0;iClu<pCluLst->GetEntriesFast();iClu++){ //clusters loop | |
102 | AliHMPIDCluster *pClu=(AliHMPIDCluster*)pCluLst->UncheckedAt(iClu); //get pointer to current cluster | |
103 | nCh = pClu->Ch(); | |
104 | fXClu[iClu] = pClu->X();fYClu[iClu] = pClu->Y(); //store x,y for fitting procedure | |
105 | fClCk[iClu] = kTRUE; //all cluster are accepted at this stage to be reconstructed | |
106 | if(pClu->Q()>qthre) fClCk[iClu] = kFALSE; // not a good photon in any case (multiple MIPs) | |
107 | if(pClu->Q()>qRef){ //searching the highest charge to select a MIP | |
108 | qRef = pClu->Q(); | |
109 | mipId=iClu; mipX=pClu->X();mipY=pClu->Y();mipQ=(Int_t)pClu->Q(); | |
110 | } | |
111 | // Printf(" n. %d x %f y %f Q %f",iClu,pClu->X(),pClu->Y(),pClu->Q()); | |
112 | }//clusters loop | |
113 | ||
114 | fNClu = pCluLst->GetEntriesFast(); | |
115 | if(qRef>qthre){ //charge compartible with MIP clusters | |
116 | fIdxMip = mipId; | |
117 | fClCk[mipId] = kFALSE; | |
118 | fMipX = mipX; fMipY=mipY; fMipQ = qRef; | |
119 | // Printf(" mipId %d x %f y %f Q %f",fIdxMip,fMipX,fMipY,fMipQ); | |
120 | if(!DoRecHiddenTrk()) { | |
121 | pTrk->SetHMPIDsignal(pParam->kNoPhotAccept); | |
122 | return kFALSE; | |
123 | } //Do track and ring reconstruction,if problems returns 1 | |
124 | pTrk->SetHMPIDtrk(fRadX,fRadY,fThTrkFit,fPhTrkFit); //store track intersection info | |
125 | pTrk->SetHMPIDmip(fMipX,fMipY,(Int_t)fMipQ,fNClu); //store mip info | |
126 | pTrk->SetHMPIDcluIdx(nCh,fIdxMip); //set cham number and index of cluster | |
127 | pTrk->SetHMPIDsignal(fCkovFit); //find best Theta ckov for ring i.e. track | |
128 | pTrk->SetHMPIDchi2(fCkovSig2); //errors squared | |
129 | // Printf(" n clusters tot %i accepted %i",pCluLst->GetEntriesFast(),fNClu); | |
130 | // Printf("CkovHiddenTrk: thetaC %f th %f ph %f",fCkovFit,fThTrkFit,fPhTrkFit); | |
131 | return kTRUE; | |
132 | } | |
133 | ||
134 | return kFALSE; | |
135 | }//CkovHiddenTrk() | |
136 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
137 | Bool_t AliHMPIDReconHTA::DoRecHiddenTrk() | |
138 | { | |
139 | // Pattern recognition method without any infos from tracking... | |
140 | // First a preclustering filter to avoid part of the noise | |
141 | // Then only ellipsed-rings are fitted (no possibility, | |
142 | // for the moment, to reconstruct very inclined tracks) | |
143 | // Finally a fitting with (th,ph) free, starting by very close values | |
144 | // previously evaluated. | |
145 | // Arguments: none | |
146 | // Returns: none | |
147 | Double_t thTrkRec,phiTrkRec,thetaCRec; | |
148 | ||
149 | if(!FindShape(thTrkRec,phiTrkRec,thetaCRec)) return kFALSE; | |
150 | ||
151 | if(!FitFree(thTrkRec,phiTrkRec)) {return kFALSE;} | |
152 | ||
153 | return kTRUE; | |
154 | }//DoRecHiddenTrk() | |
155 | /* | |
156 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
157 | Bool_t AliHMPIDReconHTA::DoRecHiddenTrk(TClonesArray *pCluLst) | |
158 | { | |
159 | // Pattern recognition method without any infos from tracking... | |
160 | // First a preclustering filter to avoid part of the noise | |
161 | // Then only ellipsed-rings are fitted (no possibility, | |
162 | // for the moment, to reconstruct very inclined tracks) | |
163 | // Finally a fitting with (th,ph) free, starting by very close values | |
164 | // previously evaluated. | |
165 | // Arguments: none | |
166 | // Returns: none | |
167 | Double_t thTrkRec,phiTrkRec,thetaCRec; | |
168 | ||
169 | if(!FindShape(thTrkRec,phiTrkRec,thetaCRec)) return kFALSE; | |
170 | ||
171 | // Printf("thTrkRec %f phiTrkRec %f ThetaCRec %f",thTrkRec*TMath::RadToDeg(),phiTrkRec*TMath::RadToDeg(),thetaCRec*TMath::RadToDeg()); | |
172 | Int_t nClTmp1 = pCluLst->GetEntriesFast()-1; //minus MIP... | |
173 | Int_t nClTmp2 = 0; | |
174 | ||
175 | while(nClTmp1 != nClTmp2){ | |
176 | SetNClu(pCluLst->GetEntriesFast()); | |
177 | if(!FitFree(thTrkRec,phiTrkRec)) {return kFALSE;} | |
178 | nClTmp2 = NClu(); | |
179 | if(nClTmp2!=nClTmp1) {nClTmp1=nClTmp2;nClTmp2=0;} | |
180 | } | |
181 | ||
182 | fNClu = nClTmp2; | |
183 | return kTRUE; | |
184 | }//DoRecHiddenTrk() | |
185 | */ | |
186 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
187 | Bool_t AliHMPIDReconHTA::CluPreFilter(TClonesArray *pCluLst) | |
188 | { | |
aa00f952 | 189 | // Pre-filter of bkg clusters |
190 | // Arguments: pSluLst - List of the clusters for a given chamber | |
191 | // Returns: status - TRUE if filtering leaves enough photons, FALSE if not | |
5a3482a0 | 192 | // |
193 | Int_t nClusTot = pCluLst->GetEntriesFast(); | |
194 | if(nClusTot<4||nClusTot>100) { | |
195 | return kFALSE; | |
196 | } else { | |
197 | InitVars(nClusTot); | |
198 | return kTRUE; | |
199 | } | |
200 | } | |
201 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
202 | Bool_t AliHMPIDReconHTA::FindShape(Double_t &thTrkRec,Double_t &phiTrkRec,Double_t &thetaCRec) | |
203 | { | |
aa00f952 | 204 | // Finds the estimates for phi and theta of the track and the ThetaCerenkov |
205 | // by using a database of the shapes of the rings | |
206 | // Arguments: none | |
207 | // Returns: thTrkRec - estimate of theta track | |
208 | // phiTrkRec - estimate of phi track | |
209 | // thetaCRec - estimate of ThetaCerenkov | |
210 | // status - TRUE if a good solution is found, FALSE if not | |
211 | ||
5a3482a0 | 212 | Double_t *phiphot = new Double_t[fNClu]; |
213 | Double_t *dist = new Double_t[fNClu]; | |
214 | Int_t *indphi = new Int_t[fNClu]; | |
215 | ||
216 | Bool_t status; | |
217 | ||
218 | // Sort in phi angle... | |
219 | for(Int_t i=0;i<fNClu;i++) { | |
220 | phiphot[i] = (TMath::ATan2(fMipY-fYClu[i],fMipX-fXClu[i])+TMath::Pi())*TMath::RadToDeg(); | |
221 | dist[i]=TMath::Sqrt((fMipX-fXClu[i])*(fMipX-fXClu[i])+(fMipY-fYClu[i])*(fMipY-fYClu[i])); | |
222 | } | |
223 | ||
224 | TMath::Sort(fNClu,phiphot,indphi,kFALSE); | |
225 | ||
226 | // Purify with a truncated mean; | |
227 | Int_t np=0; | |
228 | Double_t dMean = 0; | |
229 | Double_t dMean2 = 0; | |
230 | for(Int_t i=0;i<fNClu;i++) { | |
231 | if(!fClCk[indphi[i]]) continue; // Check if a good photon candidate or not | |
232 | dMean +=dist[indphi[i]]; | |
233 | dMean2+=dist[indphi[i]]*dist[indphi[i]]; | |
234 | np++; | |
235 | } | |
236 | ||
237 | dMean /=(Double_t)np; | |
238 | dMean2 /=(Double_t)np; | |
239 | Double_t rms = TMath::Sqrt(dMean2 - dMean*dMean); | |
240 | ||
241 | for(Int_t i=0;i<fNClu;i++) { | |
242 | if(!fClCk[indphi[i]]) continue; // Check if a good photon candidate or not | |
243 | if(TMath::Abs(dMean-dist[indphi[i]]) > 2*rms) { | |
244 | fClCk[indphi[i]] = kFALSE; | |
245 | continue; | |
246 | } | |
247 | } | |
248 | ||
249 | // | |
250 | // purify vectors for good photon candidates | |
251 | // | |
252 | Int_t npeff=0; | |
253 | Double_t *phiphotP = new Double_t[fNClu+1]; | |
254 | Double_t *distP = new Double_t[fNClu+1]; | |
255 | for(Int_t i=0;i<fNClu;i++) { | |
256 | if(!fClCk[indphi[i]]) continue; // Check if a good photon candidate or not | |
257 | phiphotP[npeff] = phiphot[indphi[i]]; | |
258 | distP[npeff] = dist[indphi[i]]; | |
259 | npeff++; | |
260 | } | |
261 | ||
262 | delete [] phiphot; | |
263 | delete [] dist; | |
264 | delete [] indphi; | |
265 | ||
266 | if(npeff<3) { | |
267 | delete [] phiphotP; | |
268 | delete [] distP; | |
269 | return kFALSE; | |
270 | } | |
271 | ||
272 | // for(Int_t i=0;i<npeff;i++) {Printf(" n. %d phiphot %f dist %f",i,phiphotP[i],distP[i]);} | |
273 | ||
274 | Double_t xA,xB; | |
275 | if(ShapeModel(npeff,phiphotP,distP,xA,xB,phiTrkRec)) { | |
276 | ||
277 | // Printf("FindShape: phi start %f",phiTrkRec*TMath::RadToDeg()); | |
278 | ||
aa00f952 | 279 | Int_t bin = fgDatabase->FindBin(xA,xB); |
280 | Int_t compact = (Int_t)fgDatabase->GetBinContent(bin); | |
5a3482a0 | 281 | thetaCRec = (Double_t)(compact%1000); |
282 | thTrkRec = (Double_t)(compact/1000); | |
283 | ||
284 | thTrkRec *= TMath::DegToRad(); | |
285 | thetaCRec *= TMath::DegToRad(); | |
286 | ||
287 | // Printf("FindShape: xA %f xB %f compact %d thTrk %f thC %f",xA,xB,compact,thTrkRec*TMath::RadToDeg(),thetaCRec*TMath::RadToDeg()); | |
288 | ||
289 | status = kTRUE; | |
290 | ||
291 | } else { | |
292 | ||
293 | status = kFALSE; | |
294 | ||
295 | } | |
296 | ||
297 | delete [] phiphotP; | |
298 | delete [] distP; | |
299 | ||
300 | return status; | |
301 | } | |
302 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
303 | Bool_t AliHMPIDReconHTA::ShapeModel(Int_t np,Double_t *phiphot,Double_t *dist,Double_t &xA,Double_t &xB,Double_t &phiStart) | |
304 | { | |
aa00f952 | 305 | // Find a Spline curve to define dist. vs. phi angle |
306 | // in order to estimate the phi of the track | |
307 | // Arguments: np - # points corresponding to # photon candidates | |
308 | // dist - distance of each photon from MIP | |
309 | // phiphot - phi of the photon in the DRS | |
310 | // Returns: xA - min. distance from MIP | |
311 | // xB - dist. from mip perpedicular to the major axis | |
312 | // phiStart- estimate of the track phi | |
313 | ||
5a3482a0 | 314 | TGraph *phigr = new TGraph(np,phiphot,dist); |
315 | TSpline3 *sphi = new TSpline3("sphi",phigr); | |
316 | if(!sphi) {Printf("Spline not created!Bye.");return kFALSE;} | |
317 | ||
318 | Int_t locMin = TMath::LocMin(np,dist); | |
319 | Int_t locMax = TMath::LocMax(np,dist); | |
320 | ||
321 | Double_t minX = phiphot[locMin]; | |
322 | // Double_t minY = dist[locMin]; | |
323 | Double_t maxX = phiphot[locMax]; | |
324 | // Double_t maxY = dist[locMax]; | |
325 | ||
326 | Int_t ip[3] = {-1,0,1}; | |
327 | if(locMin==0 ) {ip[0]= 0;ip[1]= 1;ip[2]= 2;} | |
328 | if(locMin==np-1) {ip[0]=-2;ip[1]=-1;ip[2]= 0;} | |
329 | ||
330 | Double_t minXf = VertParab(phiphot[locMin+ip[0]],dist[locMin+ip[0]], | |
331 | phiphot[locMin+ip[1]],dist[locMin+ip[1]], | |
332 | phiphot[locMin+ip[2]],dist[locMin+ip[2]]); | |
333 | if(minXf< phiphot[locMin+ip[0]] || minXf > phiphot[locMin+ip[2]]) minXf = minX; | |
334 | ||
335 | ip[0]=-1;ip[1]=0;ip[2]=1; | |
336 | if(locMax==0 ) {ip[0]= 0;ip[1]= 1;ip[2]= 2;} | |
337 | if(locMax==np-1) {ip[0]=-2;ip[1]=-1;ip[2]= 0;} | |
338 | ||
339 | Double_t maxXf = VertParab(phiphot[locMax+ip[0]],dist[locMax+ip[0]], | |
340 | phiphot[locMax+ip[1]],dist[locMax+ip[1]], | |
341 | phiphot[locMax+ip[2]],dist[locMax+ip[2]]); | |
342 | if(maxXf< phiphot[locMax+ip[0]] || maxXf > phiphot[locMax+ip[2]]) maxXf = maxX; | |
343 | ||
344 | // Printf(" phi at mindist %f and found %f",minX,minXf); | |
345 | // Printf(" phi at maxdist %f and found %f",maxX,maxXf); | |
346 | // | |
347 | if(TMath::Abs(maxXf-minXf)>30) { | |
348 | xA = sphi->Eval(minXf); | |
349 | if(xA < 0) xA = dist[sphi->FindX(xA)]; | |
350 | xB = sphi->Eval(minXf-90); | |
351 | if(xB < 0) xB = dist[sphi->FindX(xB)]; | |
352 | phiStart = minXf-180; //open ring or acceptance effect...so believe to min phi angle! | |
353 | } else { | |
354 | phiStart = 0.5*(maxXf-180+minXf); | |
355 | xA = sphi->Eval(phiStart+180); | |
356 | if(xA < 0) xA = dist[sphi->FindX(xA)]; | |
357 | xB = sphi->Eval(phiStart+90); | |
358 | if(xB < 0) xB = dist[sphi->FindX(xB)]; | |
359 | } | |
360 | // | |
361 | // Printf("ShapeModel: phiStart %f xA %f xB %f",phiStart,xA,xB); | |
362 | ||
363 | phiStart*=TMath::DegToRad(); | |
364 | return kTRUE; | |
365 | } | |
366 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
aa00f952 | 367 | Double_t AliHMPIDReconHTA::VertParab(Double_t x1,Double_t y1,Double_t x2, Double_t y2, Double_t x3, Double_t y3)const |
5a3482a0 | 368 | { |
aa00f952 | 369 | // It uses parabola from 3 points to evaluate the x-coord of the parab |
370 | // Arguments: xi,yi - points | |
371 | // Returns: x-coord of the vertex | |
372 | ||
5a3482a0 | 373 | Double_t a = ((x1-x3)*(y1-y2)-(x1-x2)*(y1-y3))/((x1*x1-x2*x2)*(x1-x3)-(x1*x1-x3*x3)*(x1-x2)); |
374 | Double_t b = (y1-y2 - a*(x1*x1-x2*x2))/(x1-x2); | |
375 | // Double_t c = y1 - a*x1*x1-b*x1; | |
376 | return -0.5*b/a; | |
377 | } | |
378 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
379 | Bool_t AliHMPIDReconHTA::FitFree(Double_t thTrkRec,Double_t phiTrkRec) | |
380 | { | |
381 | // Fit performed by minimizing RMS/sqrt(n) of the | |
382 | // photons reconstructed. First phi is fixed and theta | |
383 | // is fouond, then (th,ph) of the track | |
384 | // as free parameters | |
385 | // Arguments: PhiRec phi of the track | |
386 | // Returns: none | |
387 | ||
5a3482a0 | 388 | TMinuit *pMinuit = new TMinuit(2); |
389 | pMinuit->mncler(); | |
390 | gMinuit->SetObjectFit((TObject*)this); gMinuit->SetFCN(AliHMPIDReconHTA::FunMinPhot); //set fit function | |
391 | Double_t aArg=-1,parStep,parLow,parHigh; Int_t iErrFlg; //tmp vars for TMinuit | |
392 | Double_t d1,d2,d3; | |
393 | TString sName; | |
394 | Double_t th,ph; | |
395 | ||
396 | pMinuit->mnexcm("SET PRI",&aArg,1,iErrFlg); //suspend all printout from TMinuit | |
397 | pMinuit->mnexcm("SET NOW",&aArg,0,iErrFlg); | |
398 | ||
399 | if(thTrkRec==0) thTrkRec = 3.*TMath::DegToRad(); // not to start from the edge... | |
400 | ||
401 | pMinuit->mnparm(0," thTrk ",thTrkRec ,parStep=0.01,parLow=0,parHigh=TMath::PiOver4(),iErrFlg); | |
402 | pMinuit->mnparm(1," phiTrk ",phiTrkRec,parStep=0.01,parLow=0,parHigh=TMath::TwoPi(),iErrFlg); | |
403 | ||
404 | pMinuit->FixParameter(1); | |
405 | pMinuit->mnexcm("SIMPLEX" ,&aArg,0,iErrFlg); | |
406 | pMinuit->mnexcm("MIGRAD" ,&aArg,0,iErrFlg); | |
407 | pMinuit->Release(1); | |
408 | pMinuit->mnexcm("MIGRAD" ,&aArg,0,iErrFlg); | |
409 | ||
410 | pMinuit->mnpout(0,sName,th,d1,d2,d3,iErrFlg); | |
411 | pMinuit->mnpout(1,sName,ph,d1,d2,d3,iErrFlg); | |
412 | ||
413 | Double_t f,par[2]; | |
414 | Double_t *grad=0x0; | |
415 | par[0] = th;par[1] = ph; | |
416 | pMinuit->Eval(2,grad,f,par,3); | |
417 | ||
418 | // Printf("FitFree: theta %f phi %f",th,ph); | |
419 | ||
420 | SetTrkFit(th,ph); | |
421 | return kTRUE; | |
422 | } | |
423 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
424 | void AliHMPIDReconHTA::FunMinPhot(Int_t &/* */,Double_t* /* */,Double_t &f,Double_t *par,Int_t iflag) | |
425 | { | |
426 | // Minimization function to find best track and thetaC parameters | |
427 | // Arguments: f = function value to minimize | |
428 | // par = list of parameter to find | |
429 | // iflag = flag status. See Minuit instructions | |
430 | // Returns: none | |
431 | // | |
432 | // Note: it is necessary to call an instance of AlihMPIDParam. Not possible to use fParam | |
433 | // because of the static instantiation of the function in Minuit | |
434 | ||
435 | AliHMPIDParam *pParam=AliHMPIDParam::Instance(); | |
436 | AliHMPIDReconHTA *pRecHTA=(AliHMPIDReconHTA*)gMinuit->GetObjectFit(); | |
437 | AliHMPIDRecon pRec; | |
438 | Double_t sizeCh = 0.5*pParam->RadThick()+pParam->WinThick()+pParam->GapThick(); | |
439 | Double_t thTrk = par[0]; | |
440 | Double_t phTrk = par[1]; | |
441 | Double_t xrad = pRecHTA->MipX() - sizeCh*TMath::Tan(thTrk)*TMath::Cos(phTrk); | |
442 | Double_t yrad = pRecHTA->MipY() - sizeCh*TMath::Tan(thTrk)*TMath::Sin(phTrk); | |
443 | pRecHTA->SetRadXY(xrad,yrad); | |
444 | pRec.SetTrack(xrad,yrad,thTrk,phTrk); | |
445 | ||
446 | Double_t meanCkov =0; | |
447 | Double_t meanCkov2=0; | |
448 | Double_t thetaCer,phiCer; | |
449 | Int_t nClAcc = 0; | |
450 | Int_t nClTot=pRecHTA->NClu(); | |
451 | ||
452 | for(Int_t i=0;i<nClTot;i++) { | |
453 | if(!(pRecHTA->ClCk(i))) continue; | |
454 | pRec.FindPhotCkov(pRecHTA->XClu(i),pRecHTA->YClu(i),thetaCer,phiCer); | |
455 | meanCkov += thetaCer; | |
456 | meanCkov2 += thetaCer*thetaCer; | |
457 | nClAcc++; | |
458 | } | |
459 | if(nClAcc==0) {f=999;return;} | |
460 | meanCkov /=(Double_t)nClAcc; | |
461 | meanCkov2 /=(Double_t)nClAcc; | |
462 | Double_t rms = TMath::Sqrt(TMath::Abs(meanCkov2 - meanCkov*meanCkov)); | |
463 | f = rms/TMath::Sqrt((Double_t)nClAcc); | |
464 | ||
465 | if(iflag==3) { | |
466 | /* | |
467 | Printf("FunMinPhot before: photons candidates %i used %i",nClTot,nClAcc); | |
468 | nClAcc = 0; | |
469 | Double_t meanCkov1=0; | |
470 | Double_t meanCkov2=0; | |
471 | for(Int_t i=0;i<nClTot;i++) { | |
472 | if(!(pRec->ClCk(i))) continue; | |
473 | pRec->FindPhotCkov(pRec->XClu(i),pRec->YClu(i),thetaCer,phiCer); | |
474 | if(TMath::Abs(thetaCer-meanCkov)<2*rms) { | |
475 | meanCkov1 += thetaCer; | |
476 | meanCkov2 += thetaCer*thetaCer; | |
477 | nClAcc++; | |
478 | } else pRec->SetClCk(i,kFALSE); | |
479 | } | |
480 | meanCkov1/=nClAcc; | |
481 | Double_t rms2 = (meanCkov2 - meanCkov*meanCkov*nClAcc)/nClAcc; | |
482 | Printf("FunMinPhot after: photons candidates %i used %i thetaCer %f",nClTot,nClAcc,meanCkov1); | |
483 | pRec->SetCkovFit(meanCkov1); | |
484 | pRec->SetCkovSig2(rms2); | |
485 | pRec->SetNClu(nClAcc); | |
486 | */ | |
487 | // Printf("FunMinPhot: reconstructed theta Cerenkov %f with %d photons",meanCkov,nClAcc); | |
488 | pRecHTA->SetCkovFit(meanCkov); | |
489 | pRecHTA->SetCkovSig2(rms*rms); | |
490 | pRecHTA->SetNClu(nClAcc); | |
491 | } | |
492 | }//FunMinPhot() | |
493 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
494 | void AliHMPIDReconHTA::InitDatabase() | |
495 | { | |
aa00f952 | 496 | // Construction a database of ring shapes on fly |
497 | // Arguments: none | |
498 | // Returns : none | |
499 | // N.B. fgDatabase points to a TH2I with x-min dist from MIP | |
500 | // y-dist from the ring of the MIP perpendicular to major axis | |
501 | // The content is the packed info of track theta and thetaC in degrees | |
502 | // thetaC+1000*thTrk | |
503 | // | |
5a3482a0 | 504 | TStopwatch timer; |
505 | timer.Start(); | |
506 | ||
507 | AliInfo(Form("database HTA is being built.Please, wait...")); | |
508 | // | |
509 | Double_t x[3],y[3]; | |
510 | ||
511 | AliHMPIDRecon rec; | |
512 | ||
513 | if(!fParam) fParam=AliHMPIDParam::Instance(); | |
514 | Double_t thetaMax = TMath::ACos(1./fParam->MeanIdxRad()); | |
515 | Double_t thTrkMax = 1./TMath::ASin(fParam->MeanIdxRad()); | |
516 | ||
517 | Int_t nstepx = 1000; | |
518 | Int_t nstepy = 1000; | |
519 | ||
520 | TH2I *deconv = new TH2I("deconv","database;d1;d2;thC+1000*thTrk",500,0,50,150,0,15); | |
521 | // | |
522 | Double_t xrad = 0; | |
523 | Double_t yrad = 0; | |
524 | Double_t phTrk = 0; | |
525 | ||
526 | for(Int_t i=0;i<nstepx;i++) { //loop on thetaC | |
527 | for(Int_t j=0;j<nstepy;j++) { //loop on theta particle | |
528 | Double_t thetaC = thetaMax/nstepx*((Double_t)i+0.5); | |
529 | Double_t thTrk = thTrkMax/nstepy*((Double_t)j+0.5); | |
530 | // | |
531 | //mip position | |
532 | // | |
533 | Double_t sizeCh = fParam->RadThick()+fParam->WinThick()+fParam->GapThick(); | |
534 | Double_t xmip = xrad + sizeCh*TMath::Tan(thTrk)*TMath::Cos(phTrk); | |
535 | Double_t ymip = yrad + sizeCh*TMath::Tan(thTrk)*TMath::Sin(phTrk); | |
536 | ||
537 | Double_t dist1,dist2; | |
538 | // | |
539 | //first point at phi=0 | |
540 | // | |
541 | rec.SetTrack(xrad,yrad,thTrk,phTrk); | |
542 | TVector2 pos; | |
543 | pos=rec.TracePhot(thetaC,0); | |
544 | ||
545 | if(pos.X()==-999) { | |
546 | dist1 = 0; //open ring...anly the distance btw mip and point at 180 will be considered | |
547 | } else { | |
548 | x[0] = pos.X(); y[0] = pos.Y(); | |
549 | dist1 = TMath::Sqrt((x[0]-xmip)*(x[0]-xmip)+(y[0]-ymip)*(y[0]-ymip)); | |
550 | } | |
551 | // | |
552 | //second point at phi=180 | |
553 | // | |
554 | rec.SetTrack(xrad,yrad,thTrk,phTrk); | |
555 | pos=rec.TracePhot(thetaC,TMath::Pi()); | |
556 | ||
557 | if(pos.X()==-999) {Printf("it should not happens!Bye");return;} | |
558 | x[1] = pos.X(); y[1] = pos.Y(); | |
559 | if((x[1]-xmip)*(x[0]-xmip)>0) continue; // to avoid circles out mips (for very low ThetaC) | |
560 | dist2 = TMath::Sqrt((x[1]-xmip)*(x[1]-xmip)+(y[1]-ymip)*(y[1]-ymip)); | |
561 | ||
562 | // Double_t distA = dist1+dist2; | |
563 | Double_t distA = dist2; // only the minimum: problem of acceptance | |
564 | // | |
565 | //second point at phi=90 | |
566 | // | |
567 | rec.SetTrack(xrad,yrad,thTrk,phTrk); | |
568 | pos=rec.TracePhot(thetaC,TMath::PiOver2()); | |
569 | ||
570 | if(pos.X()==-999) continue; | |
571 | x[2] = pos.X(); y[2] = pos.Y(); | |
572 | Double_t distB = TMath::Sqrt((x[2]-xmip)*(x[2]-xmip)+(y[2]-ymip)*(y[2]-ymip)); | |
573 | // compact the infos... | |
574 | Int_t compact = (Int_t)(thetaC*TMath::RadToDeg())+1000*(Int_t)(thTrk*TMath::RadToDeg()); | |
575 | Int_t bin = deconv->FindBin(distA,distB); | |
576 | if(deconv->GetBinContent(bin)==0) deconv->Fill(distA,distB,compact); | |
577 | } | |
578 | } | |
579 | ||
580 | FillZeroChan(deconv); | |
aa00f952 | 581 | fgDatabase = deconv; |
5a3482a0 | 582 | |
583 | timer.Stop(); | |
584 | Double_t nSecs = timer.CpuTime(); | |
585 | AliInfo(Form("database HTA successfully open in %3.1f sec.(CPU). Reconstruction is started.",nSecs)); | |
586 | } | |
587 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
aa00f952 | 588 | void AliHMPIDReconHTA::FillZeroChan(TH2I *deconv)const |
5a3482a0 | 589 | { |
590 | Int_t nbinx = deconv->GetNbinsX(); | |
591 | Int_t nbiny = deconv->GetNbinsY(); | |
592 | for(Int_t i = 0;i<nbinx;i++) { | |
593 | for(Int_t j = 0;j<nbiny;j++) { | |
594 | if(deconv->GetBinContent(i,j) == 0) { | |
595 | Int_t nXmin = i-1; Int_t nXmax=i+1; | |
596 | Int_t nYmin = j-1; Int_t nYmax=j+1; | |
597 | Int_t nc = 0; | |
598 | Double_t meanC =0; | |
599 | Double_t meanTrk =0; | |
600 | for(Int_t ix=nXmin;ix<=nXmax;ix++) { | |
601 | if(ix<0||ix>nbinx) continue; | |
602 | for(Int_t iy=nYmin;iy<=nYmax;iy++) { | |
603 | if(iy<0||iy>nbiny) continue; | |
604 | meanC += (Int_t)deconv->GetBinContent(ix,iy)%1000; | |
605 | meanTrk+= (Int_t)deconv->GetBinContent(ix,iy)/1000; | |
606 | nc++; | |
607 | } | |
608 | meanC/=nc; meanTrk/=nc; | |
609 | Int_t compact = (Int_t)meanC+1000*(Int_t)meanTrk; | |
610 | if(compact>0)deconv->SetCellContent(i,j,compact); | |
611 | } | |
612 | } | |
613 | } | |
614 | } | |
615 | } | |
616 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |