]>
Commit | Line | Data |
---|---|---|
5cb4dfc3 | 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 | ||
88dd9ad4 | 16 | ////////////////////////////////////////////////////////////////////////// |
17 | // // | |
18 | // AliRICHRecon // | |
19 | // // | |
20 | // RICH class to perfom pattern recognition based on Hough transfrom // | |
21 | // // | |
22 | ////////////////////////////////////////////////////////////////////////// | |
23 | ||
5cb4dfc3 | 24 | #include <Riostream.h> |
024a7e64 | 25 | #include <TCanvas.h> |
26 | #include <TFile.h> | |
5cb4dfc3 | 27 | #include <TH2.h> |
28 | #include <TMath.h> | |
024a7e64 | 29 | #include <TMath.h> |
5cb4dfc3 | 30 | #include <TMinuit.h> |
31 | #include <TNtuple.h> | |
024a7e64 | 32 | #include <TParticle.h> |
33 | #include <TRandom.h> | |
5cb4dfc3 | 34 | #include <TRotation.h> |
5cb4dfc3 | 35 | #include <TVector3.h> |
024a7e64 | 36 | |
37 | #include "AliLoader.h" | |
38 | #include "AliRICH.h" | |
39 | #include "AliRICHChamber.h" | |
40 | #include "AliRICHParam.h" | |
41 | #include "AliRICHRecon.h" | |
42 | #include "AliRun.h" | |
43 | #include "AliStack.h" | |
5cb4dfc3 | 44 | |
45 | #define NPointsOfRing 201 | |
46 | ||
88dd9ad4 | 47 | TMinuit *gAliRICHminuit ; |
5cb4dfc3 | 48 | |
8cc78f7f | 49 | void fcnrecon(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag); |
b068561d | 50 | //__________________________________________________________________________________________________ |
51 | AliRICHRecon::AliRICHRecon(const char*name, const char*title) | |
52 | :TTask(name,title) | |
5cb4dfc3 | 53 | { |
88dd9ad4 | 54 | // main ctor |
b068561d | 55 | fThetaBin=750; fThetaMin = 0.0; fThetaMax = 0.75; |
88dd9ad4 | 56 | fDTheta = 0.001; fWindowWidth = 0.060; |
57 | fNpadX = AliRICHParam::NpadsY(); | |
58 | fNpadY = AliRICHParam::NpadsX(); | |
59 | fPadSizeX = AliRICHParam::PadSizeY(); | |
60 | fPadSizeY = AliRICHParam::PadSizeX(); | |
61 | fRadiatorWidth = AliRICHParam::FreonThickness(); | |
62 | fQuartzWidth = AliRICHParam::QuartzThickness(); | |
63 | fGapWidth = AliRICHParam::RadiatorToPads(); | |
64 | fXmin = -AliRICHParam::PcSizeY()/2.; | |
65 | fXmax = AliRICHParam::PcSizeY()/2.; | |
66 | fYmin = -AliRICHParam::PcSizeX()/2.; | |
67 | fYmax = AliRICHParam::PcSizeX()/2.; | |
5cb4dfc3 | 68 | fRich = (AliRICH*)gAlice->GetDetector("RICH"); |
b068561d | 69 | fOutFile=new TFile("Anal.root","RECREATE","My Analysis histos"); |
70 | if(fIsDISPLAY) fDisplay = new TCanvas("Display","RICH Display",0,0,1200,750); | |
71 | fNtuple=new TNtuple("hn","ntuple", | |
88dd9ad4 | 72 | "Run:Trig:VertZ:Pmod:Pt:Eta:TrackTheta:TrackPhi:TrackThetaFit:TrackPhiFit:Charge::NPhotons:NPhotonsFit:InRing:MassOfParticle:HoughArea:Multiplicity:TPCLastZ"); |
5cb4dfc3 | 73 | } |
b068561d | 74 | //__________________________________________________________________________________________________ |
5cb4dfc3 | 75 | void AliRICHRecon::StartProcessEvent() |
76 | { | |
88dd9ad4 | 77 | //start to process for pattern recognition |
5cb4dfc3 | 78 | |
88dd9ad4 | 79 | Float_t trackThetaStored = 0; |
80 | Float_t trackPhiStored = 0; | |
81 | Float_t thetaCerenkovStored = 0; | |
82 | Int_t houghPhotonsStored = 0; | |
5cb4dfc3 | 83 | |
84 | SetFreonScaleFactor(0.994); | |
85 | ||
b068561d | 86 | if(fIsDISPLAY) |
5cb4dfc3 | 87 | { |
88 | DrawEvent(0); | |
b068561d | 89 | // Waiting(); |
5cb4dfc3 | 90 | } |
91 | ||
f540341d | 92 | AliLoader * richLoader = Rich()->GetLoader(); |
93 | AliRunLoader * runLoader = richLoader->GetRunLoader(); | |
94 | ||
95 | if (richLoader->TreeH() == 0x0) richLoader->LoadHits(); | |
96 | if (richLoader->TreeR() == 0x0) richLoader->LoadRecPoints(); | |
97 | if (richLoader->TreeD() == 0x0) richLoader->LoadDigits(); | |
98 | if (runLoader->TreeE() == 0x0) runLoader->LoadHeader(); | |
99 | if (runLoader->TreeK() == 0x0) runLoader->LoadKinematics(); | |
100 | ||
101 | richLoader->TreeR()->GetEntry(0); | |
5cb4dfc3 | 102 | |
f540341d | 103 | Float_t clusX[7][500],clusY[7][500]; |
104 | Int_t clusQ[7][500],clusMul[7][500]; | |
105 | Int_t nClusters[7]; | |
5cb4dfc3 | 106 | |
f540341d | 107 | for (Int_t ich=0;ich<7;ich++) { |
108 | nClusters[ich] = Rich()->Clusters(ich+1)->GetEntries(); | |
109 | for(Int_t k=0;k<nClusters[ich];k++) { | |
110 | AliRICHcluster *pCluster = (AliRICHcluster *)Rich()->Clusters(ich+1)->At(k); | |
111 | clusX[ich][k] = pCluster->X(); | |
112 | clusY[ich][k] = pCluster->Y(); | |
113 | clusQ[ich][k] = pCluster->Q(); | |
114 | clusMul[ich][k] = pCluster->Size(); | |
115 | // pCluster->Print(); | |
5cb4dfc3 | 116 | } |
f540341d | 117 | } |
118 | ||
119 | Int_t nPrimaries = (Int_t)richLoader->TreeH()->GetEntries(); | |
120 | ||
121 | cout << " N. primaries " << nPrimaries << endl; | |
122 | ||
123 | for(Int_t i=0;i<nPrimaries;i++){ | |
124 | ||
125 | richLoader->TreeH()->GetEntry(i); | |
126 | ||
127 | // Rich()->Hits()->Print(); | |
128 | Int_t iPrim = 0; | |
129 | ||
130 | AliRICHhit* pHit=0; | |
131 | ||
132 | for(Int_t j=0;j<Rich()->Hits()->GetEntries();j++) { | |
5cb4dfc3 | 133 | |
f540341d | 134 | pHit = (AliRICHhit*)Rich()->Hits()->At(j); |
135 | if(pHit->GetTrack() < nPrimaries) break; | |
136 | iPrim++; | |
137 | } | |
138 | ||
139 | cout << " iPrim " << iPrim << " pHit " << pHit << endl; | |
140 | ||
141 | if (!pHit) return; | |
142 | ||
143 | // pHit->Print(); | |
144 | ||
145 | TParticle *pParticle = runLoader->Stack()->Particle(pHit->GetTrack()); | |
146 | Float_t pmod = pParticle->P(); | |
147 | Float_t pt = pParticle->Pt(); | |
148 | Float_t trackEta = pParticle->Eta(); | |
149 | Int_t q = (Int_t)TMath::Sign(1.,pParticle->GetPDG()->Charge()); | |
150 | ||
151 | // pParticle->Print(); | |
152 | ||
153 | cout << " pmod " << pmod << " pt " << pt << " Eta " << trackEta << " charge " << q << endl; | |
154 | ||
155 | SetTrackMomentum(pmod); | |
156 | SetTrackPt(pt); | |
157 | SetTrackEta(trackEta); | |
158 | SetTrackCharge(q); | |
159 | ||
160 | TVector3 pLocal(0,0,0);//????? | |
161 | ||
162 | TVector2 primLocal =Rich()->C(pHit->C())->Glob2Loc(pHit->InX3()); | |
163 | ||
164 | // Float_t pmodFreo = pLocal.Mag(); | |
165 | Float_t trackTheta = pLocal.Theta(); | |
166 | Float_t trackPhi = pLocal.Phi(); | |
167 | ||
168 | // cout << " trackTheta " << trackTheta << " trackPhi " << trackPhi << endl; | |
169 | ||
170 | SetTrackTheta(trackTheta); | |
171 | SetTrackPhi(trackPhi); | |
172 | ||
173 | Int_t maxInd = 0; | |
174 | Float_t minDist = 999.; | |
175 | ||
176 | // cout << " n Clusters " << nClusters[pHit->Chamber()-1] << " for chamber n. " << pHit->Chamber() << endl; | |
177 | ||
178 | for(Int_t j=0;j<nClusters[pHit->Chamber()-1];j++) | |
179 | { | |
180 | Float_t diffx = primLocal.X() - clusX[pHit->Chamber()-1][j]; | |
181 | Float_t diffy = primLocal.Y() - clusY[pHit->Chamber()-1][j]; | |
182 | ||
183 | ||
184 | Float_t diff = sqrt(diffx*diffx + diffy*diffy); | |
185 | ||
186 | if(diff < minDist) | |
187 | { | |
188 | minDist = diff; | |
189 | maxInd = j; | |
190 | } | |
191 | ||
5cb4dfc3 | 192 | } |
f540341d | 193 | |
194 | Float_t diffx = primLocal.X() - clusX[pHit->Chamber()-1][maxInd]; | |
195 | Float_t diffy = primLocal.Y() - clusY[pHit->Chamber()-1][maxInd]; | |
196 | ||
197 | cout << " diffx " << diffx << " diffy " << diffy << endl; | |
198 | ||
199 | ||
200 | SetMipIndex(maxInd); | |
201 | SetTrackIndex(i); | |
202 | ||
203 | Float_t shiftX = 0;//primLocal.X()/primLocal.Z()*(fRadiatorWidth+fQuartzWidth+fGapWidth) + primLocal.X(); ????? | |
204 | Float_t shiftY = 0;//primLocal.Y()/primLocal.Z()*(fRadiatorWidth+fQuartzWidth+fGapWidth) + primLocal.Y(); ????? | |
205 | ||
206 | SetShiftX(shiftX); | |
207 | SetShiftY(shiftY); | |
208 | ||
209 | Float_t *pclusX = &clusX[pHit->Chamber()-1][0]; | |
210 | Float_t *pclusY = &clusY[pHit->Chamber()-1][0]; | |
211 | ||
212 | SetCandidatePhotonX(pclusX); | |
213 | SetCandidatePhotonY(pclusY); | |
214 | SetCandidatePhotonsNumber(nClusters[pHit->Chamber()-1]); | |
215 | ||
216 | Int_t qch = clusQ[pHit->Chamber()-1][maxInd]; | |
217 | ||
218 | ||
219 | if(minDist < 3.0 && qch > 120 && maxInd !=0) | |
220 | { | |
221 | ||
222 | if(fIsBACKGROUND) | |
223 | { | |
224 | ||
225 | Float_t xrndm = fXmin + (fXmax-fXmin)*gRandom->Rndm(280964); | |
226 | Float_t yrndm = fYmin + (fYmax-fYmin)*gRandom->Rndm(280964); | |
227 | SetShiftX(xrndm); | |
228 | SetShiftY(yrndm); | |
229 | ||
230 | } | |
231 | ||
232 | PatRec(); | |
233 | ||
234 | trackThetaStored = GetTrackTheta(); | |
235 | trackPhiStored = GetTrackPhi(); | |
236 | thetaCerenkovStored = GetThetaCerenkov(); | |
237 | houghPhotonsStored = GetHoughPhotons(); | |
238 | ||
239 | Int_t diffNPhotons = 999; | |
240 | Int_t nsteps = 0; | |
241 | Float_t diffTrackTheta = 999.; | |
242 | Float_t diffTrackPhi = 999.; | |
243 | ||
244 | while(fIsMINIMIZER && GetHoughPhotons() > 2 | |
245 | && diffNPhotons !=0 | |
246 | && diffTrackTheta > 0.0001 | |
247 | && nsteps < 10) | |
248 | { | |
249 | ||
250 | Int_t houghPhotonsBefore = GetHoughPhotons(); | |
251 | ||
252 | Float_t trackThetaBefore = GetTrackTheta(); | |
253 | Float_t trackPhiBefore = GetTrackPhi(); | |
254 | ||
255 | Minimization(); | |
256 | ||
257 | PatRec(); | |
258 | ||
259 | diffNPhotons = TMath::Abs(houghPhotonsBefore - GetHoughPhotons()); | |
260 | ||
261 | Float_t trackThetaAfter = GetTrackTheta(); | |
262 | Float_t trackPhiAfter = GetTrackPhi(); | |
263 | ||
264 | diffTrackTheta = TMath::Abs(trackThetaAfter - trackThetaBefore); | |
265 | diffTrackPhi = TMath::Abs(trackPhiAfter - trackPhiBefore); | |
266 | ||
267 | if(fDebug) | |
88dd9ad4 | 268 | cout << " houghPhotonsBefore " << houghPhotonsBefore |
5cb4dfc3 | 269 | << " GetHoughPhotons() " << GetHoughPhotons(); |
f540341d | 270 | |
271 | nsteps++; | |
272 | } | |
273 | ||
274 | SetFittedThetaCerenkov(GetThetaCerenkov()); | |
275 | SetFittedHoughPhotons(GetHoughPhotons()); | |
276 | ||
277 | SetTrackTheta(trackThetaStored); | |
278 | SetTrackPhi(trackPhiStored); | |
279 | SetThetaCerenkov(thetaCerenkovStored); | |
280 | SetHoughPhotons(houghPhotonsStored); | |
281 | ||
282 | SetMinDist(minDist); | |
283 | ||
284 | FillHistograms(); | |
285 | ||
286 | if(fIsDISPLAY) DrawEvent(1); | |
287 | ||
288 | Waiting(); | |
289 | ||
290 | } | |
291 | } | |
b068561d | 292 | if(fIsDISPLAY) fDisplay->Print("display.ps"); |
293 | }//StartProcessEvent() | |
294 | //__________________________________________________________________________________________________ | |
5cb4dfc3 | 295 | void AliRICHRecon::EndProcessEvent() |
88dd9ad4 | 296 | { |
297 | // function called at the end of the event loop | |
5cb4dfc3 | 298 | |
b068561d | 299 | fOutFile->Write(); |
300 | fOutFile->Close(); | |
5cb4dfc3 | 301 | } |
b068561d | 302 | //__________________________________________________________________________________________________ |
5cb4dfc3 | 303 | void AliRICHRecon::PatRec() |
304 | { | |
88dd9ad4 | 305 | //pattern recognition method based on Hough transform |
5cb4dfc3 | 306 | |
9eb3dcf5 | 307 | |
88dd9ad4 | 308 | Float_t trackTheta = GetTrackTheta(); |
309 | Float_t trackPhi = GetTrackPhi(); | |
5cb4dfc3 | 310 | Float_t pmod = GetTrackMomentum(); |
88dd9ad4 | 311 | Int_t iMipIndex = GetMipIndex(); |
5cb4dfc3 | 312 | |
313 | Bool_t kPatRec = kFALSE; | |
314 | ||
88dd9ad4 | 315 | Int_t candidatePhotons = 0; |
5cb4dfc3 | 316 | |
88dd9ad4 | 317 | Float_t shiftX = GetShiftX(); |
318 | Float_t shiftY = GetShiftY(); | |
5cb4dfc3 | 319 | |
88dd9ad4 | 320 | Float_t* candidatePhotonX = GetCandidatePhotonX(); |
321 | Float_t* candidatePhotonY = GetCandidatePhotonY(); | |
5cb4dfc3 | 322 | |
88dd9ad4 | 323 | Int_t candidatePhotonsNumber = GetCandidatePhotonsNumber(); |
5cb4dfc3 | 324 | |
88dd9ad4 | 325 | if(fDebug) cout << " n " << candidatePhotonsNumber << endl; |
5cb4dfc3 | 326 | |
327 | SetThetaCerenkov(999.); | |
328 | SetHoughPhotons(0); | |
329 | SetHoughPhotonsNorm(0); | |
330 | SetHoughRMS(999.); | |
331 | ||
88dd9ad4 | 332 | for (Int_t j=0; j < candidatePhotonsNumber; j++) |
5cb4dfc3 | 333 | { |
334 | ||
335 | SetPhotonIndex(j); | |
336 | ||
337 | SetPhotonFlag(0); | |
338 | SetPhotonEta(-999.); | |
339 | SetPhotonWeight(0.); | |
340 | ||
88dd9ad4 | 341 | if (j == iMipIndex) continue; |
5cb4dfc3 | 342 | |
5cb4dfc3 | 343 | |
88dd9ad4 | 344 | if(candidatePhotonX[j] < -64.) continue; /* avoid artificial clusters from edge uesd by Yale.... */ |
5cb4dfc3 | 345 | |
88dd9ad4 | 346 | Float_t xtoentr = candidatePhotonX[j] - shiftX; |
347 | Float_t ytoentr = candidatePhotonY[j] - shiftY; | |
5cb4dfc3 | 348 | |
349 | // Float_t chargehit = fHits_charge[j]; | |
350 | // if(chargehit > 150) continue; | |
351 | ||
88dd9ad4 | 352 | SetEntranceX(xtoentr); |
353 | SetEntranceY(ytoentr); | |
5cb4dfc3 | 354 | |
355 | FindPhiPoint(); | |
356 | ||
88dd9ad4 | 357 | Int_t photonStatus = PhotonInBand(); |
5cb4dfc3 | 358 | |
359 | if(fDebug) | |
360 | { | |
88dd9ad4 | 361 | cout << " Photon n. " << j << " Status " << photonStatus << " accepted " << endl; |
362 | cout << " CandidatePhotonX[j] " << candidatePhotonX[j] << " CandidatePhotonY[j] " << candidatePhotonY[j] << endl; | |
5cb4dfc3 | 363 | } |
364 | ||
88dd9ad4 | 365 | if(photonStatus == 0) continue; |
5cb4dfc3 | 366 | |
367 | SetPhotonFlag(1); | |
368 | ||
369 | FindThetaPhotonCerenkov(); | |
370 | ||
88dd9ad4 | 371 | Float_t thetaPhotonCerenkov = GetThetaPhotonCerenkov(); |
5cb4dfc3 | 372 | |
88dd9ad4 | 373 | if(fDebug) cout << " theta photon " << thetaPhotonCerenkov << endl; |
5cb4dfc3 | 374 | |
88dd9ad4 | 375 | SetPhotonEta(thetaPhotonCerenkov); |
5cb4dfc3 | 376 | |
88dd9ad4 | 377 | candidatePhotons++; |
5cb4dfc3 | 378 | |
5cb4dfc3 | 379 | |
380 | } | |
381 | ||
88dd9ad4 | 382 | if(candidatePhotons >= 1) kPatRec = kTRUE; |
5cb4dfc3 | 383 | |
384 | if(!kPatRec) return; | |
385 | { | |
386 | SetThetaCerenkov(999.); | |
387 | SetHoughPhotons(0); | |
388 | } | |
88dd9ad4 | 389 | SetPhotonsNumber(candidatePhotonsNumber); |
5cb4dfc3 | 390 | |
391 | HoughResponse(); | |
392 | ||
b068561d | 393 | fNrings++; |
5cb4dfc3 | 394 | |
395 | FlagPhotons(); | |
88dd9ad4 | 396 | Int_t nPhotonHough = GetHoughPhotons(); |
5cb4dfc3 | 397 | |
88dd9ad4 | 398 | if(nPhotonHough < 1) |
5cb4dfc3 | 399 | { |
400 | SetThetaCerenkov(999.); | |
401 | SetHoughPhotonsNorm(0.); | |
402 | return; | |
403 | } | |
404 | ||
b068561d | 405 | if(fIsWEIGHT) FindWeightThetaCerenkov(); |
5cb4dfc3 | 406 | |
88dd9ad4 | 407 | Float_t thetaCerenkov = GetThetaCerenkov(); |
5cb4dfc3 | 408 | |
88dd9ad4 | 409 | SetThetaOfRing(thetaCerenkov); |
5cb4dfc3 | 410 | FindAreaAndPortionOfRing(); |
411 | ||
88dd9ad4 | 412 | Float_t nPhotonHoughNorm = ((Float_t)nPhotonHough)/GetPortionOfRing(); |
413 | SetHoughPhotonsNorm(nPhotonHoughNorm); | |
5cb4dfc3 | 414 | |
415 | // Calculate the area where the photon are accepted... | |
416 | ||
88dd9ad4 | 417 | Float_t thetaInternal = thetaCerenkov - 0.5*fWindowWidth; |
418 | SetThetaOfRing(thetaInternal); | |
5cb4dfc3 | 419 | FindAreaAndPortionOfRing(); |
88dd9ad4 | 420 | Float_t internalArea = GetAreaOfRing(); |
5cb4dfc3 | 421 | |
88dd9ad4 | 422 | Float_t thetaExternal = thetaCerenkov + 0.5*fWindowWidth; |
423 | SetThetaOfRing(thetaExternal); | |
5cb4dfc3 | 424 | FindAreaAndPortionOfRing(); |
88dd9ad4 | 425 | Float_t externalArea = GetAreaOfRing(); |
5cb4dfc3 | 426 | |
88dd9ad4 | 427 | Float_t houghArea = externalArea - internalArea; |
5cb4dfc3 | 428 | |
88dd9ad4 | 429 | SetHoughArea(houghArea); |
5cb4dfc3 | 430 | |
431 | if(fDebug) | |
432 | { | |
433 | cout << " ----- SUMMARY OF RECONSTRUCTION ----- " << endl; | |
88dd9ad4 | 434 | cout << " Rings found " << fNrings << " with thetac " << thetaCerenkov << endl; |
5cb4dfc3 | 435 | |
5cb4dfc3 | 436 | |
437 | cout << " Nphotons " << GetPhotonsNumber() | |
88dd9ad4 | 438 | << " Hough " << nPhotonHough |
439 | << " norm " << nPhotonHoughNorm << endl; | |
5cb4dfc3 | 440 | |
88dd9ad4 | 441 | cout << " In PatRec:p " << pmod << " theta " << trackTheta << " phi " << trackPhi << endl; |
5cb4dfc3 | 442 | cout << " ------------------------------------- " << endl; |
443 | } | |
444 | ||
88dd9ad4 | 445 | Int_t nPhotons = GetPhotonsNumber(); |
5cb4dfc3 | 446 | |
447 | Float_t xmean = 0.; | |
448 | Float_t x2mean = 0.; | |
449 | Int_t nev = 0; | |
450 | ||
88dd9ad4 | 451 | for (Int_t j=0; j < nPhotons;j++) |
5cb4dfc3 | 452 | { |
453 | SetPhotonIndex(j); | |
454 | ||
455 | Float_t eta = GetPhotonEta(); | |
456 | ||
457 | if(eta != -999.) | |
458 | { | |
459 | if(GetPhotonFlag() == 2) | |
460 | { | |
461 | ||
5cb4dfc3 | 462 | |
463 | xmean += eta; | |
464 | x2mean += eta*eta; | |
465 | nev++; | |
466 | } | |
467 | } | |
468 | } | |
469 | ||
470 | if(nev > 0) | |
471 | { | |
472 | xmean /=(Float_t)nev; | |
473 | x2mean /=(Float_t)nev; | |
474 | } else { | |
475 | xmean = 0.; | |
476 | x2mean = 0.; | |
477 | } | |
478 | ||
88dd9ad4 | 479 | Float_t vRMS = sqrt(x2mean - xmean*xmean); |
5cb4dfc3 | 480 | |
88dd9ad4 | 481 | SetHoughRMS(vRMS); |
5cb4dfc3 | 482 | |
88dd9ad4 | 483 | if(fDebug) cout << " RMS " << vRMS << endl; |
5cb4dfc3 | 484 | |
485 | } | |
486 | ||
487 | void AliRICHRecon::FindEmissionPoint() | |
88dd9ad4 | 488 | { |
489 | //estimate the emission point in radiator | |
5cb4dfc3 | 490 | |
491 | // Find emission point | |
492 | ||
88dd9ad4 | 493 | Float_t absorbtionLenght=7.83*fRadiatorWidth; //absorption length in the freon (cm) |
5cb4dfc3 | 494 | // 7.83 = -1/ln(T0) where |
495 | // T0->Trasmission freon at 180nm = 0.88 (Eph=6.85eV) | |
496 | Float_t photonLenght, photonLenghtMin, photonLenghtMax; | |
497 | ||
88dd9ad4 | 498 | photonLenght=exp(-fRadiatorWidth/(absorbtionLenght*cos(fCerenkovAnglePad))); |
499 | photonLenghtMin=fRadiatorWidth*photonLenght/(1.-photonLenght); | |
5cb4dfc3 | 500 | photonLenghtMax=absorbtionLenght*cos(fCerenkovAnglePad); |
88dd9ad4 | 501 | Float_t emissionPoint = fRadiatorWidth + photonLenghtMin - photonLenghtMax; |
5cb4dfc3 | 502 | |
88dd9ad4 | 503 | SetEmissionPoint(emissionPoint); |
5cb4dfc3 | 504 | } |
505 | ||
506 | ||
507 | Int_t AliRICHRecon::PhotonInBand() | |
508 | { | |
88dd9ad4 | 509 | //search band fro photon candidates |
5cb4dfc3 | 510 | |
88dd9ad4 | 511 | // Float_t massOfParticle; |
5cb4dfc3 | 512 | Float_t beta; |
513 | Float_t nfreon; | |
514 | ||
515 | Float_t thetacer; | |
516 | ||
88dd9ad4 | 517 | Float_t xtoentr = GetEntranceX(); |
518 | Float_t ytoentr = GetEntranceY(); | |
5cb4dfc3 | 519 | |
88dd9ad4 | 520 | Float_t innerRadius; |
521 | Float_t outerRadius; | |
5cb4dfc3 | 522 | |
523 | Float_t phpad = GetPhiPoint(); | |
524 | ||
525 | // Float_t pmod = GetTrackMomentum(); | |
88dd9ad4 | 526 | // Float_t trackTheta = GetTrackTheta(); |
527 | // Float_t trackPhi = GetTrackPhi(); | |
5cb4dfc3 | 528 | |
529 | // inner radius // | |
530 | SetPhotonEnergy(5.6); | |
88dd9ad4 | 531 | SetEmissionPoint(fRadiatorWidth -0.0001); |
5cb4dfc3 | 532 | SetMassHypotesis(0.93828); |
533 | ||
534 | SetBetaOfParticle(); | |
535 | SetFreonRefractiveIndex(); | |
536 | ||
537 | beta = GetBetaOfParticle(); | |
538 | nfreon = GetFreonRefractiveIndex(); | |
539 | ||
540 | thetacer = Cerenkovangle(nfreon,beta); | |
541 | ||
542 | thetacer = 0.; | |
543 | ||
544 | if(fDebug) cout << " thetacer in photoninband min " << thetacer << endl; | |
545 | ||
546 | FindThetaAtQuartz(thetacer); | |
547 | ||
548 | if(thetacer == 999. || GetThetaAtQuartz() == 999.) | |
549 | { | |
88dd9ad4 | 550 | innerRadius = -999.; |
5cb4dfc3 | 551 | SetXInnerRing(-999.); |
552 | SetYInnerRing(-999.); | |
553 | SetRadiusInnerRing(-999.); | |
554 | } | |
555 | else | |
556 | { | |
557 | SetThetaPhotonInDRS(GetThetaAtQuartz()); | |
558 | SetPhiPhotonInDRS(phpad); | |
559 | ||
88dd9ad4 | 560 | innerRadius = FromEmissionToCathode(); |
561 | if(innerRadius == 999.) innerRadius = -999.; | |
5cb4dfc3 | 562 | |
563 | SetXInnerRing(GetXPointOnCathode()); | |
564 | SetYInnerRing(GetYPointOnCathode()); | |
88dd9ad4 | 565 | SetRadiusInnerRing(innerRadius); |
5cb4dfc3 | 566 | } |
567 | ||
568 | // outer radius // | |
569 | SetPhotonEnergy(7.7); | |
570 | SetEmissionPoint(0.); | |
571 | // SetMassHypotesis(0.139567); | |
572 | SetMassHypotesis(0.); | |
573 | ||
574 | SetBetaOfParticle(); | |
575 | SetFreonRefractiveIndex(); | |
576 | ||
577 | beta = GetBetaOfParticle(); | |
578 | nfreon = GetFreonRefractiveIndex(); | |
579 | ||
580 | thetacer = Cerenkovangle(nfreon,beta); | |
581 | ||
582 | // thetacer = 0.75; | |
583 | ||
584 | if(fDebug) cout << " thetacer in photoninband max " << thetacer << endl; | |
585 | ||
586 | FindThetaAtQuartz(thetacer); | |
587 | ||
588 | if(thetacer == 999. || GetThetaAtQuartz() == 999.) | |
589 | { | |
88dd9ad4 | 590 | outerRadius = 999.; |
5cb4dfc3 | 591 | SetXOuterRing(999.); |
592 | SetYOuterRing(999.); | |
593 | SetRadiusOuterRing(999.); | |
594 | } | |
595 | else | |
596 | { | |
597 | SetThetaPhotonInDRS(GetThetaAtQuartz()); | |
598 | SetPhiPhotonInDRS(phpad); | |
599 | ||
88dd9ad4 | 600 | outerRadius = FromEmissionToCathode(); |
601 | // cout << " outerRadius " << outerRadius << endl; | |
5cb4dfc3 | 602 | SetXOuterRing(GetXPointOnCathode()); |
603 | SetYOuterRing(GetYPointOnCathode()); | |
88dd9ad4 | 604 | SetRadiusOuterRing(outerRadius); |
5cb4dfc3 | 605 | } |
606 | ||
88dd9ad4 | 607 | Float_t padradius = sqrt(TMath::Power(xtoentr,2)+TMath::Power(ytoentr,2)); |
5cb4dfc3 | 608 | |
88dd9ad4 | 609 | if(fDebug) printf(" rmin %f r %f rmax %f \n",innerRadius,padradius,outerRadius); |
5cb4dfc3 | 610 | |
88dd9ad4 | 611 | if(padradius>=innerRadius && padradius<=outerRadius) return 1; |
5cb4dfc3 | 612 | return 0; |
613 | } | |
614 | ||
88dd9ad4 | 615 | void AliRICHRecon::FindThetaAtQuartz(Float_t thetaCerenkov) |
5cb4dfc3 | 616 | { |
88dd9ad4 | 617 | //find the theta at the quartz plate |
5cb4dfc3 | 618 | |
88dd9ad4 | 619 | if(thetaCerenkov == 999.) |
5cb4dfc3 | 620 | { |
621 | SetThetaAtQuartz(999.); | |
622 | return; | |
623 | } | |
624 | ||
88dd9ad4 | 625 | Float_t thetaAtQuartz = 999.; |
5cb4dfc3 | 626 | |
88dd9ad4 | 627 | Float_t trackTheta = GetTrackTheta(); |
5cb4dfc3 | 628 | |
88dd9ad4 | 629 | if(trackTheta == 0) { |
5cb4dfc3 | 630 | |
88dd9ad4 | 631 | if(fDebug) cout << " Theta sol unique " << thetaCerenkov << endl; |
5cb4dfc3 | 632 | |
88dd9ad4 | 633 | thetaAtQuartz = thetaCerenkov; |
634 | SetThetaAtQuartz(thetaAtQuartz); | |
5cb4dfc3 | 635 | return; |
636 | } | |
637 | ||
88dd9ad4 | 638 | Float_t trackPhi = GetTrackPhi(); |
639 | Float_t phiPoint = GetPhiPoint(); | |
5cb4dfc3 | 640 | |
88dd9ad4 | 641 | Double_t den = TMath::Sin((Double_t)trackTheta) |
642 | *TMath::Cos((Double_t)trackPhi) | |
643 | *TMath::Cos((Double_t)phiPoint) + | |
644 | TMath::Sin((Double_t)trackTheta) | |
645 | *TMath::Sin((Double_t)trackPhi) | |
646 | *TMath::Sin((Double_t)phiPoint); | |
647 | Double_t b = TMath::Cos((Double_t)trackTheta)/den; | |
648 | Double_t c = -TMath::Cos((Double_t)thetaCerenkov)/den; | |
5cb4dfc3 | 649 | |
88dd9ad4 | 650 | Double_t underSqrt = 1 + b*b - c*c; |
5cb4dfc3 | 651 | |
652 | if(fDebug) | |
653 | { | |
88dd9ad4 | 654 | cout << " trackTheta " << trackTheta << endl; |
655 | cout << " TrackPhi " << trackPhi << endl; | |
656 | cout << " PhiPoint " << phiPoint << endl; | |
657 | cout << " ThetaCerenkov " << thetaCerenkov << endl; | |
5cb4dfc3 | 658 | cout << " den b c " << den << " b " << b << " c " << c << endl; |
659 | } | |
660 | ||
88dd9ad4 | 661 | if(underSqrt < 0) { |
662 | if(fDebug) cout << " sqrt negative !!!!" << underSqrt << endl; | |
5cb4dfc3 | 663 | SetThetaAtQuartz(999.); |
664 | return; | |
665 | } | |
666 | ||
88dd9ad4 | 667 | Double_t sol1 = (1+TMath::Sqrt(underSqrt))/(b-c); |
668 | Double_t sol2 = (1-TMath::Sqrt(underSqrt))/(b-c); | |
5cb4dfc3 | 669 | |
88dd9ad4 | 670 | Double_t thetaSol1 = 2*TMath::ATan(sol1); |
671 | Double_t thetaSol2 = 2*TMath::ATan(sol2); | |
5cb4dfc3 | 672 | |
88dd9ad4 | 673 | if(fDebug) cout << " Theta sol 1 " << thetaSol1 |
674 | << " Theta sol 2 " << thetaSol2 << endl; | |
5cb4dfc3 | 675 | |
88dd9ad4 | 676 | if(thetaSol1>0 && thetaSol1 < TMath::Pi()) thetaAtQuartz = (Float_t)thetaSol1; |
677 | if(thetaSol2>0 && thetaSol2 < TMath::Pi()) thetaAtQuartz = (Float_t)thetaSol2; | |
5cb4dfc3 | 678 | |
88dd9ad4 | 679 | SetThetaAtQuartz(thetaAtQuartz); |
5cb4dfc3 | 680 | } |
681 | ||
682 | void AliRICHRecon::FindThetaPhotonCerenkov() | |
683 | { | |
88dd9ad4 | 684 | //find theta cerenkov of ring |
5cb4dfc3 | 685 | |
88dd9ad4 | 686 | Float_t thetaCerMin = 0.; |
687 | Float_t thetaCerMax = 0.75; | |
688 | Float_t thetaCerMean; | |
5cb4dfc3 | 689 | |
88dd9ad4 | 690 | Float_t radiusMin, radiusMax, radiusMean; |
5cb4dfc3 | 691 | Int_t nIteration = 0; |
692 | ||
88dd9ad4 | 693 | const Float_t kTollerance = 0.05; |
5cb4dfc3 | 694 | |
695 | // Float_t pmod = GetTrackMomentum(); | |
88dd9ad4 | 696 | // Float_t trackTheta = GetTrackTheta(); |
697 | // Float_t trackPhi = GetTrackPhi(); | |
5cb4dfc3 | 698 | |
88dd9ad4 | 699 | Float_t phiPoint = GetPhiPoint(); |
5cb4dfc3 | 700 | |
701 | SetPhotonEnergy(6.85); | |
88dd9ad4 | 702 | SetEmissionPoint(fRadiatorWidth/2); |
5cb4dfc3 | 703 | |
88dd9ad4 | 704 | Float_t xPoint = GetEntranceX(); |
705 | Float_t yPoint = GetEntranceY(); | |
706 | Float_t distPoint = sqrt(xPoint*xPoint + yPoint*yPoint); | |
5cb4dfc3 | 707 | |
88dd9ad4 | 708 | if(fDebug) cout << " DistPoint " << distPoint << endl; |
5cb4dfc3 | 709 | |
710 | // Star minimization... | |
711 | ||
712 | // First value... | |
713 | ||
88dd9ad4 | 714 | FindThetaAtQuartz(thetaCerMin); |
5cb4dfc3 | 715 | |
716 | if(GetThetaAtQuartz() == 999.) | |
717 | { | |
88dd9ad4 | 718 | radiusMin = -999.; |
5cb4dfc3 | 719 | } |
720 | else | |
721 | { | |
722 | SetThetaPhotonInDRS(GetThetaAtQuartz()); | |
88dd9ad4 | 723 | SetPhiPhotonInDRS(phiPoint); |
5cb4dfc3 | 724 | |
88dd9ad4 | 725 | radiusMin = FromEmissionToCathode(); |
5cb4dfc3 | 726 | } |
727 | ||
728 | // Second value... | |
729 | ||
88dd9ad4 | 730 | FindThetaAtQuartz(thetaCerMax); |
5cb4dfc3 | 731 | if(GetThetaAtQuartz() == 999.) |
732 | { | |
88dd9ad4 | 733 | radiusMax = 999.; |
5cb4dfc3 | 734 | } |
735 | else | |
736 | { | |
737 | SetThetaPhotonInDRS(GetThetaAtQuartz()); | |
88dd9ad4 | 738 | SetPhiPhotonInDRS(phiPoint); |
5cb4dfc3 | 739 | |
88dd9ad4 | 740 | radiusMax = FromEmissionToCathode(); |
5cb4dfc3 | 741 | } |
742 | // Mean value... | |
743 | ||
88dd9ad4 | 744 | thetaCerMean = (thetaCerMax + thetaCerMin)/2; |
5cb4dfc3 | 745 | |
88dd9ad4 | 746 | FindThetaAtQuartz(thetaCerMean); |
5cb4dfc3 | 747 | if(GetThetaAtQuartz() == 999.) |
748 | { | |
88dd9ad4 | 749 | radiusMean = 999.; |
5cb4dfc3 | 750 | } |
751 | else | |
752 | { | |
753 | SetThetaPhotonInDRS(GetThetaAtQuartz()); | |
88dd9ad4 | 754 | SetPhiPhotonInDRS(phiPoint); |
5cb4dfc3 | 755 | |
88dd9ad4 | 756 | radiusMean = FromEmissionToCathode(); |
5cb4dfc3 | 757 | } |
758 | ||
88dd9ad4 | 759 | if(fDebug) cout << " r1 " << radiusMin << " rmean " |
760 | << radiusMean << " r2 " << radiusMax << endl; | |
5cb4dfc3 | 761 | |
88dd9ad4 | 762 | while (TMath::Abs(radiusMean-distPoint) > kTollerance) |
5cb4dfc3 | 763 | { |
764 | ||
88dd9ad4 | 765 | if((radiusMin-distPoint)*(radiusMean-distPoint) < 0) thetaCerMax = thetaCerMean; |
766 | if((radiusMin-distPoint)*(radiusMean-distPoint) > 0) { | |
5cb4dfc3 | 767 | |
88dd9ad4 | 768 | thetaCerMin = thetaCerMean; |
5cb4dfc3 | 769 | |
88dd9ad4 | 770 | FindThetaAtQuartz(thetaCerMin); |
5cb4dfc3 | 771 | SetThetaPhotonInDRS(GetThetaAtQuartz()); |
88dd9ad4 | 772 | SetPhiPhotonInDRS(phiPoint); |
5cb4dfc3 | 773 | |
88dd9ad4 | 774 | radiusMin =FromEmissionToCathode(); |
5cb4dfc3 | 775 | } |
776 | ||
88dd9ad4 | 777 | thetaCerMean = (thetaCerMax + thetaCerMin)/2; |
5cb4dfc3 | 778 | |
88dd9ad4 | 779 | FindThetaAtQuartz(thetaCerMean); |
5cb4dfc3 | 780 | SetThetaPhotonInDRS(GetThetaAtQuartz()); |
88dd9ad4 | 781 | SetPhiPhotonInDRS(phiPoint); |
5cb4dfc3 | 782 | |
88dd9ad4 | 783 | radiusMean = FromEmissionToCathode(); |
5cb4dfc3 | 784 | |
785 | nIteration++; | |
786 | if(nIteration>=50) { | |
787 | if(fDebug) printf(" max iterations in FindPhotonCerenkov\n"); | |
788 | SetThetaPhotonCerenkov(999.); | |
789 | return; | |
790 | } | |
791 | } | |
792 | ||
88dd9ad4 | 793 | SetThetaPhotonCerenkov(thetaCerMean); |
5cb4dfc3 | 794 | |
795 | } | |
796 | ||
797 | void AliRICHRecon::FindAreaAndPortionOfRing() | |
798 | { | |
88dd9ad4 | 799 | //find fraction of the ring accepted by the RICH |
5cb4dfc3 | 800 | |
88dd9ad4 | 801 | Float_t xPoint[NPointsOfRing], yPoint[NPointsOfRing]; |
5cb4dfc3 | 802 | |
88dd9ad4 | 803 | // Float_t xtoentr = GetEntranceX(); |
804 | // Float_t ytoentr = GetEntranceY(); | |
805 | Float_t shiftX = GetShiftX(); | |
806 | Float_t shiftY = GetShiftY(); | |
5cb4dfc3 | 807 | |
88dd9ad4 | 808 | Float_t xemiss = GetXCoordOfEmission(); |
809 | Float_t yemiss = GetYCoordOfEmission(); | |
5cb4dfc3 | 810 | |
88dd9ad4 | 811 | Float_t x0 = xemiss + shiftX; |
812 | Float_t y0 = yemiss + shiftY; | |
5cb4dfc3 | 813 | |
814 | // Float_t pmod = GetTrackMomentum(); | |
88dd9ad4 | 815 | // Float_t trackTheta = GetTrackTheta(); |
816 | // Float_t trackPhi = GetTrackPhi(); | |
5cb4dfc3 | 817 | |
818 | SetPhotonEnergy(6.85); | |
819 | SetFreonRefractiveIndex(); | |
820 | ||
88dd9ad4 | 821 | SetEmissionPoint(fRadiatorWidth/2.); |
5cb4dfc3 | 822 | |
88dd9ad4 | 823 | Float_t theta = GetThetaOfRing(); |
5cb4dfc3 | 824 | |
825 | Int_t nPoints = 0; | |
88dd9ad4 | 826 | Int_t nPsiAccepted = 0; |
827 | Int_t nPsiTotal = 0; | |
5cb4dfc3 | 828 | |
829 | for(Int_t i=0;i<NPointsOfRing-1;i++) | |
830 | { | |
831 | ||
88dd9ad4 | 832 | Float_t psi = 2*TMath::Pi()*i/NPointsOfRing; |
5cb4dfc3 | 833 | |
88dd9ad4 | 834 | SetThetaPhotonInTRS(theta); |
835 | SetPhiPhotonInTRS(psi); | |
5cb4dfc3 | 836 | FindPhotonAnglesInDRS(); |
837 | ||
88dd9ad4 | 838 | Float_t radius = FromEmissionToCathode(); |
839 | if (radius == 999.) continue; | |
5cb4dfc3 | 840 | |
88dd9ad4 | 841 | nPsiTotal++; |
5cb4dfc3 | 842 | |
88dd9ad4 | 843 | Float_t xPointRing = GetXPointOnCathode() + shiftX; |
844 | Float_t yPointRing = GetYPointOnCathode() + shiftY; | |
5cb4dfc3 | 845 | |
88dd9ad4 | 846 | SetDetectorWhereX(xPointRing); |
847 | SetDetectorWhereY(yPointRing); | |
5cb4dfc3 | 848 | |
88dd9ad4 | 849 | Int_t zone = CheckDetectorAcceptance(); |
5cb4dfc3 | 850 | |
5cb4dfc3 | 851 | |
88dd9ad4 | 852 | if (zone != 0) |
5cb4dfc3 | 853 | { |
854 | FindIntersectionWithDetector(); | |
88dd9ad4 | 855 | xPoint[nPoints] = GetIntersectionX(); |
856 | yPoint[nPoints] = GetIntersectionY(); | |
5cb4dfc3 | 857 | } |
858 | else | |
859 | { | |
88dd9ad4 | 860 | xPoint[nPoints] = xPointRing; |
861 | yPoint[nPoints] = yPointRing; | |
862 | nPsiAccepted++; | |
5cb4dfc3 | 863 | } |
864 | ||
865 | nPoints++; | |
866 | ||
867 | } | |
868 | ||
88dd9ad4 | 869 | xPoint[nPoints] = xPoint[0]; |
870 | yPoint[nPoints] = yPoint[0]; | |
5cb4dfc3 | 871 | |
872 | // find area... | |
873 | ||
88dd9ad4 | 874 | Float_t area = 0; |
5cb4dfc3 | 875 | |
876 | for (Int_t i = 0; i < nPoints; i++) | |
877 | { | |
88dd9ad4 | 878 | area += TMath::Abs((xPoint[i]-x0)*(yPoint[i+1]-y0) - (xPoint[i+1]-x0)*(yPoint[i]-y0)); |
5cb4dfc3 | 879 | } |
880 | ||
88dd9ad4 | 881 | area *= 0.5; |
5cb4dfc3 | 882 | |
88dd9ad4 | 883 | Float_t portionOfRing = ((Float_t)nPsiAccepted)/((Float_t)(nPsiTotal)); |
5cb4dfc3 | 884 | |
5cb4dfc3 | 885 | |
88dd9ad4 | 886 | SetAreaOfRing(area); |
887 | SetPortionOfRing(portionOfRing); | |
5cb4dfc3 | 888 | } |
889 | ||
890 | void AliRICHRecon::FindIntersectionWithDetector() | |
891 | { | |
88dd9ad4 | 892 | // find ring intersection with CsI edges |
5cb4dfc3 | 893 | |
88dd9ad4 | 894 | Float_t xIntersect, yIntersect; |
5cb4dfc3 | 895 | Float_t x1, x2, y1, y2; |
896 | ||
88dd9ad4 | 897 | Float_t shiftX = GetShiftX(); |
898 | Float_t shiftY = GetShiftY(); | |
5cb4dfc3 | 899 | |
88dd9ad4 | 900 | Float_t xPoint = GetXPointOnCathode() + shiftX; |
901 | Float_t yPoint = GetYPointOnCathode() + shiftY; | |
5cb4dfc3 | 902 | |
88dd9ad4 | 903 | Float_t xemiss = GetXCoordOfEmission(); |
904 | Float_t yemiss = GetYCoordOfEmission(); | |
5cb4dfc3 | 905 | |
88dd9ad4 | 906 | Float_t phi = GetPhiPhotonInDRS(); |
907 | Float_t m = tan(phi); | |
5cb4dfc3 | 908 | |
88dd9ad4 | 909 | Float_t x0 = xemiss + shiftX; |
910 | Float_t y0 = yemiss + shiftY; | |
5cb4dfc3 | 911 | |
88dd9ad4 | 912 | if(xPoint > x0) |
5cb4dfc3 | 913 | { |
914 | x1 = x0; | |
88dd9ad4 | 915 | x2 = xPoint; |
5cb4dfc3 | 916 | } |
917 | else | |
918 | { | |
919 | x2 = x0; | |
88dd9ad4 | 920 | x1 = xPoint; |
5cb4dfc3 | 921 | } |
88dd9ad4 | 922 | if(yPoint > y0) |
5cb4dfc3 | 923 | { |
924 | y1 = y0; | |
88dd9ad4 | 925 | y2 = yPoint; |
5cb4dfc3 | 926 | } |
927 | else | |
928 | { | |
929 | y2 = y0; | |
88dd9ad4 | 930 | y1 = yPoint; |
5cb4dfc3 | 931 | } |
932 | // | |
88dd9ad4 | 933 | xIntersect = fXmax; |
934 | yIntersect = m*(xIntersect - x0) + y0; | |
935 | if (yIntersect >= fYmin && yIntersect <= fYmax && xIntersect >= x1 && xIntersect <= x2) | |
5cb4dfc3 | 936 | { |
88dd9ad4 | 937 | SetIntersectionX(xIntersect); |
938 | SetIntersectionY(yIntersect); | |
5cb4dfc3 | 939 | return; |
940 | } | |
941 | // | |
88dd9ad4 | 942 | xIntersect = fXmin; |
943 | yIntersect = m*(xIntersect - x0) + y0; | |
944 | if (yIntersect >= fYmin && yIntersect <= fYmax && xIntersect >= x1 && xIntersect <= x2) | |
5cb4dfc3 | 945 | { |
88dd9ad4 | 946 | SetIntersectionX(xIntersect); |
947 | SetIntersectionY(yIntersect); | |
5cb4dfc3 | 948 | return; |
949 | } | |
950 | // | |
88dd9ad4 | 951 | yIntersect = fYmax; |
952 | xIntersect = (yIntersect - y0)/m + x0; | |
953 | if (xIntersect >= fXmin && xIntersect <= fXmax && yIntersect >= y1 && yIntersect <= y2) | |
5cb4dfc3 | 954 | { |
88dd9ad4 | 955 | SetIntersectionX(xIntersect); |
956 | SetIntersectionY(yIntersect); | |
5cb4dfc3 | 957 | return; |
958 | } | |
959 | // | |
88dd9ad4 | 960 | yIntersect = fYmin; |
961 | xIntersect = (yIntersect - y0)/m + x0; | |
962 | if (xIntersect >= fXmin && xIntersect <= fXmax && yIntersect >= y1 && yIntersect <= y2) | |
5cb4dfc3 | 963 | { |
88dd9ad4 | 964 | SetIntersectionX(xIntersect); |
965 | SetIntersectionY(yIntersect); | |
5cb4dfc3 | 966 | return; |
967 | } | |
968 | ||
969 | cout << " sono fuori!!!!!!" << endl; | |
5cb4dfc3 | 970 | |
971 | } | |
b068561d | 972 | //__________________________________________________________________________________________________ |
88dd9ad4 | 973 | Int_t AliRICHRecon::CheckDetectorAcceptance() const |
5cb4dfc3 | 974 | { |
88dd9ad4 | 975 | // check for the acceptance |
5cb4dfc3 | 976 | |
977 | // crosses X -2.6 2.6 cm | |
978 | // crosses Y -1 1 cm | |
979 | ||
88dd9ad4 | 980 | Float_t xcoord = GetDetectorWhereX(); |
981 | Float_t ycoord = GetDetectorWhereY(); | |
5cb4dfc3 | 982 | |
88dd9ad4 | 983 | if(xcoord > fXmax) |
5cb4dfc3 | 984 | { |
88dd9ad4 | 985 | if(ycoord > fYmax) return 2; |
986 | if(ycoord > fYmin && ycoord < fYmax) return 3; | |
987 | if(ycoord < fYmin) return 4; | |
5cb4dfc3 | 988 | } |
88dd9ad4 | 989 | if(xcoord < fXmin) |
5cb4dfc3 | 990 | { |
88dd9ad4 | 991 | if(ycoord > fYmax) return 8; |
992 | if(ycoord > fYmin && ycoord < fYmax) return 7; | |
993 | if(ycoord < fYmin) return 6; | |
5cb4dfc3 | 994 | } |
88dd9ad4 | 995 | if(xcoord > fXmin && xcoord < fXmax) |
5cb4dfc3 | 996 | { |
88dd9ad4 | 997 | if(ycoord > fYmax) return 1; |
998 | if(ycoord > fYmin && ycoord < fYmax) return 0; | |
999 | if(ycoord < fYmin) return 5; | |
5cb4dfc3 | 1000 | } |
1001 | return 999; | |
1002 | } | |
b068561d | 1003 | //__________________________________________________________________________________________________ |
5cb4dfc3 | 1004 | Float_t AliRICHRecon::PhotonPositionOnCathode() |
1005 | { | |
88dd9ad4 | 1006 | // find the photon position on the CsI |
1007 | // Float_t massOfParticle; | |
5cb4dfc3 | 1008 | Float_t beta; |
1009 | Float_t nfreon; | |
1010 | ||
5cb4dfc3 | 1011 | |
1012 | SetPhotonEnergy(6.85); | |
88dd9ad4 | 1013 | SetEmissionPoint(fRadiatorWidth/2.); |
5cb4dfc3 | 1014 | SetMassHypotesis(0.139567); |
1015 | ||
1016 | SetBetaOfParticle(); | |
1017 | SetFreonRefractiveIndex(); | |
1018 | ||
1019 | beta = GetBetaOfParticle(); | |
1020 | nfreon = GetFreonRefractiveIndex(); | |
1021 | ||
5cb4dfc3 | 1022 | |
88dd9ad4 | 1023 | Float_t radius = FromEmissionToCathode(); |
1024 | if (radius == 999.) return 999.; | |
5cb4dfc3 | 1025 | |
5cb4dfc3 | 1026 | return 0; |
1027 | } | |
1028 | ||
1029 | void AliRICHRecon::FindPhotonAnglesInDRS() | |
1030 | { | |
1031 | // Setup the rotation matrix of the track... | |
1032 | ||
88dd9ad4 | 1033 | TRotation mtheta; |
1034 | TRotation mphi; | |
1035 | TRotation minv; | |
1036 | TRotation mrot; | |
5cb4dfc3 | 1037 | |
88dd9ad4 | 1038 | Float_t trackTheta = GetTrackTheta(); |
1039 | Float_t trackPhi = GetTrackPhi(); | |
5cb4dfc3 | 1040 | |
88dd9ad4 | 1041 | mtheta.RotateY(trackTheta); |
1042 | mphi.RotateZ(trackPhi); | |
5cb4dfc3 | 1043 | |
88dd9ad4 | 1044 | mrot = mphi * mtheta; |
1045 | // minv = mrot.Inverse(); | |
5cb4dfc3 | 1046 | |
88dd9ad4 | 1047 | TVector3 photonInRadiator(1,1,1); |
5cb4dfc3 | 1048 | |
88dd9ad4 | 1049 | Float_t thetaCerenkov = GetThetaPhotonInTRS(); |
1050 | Float_t phiCerenkov = GetPhiPhotonInTRS(); | |
5cb4dfc3 | 1051 | |
88dd9ad4 | 1052 | photonInRadiator.SetTheta(thetaCerenkov); |
1053 | photonInRadiator.SetPhi(phiCerenkov); | |
1054 | photonInRadiator = mrot * photonInRadiator; | |
1055 | Float_t theta = photonInRadiator.Theta(); | |
1056 | Float_t phi = photonInRadiator.Phi(); | |
1057 | SetThetaPhotonInDRS(theta); | |
1058 | SetPhiPhotonInDRS(phi); | |
5cb4dfc3 | 1059 | |
1060 | } | |
1061 | ||
1062 | Float_t AliRICHRecon::FromEmissionToCathode() | |
1063 | { | |
88dd9ad4 | 1064 | // trace from emission point to cathode |
5cb4dfc3 | 1065 | |
1066 | Float_t nfreon, nquartz, ngas; | |
1067 | ||
1068 | SetFreonRefractiveIndex(); | |
1069 | SetQuartzRefractiveIndex(); | |
1070 | SetGasRefractiveIndex(); | |
1071 | ||
1072 | nfreon = GetFreonRefractiveIndex(); | |
1073 | nquartz = GetQuartzRefractiveIndex(); | |
1074 | ngas = GetGasRefractiveIndex(); | |
1075 | ||
88dd9ad4 | 1076 | Float_t trackTheta = GetTrackTheta(); |
1077 | Float_t trackPhi = GetTrackPhi(); | |
1078 | Float_t lengthOfEmissionPoint = GetEmissionPoint(); | |
5cb4dfc3 | 1079 | |
88dd9ad4 | 1080 | Float_t theta = GetThetaPhotonInDRS(); |
1081 | Float_t phi = GetPhiPhotonInDRS(); | |
5cb4dfc3 | 1082 | |
1083 | // cout << " Theta " << Theta << " Phi " << Phi << endl; | |
1084 | ||
88dd9ad4 | 1085 | Float_t xemiss = lengthOfEmissionPoint*tan(trackTheta)*cos(trackPhi); |
1086 | Float_t yemiss = lengthOfEmissionPoint*tan(trackTheta)*sin(trackPhi); | |
5cb4dfc3 | 1087 | |
88dd9ad4 | 1088 | SetXCoordOfEmission(xemiss); |
1089 | SetYCoordOfEmission(yemiss); | |
5cb4dfc3 | 1090 | |
88dd9ad4 | 1091 | Float_t thetaquar = SnellAngle(nfreon, nquartz, theta); |
5cb4dfc3 | 1092 | |
1093 | if(thetaquar == 999.) | |
1094 | { | |
1095 | SetXPointOnCathode(999.); | |
1096 | SetYPointOnCathode(999.); | |
1097 | return thetaquar; | |
1098 | } | |
1099 | ||
1100 | Float_t thetagap = SnellAngle( nquartz, ngas, thetaquar); | |
1101 | ||
1102 | if(thetagap == 999.) | |
1103 | { | |
1104 | SetXPointOnCathode(999.); | |
1105 | SetYPointOnCathode(999.); | |
1106 | return thetagap; | |
1107 | } | |
1108 | ||
88dd9ad4 | 1109 | Float_t xw = (fRadiatorWidth - lengthOfEmissionPoint)*cos(phi)*tan(theta); |
1110 | Float_t xq = fQuartzWidth*cos(phi)*tan(thetaquar); | |
1111 | Float_t xg = fGapWidth*cos(phi)*tan(thetagap); | |
1112 | Float_t yw = (fRadiatorWidth - lengthOfEmissionPoint)*sin(phi)*tan(theta); | |
1113 | Float_t yq = fQuartzWidth*sin(phi)*tan(thetaquar); | |
1114 | Float_t yg = fGapWidth*sin(phi)*tan(thetagap); | |
5cb4dfc3 | 1115 | |
5cb4dfc3 | 1116 | |
88dd9ad4 | 1117 | Float_t xtot = xemiss + xw + xq + xg; |
1118 | Float_t ytot = yemiss + yw + yq + yg; | |
5cb4dfc3 | 1119 | |
1120 | SetXPointOnCathode(xtot); | |
1121 | SetYPointOnCathode(ytot); | |
1122 | ||
5cb4dfc3 | 1123 | |
88dd9ad4 | 1124 | Float_t distanceFromEntrance = sqrt(TMath::Power(fPhotonLimitX,2) |
5cb4dfc3 | 1125 | +TMath::Power(fPhotonLimitY,2)); |
1126 | ||
88dd9ad4 | 1127 | return distanceFromEntrance; |
5cb4dfc3 | 1128 | |
1129 | } | |
1130 | ||
1131 | ||
1132 | void AliRICHRecon::FindPhiPoint() | |
1133 | { | |
88dd9ad4 | 1134 | //find phi of generated point |
5cb4dfc3 | 1135 | |
88dd9ad4 | 1136 | Float_t xtoentr = GetEntranceX(); |
1137 | Float_t ytoentr = GetEntranceY(); | |
5cb4dfc3 | 1138 | |
88dd9ad4 | 1139 | Float_t trackTheta = GetTrackTheta(); |
1140 | Float_t trackPhi = GetTrackPhi(); | |
5cb4dfc3 | 1141 | |
88dd9ad4 | 1142 | Float_t emissionPoint = GetEmissionPoint(); |
5cb4dfc3 | 1143 | |
88dd9ad4 | 1144 | Float_t argY = ytoentr - emissionPoint*tan(trackTheta)*sin(trackPhi); |
1145 | Float_t argX = xtoentr - emissionPoint*tan(trackTheta)*cos(trackPhi); | |
5cb4dfc3 | 1146 | Float_t phipad = atan2(argY,argX); |
1147 | ||
1148 | SetPhiPoint(phipad); | |
1149 | ||
1150 | } | |
1151 | ||
1152 | Float_t AliRICHRecon::Cerenkovangle(Float_t n, Float_t beta) | |
1153 | { | |
88dd9ad4 | 1154 | // cerenkov angle from n and beta |
5cb4dfc3 | 1155 | |
1156 | // Compute the cerenkov angle | |
1157 | ||
1158 | Float_t thetacer; | |
1159 | ||
1160 | if((n*beta)<1.) { | |
1161 | thetacer = 999.; | |
1162 | // cout << " warning in Cerenkoangle !!!!!! " << endl; | |
1163 | return thetacer; | |
1164 | } | |
1165 | ||
1166 | thetacer = acos (1./(n*beta)); | |
1167 | return thetacer; | |
1168 | } | |
1169 | ||
1170 | Float_t AliRICHRecon::SnellAngle(Float_t n1, Float_t n2, Float_t theta1) | |
88dd9ad4 | 1171 | { |
1172 | // Snell law | |
5cb4dfc3 | 1173 | |
1174 | // Compute the Snell angle | |
1175 | ||
1176 | Float_t sinrefractangle; | |
1177 | Float_t refractangle; | |
1178 | ||
1179 | sinrefractangle = (n1/n2)*sin(theta1); | |
1180 | ||
1181 | if(sinrefractangle>1.) { | |
1182 | // cout << " PROBLEMS IN SNELL ANGLE !!!!! " << endl; | |
1183 | refractangle = 999.; | |
1184 | return refractangle; | |
1185 | } | |
1186 | ||
1187 | refractangle = asin(sinrefractangle); | |
1188 | return refractangle; | |
1189 | } | |
1190 | ||
1191 | ||
1192 | void AliRICHRecon::HoughResponse() | |
88dd9ad4 | 1193 | { |
1194 | //Hough response | |
5cb4dfc3 | 1195 | |
1196 | // Implement Hough response pat. rec. method | |
1197 | ||
88dd9ad4 | 1198 | Float_t *hCSspace; |
5cb4dfc3 | 1199 | |
1200 | int bin=0; | |
1201 | int bin1=0; | |
1202 | int bin2=0; | |
1203 | int i, j, k, nCorrBand; | |
1204 | float hcs[750],hcsw[750]; | |
1205 | float angle, weight; | |
1206 | float lowerlimit,upperlimit; | |
1207 | ||
1208 | float etaPeak[100]; | |
1209 | ||
1210 | int nBin; | |
1211 | ||
1212 | float etaPeakPos = -1; | |
1213 | ||
1214 | Int_t etaPeakCount = -1; | |
1215 | ||
88dd9ad4 | 1216 | Float_t thetaCerenkov = 0.; |
5cb4dfc3 | 1217 | |
1218 | nBin = (int)(0.5+fThetaMax/(fDTheta)); | |
1219 | nCorrBand = (int)(0.5+ fWindowWidth/(2 * fDTheta)); | |
1220 | ||
1221 | memset ((void *)hcs, 0, fThetaBin*sizeof(float)); | |
1222 | memset ((void *)hcsw, 0, fThetaBin*sizeof(float)); | |
1223 | ||
88dd9ad4 | 1224 | Int_t nPhotons = GetPhotonsNumber(); |
5cb4dfc3 | 1225 | |
88dd9ad4 | 1226 | Int_t weightFlag = 0; |
5cb4dfc3 | 1227 | |
88dd9ad4 | 1228 | for (k=0; k< nPhotons; k++) { |
5cb4dfc3 | 1229 | |
1230 | SetPhotonIndex(k); | |
1231 | ||
1232 | angle = GetPhotonEta(); | |
1233 | ||
1234 | if(angle == -999.) continue; | |
1235 | ||
1236 | if (angle>=fThetaMin && angle<= fThetaMax) | |
1237 | ||
1238 | { | |
1239 | ||
1240 | bin = (int)(0.5+angle/(fDTheta)); | |
1241 | ||
1242 | bin1= bin-nCorrBand; | |
1243 | bin2= bin+nCorrBand; | |
1244 | ||
1245 | // calculate weights | |
1246 | ||
b068561d | 1247 | if(fIsWEIGHT) |
5cb4dfc3 | 1248 | { |
1249 | lowerlimit = ((Float_t)bin1)*fDTheta + 0.5*fDTheta; | |
1250 | SetThetaOfRing(lowerlimit); | |
1251 | FindAreaAndPortionOfRing(); | |
1252 | Float_t area1 = GetAreaOfRing(); | |
1253 | ||
1254 | upperlimit = ((Float_t)bin2)*fDTheta + 0.5*fDTheta; | |
1255 | SetThetaOfRing(upperlimit); | |
1256 | FindAreaAndPortionOfRing(); | |
1257 | Float_t area2 = GetAreaOfRing(); | |
1258 | ||
1259 | // cout << "lowerlimit" << lowerlimit << "upperlimit " << upperlimit << endl; | |
1260 | Float_t diffarea = area2 - area1; | |
1261 | ||
1262 | if(diffarea>0) | |
1263 | { | |
1264 | weight = 1./(area2-area1); | |
1265 | } | |
1266 | else | |
1267 | { | |
88dd9ad4 | 1268 | weightFlag = 1; |
5cb4dfc3 | 1269 | weight = 1.; |
1270 | } | |
1271 | ||
1272 | // cout <<" low "<< lowerlimit << " up " << upperlimit << | |
1273 | // " area1 " << area1 << " area2 " << area2 << " weight " << weight << endl; | |
1274 | ||
1275 | } | |
1276 | else | |
1277 | { | |
1278 | weight = 1.; | |
1279 | } | |
1280 | ||
1281 | SetPhotonWeight(weight); | |
1282 | ||
1283 | // cout << "weight..." << weight << endl; | |
1284 | ||
5cb4dfc3 | 1285 | |
1286 | if (bin1<0) bin1=0; | |
1287 | if (bin2>nBin) bin2=nBin; | |
1288 | ||
1289 | for (j=bin1; j<bin2; j++) | |
1290 | { | |
1291 | hcs[j] += 1; | |
1292 | hcsw[j] += weight; | |
1293 | } | |
1294 | } | |
1295 | } | |
1296 | ||
5cb4dfc3 | 1297 | |
88dd9ad4 | 1298 | if(weightFlag == 0) |
5cb4dfc3 | 1299 | { |
88dd9ad4 | 1300 | hCSspace = hcsw; |
5cb4dfc3 | 1301 | } |
1302 | else | |
1303 | { | |
88dd9ad4 | 1304 | hCSspace = hcs; |
5cb4dfc3 | 1305 | // cout << " probems with weight...normal procedure adopted " << endl; |
1306 | } | |
1307 | ||
88dd9ad4 | 1308 | HoughFiltering(hCSspace); |
5cb4dfc3 | 1309 | |
1310 | for (bin=0; bin <nBin; bin++) { | |
1311 | angle = (bin+0.5) * (fDTheta); | |
88dd9ad4 | 1312 | if (hCSspace[bin] && hCSspace[bin] > etaPeakPos) { |
5cb4dfc3 | 1313 | etaPeakCount = 0; |
88dd9ad4 | 1314 | etaPeakPos = hCSspace[bin]; |
5cb4dfc3 | 1315 | etaPeak[0]=angle; |
1316 | } | |
1317 | else { | |
88dd9ad4 | 1318 | if (hCSspace[bin] == etaPeakPos) { |
5cb4dfc3 | 1319 | etaPeak[++etaPeakCount] = angle; |
1320 | } | |
1321 | } | |
1322 | } | |
1323 | ||
1324 | for (i=0; i<etaPeakCount+1; i++) { | |
88dd9ad4 | 1325 | thetaCerenkov += etaPeak[i]; |
5cb4dfc3 | 1326 | } |
1327 | if (etaPeakCount>=0) { | |
88dd9ad4 | 1328 | thetaCerenkov /= etaPeakCount+1; |
5cb4dfc3 | 1329 | fThetaPeakPos = etaPeakPos; |
1330 | } | |
1331 | ||
88dd9ad4 | 1332 | SetThetaCerenkov(thetaCerenkov); |
5cb4dfc3 | 1333 | } |
1334 | ||
1335 | ||
1336 | void AliRICHRecon::HoughFiltering(float hcs[]) | |
1337 | { | |
88dd9ad4 | 1338 | // filter for Hough |
5cb4dfc3 | 1339 | |
1340 | // hough filtering | |
1341 | ||
1342 | float hcsFilt[750]; | |
1343 | float k[5] = {0.05, 0.25, 0.4, 0.25, 0.05}; | |
1344 | int nx, i, nxDx; | |
1345 | int sizeHCS; | |
1346 | int nBin; | |
1347 | ||
1348 | nBin = (int)(1+fThetaMax/fDTheta); | |
1349 | sizeHCS = fThetaBin*sizeof(float); | |
1350 | ||
1351 | memset ((void *)hcsFilt, 0, sizeHCS); | |
1352 | ||
1353 | for (nx = 0; nx < nBin; nx++) { | |
1354 | for (i = 0; i < 5; i++) { | |
1355 | nxDx = nx + (i-2); | |
1356 | if (nxDx> -1 && nxDx<nBin) | |
1357 | hcsFilt[nx] += hcs[nxDx] * k[i]; | |
1358 | } | |
1359 | } | |
1360 | ||
1361 | for (nx = 0; nx < nBin; nx++) { | |
1362 | hcs[nx] = hcsFilt[nx]; | |
1363 | } | |
1364 | } | |
1365 | ||
1366 | void AliRICHRecon::FindWeightThetaCerenkov() | |
1367 | { | |
88dd9ad4 | 1368 | // manage with weight for photons |
5cb4dfc3 | 1369 | |
1370 | Float_t wei = 0.; | |
88dd9ad4 | 1371 | Float_t weightThetaCerenkov = 0.; |
5cb4dfc3 | 1372 | |
88dd9ad4 | 1373 | Int_t nPhotons = GetPhotonsNumber(); |
1374 | for(Int_t i=0;i<nPhotons;i++) | |
5cb4dfc3 | 1375 | { |
1376 | SetPhotonIndex(i); | |
1377 | ||
1378 | if(GetPhotonFlag() == 2) | |
1379 | { | |
88dd9ad4 | 1380 | Float_t photonEta = GetPhotonEta(); |
1381 | Float_t photonWeight = GetPhotonWeight(); | |
1382 | weightThetaCerenkov += photonEta*photonWeight; | |
1383 | wei += photonWeight; | |
5cb4dfc3 | 1384 | } |
1385 | } | |
1386 | ||
1387 | if(wei != 0.) | |
1388 | { | |
88dd9ad4 | 1389 | weightThetaCerenkov /= wei; |
5cb4dfc3 | 1390 | } |
1391 | else | |
1392 | { | |
88dd9ad4 | 1393 | weightThetaCerenkov = 0.; |
5cb4dfc3 | 1394 | } |
1395 | ||
88dd9ad4 | 1396 | SetThetaCerenkov(weightThetaCerenkov); |
5cb4dfc3 | 1397 | |
88dd9ad4 | 1398 | cout << " thetac weighted -> " << weightThetaCerenkov << endl; |
5cb4dfc3 | 1399 | } |
1400 | ||
1401 | ||
1402 | void AliRICHRecon::FlagPhotons() | |
1403 | { | |
88dd9ad4 | 1404 | // flag photons |
5cb4dfc3 | 1405 | |
88dd9ad4 | 1406 | Int_t nPhotonHough = 0; |
5cb4dfc3 | 1407 | |
88dd9ad4 | 1408 | Float_t thetaCerenkov = GetThetaCerenkov(); |
1409 | if(fDebug) cout << " fThetaCerenkov " << thetaCerenkov << endl; | |
5cb4dfc3 | 1410 | |
88dd9ad4 | 1411 | Float_t thetaDist= thetaCerenkov - fThetaMin; |
1412 | Int_t steps = (Int_t)(thetaDist / fDTheta); | |
5cb4dfc3 | 1413 | |
1414 | Float_t tmin = fThetaMin + (Float_t)(steps - 1)*fDTheta; | |
1415 | Float_t tmax = fThetaMin + (Float_t)(steps)*fDTheta; | |
1416 | Float_t tavg = 0.5*(tmin+tmax); | |
1417 | ||
1418 | tmin = tavg - 0.5*fWindowWidth; | |
1419 | tmax = tavg + 0.5*fWindowWidth; | |
1420 | ||
1421 | if(fDebug) cout << " tmin " << tmin << " tmax " << tmax << endl; | |
88dd9ad4 | 1422 | if(fDebug) cout << " thetac " << thetaCerenkov << endl; |
5cb4dfc3 | 1423 | |
88dd9ad4 | 1424 | // Int_t candidatePhotonsNumber = GetCandidatePhotonsNumber(); |
5cb4dfc3 | 1425 | |
88dd9ad4 | 1426 | Int_t nPhotons = GetPhotonsNumber(); |
5cb4dfc3 | 1427 | |
88dd9ad4 | 1428 | // for(Int_t i=0;i<candidatePhotonsNumber;i++) |
5cb4dfc3 | 1429 | |
88dd9ad4 | 1430 | for(Int_t i=0;i<nPhotons;i++) |
5cb4dfc3 | 1431 | { |
1432 | SetPhotonIndex(i); | |
1433 | ||
88dd9ad4 | 1434 | Float_t photonEta = GetPhotonEta(); |
5cb4dfc3 | 1435 | |
88dd9ad4 | 1436 | if(photonEta == -999.) continue; |
5cb4dfc3 | 1437 | |
88dd9ad4 | 1438 | if(photonEta >= tmin && photonEta <= tmax) |
5cb4dfc3 | 1439 | { |
1440 | SetPhotonFlag(2); | |
88dd9ad4 | 1441 | nPhotonHough++; |
5cb4dfc3 | 1442 | } |
1443 | } | |
88dd9ad4 | 1444 | SetHoughPhotons(nPhotonHough); |
5cb4dfc3 | 1445 | } |
1446 | ||
88dd9ad4 | 1447 | void AliRICHRecon::DrawEvent(Int_t flag) const |
5cb4dfc3 | 1448 | { |
88dd9ad4 | 1449 | // draw event with rings |
5cb4dfc3 | 1450 | |
1451 | flag=1; // dummy to be removed... | |
5cb4dfc3 | 1452 | } |
1453 | ||
1454 | Float_t AliRICHRecon::FindMassOfParticle() | |
1455 | { | |
88dd9ad4 | 1456 | // find mass of the particle from theta cerenkov |
5cb4dfc3 | 1457 | |
1458 | Float_t pmod = GetTrackMomentum(); | |
1459 | ||
1460 | SetPhotonEnergy(6.85); | |
1461 | SetFreonRefractiveIndex(); | |
1462 | ||
88dd9ad4 | 1463 | Float_t thetaCerenkov = GetThetaCerenkov(); |
1464 | FindBetaFromTheta(thetaCerenkov); | |
5cb4dfc3 | 1465 | |
1466 | Double_t beta = (Double_t)(GetBetaOfParticle()); | |
1467 | Double_t den = 1. - beta*beta; | |
1468 | if(den<=0.) return 999.; | |
1469 | ||
1470 | Double_t gamma = 1./TMath::Sqrt(den); | |
1471 | ||
1472 | Float_t mass = pmod/(beta*(Float_t)gamma); | |
1473 | ||
1474 | return mass; | |
1475 | } | |
1476 | ||
1477 | ||
1478 | void AliRICHRecon::FillHistograms() | |
1479 | { | |
88dd9ad4 | 1480 | // fill histograms.. |
5cb4dfc3 | 1481 | |
88dd9ad4 | 1482 | Float_t fittedTrackTheta, fittedTrackPhi; |
5cb4dfc3 | 1483 | |
88dd9ad4 | 1484 | Float_t thetaCerenkov = GetThetaCerenkov(); |
1485 | if(thetaCerenkov == 999.) return; | |
5cb4dfc3 | 1486 | |
88dd9ad4 | 1487 | Float_t vertZ = GetEventVertexZ(); |
5cb4dfc3 | 1488 | |
88dd9ad4 | 1489 | Float_t trackTheta = GetTrackTheta(); |
1490 | Float_t trackPhi = GetTrackPhi(); | |
5cb4dfc3 | 1491 | Float_t pmod = GetTrackMomentum(); |
1492 | Float_t pt = GetTrackPt(); | |
88dd9ad4 | 1493 | Float_t trackEta = GetTrackEta(); |
5cb4dfc3 | 1494 | Int_t q = GetTrackCharge(); |
88dd9ad4 | 1495 | Float_t tPCLastZ = GetTrackTPCLastZ(); |
1496 | Float_t minDist = GetMinDist(); | |
5cb4dfc3 | 1497 | |
88dd9ad4 | 1498 | fittedTrackTheta = GetFittedTrackTheta(); |
1499 | fittedTrackPhi = GetFittedTrackPhi(); | |
1500 | Int_t fittednPhotonHough = GetFittedHoughPhotons(); | |
5cb4dfc3 | 1501 | |
1502 | if(fDebug) | |
1503 | { | |
88dd9ad4 | 1504 | cout << " p " << pmod << " ThetaC " << thetaCerenkov |
b068561d | 1505 | << " rings " << fNrings << endl; |
5cb4dfc3 | 1506 | } |
1507 | ||
88dd9ad4 | 1508 | Int_t nPhotonHough = GetHoughPhotons(); |
1509 | // Float_t nPhotonHoughNorm = GetHoughPhotonsNorm(); | |
1510 | Float_t inRing = GetPortionOfRing(); | |
5cb4dfc3 | 1511 | |
88dd9ad4 | 1512 | Float_t massOfParticle = FindMassOfParticle(); |
5cb4dfc3 | 1513 | |
88dd9ad4 | 1514 | Float_t houghArea = GetHoughArea(); |
1515 | Float_t multiplicity = GetEventMultiplicity(); | |
5cb4dfc3 | 1516 | |
5cb4dfc3 | 1517 | |
1518 | Float_t var[20]; | |
1519 | ||
5cb4dfc3 | 1520 | var[0] = 0; |
1521 | var[1] = 0; | |
88dd9ad4 | 1522 | var[2] = vertZ; |
5cb4dfc3 | 1523 | var[3] = pmod; |
1524 | var[4] = pt; | |
88dd9ad4 | 1525 | var[5] = trackEta; |
1526 | var[6] = trackTheta; | |
1527 | var[7] = trackPhi; | |
1528 | var[8] = fittedTrackTheta; | |
1529 | var[9] = fittedTrackPhi; | |
5cb4dfc3 | 1530 | var[10] = q; |
88dd9ad4 | 1531 | var[11] = thetaCerenkov; |
1532 | var[12] = (Float_t)nPhotonHough; | |
1533 | var[13] = (Float_t)fittednPhotonHough; | |
1534 | var[14] = inRing; | |
1535 | var[15] = massOfParticle; | |
1536 | var[16] = houghArea; | |
1537 | var[17] = multiplicity; | |
1538 | var[18] = tPCLastZ; | |
1539 | var[19] = minDist; | |
5cb4dfc3 | 1540 | |
b068561d | 1541 | fNtuple->Fill(var); |
5cb4dfc3 | 1542 | |
5cb4dfc3 | 1543 | |
88dd9ad4 | 1544 | fittedTrackTheta = GetFittedTrackTheta(); |
1545 | fittedTrackPhi = GetFittedTrackPhi(); | |
5cb4dfc3 | 1546 | |
5cb4dfc3 | 1547 | |
5cb4dfc3 | 1548 | |
88dd9ad4 | 1549 | if(thetaCerenkov > 0.505 && thetaCerenkov < 0.605) { |
5cb4dfc3 | 1550 | SetPhotonEnergy(6.85); |
1551 | SetFreonRefractiveIndex(); | |
5cb4dfc3 | 1552 | } |
1553 | ||
88dd9ad4 | 1554 | Int_t nPhotons = GetPhotonsNumber(); |
5cb4dfc3 | 1555 | |
88dd9ad4 | 1556 | for (Int_t j=0; j < nPhotons;j++) |
b068561d | 1557 | SetPhotonIndex(j); |
1558 | }//FillHistograms() | |
1559 | //__________________________________________________________________________________________________ | |
5cb4dfc3 | 1560 | void AliRICHRecon::Minimization() |
1561 | { | |
88dd9ad4 | 1562 | // minimization to find the best theta and phi of the track |
5cb4dfc3 | 1563 | |
1564 | Double_t arglist; | |
1565 | Int_t ierflag = 0; | |
1566 | ||
1567 | static Double_t vstart[2]; | |
1568 | static Double_t lower[2], upper[2]; | |
1569 | static Double_t step[2]={0.001,0.001}; | |
1570 | ||
88dd9ad4 | 1571 | Double_t trackThetaNew,trackPhiNew; |
5cb4dfc3 | 1572 | TString chname; |
1573 | Double_t eps, b1, b2; | |
1574 | Int_t ierflg; | |
1575 | ||
88dd9ad4 | 1576 | gAliRICHminuit = new TMinuit(2); |
1577 | gAliRICHminuit->SetObjectFit((TObject *)this); | |
1578 | gAliRICHminuit->SetFCN(fcnrecon); | |
1579 | gAliRICHminuit->mninit(5,10,7); | |
5cb4dfc3 | 1580 | |
1581 | vstart[0] = (Double_t)GetTrackTheta(); | |
1582 | vstart[1] = (Double_t)GetTrackPhi(); | |
1583 | ||
1584 | lower[0] = vstart[0] - 0.03; | |
1585 | if(lower[0] < 0) lower[0] = 0.; | |
1586 | upper[0] = vstart[0] + 0.03; | |
1587 | lower[1] = vstart[1] - 0.03; | |
1588 | upper[1] = vstart[1] + 0.03; | |
1589 | ||
1590 | ||
88dd9ad4 | 1591 | gAliRICHminuit->mnparm(0,"theta",vstart[0],step[0],lower[0],upper[0],ierflag); |
1592 | gAliRICHminuit->mnparm(1," phi ",vstart[1],step[1],lower[1],upper[1],ierflag); | |
5cb4dfc3 | 1593 | |
1594 | arglist = -1; | |
1595 | ||
88dd9ad4 | 1596 | // gAliRICHminuit->FixParameter(0); |
5cb4dfc3 | 1597 | |
88dd9ad4 | 1598 | gAliRICHminuit->SetPrintLevel(-1); |
1599 | // gAliRICHminuit->mnexcm("SET PRI",&arglist, 1, ierflag); | |
1600 | gAliRICHminuit->mnexcm("SET NOGR",&arglist, 1, ierflag); | |
1601 | gAliRICHminuit->mnexcm("SET NOW",&arglist, 1, ierflag); | |
5cb4dfc3 | 1602 | arglist = 1; |
88dd9ad4 | 1603 | gAliRICHminuit->mnexcm("SET ERR", &arglist, 1,ierflg); |
5cb4dfc3 | 1604 | arglist = -1; |
1605 | ||
88dd9ad4 | 1606 | // gAliRICHminuit->mnscan(); |
5cb4dfc3 | 1607 | |
88dd9ad4 | 1608 | // gAliRICHminuit->mnexcm("SIMPLEX",&arglist, 0, ierflag); |
1609 | gAliRICHminuit->mnexcm("MIGRAD",&arglist, 0, ierflag); | |
1610 | gAliRICHminuit->mnexcm("EXIT" ,&arglist, 0, ierflag); | |
5cb4dfc3 | 1611 | |
88dd9ad4 | 1612 | gAliRICHminuit->mnpout(0,chname, trackThetaNew, eps , b1, b2, ierflg); |
1613 | gAliRICHminuit->mnpout(1,chname, trackPhiNew, eps , b1, b2, ierflg); | |
5cb4dfc3 | 1614 | |
1615 | //values after the fit... | |
88dd9ad4 | 1616 | SetFittedTrackTheta((Float_t)trackThetaNew); |
1617 | SetFittedTrackPhi((Float_t)trackPhiNew); | |
5cb4dfc3 | 1618 | |
88dd9ad4 | 1619 | delete gAliRICHminuit; |
5cb4dfc3 | 1620 | |
1621 | } | |
1622 | ||
1623 | void AliRICHRecon::EstimationOfTheta() | |
1624 | { | |
88dd9ad4 | 1625 | // theta estimate |
5cb4dfc3 | 1626 | |
88dd9ad4 | 1627 | Int_t nPhotons = 0; |
5cb4dfc3 | 1628 | |
88dd9ad4 | 1629 | Float_t shiftX = GetShiftX(); |
1630 | Float_t shiftY = GetShiftY(); | |
5cb4dfc3 | 1631 | |
88dd9ad4 | 1632 | Float_t *candidatePhotonX = GetCandidatePhotonX(); |
1633 | Float_t *candidatePhotonY = GetCandidatePhotonY(); | |
5cb4dfc3 | 1634 | |
88dd9ad4 | 1635 | Int_t nPhotonsCandidates = GetCandidatePhotonsNumber(); |
5cb4dfc3 | 1636 | |
88dd9ad4 | 1637 | // cout << "MINIM: Nphotons " << nPhotonsCandidates << endl; |
5cb4dfc3 | 1638 | |
88dd9ad4 | 1639 | for (Int_t j=0; j < nPhotonsCandidates; j++) |
5cb4dfc3 | 1640 | { |
1641 | ||
1642 | SetPhotonIndex(j); | |
1643 | ||
1644 | if(!GetPhotonFlag()) continue; | |
1645 | ||
88dd9ad4 | 1646 | Float_t xtoentr = candidatePhotonX[j] - shiftX; |
1647 | Float_t ytoentr = candidatePhotonY[j] - shiftY; | |
5cb4dfc3 | 1648 | |
88dd9ad4 | 1649 | SetEntranceX(xtoentr); |
1650 | SetEntranceY(ytoentr); | |
5cb4dfc3 | 1651 | |
1652 | FindPhiPoint(); | |
1653 | ||
1654 | FindThetaPhotonCerenkov(); | |
1655 | ||
88dd9ad4 | 1656 | Float_t thetaPhotonCerenkov = GetThetaPhotonCerenkov(); |
5cb4dfc3 | 1657 | |
88dd9ad4 | 1658 | // cout << " ACCEPTED!!! " << thetaPhotonCerenkov << endl; |
5cb4dfc3 | 1659 | |
88dd9ad4 | 1660 | SetPhotonEta(thetaPhotonCerenkov); |
5cb4dfc3 | 1661 | |
88dd9ad4 | 1662 | nPhotons++; |
5cb4dfc3 | 1663 | |
1664 | } | |
1665 | ||
1666 | Float_t xmean = 0.; | |
1667 | Float_t x2mean = 0.; | |
1668 | Int_t nev = 0; | |
1669 | ||
88dd9ad4 | 1670 | for (Int_t j=0; j < nPhotonsCandidates;j++) |
5cb4dfc3 | 1671 | { |
1672 | SetPhotonIndex(j); | |
1673 | ||
1674 | Float_t eta = GetPhotonEta(); | |
1675 | ||
1676 | if(eta != -999.) | |
1677 | { | |
1678 | if(GetPhotonFlag() == 2) | |
1679 | { | |
1680 | xmean += eta; | |
1681 | x2mean += eta*eta; | |
1682 | nev++; | |
1683 | } | |
1684 | } | |
1685 | } | |
1686 | ||
1687 | if(nev > 0) | |
1688 | { | |
1689 | xmean /=(Float_t)nev; | |
1690 | x2mean /=(Float_t)nev; | |
1691 | } else { | |
1692 | xmean = 0.; | |
1693 | x2mean = 0.; | |
1694 | } | |
1695 | ||
88dd9ad4 | 1696 | Float_t vRMS = sqrt(x2mean - xmean*xmean); |
5cb4dfc3 | 1697 | |
88dd9ad4 | 1698 | // cout << " RMS " << vRMS; |
5cb4dfc3 | 1699 | |
1700 | SetEstimationOfTheta(xmean); | |
88dd9ad4 | 1701 | SetEstimationOfThetaRMS(vRMS); |
5cb4dfc3 | 1702 | } |
1703 | ||
b068561d | 1704 | void fcnrecon(Int_t& /*npar*/, Double_t* /*gin*/, Double_t &f, Double_t *par, Int_t) |
5cb4dfc3 | 1705 | { |
88dd9ad4 | 1706 | // function to be minimized |
1707 | AliRICHRecon *gMyRecon = (AliRICHRecon*)gAliRICHminuit->GetObjectFit(); | |
5cb4dfc3 | 1708 | |
1709 | Float_t p0 = (Float_t)par[0]; | |
1710 | Float_t p1 = (Float_t)par[1]; | |
1711 | ||
1712 | gMyRecon->SetTrackTheta(p0); | |
1713 | gMyRecon->SetTrackPhi(p1); | |
1714 | ||
1715 | gMyRecon->EstimationOfTheta(); | |
88dd9ad4 | 1716 | Float_t vRMS = gMyRecon->GetEstimationOfThetaRMS(); |
5cb4dfc3 | 1717 | |
88dd9ad4 | 1718 | Int_t houghPhotons = gMyRecon->GetHoughPhotons(); |
5cb4dfc3 | 1719 | |
1720 | ||
88dd9ad4 | 1721 | f = (Double_t)(1000*vRMS/(Float_t)houghPhotons); |
5cb4dfc3 | 1722 | |
b068561d | 1723 | // if(fDebug) cout << " f " << f |
1724 | // << " theta " << par[0] << " phi " << par[1] | |
88dd9ad4 | 1725 | // << " HoughPhotons " << houghPhotons << endl; |
b068561d | 1726 | // |
1727 | // if(fDebug&&iflag == 3) | |
1728 | // { | |
1729 | // cout << " --- end convergence...summary --- " << endl; | |
1730 | // cout << " theta " << par[0] << endl; | |
1731 | // cout << " phi " << par[1] << endl; | |
1732 | // } | |
5cb4dfc3 | 1733 | } |
1734 | ||
b068561d | 1735 | void AliRICHRecon::Waiting() |
5cb4dfc3 | 1736 | { |
88dd9ad4 | 1737 | // wait, wait.... |
b068561d | 1738 | if(!fIsDISPLAY) return; |
5cb4dfc3 | 1739 | cout << " Press any key to continue..."; |
1740 | ||
1741 | // gSystem->ProcessEvents(); | |
1742 | getchar(); | |
1743 | ||
1744 | cout << endl; | |
1745 | ||
1746 | return; | |
1747 | } | |
1748 |