]>
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() | |
9785d5fb | 27 | #include <TFile.h> //Database() |
5a3482a0 | 28 | #include <TMinuit.h> //FitFree() |
29 | #include <TClonesArray.h> //CkovHiddenTrk() | |
30 | #include <AliESDtrack.h> //CkovHiddenTrk() | |
9785d5fb | 31 | #include <TH2F.h> //InitDatabase() |
5a3482a0 | 32 | #include <TGraph.h> //ShapeModel() |
5a3482a0 | 33 | #include <TSpline.h> //ShapeModel() |
2e81abb0 | 34 | #include <TCanvas.h> //ShapeModel() |
5a3482a0 | 35 | #include "TStopwatch.h" // |
36 | ||
45118213 | 37 | Int_t AliHMPIDReconHTA::fgDB[500][150]={{75000*0}}; |
5a3482a0 | 38 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
c61a7285 | 39 | AliHMPIDReconHTA::AliHMPIDReconHTA(): |
40 | TTask("RichRec","RichPat"), | |
41 | fMipX(-999), | |
42 | fMipY(-999), | |
43 | fMipQ(-999), | |
44 | fRadX(-999), | |
45 | fRadY(-999), | |
46 | fIdxMip(0), | |
47 | fNClu(0), | |
48 | fXClu(0), | |
49 | fYClu(0), | |
2e81abb0 | 50 | fPhiPhot(0), |
51 | fThetaPhot(0), | |
c61a7285 | 52 | fClCk(0), |
ee2f3539 | 53 | fThTrkIn(-999), |
54 | fPhTrkIn(-999), | |
c61a7285 | 55 | fThTrkFit(-999), |
56 | fPhTrkFit(-999), | |
57 | fCkovFit(-999), | |
2e81abb0 | 58 | fNCluFit(0), |
c61a7285 | 59 | fCkovSig2(0), |
2e81abb0 | 60 | fFitStatus(0), |
c61a7285 | 61 | fParam(AliHMPIDParam::Instance()) |
5a3482a0 | 62 | { |
63 | //.. | |
64 | //hidden algorithm | |
65 | //.. | |
5a3482a0 | 66 | fParam->SetRefIdx(fParam->MeanIdxRad()); // initialization of ref index to a default one |
5565f017 | 67 | InitDatabase(); |
5a3482a0 | 68 | } |
69 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
70 | AliHMPIDReconHTA::~AliHMPIDReconHTA() | |
71 | { | |
72 | DeleteVars(); | |
73 | } | |
74 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
75 | void AliHMPIDReconHTA::InitVars(Int_t n) | |
76 | { | |
77 | //.. | |
78 | //Init some variables | |
79 | //.. | |
80 | fXClu = new Double_t[n]; | |
81 | fYClu = new Double_t[n]; | |
2e81abb0 | 82 | fPhiPhot = new Double_t[n]; |
83 | fThetaPhot = new Double_t[n]; | |
5a3482a0 | 84 | fClCk = new Bool_t[n]; |
85 | for(Int_t i=0;i<n;i++) fClCk[i] = kTRUE; | |
86 | // | |
87 | } | |
88 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
aa00f952 | 89 | void AliHMPIDReconHTA::DeleteVars()const |
5a3482a0 | 90 | { |
91 | //.. | |
92 | //Delete variables | |
93 | //.. | |
6fefa5ce | 94 | if(fXClu) delete [] fXClu; |
95 | if(fYClu) delete [] fYClu; | |
96 | if(fPhiPhot) delete [] fPhiPhot; | |
97 | if(fThetaPhot) delete [] fThetaPhot; | |
98 | if(fClCk) delete [] fClCk; | |
5a3482a0 | 99 | } |
100 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
9785d5fb | 101 | Bool_t AliHMPIDReconHTA::CkovHiddenTrk(AliESDtrack *pTrk,TClonesArray *pCluLst,Int_t index, Double_t nmean) |
5a3482a0 | 102 | { |
103 | // Pattern recognition method without any infos from tracking:HTA (Hidden Track Algorithm)... | |
104 | // The method finds in the chmuber the cluster with the highest charge | |
105 | // compatibile with a MIP, then the strategy is applied | |
106 | // Arguments: pTrk - pointer to ESD track | |
107 | // pCluLs - list of clusters for a given chamber | |
9785d5fb | 108 | // pNmean - pointer to ref. index |
109 | // pQthre - pointer to qthre | |
5a3482a0 | 110 | // Returns: - 0=ok,1=not fitted |
111 | ||
112 | AliHMPIDParam *pParam = AliHMPIDParam::Instance(); | |
5a3482a0 | 113 | |
9785d5fb | 114 | if(!CluPreFilter(pCluLst)) return kFALSE; |
115 | ||
5a3482a0 | 116 | Int_t nCh=0; |
43d3333b | 117 | Int_t sizeClu=0; |
9785d5fb | 118 | |
119 | fNClu = pCluLst->GetEntriesFast(); | |
120 | ||
121 | for (Int_t iClu=0;iClu<fNClu;iClu++){ //clusters loop | |
5a3482a0 | 122 | AliHMPIDCluster *pClu=(AliHMPIDCluster*)pCluLst->UncheckedAt(iClu); //get pointer to current cluster |
5a3482a0 | 123 | fXClu[iClu] = pClu->X();fYClu[iClu] = pClu->Y(); //store x,y for fitting procedure |
124 | fClCk[iClu] = kTRUE; //all cluster are accepted at this stage to be reconstructed | |
9785d5fb | 125 | |
126 | if(iClu == index) { | |
127 | ||
128 | fMipX = pClu->X(); | |
129 | fMipY = pClu->Y(); | |
130 | fMipQ = pClu->Q(); | |
131 | sizeClu = pClu->Size(); | |
132 | nCh = pClu->Ch(); | |
133 | fClCk[index] = kFALSE; | |
134 | fIdxMip = index; | |
2e81abb0 | 135 | AliDebug(1,Form(" MIP n. %i x %f y %f Q %f",iClu,pClu->X(),pClu->Y(),pClu->Q())); |
9785d5fb | 136 | } |
5a3482a0 | 137 | }//clusters loop |
9785d5fb | 138 | |
139 | pParam->SetRefIdx(nmean); | |
140 | ||
e56b695f | 141 | // |
142 | Float_t xra,yra,th,ph; pTrk->GetHMPIDtrk(xra,yra,th,ph); | |
2e81abb0 | 143 | AliDebug(1,Form(" simulated phiTRK %6.2f thetaTRK %6.2f",ph*TMath::RadToDeg(),th*TMath::RadToDeg())); |
e56b695f | 144 | // |
145 | ||
9785d5fb | 146 | if(!DoRecHiddenTrk()) { |
147 | pTrk->SetHMPIDsignal(pParam->kNoPhotAccept); | |
148 | return kFALSE; | |
149 | } //Do track and ring reconstruction,if problems returns 1 | |
2e81abb0 | 150 | AliDebug(1,Form(" fitted phi %6.2f ",fPhTrkFit*TMath::RadToDeg())); |
9785d5fb | 151 | |
152 | pTrk->SetHMPIDtrk(fRadX,fRadY,fThTrkFit,fPhTrkFit); //store track intersection info | |
2e81abb0 | 153 | pTrk->SetHMPIDmip(fMipX,fMipY,(Int_t)fMipQ,NCluFit()); //store mip info + n. phots of the ring |
9785d5fb | 154 | pTrk->SetHMPIDcluIdx(nCh,fIdxMip+1000*sizeClu); //set cham number, index of cluster + cluster size |
155 | pTrk->SetHMPIDsignal(fCkovFit); //find best Theta ckov for ring i.e. track | |
156 | pTrk->SetHMPIDchi2(fCkovSig2); //errors squared | |
2e81abb0 | 157 | AliDebug(1,Form(" n clusters tot %i fitted to ring %i",fNClu,NCluFit())); |
158 | for(Int_t i=0;i<fNClu;i++) { | |
159 | AliDebug(1,Form(" n.%3i ThetaCer %8.3f PhiCer %8.3f check %i",i,fThetaPhot[i],fPhiPhot[i],fClCk[i])); | |
160 | } | |
161 | AliDebug(1,Form("CkovHiddenTrk: thetaC %f th %f ph %f",fCkovFit,fThTrkFit*TMath::RadToDeg(),fPhTrkFit*TMath::RadToDeg())); | |
5a3482a0 | 162 | |
9785d5fb | 163 | return kTRUE; |
164 | ||
5a3482a0 | 165 | }//CkovHiddenTrk() |
166 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
167 | Bool_t AliHMPIDReconHTA::DoRecHiddenTrk() | |
168 | { | |
169 | // Pattern recognition method without any infos from tracking... | |
170 | // First a preclustering filter to avoid part of the noise | |
171 | // Then only ellipsed-rings are fitted (no possibility, | |
172 | // for the moment, to reconstruct very inclined tracks) | |
173 | // Finally a fitting with (th,ph) free, starting by very close values | |
174 | // previously evaluated. | |
175 | // Arguments: none | |
176 | // Returns: none | |
177 | Double_t thTrkRec,phiTrkRec,thetaCRec; | |
178 | ||
ee2f3539 | 179 | if(!FindShape(thTrkRec,phiTrkRec,thetaCRec)) { |
2e81abb0 | 180 | AliDebug(1,Form(" FindShape failed...!")); |
ee2f3539 | 181 | return kFALSE; |
182 | } | |
2e81abb0 | 183 | AliDebug(1,Form(" FindShape accepted...!")); |
5a3482a0 | 184 | |
2e81abb0 | 185 | if(!FitRing(thTrkRec,phiTrkRec)) { |
186 | AliDebug(1,Form(" FitRing failed...!")); | |
ee2f3539 | 187 | return kFALSE; |
188 | } | |
2e81abb0 | 189 | AliDebug(1,Form(" FitRing accepted...!")); |
190 | ||
191 | if(!UniformDistrib()) { | |
192 | AliDebug(1,Form(" UniformDistrib failed...!")); | |
193 | return kFALSE; | |
194 | } | |
195 | AliDebug(1,Form(" UniformDistrib accepted...!")); | |
196 | ||
5a3482a0 | 197 | return kTRUE; |
198 | }//DoRecHiddenTrk() | |
5a3482a0 | 199 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
200 | Bool_t AliHMPIDReconHTA::CluPreFilter(TClonesArray *pCluLst) | |
201 | { | |
aa00f952 | 202 | // Pre-filter of bkg clusters |
203 | // Arguments: pSluLst - List of the clusters for a given chamber | |
204 | // Returns: status - TRUE if filtering leaves enough photons, FALSE if not | |
5a3482a0 | 205 | // |
206 | Int_t nClusTot = pCluLst->GetEntriesFast(); | |
207 | if(nClusTot<4||nClusTot>100) { | |
208 | return kFALSE; | |
209 | } else { | |
210 | InitVars(nClusTot); | |
211 | return kTRUE; | |
212 | } | |
213 | } | |
214 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
215 | Bool_t AliHMPIDReconHTA::FindShape(Double_t &thTrkRec,Double_t &phiTrkRec,Double_t &thetaCRec) | |
216 | { | |
aa00f952 | 217 | // Finds the estimates for phi and theta of the track and the ThetaCerenkov |
218 | // by using a database of the shapes of the rings | |
219 | // Arguments: none | |
220 | // Returns: thTrkRec - estimate of theta track | |
221 | // phiTrkRec - estimate of phi track | |
222 | // thetaCRec - estimate of ThetaCerenkov | |
223 | // status - TRUE if a good solution is found, FALSE if not | |
224 | ||
aa8db9d6 | 225 | Double_t phiphot[1000]; |
226 | Double_t dist[1000]; | |
227 | Int_t indphi[1000]; | |
5a3482a0 | 228 | |
229 | Bool_t status; | |
230 | ||
aa8db9d6 | 231 | if(fNClu>1000) return kFALSE; // too many clusters.... |
232 | ||
9785d5fb | 233 | // Sort in phi angle... |
5a3482a0 | 234 | for(Int_t i=0;i<fNClu;i++) { |
9785d5fb | 235 | if(!fClCk[i]) { |
2e81abb0 | 236 | AliDebug(1,Form(" n.%3i xMIP %8.3f yMIP %8.3f check %i",i,fMipX,fMipY,fClCk[i])); |
9785d5fb | 237 | phiphot[i] = 999.; |
2e81abb0 | 238 | dist[i] = 999.; |
9785d5fb | 239 | continue; |
240 | } | |
5a3482a0 | 241 | phiphot[i] = (TMath::ATan2(fMipY-fYClu[i],fMipX-fXClu[i])+TMath::Pi())*TMath::RadToDeg(); |
242 | dist[i]=TMath::Sqrt((fMipX-fXClu[i])*(fMipX-fXClu[i])+(fMipY-fYClu[i])*(fMipY-fYClu[i])); | |
2e81abb0 | 243 | AliDebug(1,Form(" n.%3i phiphot %8.3f dist %8.3f check %i",i,phiphot[i],dist[i],fClCk[i])); |
5a3482a0 | 244 | } |
245 | ||
246 | TMath::Sort(fNClu,phiphot,indphi,kFALSE); | |
247 | ||
248 | // Purify with a truncated mean; | |
249 | Int_t np=0; | |
250 | Double_t dMean = 0; | |
251 | Double_t dMean2 = 0; | |
252 | for(Int_t i=0;i<fNClu;i++) { | |
253 | if(!fClCk[indphi[i]]) continue; // Check if a good photon candidate or not | |
254 | dMean +=dist[indphi[i]]; | |
255 | dMean2+=dist[indphi[i]]*dist[indphi[i]]; | |
256 | np++; | |
257 | } | |
258 | ||
259 | dMean /=(Double_t)np; | |
260 | dMean2 /=(Double_t)np; | |
261 | Double_t rms = TMath::Sqrt(dMean2 - dMean*dMean); | |
262 | ||
263 | for(Int_t i=0;i<fNClu;i++) { | |
264 | if(!fClCk[indphi[i]]) continue; // Check if a good photon candidate or not | |
2e81abb0 | 265 | if(TMath::Abs(dMean-dist[indphi[i]]) > 1.5*rms) { |
5a3482a0 | 266 | fClCk[indphi[i]] = kFALSE; |
267 | continue; | |
268 | } | |
269 | } | |
270 | ||
2e81abb0 | 271 | AliDebug(1,"Purification of photons..."); |
5a3482a0 | 272 | // |
273 | // purify vectors for good photon candidates | |
274 | // | |
275 | Int_t npeff=0; | |
276 | Double_t *phiphotP = new Double_t[fNClu+1]; | |
277 | Double_t *distP = new Double_t[fNClu+1]; | |
278 | for(Int_t i=0;i<fNClu;i++) { | |
2e81abb0 | 279 | AliDebug(1,Form(" n. %3i phiphot %8.3f dist %8.3f check %i",i,phiphot[indphi[i]],dist[indphi[i]],fClCk[indphi[i]])); |
5a3482a0 | 280 | if(!fClCk[indphi[i]]) continue; // Check if a good photon candidate or not |
281 | phiphotP[npeff] = phiphot[indphi[i]]; | |
282 | distP[npeff] = dist[indphi[i]]; | |
283 | npeff++; | |
284 | } | |
285 | ||
5a3482a0 | 286 | if(npeff<3) { |
2e81abb0 | 287 | AliDebug(1,Form("FindShape failed: no enough photons = %i...",npeff)); |
5a3482a0 | 288 | delete [] phiphotP; |
289 | delete [] distP; | |
290 | return kFALSE; | |
291 | } | |
292 | ||
5a3482a0 | 293 | Double_t xA,xB; |
9785d5fb | 294 | status = kFALSE; |
ee2f3539 | 295 | |
2e81abb0 | 296 | if(!ShapeModel(npeff,phiphotP,distP,xA,xB,phiTrkRec)) {AliDebug(1,Form("ShapeModel failed ")); return kFALSE;} |
297 | ||
298 | // if(xA > 50 || xB > 15) {AliDebug(1,Form("xA and xB failed out of range")); return kFALSE;} | |
5a3482a0 | 299 | |
170a4fc5 | 300 | Int_t binxDB,binyDB; |
2e81abb0 | 301 | Int_t compactDB=-1; |
302 | ||
303 | if(xA > xB) { //geometrically not possible, try to recover on a mean values... | |
304 | ||
305 | FindBinDB(xA,xA,binxDB,binyDB); | |
306 | if(binxDB<0 || binyDB<0) {AliDebug(1,Form("bin < 0 ! failed ")); return kFALSE;} | |
307 | Int_t compactDB1 = CompactDB(binxDB,binyDB); | |
308 | FindBinDB(xB,xB,binxDB,binyDB); | |
309 | if(binxDB<0 || binyDB<0) {AliDebug(1,Form("bin < 0 ! failed ")); return kFALSE;} | |
310 | Int_t compactDB2 = CompactDB(binxDB,binyDB); | |
311 | Double_t thetaCRec1 = (Double_t)(compactDB1%1000); | |
312 | Double_t thetaCRec2 = (Double_t)(compactDB2%1000); | |
313 | Double_t thTrkRec1 = (Double_t)(compactDB1/1000); | |
314 | Double_t thTrkRec2 = (Double_t)(compactDB2/1000); | |
315 | thetaCRec = 0.5*(thetaCRec1+thetaCRec2); | |
316 | thTrkRec = 0.5*( thTrkRec1+ thTrkRec2); | |
317 | ||
318 | } else { | |
319 | ||
320 | FindBinDB(xA,xB,binxDB,binyDB); | |
321 | if(binxDB<0 || binyDB<0) {AliDebug(1,Form("bin < 0 ! failed ")); return kFALSE;} | |
5a3482a0 | 322 | |
2e81abb0 | 323 | compactDB = CompactDB(binxDB,binyDB); |
324 | ||
325 | if(compactDB<0) {AliDebug(1,Form("compact< 0! failed ")); return kFALSE;} | |
326 | // | |
327 | // | |
328 | thetaCRec = (Double_t)(compactDB%1000); | |
329 | thTrkRec = (Double_t)(compactDB/1000); | |
330 | ||
331 | } | |
332 | ||
333 | AliDebug(1,Form(" CompactDB %i thTrkRec %8.3f thetaCRec %8.3f ",compactDB,thTrkRec,thetaCRec)); | |
334 | ||
335 | phiTrkRec *= TMath::DegToRad(); | |
ee2f3539 | 336 | thTrkRec *= TMath::DegToRad(); |
337 | thetaCRec *= TMath::DegToRad(); | |
9785d5fb | 338 | |
ee2f3539 | 339 | status = kTRUE; |
5a3482a0 | 340 | |
341 | delete [] phiphotP; | |
342 | delete [] distP; | |
343 | ||
344 | return status; | |
345 | } | |
346 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
347 | Bool_t AliHMPIDReconHTA::ShapeModel(Int_t np,Double_t *phiphot,Double_t *dist,Double_t &xA,Double_t &xB,Double_t &phiStart) | |
348 | { | |
aa00f952 | 349 | // Find a Spline curve to define dist. vs. phi angle |
350 | // in order to estimate the phi of the track | |
351 | // Arguments: np - # points corresponding to # photon candidates | |
352 | // dist - distance of each photon from MIP | |
353 | // phiphot - phi of the photon in the DRS | |
354 | // Returns: xA - min. distance from MIP | |
355 | // xB - dist. from mip perpedicular to the major axis | |
356 | // phiStart- estimate of the track phi | |
357 | ||
5a3482a0 | 358 | TGraph *phigr = new TGraph(np,phiphot,dist); |
2e81abb0 | 359 | phiStart = FindSimmPhi(); |
360 | ||
361 | Double_t phiStart1 = phiStart; | |
362 | if(phiStart1 > 360) phiStart1 -= 360; | |
363 | Double_t phiStart2 = phiStart+90; | |
364 | if(phiStart2 > 360) phiStart2 -= 360; | |
365 | xA = phigr->Eval(phiStart1); | |
366 | xB = phigr->Eval(phiStart2); | |
367 | //--- | |
368 | AliDebug(1,Form("phiStart %f phiStart1 %f phiStart2 %f ",phiStart,phiStart1,phiStart2)); | |
369 | AliDebug(1,Form("xA %f xB %f",xA,xB)); | |
5a3482a0 | 370 | |
2e81abb0 | 371 | phiStart += 180; |
372 | if(phiStart>360) phiStart-=360; | |
5a3482a0 | 373 | |
5a3482a0 | 374 | return kTRUE; |
375 | } | |
376 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
aa00f952 | 377 | Double_t AliHMPIDReconHTA::VertParab(Double_t x1,Double_t y1,Double_t x2, Double_t y2, Double_t x3, Double_t y3)const |
5a3482a0 | 378 | { |
aa00f952 | 379 | // It uses parabola from 3 points to evaluate the x-coord of the parab |
380 | // Arguments: xi,yi - points | |
381 | // Returns: x-coord of the vertex | |
382 | ||
5a3482a0 | 383 | Double_t a = ((x1-x3)*(y1-y2)-(x1-x2)*(y1-y3))/((x1*x1-x2*x2)*(x1-x3)-(x1*x1-x3*x3)*(x1-x2)); |
384 | Double_t b = (y1-y2 - a*(x1*x1-x2*x2))/(x1-x2); | |
385 | // Double_t c = y1 - a*x1*x1-b*x1; | |
386 | return -0.5*b/a; | |
387 | } | |
388 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
2e81abb0 | 389 | Bool_t AliHMPIDReconHTA::FitRing(Double_t thTrkRec,Double_t phiTrkRec) |
390 | { | |
391 | Double_t th = thTrkRec; | |
392 | Double_t ph = phiTrkRec; | |
393 | ||
394 | FitFree(th,ph); | |
395 | while(FitStatus()) { | |
396 | th = fThTrkFit; | |
397 | ph = fPhTrkFit; | |
398 | FitFree(th,ph); | |
399 | } | |
400 | return kTRUE; | |
401 | } | |
402 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
5a3482a0 | 403 | Bool_t AliHMPIDReconHTA::FitFree(Double_t thTrkRec,Double_t phiTrkRec) |
404 | { | |
405 | // Fit performed by minimizing RMS/sqrt(n) of the | |
406 | // photons reconstructed. First phi is fixed and theta | |
407 | // is fouond, then (th,ph) of the track | |
408 | // as free parameters | |
409 | // Arguments: PhiRec phi of the track | |
410 | // Returns: none | |
411 | ||
5a3482a0 | 412 | TMinuit *pMinuit = new TMinuit(2); |
413 | pMinuit->mncler(); | |
414 | gMinuit->SetObjectFit((TObject*)this); gMinuit->SetFCN(AliHMPIDReconHTA::FunMinPhot); //set fit function | |
415 | Double_t aArg=-1,parStep,parLow,parHigh; Int_t iErrFlg; //tmp vars for TMinuit | |
416 | Double_t d1,d2,d3; | |
417 | TString sName; | |
418 | Double_t th,ph; | |
419 | ||
420 | pMinuit->mnexcm("SET PRI",&aArg,1,iErrFlg); //suspend all printout from TMinuit | |
421 | pMinuit->mnexcm("SET NOW",&aArg,0,iErrFlg); | |
422 | ||
423 | if(thTrkRec==0) thTrkRec = 3.*TMath::DegToRad(); // not to start from the edge... | |
424 | ||
2e81abb0 | 425 | AliDebug(1,Form(" Minimization STARTED with phiTRK %6.2f thetaTRK %6.2f",phiTrkRec,thTrkRec)); |
426 | ||
5a3482a0 | 427 | pMinuit->mnparm(0," thTrk ",thTrkRec ,parStep=0.01,parLow=0,parHigh=TMath::PiOver4(),iErrFlg); |
428 | pMinuit->mnparm(1," phiTrk ",phiTrkRec,parStep=0.01,parLow=0,parHigh=TMath::TwoPi(),iErrFlg); | |
429 | ||
430 | pMinuit->FixParameter(1); | |
431 | pMinuit->mnexcm("SIMPLEX" ,&aArg,0,iErrFlg); | |
432 | pMinuit->mnexcm("MIGRAD" ,&aArg,0,iErrFlg); | |
433 | pMinuit->Release(1); | |
434 | pMinuit->mnexcm("MIGRAD" ,&aArg,0,iErrFlg); | |
435 | ||
436 | pMinuit->mnpout(0,sName,th,d1,d2,d3,iErrFlg); | |
437 | pMinuit->mnpout(1,sName,ph,d1,d2,d3,iErrFlg); | |
438 | ||
439 | Double_t f,par[2]; | |
440 | Double_t *grad=0x0; | |
441 | par[0] = th;par[1] = ph; | |
442 | pMinuit->Eval(2,grad,f,par,3); | |
443 | ||
5a3482a0 | 444 | SetTrkFit(th,ph); |
445 | return kTRUE; | |
446 | } | |
447 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
448 | void AliHMPIDReconHTA::FunMinPhot(Int_t &/* */,Double_t* /* */,Double_t &f,Double_t *par,Int_t iflag) | |
449 | { | |
450 | // Minimization function to find best track and thetaC parameters | |
451 | // Arguments: f = function value to minimize | |
452 | // par = list of parameter to find | |
453 | // iflag = flag status. See Minuit instructions | |
454 | // Returns: none | |
455 | // | |
456 | // Note: it is necessary to call an instance of AlihMPIDParam. Not possible to use fParam | |
457 | // because of the static instantiation of the function in Minuit | |
458 | ||
459 | AliHMPIDParam *pParam=AliHMPIDParam::Instance(); | |
460 | AliHMPIDReconHTA *pRecHTA=(AliHMPIDReconHTA*)gMinuit->GetObjectFit(); | |
461 | AliHMPIDRecon pRec; | |
462 | Double_t sizeCh = 0.5*pParam->RadThick()+pParam->WinThick()+pParam->GapThick(); | |
463 | Double_t thTrk = par[0]; | |
464 | Double_t phTrk = par[1]; | |
465 | Double_t xrad = pRecHTA->MipX() - sizeCh*TMath::Tan(thTrk)*TMath::Cos(phTrk); | |
466 | Double_t yrad = pRecHTA->MipY() - sizeCh*TMath::Tan(thTrk)*TMath::Sin(phTrk); | |
467 | pRecHTA->SetRadXY(xrad,yrad); | |
468 | pRec.SetTrack(xrad,yrad,thTrk,phTrk); | |
469 | ||
470 | Double_t meanCkov =0; | |
471 | Double_t meanCkov2=0; | |
472 | Double_t thetaCer,phiCer; | |
473 | Int_t nClAcc = 0; | |
474 | Int_t nClTot=pRecHTA->NClu(); | |
2e81abb0 | 475 | |
5a3482a0 | 476 | for(Int_t i=0;i<nClTot;i++) { |
2e81abb0 | 477 | |
478 | if(pRecHTA->IdxMip() == i) { | |
479 | pRecHTA->SetPhotAngles(i,999.,999.); | |
480 | continue; | |
481 | } | |
482 | ||
5a3482a0 | 483 | if(!(pRecHTA->ClCk(i))) continue; |
2e81abb0 | 484 | |
485 | Bool_t status = pRec.FindPhotCkov(pRecHTA->XClu(i),pRecHTA->YClu(i),thetaCer,phiCer); | |
486 | if(!status) { | |
487 | pRecHTA->SetPhotAngles(i,999.,999.); | |
488 | continue; | |
489 | } | |
490 | pRecHTA->SetPhotAngles(i,thetaCer,phiCer); | |
5a3482a0 | 491 | meanCkov += thetaCer; |
492 | meanCkov2 += thetaCer*thetaCer; | |
493 | nClAcc++; | |
494 | } | |
2e81abb0 | 495 | |
5a3482a0 | 496 | if(nClAcc==0) {f=999;return;} |
497 | meanCkov /=(Double_t)nClAcc; | |
498 | meanCkov2 /=(Double_t)nClAcc; | |
499 | Double_t rms = TMath::Sqrt(TMath::Abs(meanCkov2 - meanCkov*meanCkov)); | |
500 | f = rms/TMath::Sqrt((Double_t)nClAcc); | |
501 | ||
502 | if(iflag==3) { | |
2e81abb0 | 503 | Int_t nClAccStep1 = nClAcc; |
5a3482a0 | 504 | nClAcc = 0; |
505 | Double_t meanCkov1=0; | |
876a9d77 | 506 | Double_t meanCkov3=0; |
5a3482a0 | 507 | for(Int_t i=0;i<nClTot;i++) { |
2e81abb0 | 508 | |
509 | if(!(pRecHTA->ClCk(i))) continue; | |
510 | thetaCer = pRecHTA->PhotTheta(i); | |
5a3482a0 | 511 | if(TMath::Abs(thetaCer-meanCkov)<2*rms) { |
512 | meanCkov1 += thetaCer; | |
876a9d77 | 513 | meanCkov3 += thetaCer*thetaCer; |
5a3482a0 | 514 | nClAcc++; |
2e81abb0 | 515 | } else pRecHTA->SetClCk(i,kFALSE); |
516 | } | |
517 | ||
518 | if(nClAcc<3) { | |
519 | pRecHTA->SetFitStatus(kFALSE); | |
520 | pRecHTA->SetCkovFit(meanCkov); | |
521 | pRecHTA->SetCkovSig2(rms*rms); | |
522 | pRecHTA->SetNCluFit(nClAccStep1); | |
523 | return; | |
5a3482a0 | 524 | } |
2e81abb0 | 525 | |
5a3482a0 | 526 | meanCkov1/=nClAcc; |
876a9d77 | 527 | Double_t rms2 = (meanCkov3 - meanCkov*meanCkov*nClAcc)/nClAcc; |
6d7793cf | 528 | |
529 | // get a logger instance | |
530 | // what for?? | |
731dd0f2 | 531 | //AliLog::GetRootLogger(); |
6d7793cf | 532 | |
2e81abb0 | 533 | if(nClAcc!=nClAccStep1) pRecHTA->SetFitStatus(kTRUE); else pRecHTA->SetFitStatus(kFALSE); |
534 | ||
535 | pRecHTA->SetCkovFit(meanCkov1); | |
536 | pRecHTA->SetCkovSig2(rms2); | |
537 | pRecHTA->SetNCluFit(nClAcc); | |
5a3482a0 | 538 | } |
539 | }//FunMinPhot() | |
540 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
541 | void AliHMPIDReconHTA::InitDatabase() | |
542 | { | |
aa00f952 | 543 | // Construction a database of ring shapes on fly |
544 | // Arguments: none | |
545 | // Returns : none | |
170a4fc5 | 546 | // N.B. fgDB is the distance with x-min from MIP |
547 | // y-dist from the ring of the MIP perpendicular to major axis | |
aa00f952 | 548 | // The content is the packed info of track theta and thetaC in degrees |
549 | // thetaC+1000*thTrk | |
550 | // | |
9785d5fb | 551 | // TFile *pout = new TFile("./database.root","recreate"); |
ee2f3539 | 552 | |
170a4fc5 | 553 | static Bool_t isDone = kFALSE; |
554 | ||
5a3482a0 | 555 | TStopwatch timer; |
170a4fc5 | 556 | |
557 | if(isDone) { | |
5565f017 | 558 | return; |
559 | } | |
170a4fc5 | 560 | |
561 | if(!isDone) { | |
562 | timer.Start(); | |
563 | } | |
564 | ||
5a3482a0 | 565 | AliInfo(Form("database HTA is being built.Please, wait...")); |
566 | // | |
c61a7285 | 567 | Double_t x[3]={0,0,0},y[3]; |
5a3482a0 | 568 | |
569 | AliHMPIDRecon rec; | |
570 | ||
571 | if(!fParam) fParam=AliHMPIDParam::Instance(); | |
572 | Double_t thetaMax = TMath::ACos(1./fParam->MeanIdxRad()); | |
573 | Double_t thTrkMax = 1./TMath::ASin(fParam->MeanIdxRad()); | |
574 | ||
575 | Int_t nstepx = 1000; | |
576 | Int_t nstepy = 1000; | |
577 | ||
5a3482a0 | 578 | // |
579 | Double_t xrad = 0; | |
580 | Double_t yrad = 0; | |
581 | Double_t phTrk = 0; | |
582 | ||
583 | for(Int_t i=0;i<nstepx;i++) { //loop on thetaC | |
584 | for(Int_t j=0;j<nstepy;j++) { //loop on theta particle | |
585 | Double_t thetaC = thetaMax/nstepx*((Double_t)i+0.5); | |
586 | Double_t thTrk = thTrkMax/nstepy*((Double_t)j+0.5); | |
587 | // | |
588 | //mip position | |
589 | // | |
fcaff63d | 590 | Double_t sizeCh = 0.5*fParam->RadThick()+fParam->WinThick()+fParam->GapThick(); |
5a3482a0 | 591 | Double_t xmip = xrad + sizeCh*TMath::Tan(thTrk)*TMath::Cos(phTrk); |
592 | Double_t ymip = yrad + sizeCh*TMath::Tan(thTrk)*TMath::Sin(phTrk); | |
593 | ||
594 | Double_t dist1,dist2; | |
595 | // | |
596 | //first point at phi=0 | |
597 | // | |
598 | rec.SetTrack(xrad,yrad,thTrk,phTrk); | |
599 | TVector2 pos; | |
600 | pos=rec.TracePhot(thetaC,0); | |
601 | ||
602 | if(pos.X()==-999) { | |
ee2f3539 | 603 | dist1 = 0; //open ring...only the distance btw mip and point at 180 will be considered |
5a3482a0 | 604 | } else { |
605 | x[0] = pos.X(); y[0] = pos.Y(); | |
606 | dist1 = TMath::Sqrt((x[0]-xmip)*(x[0]-xmip)+(y[0]-ymip)*(y[0]-ymip)); | |
607 | } | |
608 | // | |
609 | //second point at phi=180 | |
610 | // | |
611 | rec.SetTrack(xrad,yrad,thTrk,phTrk); | |
612 | pos=rec.TracePhot(thetaC,TMath::Pi()); | |
613 | ||
2e81abb0 | 614 | if(pos.X()==-999) {AliDebug(1,Form("it should not happens!Bye"));return;} |
5a3482a0 | 615 | x[1] = pos.X(); y[1] = pos.Y(); |
616 | if((x[1]-xmip)*(x[0]-xmip)>0) continue; // to avoid circles out mips (for very low ThetaC) | |
617 | dist2 = TMath::Sqrt((x[1]-xmip)*(x[1]-xmip)+(y[1]-ymip)*(y[1]-ymip)); | |
618 | ||
619 | // Double_t distA = dist1+dist2; | |
620 | Double_t distA = dist2; // only the minimum: problem of acceptance | |
621 | // | |
622 | //second point at phi=90 | |
623 | // | |
624 | rec.SetTrack(xrad,yrad,thTrk,phTrk); | |
625 | pos=rec.TracePhot(thetaC,TMath::PiOver2()); | |
626 | ||
627 | if(pos.X()==-999) continue; | |
628 | x[2] = pos.X(); y[2] = pos.Y(); | |
629 | Double_t distB = TMath::Sqrt((x[2]-xmip)*(x[2]-xmip)+(y[2]-ymip)*(y[2]-ymip)); | |
630 | // compact the infos... | |
631 | Int_t compact = (Int_t)(thetaC*TMath::RadToDeg())+1000*(Int_t)(thTrk*TMath::RadToDeg()); | |
ee2f3539 | 632 | Int_t binxDB,binyDB; |
633 | FindBinDB(distA,distB,binxDB,binyDB); | |
170a4fc5 | 634 | if(fgDB[binxDB][binyDB]==0) fgDB[binxDB][binyDB] = compact; |
5a3482a0 | 635 | } |
636 | } | |
637 | ||
9785d5fb | 638 | FillZeroChan(); |
5a3482a0 | 639 | |
170a4fc5 | 640 | if(!isDone) { |
641 | timer.Stop(); | |
642 | Double_t nSecs = timer.CpuTime(); | |
643 | AliInfo(Form("database HTA successfully open in %3.1f sec.(CPU). Reconstruction is started.",nSecs)); | |
644 | isDone = kTRUE; | |
645 | } | |
9785d5fb | 646 | |
647 | // pout->Write(); | |
648 | // pout->Close(); | |
649 | ||
9a573d52 | 650 | }//InitDatabase() |
5a3482a0 | 651 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
9785d5fb | 652 | void AliHMPIDReconHTA::FillZeroChan()const |
5a3482a0 | 653 | { |
606697a8 | 654 | //If fills eventually channel without entries |
655 | //inthe histo "database" jyst interpolating the neighboring cells | |
656 | // Arguments: histogram pointer of the database | |
657 | // Returns: none | |
658 | // | |
170a4fc5 | 659 | const Int_t nxDB = 500; |
660 | const Int_t nyDB = 150; | |
661 | ||
662 | for(Int_t i = 0;i<nxDB;i++) { | |
663 | for(Int_t j = 0;j<nyDB;j++) { | |
664 | if(fgDB[i][j] == 0) { | |
665 | fgDB[i][j] = -1; | |
5a3482a0 | 666 | Int_t nXmin = i-1; Int_t nXmax=i+1; |
667 | Int_t nYmin = j-1; Int_t nYmax=j+1; | |
668 | Int_t nc = 0; | |
669 | Double_t meanC =0; | |
670 | Double_t meanTrk =0; | |
671 | for(Int_t ix=nXmin;ix<=nXmax;ix++) { | |
170a4fc5 | 672 | if(ix<0||ix>=nxDB) continue; |
5a3482a0 | 673 | for(Int_t iy=nYmin;iy<=nYmax;iy++) { |
170a4fc5 | 674 | if(iy<0||iy>=nyDB) continue; |
675 | meanC += (Int_t)(fgDB[ix][iy]%1000); | |
676 | meanTrk+= (Int_t)(fgDB[ix][iy]/1000); | |
5a3482a0 | 677 | nc++; |
678 | } | |
679 | meanC/=nc; meanTrk/=nc; | |
680 | Int_t compact = (Int_t)meanC+1000*(Int_t)meanTrk; | |
170a4fc5 | 681 | if(compact>0)fgDB[i][j] = compact; |
5a3482a0 | 682 | } |
683 | } | |
684 | } | |
685 | } | |
9a573d52 | 686 | }//FillZeroChan() |
5a3482a0 | 687 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
ee2f3539 | 688 | Int_t AliHMPIDReconHTA::r2(Double_t *coef, Double_t &x1, Double_t &x2) |
689 | { | |
690 | //2nd deg. equation | |
691 | //solution | |
692 | // Arguments: coef 2 1 0: ax^2+bx+c=0 | |
693 | // Returns: n. of solutions | |
694 | // x1= 1st sol | |
695 | // x2= 2nd sol | |
696 | Double_t a,b,c; | |
697 | a = coef[2]; | |
698 | b = coef[1]; | |
699 | c = coef[0]; | |
700 | Double_t delta = b*b-4*a*c; | |
701 | if(delta<0) {return 0;} | |
702 | if(delta==0) { | |
703 | x1=x2=-b/(2*a); | |
704 | return 1; | |
705 | } | |
3d1548fd | 706 | if(a==0) { |
707 | x1 = -c/b; | |
708 | return 1; | |
709 | } | |
ee2f3539 | 710 | // delta>0 |
711 | x1 = (-b+TMath::Sqrt(delta))/(2*a); | |
712 | x2 = (-b-TMath::Sqrt(delta))/(2*a); | |
713 | return 2; | |
714 | }//r2() | |
715 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
9a573d52 | 716 | |
717 | Double_t AliHMPIDReconHTA::FindSimmPhi() | |
ee2f3539 | 718 | { |
719 | //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
2e81abb0 | 720 | // Reconstruction of phiTRK angle with two methods (in switching) |
721 | // | |
722 | // - least square method (for closed rings) | |
723 | // - by minimum distance mip-photon (for open rings) | |
ee2f3539 | 724 | |
725 | Float_t coeff1ord=0; Float_t coeff2ord=0; Float_t coeff0ord=0; | |
9a573d52 | 726 | Float_t xrotsumm =0; Float_t yrotsumm =0; Float_t xx =0; |
2e81abb0 | 727 | Float_t yy =0; Float_t xy =0; |
ee2f3539 | 728 | Double_t xmin=0; |
729 | Double_t ymin=0; | |
9a573d52 | 730 | |
731 | Int_t np=0; | |
ee2f3539 | 732 | |
733 | Double_t distMin = 999.; | |
734 | ||
9a573d52 | 735 | for(Int_t i=0;i<fNClu;i++) { |
736 | if(!fClCk[i]) continue; | |
737 | np++; | |
738 | xrotsumm+=fXClu[i]; // summ xi | |
739 | yrotsumm+=fYClu[i]; // summ yi | |
740 | xx+=fXClu[i]*fXClu[i]; // summ xixi | |
741 | yy+=fYClu[i]*fYClu[i]; // summ yiyi | |
742 | xy+=fXClu[i]*fYClu[i]; // summ yixi | |
ee2f3539 | 743 | Double_t dist2= (fXClu[i]-fMipX)*(fXClu[i]-fMipX)+(fYClu[i]-fMipY)*(fYClu[i]-fMipY); |
744 | if(dist2<distMin) { | |
745 | distMin = dist2; | |
746 | xmin = fXClu[i]; | |
747 | ymin = fYClu[i]; | |
748 | } | |
9a573d52 | 749 | } |
ee2f3539 | 750 | |
2e81abb0 | 751 | Double_t AngM = TMath::ATan2(ymin-fMipY,(xmin-fMipX))*TMath::RadToDeg(); |
752 | if (AngM<0) AngM+=360; | |
753 | ||
754 | AliDebug(1,Form(" Simple angle prediction with MIN phi = %f",AngM)); | |
9a573d52 | 755 | |
ee2f3539 | 756 | //_____calc. met min quadr using effective distance _________________________________________________ |
9a573d52 | 757 | |
2e81abb0 | 758 | coeff2ord = xy-xrotsumm*yrotsumm/np; |
759 | coeff1ord = yrotsumm*yrotsumm/np - xrotsumm*xrotsumm/np - yy + xx; | |
760 | coeff0ord = -coeff2ord; | |
761 | ||
762 | AliDebug(1,Form(" a = %f b = %f c = %f",coeff2ord,coeff1ord,coeff0ord)); | |
763 | ||
ee2f3539 | 764 | Double_t m1=0, m2=0; Double_t n1=0, n2=0; |
765 | // c // b // a | |
766 | Double_t coeff[3]={coeff0ord,coeff1ord,coeff2ord}; | |
9a573d52 | 767 | |
768 | r2(coeff,m1,m2); | |
769 | ||
770 | n1=(yrotsumm-m1*xrotsumm)/np; | |
771 | n2=(yrotsumm-m2*xrotsumm)/np; | |
ee2f3539 | 772 | // 2 solutions................. |
773 | ||
ee2f3539 | 774 | // negative angles solved... |
775 | ||
9a573d52 | 776 | Double_t d1 =(1/(m1*m1+1))*(yy+m1*m1*xx+np*n1*n1-2*m1*xy-2*n1*yrotsumm+2*m1*n1*xrotsumm); |
777 | Double_t d2 =(1/(m2*m2+1))*(yy+m2*m2*xx+np*n2*n2-2*m2*xy-2*n2*yrotsumm+2*m2*n2*xrotsumm); | |
2e81abb0 | 778 | |
779 | AliDebug(1,Form(" predicted distance d1 %f for angle %6.2f d2 %f for angle %6.2f",d1,TMath::ATan(m1)*TMath::RadToDeg(), | |
780 | d2,TMath::ATan(m2)*TMath::RadToDeg())); | |
9a573d52 | 781 | Double_t mMin; |
782 | if(d1 > d2) mMin = m2; else mMin = m1; | |
783 | ||
2e81abb0 | 784 | Double_t phiTrk=0; |
ee2f3539 | 785 | // |
2e81abb0 | 786 | if(ymin < fMipY && xmin > fMipX) {phiTrk = TMath::ATan(mMin)*TMath::RadToDeg()+180;} |
787 | if(ymin > fMipY && xmin < fMipX) {phiTrk = TMath::ATan(mMin)*TMath::RadToDeg()+360;} | |
788 | if(ymin > fMipY && xmin > fMipX) {phiTrk = TMath::ATan(mMin)*TMath::RadToDeg()+180;} | |
789 | if(ymin < fMipY && xmin < fMipX) {phiTrk = TMath::ATan(mMin)*TMath::RadToDeg();} | |
790 | if(ymin == fMipY && xmin > fMipX) {phiTrk = TMath::ATan(mMin)*TMath::RadToDeg()+180;} | |
791 | if(ymin == fMipY && xmin < fMipX) {phiTrk = TMath::ATan(mMin)*TMath::RadToDeg();} | |
792 | if(ymin < fMipY && xmin == fMipX) {phiTrk = 90;} | |
793 | if(ymin > fMipY && xmin == fMipX) {phiTrk = 270;} | |
9a573d52 | 794 | |
ee2f3539 | 795 | // ------------------------- choose the best----------------------- |
9a573d52 | 796 | |
ee2f3539 | 797 | |
2e81abb0 | 798 | if( AngM-40 <= phiTrk && AngM+40 >= phiTrk) return phiTrk; else return AngM; |
9a573d52 | 799 | } |
9a573d52 | 800 | //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
ee2f3539 | 801 | void AliHMPIDReconHTA::FindBinDB(Double_t x,Double_t y,Int_t &binX,Int_t &binY) |
9a573d52 | 802 | { |
170a4fc5 | 803 | const Int_t nxDB = 500; |
804 | const Int_t nyDB = 150; | |
ee2f3539 | 805 | const Double_t xlowDB = 0; |
806 | const Double_t xhigDB = 50; | |
807 | const Double_t ylowDB = 0; | |
2e81abb0 | 808 | const Double_t yhigDB = 30; |
ee2f3539 | 809 | |
810 | binX = -1; | |
811 | binY = -1; | |
170a4fc5 | 812 | if(x<xlowDB && x>xhigDB && |
813 | y<ylowDB && y>yhigDB) return; | |
814 | binX = Int_t((x-xlowDB)/(xhigDB-xlowDB)*nxDB); | |
815 | binY = Int_t((y-ylowDB)/(yhigDB-ylowDB)*nyDB); | |
816 | if(binX>=nxDB || binY>=nyDB) { | |
817 | binX = -1; | |
818 | binY = -1; | |
819 | } | |
820 | ||
ee2f3539 | 821 | } |
2e81abb0 | 822 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
823 | Bool_t AliHMPIDReconHTA::UniformDistrib() | |
824 | { | |
825 | AliHMPIDParam *pParam=AliHMPIDParam::Instance(); | |
826 | AliHMPIDRecon pRec; | |
827 | ||
828 | Double_t sizeCh = 0.5*pParam->RadThick()+pParam->WinThick()+pParam->GapThick(); | |
829 | Double_t xrad = MipX() - sizeCh*TMath::Tan(fThTrkFit)*TMath::Cos(fPhTrkFit); | |
830 | Double_t yrad = MipY() - sizeCh*TMath::Tan(fThTrkFit)*TMath::Sin(fPhTrkFit); | |
831 | ||
832 | pRec.SetTrack(xrad,yrad,fThTrkFit,fPhTrkFit); | |
833 | ||
834 | Int_t npeff=0; | |
835 | Int_t nPhotPerBin = 4; | |
836 | for(Int_t i=0;i<fNClu;i++) { | |
837 | if(!ClCk(i)) continue; | |
838 | npeff++; | |
839 | } | |
840 | ||
841 | Int_t nPhiBins = npeff/nPhotPerBin+1; | |
842 | if(nPhiBins<=1) nPhiBins = 2; | |
843 | ||
844 | Double_t *iPhiBin; | |
845 | iPhiBin = new Double_t[nPhiBins]; | |
846 | ||
847 | for(Int_t i=0;i<nPhiBins;i++) iPhiBin[i] =0; | |
848 | ||
849 | for(Int_t i=0;i<fNClu;i++) { | |
850 | if(!ClCk(i)) continue; | |
851 | Double_t phiCer = PhotPhi(i); | |
852 | ||
853 | Double_t PhiProva = phiCer*TMath::RadToDeg(); | |
854 | if(PhiProva<0) PhiProva+= 360; | |
855 | Int_t index = (Int_t)((Float_t)nPhiBins*PhiProva/360.); | |
856 | iPhiBin[index]++; | |
857 | } | |
858 | ||
859 | Double_t chi2 = 0; | |
860 | for(Int_t i=0;i<nPhiBins;i++) { | |
861 | Double_t theo = (Double_t)npeff/(Double_t)nPhiBins; | |
3d1548fd | 862 | if(theo==0) continue; |
2e81abb0 | 863 | chi2+= (iPhiBin[i] - theo)*(iPhiBin[i] - theo)/theo; |
864 | } | |
865 | ||
4ce766eb | 866 | delete [] iPhiBin; |
2e81abb0 | 867 | |
868 | Double_t prob = TMath::Prob(chi2, nPhiBins-1); | |
869 | AliDebug(1,Form(" Probability for uniform distrib: %6f.3 %s",prob,(prob<0.05) ? "rejected" : "accepted")); | |
870 | if(prob<0.05) return kFALSE; | |
871 | return kTRUE; | |
872 | ||
873 | } | |
874 | //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |