]>
Commit | Line | Data |
---|---|---|
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 | // AliRICHRecon // | |
19 | // // | |
20 | // RICH class to perfom pattern recognition based on Hough transfrom // | |
21 | // // | |
22 | ////////////////////////////////////////////////////////////////////////// | |
23 | ||
24 | #include <Riostream.h> | |
25 | #include <TCanvas.h> | |
26 | #include <TFile.h> | |
27 | #include <TH2.h> | |
28 | #include <TMath.h> | |
29 | #include <TMath.h> | |
30 | #include <TMinuit.h> | |
31 | #include <TNtuple.h> | |
32 | #include <TParticle.h> | |
33 | #include <TRandom.h> | |
34 | #include <TRotation.h> | |
35 | #include <TVector3.h> | |
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" | |
44 | ||
45 | #define NPointsOfRing 201 | |
46 | ||
47 | TMinuit *gAliRICHminuit ; | |
48 | ||
49 | void fcnrecon(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag); | |
50 | //__________________________________________________________________________________________________ | |
51 | AliRICHRecon::AliRICHRecon(const char*name, const char*title) | |
52 | :TTask(name,title) | |
53 | { | |
54 | // main ctor | |
55 | fThetaBin=750; fThetaMin = 0.0; fThetaMax = 0.75; | |
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.; | |
68 | fRich = (AliRICH*)gAlice->GetDetector("RICH"); | |
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", | |
72 | "Run:Trig:VertZ:Pmod:Pt:Eta:TrackTheta:TrackPhi:TrackThetaFit:TrackPhiFit:Charge::NPhotons:NPhotonsFit:InRing:MassOfParticle:HoughArea:Multiplicity:TPCLastZ"); | |
73 | } | |
74 | //__________________________________________________________________________________________________ | |
75 | void AliRICHRecon::StartProcessEvent() | |
76 | { | |
77 | //start to process for pattern recognition | |
78 | ||
79 | Float_t trackThetaStored = 0; | |
80 | Float_t trackPhiStored = 0; | |
81 | Float_t thetaCerenkovStored = 0; | |
82 | Int_t houghPhotonsStored = 0; | |
83 | ||
84 | SetFreonScaleFactor(0.994); | |
85 | ||
86 | if(fIsDISPLAY) | |
87 | { | |
88 | DrawEvent(0); | |
89 | // Waiting(); | |
90 | } | |
91 | ||
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); | |
102 | ||
103 | Float_t clusX[7][500],clusY[7][500]; | |
104 | Int_t clusQ[7][500],clusMul[7][500]; | |
105 | Int_t nClusters[7]; | |
106 | ||
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(); | |
116 | } | |
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++) { | |
133 | ||
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 | ||
192 | } | |
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) | |
268 | cout << " houghPhotonsBefore " << houghPhotonsBefore | |
269 | << " GetHoughPhotons() " << GetHoughPhotons(); | |
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 | } | |
292 | if(fIsDISPLAY) fDisplay->Print("display.ps"); | |
293 | }//StartProcessEvent() | |
294 | //__________________________________________________________________________________________________ | |
295 | void AliRICHRecon::EndProcessEvent() | |
296 | { | |
297 | // function called at the end of the event loop | |
298 | ||
299 | fOutFile->Write(); | |
300 | fOutFile->Close(); | |
301 | } | |
302 | //__________________________________________________________________________________________________ | |
303 | void AliRICHRecon::PatRec() | |
304 | { | |
305 | //pattern recognition method based on Hough transform | |
306 | ||
307 | ||
308 | Float_t trackTheta = GetTrackTheta(); | |
309 | Float_t trackPhi = GetTrackPhi(); | |
310 | Float_t pmod = GetTrackMomentum(); | |
311 | Int_t iMipIndex = GetMipIndex(); | |
312 | ||
313 | Bool_t kPatRec = kFALSE; | |
314 | ||
315 | Int_t candidatePhotons = 0; | |
316 | ||
317 | Float_t shiftX = GetShiftX(); | |
318 | Float_t shiftY = GetShiftY(); | |
319 | ||
320 | Float_t* candidatePhotonX = GetCandidatePhotonX(); | |
321 | Float_t* candidatePhotonY = GetCandidatePhotonY(); | |
322 | ||
323 | Int_t candidatePhotonsNumber = GetCandidatePhotonsNumber(); | |
324 | ||
325 | if(fDebug) cout << " n " << candidatePhotonsNumber << endl; | |
326 | ||
327 | SetThetaCerenkov(999.); | |
328 | SetHoughPhotons(0); | |
329 | SetHoughPhotonsNorm(0); | |
330 | SetHoughRMS(999.); | |
331 | ||
332 | for (Int_t j=0; j < candidatePhotonsNumber; j++) | |
333 | { | |
334 | ||
335 | SetPhotonIndex(j); | |
336 | ||
337 | SetPhotonFlag(0); | |
338 | SetPhotonEta(-999.); | |
339 | SetPhotonWeight(0.); | |
340 | ||
341 | if (j == iMipIndex) continue; | |
342 | ||
343 | ||
344 | if(candidatePhotonX[j] < -64.) continue; /* avoid artificial clusters from edge uesd by Yale.... */ | |
345 | ||
346 | Float_t xtoentr = candidatePhotonX[j] - shiftX; | |
347 | Float_t ytoentr = candidatePhotonY[j] - shiftY; | |
348 | ||
349 | // Float_t chargehit = fHits_charge[j]; | |
350 | // if(chargehit > 150) continue; | |
351 | ||
352 | SetEntranceX(xtoentr); | |
353 | SetEntranceY(ytoentr); | |
354 | ||
355 | FindPhiPoint(); | |
356 | ||
357 | Int_t photonStatus = PhotonInBand(); | |
358 | ||
359 | if(fDebug) | |
360 | { | |
361 | cout << " Photon n. " << j << " Status " << photonStatus << " accepted " << endl; | |
362 | cout << " CandidatePhotonX[j] " << candidatePhotonX[j] << " CandidatePhotonY[j] " << candidatePhotonY[j] << endl; | |
363 | } | |
364 | ||
365 | if(photonStatus == 0) continue; | |
366 | ||
367 | SetPhotonFlag(1); | |
368 | ||
369 | FindThetaPhotonCerenkov(); | |
370 | ||
371 | Float_t thetaPhotonCerenkov = GetThetaPhotonCerenkov(); | |
372 | ||
373 | if(fDebug) cout << " theta photon " << thetaPhotonCerenkov << endl; | |
374 | ||
375 | SetPhotonEta(thetaPhotonCerenkov); | |
376 | ||
377 | candidatePhotons++; | |
378 | ||
379 | ||
380 | } | |
381 | ||
382 | if(candidatePhotons >= 1) kPatRec = kTRUE; | |
383 | ||
384 | if(!kPatRec) return; | |
385 | { | |
386 | SetThetaCerenkov(999.); | |
387 | SetHoughPhotons(0); | |
388 | } | |
389 | SetPhotonsNumber(candidatePhotonsNumber); | |
390 | ||
391 | HoughResponse(); | |
392 | ||
393 | fNrings++; | |
394 | ||
395 | FlagPhotons(); | |
396 | Int_t nPhotonHough = GetHoughPhotons(); | |
397 | ||
398 | if(nPhotonHough < 1) | |
399 | { | |
400 | SetThetaCerenkov(999.); | |
401 | SetHoughPhotonsNorm(0.); | |
402 | return; | |
403 | } | |
404 | ||
405 | if(fIsWEIGHT) FindWeightThetaCerenkov(); | |
406 | ||
407 | Float_t thetaCerenkov = GetThetaCerenkov(); | |
408 | ||
409 | SetThetaOfRing(thetaCerenkov); | |
410 | FindAreaAndPortionOfRing(); | |
411 | ||
412 | Float_t nPhotonHoughNorm = ((Float_t)nPhotonHough)/GetPortionOfRing(); | |
413 | SetHoughPhotonsNorm(nPhotonHoughNorm); | |
414 | ||
415 | // Calculate the area where the photon are accepted... | |
416 | ||
417 | Float_t thetaInternal = thetaCerenkov - 0.5*fWindowWidth; | |
418 | SetThetaOfRing(thetaInternal); | |
419 | FindAreaAndPortionOfRing(); | |
420 | Float_t internalArea = GetAreaOfRing(); | |
421 | ||
422 | Float_t thetaExternal = thetaCerenkov + 0.5*fWindowWidth; | |
423 | SetThetaOfRing(thetaExternal); | |
424 | FindAreaAndPortionOfRing(); | |
425 | Float_t externalArea = GetAreaOfRing(); | |
426 | ||
427 | Float_t houghArea = externalArea - internalArea; | |
428 | ||
429 | SetHoughArea(houghArea); | |
430 | ||
431 | if(fDebug) | |
432 | { | |
433 | cout << " ----- SUMMARY OF RECONSTRUCTION ----- " << endl; | |
434 | cout << " Rings found " << fNrings << " with thetac " << thetaCerenkov << endl; | |
435 | ||
436 | ||
437 | cout << " Nphotons " << GetPhotonsNumber() | |
438 | << " Hough " << nPhotonHough | |
439 | << " norm " << nPhotonHoughNorm << endl; | |
440 | ||
441 | cout << " In PatRec:p " << pmod << " theta " << trackTheta << " phi " << trackPhi << endl; | |
442 | cout << " ------------------------------------- " << endl; | |
443 | } | |
444 | ||
445 | Int_t nPhotons = GetPhotonsNumber(); | |
446 | ||
447 | Float_t xmean = 0.; | |
448 | Float_t x2mean = 0.; | |
449 | Int_t nev = 0; | |
450 | ||
451 | for (Int_t j=0; j < nPhotons;j++) | |
452 | { | |
453 | SetPhotonIndex(j); | |
454 | ||
455 | Float_t eta = GetPhotonEta(); | |
456 | ||
457 | if(eta != -999.) | |
458 | { | |
459 | if(GetPhotonFlag() == 2) | |
460 | { | |
461 | ||
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 | ||
479 | Float_t vRMS = sqrt(x2mean - xmean*xmean); | |
480 | ||
481 | SetHoughRMS(vRMS); | |
482 | ||
483 | if(fDebug) cout << " RMS " << vRMS << endl; | |
484 | ||
485 | } | |
486 | ||
487 | void AliRICHRecon::FindEmissionPoint() | |
488 | { | |
489 | //estimate the emission point in radiator | |
490 | ||
491 | // Find emission point | |
492 | ||
493 | Float_t absorbtionLenght=7.83*fRadiatorWidth; //absorption length in the freon (cm) | |
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 | ||
498 | photonLenght=exp(-fRadiatorWidth/(absorbtionLenght*cos(fCerenkovAnglePad))); | |
499 | photonLenghtMin=fRadiatorWidth*photonLenght/(1.-photonLenght); | |
500 | photonLenghtMax=absorbtionLenght*cos(fCerenkovAnglePad); | |
501 | Float_t emissionPoint = fRadiatorWidth + photonLenghtMin - photonLenghtMax; | |
502 | ||
503 | SetEmissionPoint(emissionPoint); | |
504 | } | |
505 | ||
506 | ||
507 | Int_t AliRICHRecon::PhotonInBand() | |
508 | { | |
509 | //search band fro photon candidates | |
510 | ||
511 | // Float_t massOfParticle; | |
512 | Float_t beta; | |
513 | Float_t nfreon; | |
514 | ||
515 | Float_t thetacer; | |
516 | ||
517 | Float_t xtoentr = GetEntranceX(); | |
518 | Float_t ytoentr = GetEntranceY(); | |
519 | ||
520 | Float_t innerRadius; | |
521 | Float_t outerRadius; | |
522 | ||
523 | Float_t phpad = GetPhiPoint(); | |
524 | ||
525 | // Float_t pmod = GetTrackMomentum(); | |
526 | // Float_t trackTheta = GetTrackTheta(); | |
527 | // Float_t trackPhi = GetTrackPhi(); | |
528 | ||
529 | // inner radius // | |
530 | SetPhotonEnergy(5.6); | |
531 | SetEmissionPoint(fRadiatorWidth -0.0001); | |
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 | { | |
550 | innerRadius = -999.; | |
551 | SetXInnerRing(-999.); | |
552 | SetYInnerRing(-999.); | |
553 | SetRadiusInnerRing(-999.); | |
554 | } | |
555 | else | |
556 | { | |
557 | SetThetaPhotonInDRS(GetThetaAtQuartz()); | |
558 | SetPhiPhotonInDRS(phpad); | |
559 | ||
560 | innerRadius = FromEmissionToCathode(); | |
561 | if(innerRadius == 999.) innerRadius = -999.; | |
562 | ||
563 | SetXInnerRing(GetXPointOnCathode()); | |
564 | SetYInnerRing(GetYPointOnCathode()); | |
565 | SetRadiusInnerRing(innerRadius); | |
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 | { | |
590 | outerRadius = 999.; | |
591 | SetXOuterRing(999.); | |
592 | SetYOuterRing(999.); | |
593 | SetRadiusOuterRing(999.); | |
594 | } | |
595 | else | |
596 | { | |
597 | SetThetaPhotonInDRS(GetThetaAtQuartz()); | |
598 | SetPhiPhotonInDRS(phpad); | |
599 | ||
600 | outerRadius = FromEmissionToCathode(); | |
601 | // cout << " outerRadius " << outerRadius << endl; | |
602 | SetXOuterRing(GetXPointOnCathode()); | |
603 | SetYOuterRing(GetYPointOnCathode()); | |
604 | SetRadiusOuterRing(outerRadius); | |
605 | } | |
606 | ||
607 | Float_t padradius = sqrt(TMath::Power(xtoentr,2)+TMath::Power(ytoentr,2)); | |
608 | ||
609 | if(fDebug) printf(" rmin %f r %f rmax %f \n",innerRadius,padradius,outerRadius); | |
610 | ||
611 | if(padradius>=innerRadius && padradius<=outerRadius) return 1; | |
612 | return 0; | |
613 | } | |
614 | ||
615 | void AliRICHRecon::FindThetaAtQuartz(Float_t thetaCerenkov) | |
616 | { | |
617 | //find the theta at the quartz plate | |
618 | ||
619 | if(thetaCerenkov == 999.) | |
620 | { | |
621 | SetThetaAtQuartz(999.); | |
622 | return; | |
623 | } | |
624 | ||
625 | Float_t thetaAtQuartz = 999.; | |
626 | ||
627 | Float_t trackTheta = GetTrackTheta(); | |
628 | ||
629 | if(trackTheta == 0) { | |
630 | ||
631 | if(fDebug) cout << " Theta sol unique " << thetaCerenkov << endl; | |
632 | ||
633 | thetaAtQuartz = thetaCerenkov; | |
634 | SetThetaAtQuartz(thetaAtQuartz); | |
635 | return; | |
636 | } | |
637 | ||
638 | Float_t trackPhi = GetTrackPhi(); | |
639 | Float_t phiPoint = GetPhiPoint(); | |
640 | ||
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; | |
649 | ||
650 | Double_t underSqrt = 1 + b*b - c*c; | |
651 | ||
652 | if(fDebug) | |
653 | { | |
654 | cout << " trackTheta " << trackTheta << endl; | |
655 | cout << " TrackPhi " << trackPhi << endl; | |
656 | cout << " PhiPoint " << phiPoint << endl; | |
657 | cout << " ThetaCerenkov " << thetaCerenkov << endl; | |
658 | cout << " den b c " << den << " b " << b << " c " << c << endl; | |
659 | } | |
660 | ||
661 | if(underSqrt < 0) { | |
662 | if(fDebug) cout << " sqrt negative !!!!" << underSqrt << endl; | |
663 | SetThetaAtQuartz(999.); | |
664 | return; | |
665 | } | |
666 | ||
667 | Double_t sol1 = (1+TMath::Sqrt(underSqrt))/(b-c); | |
668 | Double_t sol2 = (1-TMath::Sqrt(underSqrt))/(b-c); | |
669 | ||
670 | Double_t thetaSol1 = 2*TMath::ATan(sol1); | |
671 | Double_t thetaSol2 = 2*TMath::ATan(sol2); | |
672 | ||
673 | if(fDebug) cout << " Theta sol 1 " << thetaSol1 | |
674 | << " Theta sol 2 " << thetaSol2 << endl; | |
675 | ||
676 | if(thetaSol1>0 && thetaSol1 < TMath::Pi()) thetaAtQuartz = (Float_t)thetaSol1; | |
677 | if(thetaSol2>0 && thetaSol2 < TMath::Pi()) thetaAtQuartz = (Float_t)thetaSol2; | |
678 | ||
679 | SetThetaAtQuartz(thetaAtQuartz); | |
680 | } | |
681 | ||
682 | void AliRICHRecon::FindThetaPhotonCerenkov() | |
683 | { | |
684 | //find theta cerenkov of ring | |
685 | ||
686 | Float_t thetaCerMin = 0.; | |
687 | Float_t thetaCerMax = 0.75; | |
688 | Float_t thetaCerMean; | |
689 | ||
690 | Float_t radiusMin, radiusMax, radiusMean; | |
691 | Int_t nIteration = 0; | |
692 | ||
693 | const Float_t kTollerance = 0.05; | |
694 | ||
695 | // Float_t pmod = GetTrackMomentum(); | |
696 | // Float_t trackTheta = GetTrackTheta(); | |
697 | // Float_t trackPhi = GetTrackPhi(); | |
698 | ||
699 | Float_t phiPoint = GetPhiPoint(); | |
700 | ||
701 | SetPhotonEnergy(6.85); | |
702 | SetEmissionPoint(fRadiatorWidth/2); | |
703 | ||
704 | Float_t xPoint = GetEntranceX(); | |
705 | Float_t yPoint = GetEntranceY(); | |
706 | Float_t distPoint = sqrt(xPoint*xPoint + yPoint*yPoint); | |
707 | ||
708 | if(fDebug) cout << " DistPoint " << distPoint << endl; | |
709 | ||
710 | // Star minimization... | |
711 | ||
712 | // First value... | |
713 | ||
714 | FindThetaAtQuartz(thetaCerMin); | |
715 | ||
716 | if(GetThetaAtQuartz() == 999.) | |
717 | { | |
718 | radiusMin = -999.; | |
719 | } | |
720 | else | |
721 | { | |
722 | SetThetaPhotonInDRS(GetThetaAtQuartz()); | |
723 | SetPhiPhotonInDRS(phiPoint); | |
724 | ||
725 | radiusMin = FromEmissionToCathode(); | |
726 | } | |
727 | ||
728 | // Second value... | |
729 | ||
730 | FindThetaAtQuartz(thetaCerMax); | |
731 | if(GetThetaAtQuartz() == 999.) | |
732 | { | |
733 | radiusMax = 999.; | |
734 | } | |
735 | else | |
736 | { | |
737 | SetThetaPhotonInDRS(GetThetaAtQuartz()); | |
738 | SetPhiPhotonInDRS(phiPoint); | |
739 | ||
740 | radiusMax = FromEmissionToCathode(); | |
741 | } | |
742 | // Mean value... | |
743 | ||
744 | thetaCerMean = (thetaCerMax + thetaCerMin)/2; | |
745 | ||
746 | FindThetaAtQuartz(thetaCerMean); | |
747 | if(GetThetaAtQuartz() == 999.) | |
748 | { | |
749 | radiusMean = 999.; | |
750 | } | |
751 | else | |
752 | { | |
753 | SetThetaPhotonInDRS(GetThetaAtQuartz()); | |
754 | SetPhiPhotonInDRS(phiPoint); | |
755 | ||
756 | radiusMean = FromEmissionToCathode(); | |
757 | } | |
758 | ||
759 | if(fDebug) cout << " r1 " << radiusMin << " rmean " | |
760 | << radiusMean << " r2 " << radiusMax << endl; | |
761 | ||
762 | while (TMath::Abs(radiusMean-distPoint) > kTollerance) | |
763 | { | |
764 | ||
765 | if((radiusMin-distPoint)*(radiusMean-distPoint) < 0) thetaCerMax = thetaCerMean; | |
766 | if((radiusMin-distPoint)*(radiusMean-distPoint) > 0) { | |
767 | ||
768 | thetaCerMin = thetaCerMean; | |
769 | ||
770 | FindThetaAtQuartz(thetaCerMin); | |
771 | SetThetaPhotonInDRS(GetThetaAtQuartz()); | |
772 | SetPhiPhotonInDRS(phiPoint); | |
773 | ||
774 | radiusMin =FromEmissionToCathode(); | |
775 | } | |
776 | ||
777 | thetaCerMean = (thetaCerMax + thetaCerMin)/2; | |
778 | ||
779 | FindThetaAtQuartz(thetaCerMean); | |
780 | SetThetaPhotonInDRS(GetThetaAtQuartz()); | |
781 | SetPhiPhotonInDRS(phiPoint); | |
782 | ||
783 | radiusMean = FromEmissionToCathode(); | |
784 | ||
785 | nIteration++; | |
786 | if(nIteration>=50) { | |
787 | if(fDebug) printf(" max iterations in FindPhotonCerenkov\n"); | |
788 | SetThetaPhotonCerenkov(999.); | |
789 | return; | |
790 | } | |
791 | } | |
792 | ||
793 | SetThetaPhotonCerenkov(thetaCerMean); | |
794 | ||
795 | } | |
796 | ||
797 | void AliRICHRecon::FindAreaAndPortionOfRing() | |
798 | { | |
799 | //find fraction of the ring accepted by the RICH | |
800 | ||
801 | Float_t xPoint[NPointsOfRing], yPoint[NPointsOfRing]; | |
802 | ||
803 | // Float_t xtoentr = GetEntranceX(); | |
804 | // Float_t ytoentr = GetEntranceY(); | |
805 | Float_t shiftX = GetShiftX(); | |
806 | Float_t shiftY = GetShiftY(); | |
807 | ||
808 | Float_t xemiss = GetXCoordOfEmission(); | |
809 | Float_t yemiss = GetYCoordOfEmission(); | |
810 | ||
811 | Float_t x0 = xemiss + shiftX; | |
812 | Float_t y0 = yemiss + shiftY; | |
813 | ||
814 | // Float_t pmod = GetTrackMomentum(); | |
815 | // Float_t trackTheta = GetTrackTheta(); | |
816 | // Float_t trackPhi = GetTrackPhi(); | |
817 | ||
818 | SetPhotonEnergy(6.85); | |
819 | SetFreonRefractiveIndex(); | |
820 | ||
821 | SetEmissionPoint(fRadiatorWidth/2.); | |
822 | ||
823 | Float_t theta = GetThetaOfRing(); | |
824 | ||
825 | Int_t nPoints = 0; | |
826 | Int_t nPsiAccepted = 0; | |
827 | Int_t nPsiTotal = 0; | |
828 | ||
829 | for(Int_t i=0;i<NPointsOfRing-1;i++) | |
830 | { | |
831 | ||
832 | Float_t psi = 2*TMath::Pi()*i/NPointsOfRing; | |
833 | ||
834 | SetThetaPhotonInTRS(theta); | |
835 | SetPhiPhotonInTRS(psi); | |
836 | FindPhotonAnglesInDRS(); | |
837 | ||
838 | Float_t radius = FromEmissionToCathode(); | |
839 | if (radius == 999.) continue; | |
840 | ||
841 | nPsiTotal++; | |
842 | ||
843 | Float_t xPointRing = GetXPointOnCathode() + shiftX; | |
844 | Float_t yPointRing = GetYPointOnCathode() + shiftY; | |
845 | ||
846 | SetDetectorWhereX(xPointRing); | |
847 | SetDetectorWhereY(yPointRing); | |
848 | ||
849 | Int_t zone = CheckDetectorAcceptance(); | |
850 | ||
851 | ||
852 | if (zone != 0) | |
853 | { | |
854 | FindIntersectionWithDetector(); | |
855 | xPoint[nPoints] = GetIntersectionX(); | |
856 | yPoint[nPoints] = GetIntersectionY(); | |
857 | } | |
858 | else | |
859 | { | |
860 | xPoint[nPoints] = xPointRing; | |
861 | yPoint[nPoints] = yPointRing; | |
862 | nPsiAccepted++; | |
863 | } | |
864 | ||
865 | nPoints++; | |
866 | ||
867 | } | |
868 | ||
869 | xPoint[nPoints] = xPoint[0]; | |
870 | yPoint[nPoints] = yPoint[0]; | |
871 | ||
872 | // find area... | |
873 | ||
874 | Float_t area = 0; | |
875 | ||
876 | for (Int_t i = 0; i < nPoints; i++) | |
877 | { | |
878 | area += TMath::Abs((xPoint[i]-x0)*(yPoint[i+1]-y0) - (xPoint[i+1]-x0)*(yPoint[i]-y0)); | |
879 | } | |
880 | ||
881 | area *= 0.5; | |
882 | ||
883 | Float_t portionOfRing = ((Float_t)nPsiAccepted)/((Float_t)(nPsiTotal)); | |
884 | ||
885 | ||
886 | SetAreaOfRing(area); | |
887 | SetPortionOfRing(portionOfRing); | |
888 | } | |
889 | ||
890 | void AliRICHRecon::FindIntersectionWithDetector() | |
891 | { | |
892 | // find ring intersection with CsI edges | |
893 | ||
894 | Float_t xIntersect, yIntersect; | |
895 | Float_t x1, x2, y1, y2; | |
896 | ||
897 | Float_t shiftX = GetShiftX(); | |
898 | Float_t shiftY = GetShiftY(); | |
899 | ||
900 | Float_t xPoint = GetXPointOnCathode() + shiftX; | |
901 | Float_t yPoint = GetYPointOnCathode() + shiftY; | |
902 | ||
903 | Float_t xemiss = GetXCoordOfEmission(); | |
904 | Float_t yemiss = GetYCoordOfEmission(); | |
905 | ||
906 | Float_t phi = GetPhiPhotonInDRS(); | |
907 | Float_t m = tan(phi); | |
908 | ||
909 | Float_t x0 = xemiss + shiftX; | |
910 | Float_t y0 = yemiss + shiftY; | |
911 | ||
912 | if(xPoint > x0) | |
913 | { | |
914 | x1 = x0; | |
915 | x2 = xPoint; | |
916 | } | |
917 | else | |
918 | { | |
919 | x2 = x0; | |
920 | x1 = xPoint; | |
921 | } | |
922 | if(yPoint > y0) | |
923 | { | |
924 | y1 = y0; | |
925 | y2 = yPoint; | |
926 | } | |
927 | else | |
928 | { | |
929 | y2 = y0; | |
930 | y1 = yPoint; | |
931 | } | |
932 | // | |
933 | xIntersect = fXmax; | |
934 | yIntersect = m*(xIntersect - x0) + y0; | |
935 | if (yIntersect >= fYmin && yIntersect <= fYmax && xIntersect >= x1 && xIntersect <= x2) | |
936 | { | |
937 | SetIntersectionX(xIntersect); | |
938 | SetIntersectionY(yIntersect); | |
939 | return; | |
940 | } | |
941 | // | |
942 | xIntersect = fXmin; | |
943 | yIntersect = m*(xIntersect - x0) + y0; | |
944 | if (yIntersect >= fYmin && yIntersect <= fYmax && xIntersect >= x1 && xIntersect <= x2) | |
945 | { | |
946 | SetIntersectionX(xIntersect); | |
947 | SetIntersectionY(yIntersect); | |
948 | return; | |
949 | } | |
950 | // | |
951 | yIntersect = fYmax; | |
952 | xIntersect = (yIntersect - y0)/m + x0; | |
953 | if (xIntersect >= fXmin && xIntersect <= fXmax && yIntersect >= y1 && yIntersect <= y2) | |
954 | { | |
955 | SetIntersectionX(xIntersect); | |
956 | SetIntersectionY(yIntersect); | |
957 | return; | |
958 | } | |
959 | // | |
960 | yIntersect = fYmin; | |
961 | xIntersect = (yIntersect - y0)/m + x0; | |
962 | if (xIntersect >= fXmin && xIntersect <= fXmax && yIntersect >= y1 && yIntersect <= y2) | |
963 | { | |
964 | SetIntersectionX(xIntersect); | |
965 | SetIntersectionY(yIntersect); | |
966 | return; | |
967 | } | |
968 | ||
969 | cout << " sono fuori!!!!!!" << endl; | |
970 | ||
971 | } | |
972 | //__________________________________________________________________________________________________ | |
973 | Int_t AliRICHRecon::CheckDetectorAcceptance() const | |
974 | { | |
975 | // check for the acceptance | |
976 | ||
977 | // crosses X -2.6 2.6 cm | |
978 | // crosses Y -1 1 cm | |
979 | ||
980 | Float_t xcoord = GetDetectorWhereX(); | |
981 | Float_t ycoord = GetDetectorWhereY(); | |
982 | ||
983 | if(xcoord > fXmax) | |
984 | { | |
985 | if(ycoord > fYmax) return 2; | |
986 | if(ycoord > fYmin && ycoord < fYmax) return 3; | |
987 | if(ycoord < fYmin) return 4; | |
988 | } | |
989 | if(xcoord < fXmin) | |
990 | { | |
991 | if(ycoord > fYmax) return 8; | |
992 | if(ycoord > fYmin && ycoord < fYmax) return 7; | |
993 | if(ycoord < fYmin) return 6; | |
994 | } | |
995 | if(xcoord > fXmin && xcoord < fXmax) | |
996 | { | |
997 | if(ycoord > fYmax) return 1; | |
998 | if(ycoord > fYmin && ycoord < fYmax) return 0; | |
999 | if(ycoord < fYmin) return 5; | |
1000 | } | |
1001 | return 999; | |
1002 | } | |
1003 | //__________________________________________________________________________________________________ | |
1004 | Float_t AliRICHRecon::PhotonPositionOnCathode() | |
1005 | { | |
1006 | // find the photon position on the CsI | |
1007 | // Float_t massOfParticle; | |
1008 | Float_t beta; | |
1009 | Float_t nfreon; | |
1010 | ||
1011 | ||
1012 | SetPhotonEnergy(6.85); | |
1013 | SetEmissionPoint(fRadiatorWidth/2.); | |
1014 | SetMassHypotesis(0.139567); | |
1015 | ||
1016 | SetBetaOfParticle(); | |
1017 | SetFreonRefractiveIndex(); | |
1018 | ||
1019 | beta = GetBetaOfParticle(); | |
1020 | nfreon = GetFreonRefractiveIndex(); | |
1021 | ||
1022 | ||
1023 | Float_t radius = FromEmissionToCathode(); | |
1024 | if (radius == 999.) return 999.; | |
1025 | ||
1026 | return 0; | |
1027 | } | |
1028 | ||
1029 | void AliRICHRecon::FindPhotonAnglesInDRS() | |
1030 | { | |
1031 | // Setup the rotation matrix of the track... | |
1032 | ||
1033 | TRotation mtheta; | |
1034 | TRotation mphi; | |
1035 | TRotation minv; | |
1036 | TRotation mrot; | |
1037 | ||
1038 | Float_t trackTheta = GetTrackTheta(); | |
1039 | Float_t trackPhi = GetTrackPhi(); | |
1040 | ||
1041 | mtheta.RotateY(trackTheta); | |
1042 | mphi.RotateZ(trackPhi); | |
1043 | ||
1044 | mrot = mphi * mtheta; | |
1045 | // minv = mrot.Inverse(); | |
1046 | ||
1047 | TVector3 photonInRadiator(1,1,1); | |
1048 | ||
1049 | Float_t thetaCerenkov = GetThetaPhotonInTRS(); | |
1050 | Float_t phiCerenkov = GetPhiPhotonInTRS(); | |
1051 | ||
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); | |
1059 | ||
1060 | } | |
1061 | ||
1062 | Float_t AliRICHRecon::FromEmissionToCathode() | |
1063 | { | |
1064 | // trace from emission point to cathode | |
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 | ||
1076 | Float_t trackTheta = GetTrackTheta(); | |
1077 | Float_t trackPhi = GetTrackPhi(); | |
1078 | Float_t lengthOfEmissionPoint = GetEmissionPoint(); | |
1079 | ||
1080 | Float_t theta = GetThetaPhotonInDRS(); | |
1081 | Float_t phi = GetPhiPhotonInDRS(); | |
1082 | ||
1083 | // cout << " Theta " << Theta << " Phi " << Phi << endl; | |
1084 | ||
1085 | Float_t xemiss = lengthOfEmissionPoint*tan(trackTheta)*cos(trackPhi); | |
1086 | Float_t yemiss = lengthOfEmissionPoint*tan(trackTheta)*sin(trackPhi); | |
1087 | ||
1088 | SetXCoordOfEmission(xemiss); | |
1089 | SetYCoordOfEmission(yemiss); | |
1090 | ||
1091 | Float_t thetaquar = SnellAngle(nfreon, nquartz, theta); | |
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 | ||
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); | |
1115 | ||
1116 | ||
1117 | Float_t xtot = xemiss + xw + xq + xg; | |
1118 | Float_t ytot = yemiss + yw + yq + yg; | |
1119 | ||
1120 | SetXPointOnCathode(xtot); | |
1121 | SetYPointOnCathode(ytot); | |
1122 | ||
1123 | ||
1124 | Float_t distanceFromEntrance = sqrt(TMath::Power(fPhotonLimitX,2) | |
1125 | +TMath::Power(fPhotonLimitY,2)); | |
1126 | ||
1127 | return distanceFromEntrance; | |
1128 | ||
1129 | } | |
1130 | ||
1131 | ||
1132 | void AliRICHRecon::FindPhiPoint() | |
1133 | { | |
1134 | //find phi of generated point | |
1135 | ||
1136 | Float_t xtoentr = GetEntranceX(); | |
1137 | Float_t ytoentr = GetEntranceY(); | |
1138 | ||
1139 | Float_t trackTheta = GetTrackTheta(); | |
1140 | Float_t trackPhi = GetTrackPhi(); | |
1141 | ||
1142 | Float_t emissionPoint = GetEmissionPoint(); | |
1143 | ||
1144 | Float_t argY = ytoentr - emissionPoint*tan(trackTheta)*sin(trackPhi); | |
1145 | Float_t argX = xtoentr - emissionPoint*tan(trackTheta)*cos(trackPhi); | |
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 | { | |
1154 | // cerenkov angle from n and beta | |
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) | |
1171 | { | |
1172 | // Snell law | |
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() | |
1193 | { | |
1194 | //Hough response | |
1195 | ||
1196 | // Implement Hough response pat. rec. method | |
1197 | ||
1198 | Float_t *hCSspace; | |
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 | ||
1216 | Float_t thetaCerenkov = 0.; | |
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 | ||
1224 | Int_t nPhotons = GetPhotonsNumber(); | |
1225 | ||
1226 | Int_t weightFlag = 0; | |
1227 | ||
1228 | for (k=0; k< nPhotons; k++) { | |
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 | ||
1247 | if(fIsWEIGHT) | |
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 | { | |
1268 | weightFlag = 1; | |
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 | ||
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 | ||
1297 | ||
1298 | if(weightFlag == 0) | |
1299 | { | |
1300 | hCSspace = hcsw; | |
1301 | } | |
1302 | else | |
1303 | { | |
1304 | hCSspace = hcs; | |
1305 | // cout << " probems with weight...normal procedure adopted " << endl; | |
1306 | } | |
1307 | ||
1308 | HoughFiltering(hCSspace); | |
1309 | ||
1310 | for (bin=0; bin <nBin; bin++) { | |
1311 | angle = (bin+0.5) * (fDTheta); | |
1312 | if (hCSspace[bin] && hCSspace[bin] > etaPeakPos) { | |
1313 | etaPeakCount = 0; | |
1314 | etaPeakPos = hCSspace[bin]; | |
1315 | etaPeak[0]=angle; | |
1316 | } | |
1317 | else { | |
1318 | if (hCSspace[bin] == etaPeakPos) { | |
1319 | etaPeak[++etaPeakCount] = angle; | |
1320 | } | |
1321 | } | |
1322 | } | |
1323 | ||
1324 | for (i=0; i<etaPeakCount+1; i++) { | |
1325 | thetaCerenkov += etaPeak[i]; | |
1326 | } | |
1327 | if (etaPeakCount>=0) { | |
1328 | thetaCerenkov /= etaPeakCount+1; | |
1329 | fThetaPeakPos = etaPeakPos; | |
1330 | } | |
1331 | ||
1332 | SetThetaCerenkov(thetaCerenkov); | |
1333 | } | |
1334 | ||
1335 | ||
1336 | void AliRICHRecon::HoughFiltering(float hcs[]) | |
1337 | { | |
1338 | // filter for Hough | |
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 | { | |
1368 | // manage with weight for photons | |
1369 | ||
1370 | Float_t wei = 0.; | |
1371 | Float_t weightThetaCerenkov = 0.; | |
1372 | ||
1373 | Int_t nPhotons = GetPhotonsNumber(); | |
1374 | for(Int_t i=0;i<nPhotons;i++) | |
1375 | { | |
1376 | SetPhotonIndex(i); | |
1377 | ||
1378 | if(GetPhotonFlag() == 2) | |
1379 | { | |
1380 | Float_t photonEta = GetPhotonEta(); | |
1381 | Float_t photonWeight = GetPhotonWeight(); | |
1382 | weightThetaCerenkov += photonEta*photonWeight; | |
1383 | wei += photonWeight; | |
1384 | } | |
1385 | } | |
1386 | ||
1387 | if(wei != 0.) | |
1388 | { | |
1389 | weightThetaCerenkov /= wei; | |
1390 | } | |
1391 | else | |
1392 | { | |
1393 | weightThetaCerenkov = 0.; | |
1394 | } | |
1395 | ||
1396 | SetThetaCerenkov(weightThetaCerenkov); | |
1397 | ||
1398 | cout << " thetac weighted -> " << weightThetaCerenkov << endl; | |
1399 | } | |
1400 | ||
1401 | ||
1402 | void AliRICHRecon::FlagPhotons() | |
1403 | { | |
1404 | // flag photons | |
1405 | ||
1406 | Int_t nPhotonHough = 0; | |
1407 | ||
1408 | Float_t thetaCerenkov = GetThetaCerenkov(); | |
1409 | if(fDebug) cout << " fThetaCerenkov " << thetaCerenkov << endl; | |
1410 | ||
1411 | Float_t thetaDist= thetaCerenkov - fThetaMin; | |
1412 | Int_t steps = (Int_t)(thetaDist / fDTheta); | |
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; | |
1422 | if(fDebug) cout << " thetac " << thetaCerenkov << endl; | |
1423 | ||
1424 | // Int_t candidatePhotonsNumber = GetCandidatePhotonsNumber(); | |
1425 | ||
1426 | Int_t nPhotons = GetPhotonsNumber(); | |
1427 | ||
1428 | // for(Int_t i=0;i<candidatePhotonsNumber;i++) | |
1429 | ||
1430 | for(Int_t i=0;i<nPhotons;i++) | |
1431 | { | |
1432 | SetPhotonIndex(i); | |
1433 | ||
1434 | Float_t photonEta = GetPhotonEta(); | |
1435 | ||
1436 | if(photonEta == -999.) continue; | |
1437 | ||
1438 | if(photonEta >= tmin && photonEta <= tmax) | |
1439 | { | |
1440 | SetPhotonFlag(2); | |
1441 | nPhotonHough++; | |
1442 | } | |
1443 | } | |
1444 | SetHoughPhotons(nPhotonHough); | |
1445 | } | |
1446 | ||
1447 | void AliRICHRecon::DrawEvent(Int_t flag) const | |
1448 | { | |
1449 | // draw event with rings | |
1450 | ||
1451 | flag=1; // dummy to be removed... | |
1452 | } | |
1453 | ||
1454 | Float_t AliRICHRecon::FindMassOfParticle() | |
1455 | { | |
1456 | // find mass of the particle from theta cerenkov | |
1457 | ||
1458 | Float_t pmod = GetTrackMomentum(); | |
1459 | ||
1460 | SetPhotonEnergy(6.85); | |
1461 | SetFreonRefractiveIndex(); | |
1462 | ||
1463 | Float_t thetaCerenkov = GetThetaCerenkov(); | |
1464 | FindBetaFromTheta(thetaCerenkov); | |
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 | { | |
1480 | // fill histograms.. | |
1481 | ||
1482 | Float_t fittedTrackTheta, fittedTrackPhi; | |
1483 | ||
1484 | Float_t thetaCerenkov = GetThetaCerenkov(); | |
1485 | if(thetaCerenkov == 999.) return; | |
1486 | ||
1487 | Float_t vertZ = GetEventVertexZ(); | |
1488 | ||
1489 | Float_t trackTheta = GetTrackTheta(); | |
1490 | Float_t trackPhi = GetTrackPhi(); | |
1491 | Float_t pmod = GetTrackMomentum(); | |
1492 | Float_t pt = GetTrackPt(); | |
1493 | Float_t trackEta = GetTrackEta(); | |
1494 | Int_t q = GetTrackCharge(); | |
1495 | Float_t tPCLastZ = GetTrackTPCLastZ(); | |
1496 | Float_t minDist = GetMinDist(); | |
1497 | ||
1498 | fittedTrackTheta = GetFittedTrackTheta(); | |
1499 | fittedTrackPhi = GetFittedTrackPhi(); | |
1500 | Int_t fittednPhotonHough = GetFittedHoughPhotons(); | |
1501 | ||
1502 | if(fDebug) | |
1503 | { | |
1504 | cout << " p " << pmod << " ThetaC " << thetaCerenkov | |
1505 | << " rings " << fNrings << endl; | |
1506 | } | |
1507 | ||
1508 | Int_t nPhotonHough = GetHoughPhotons(); | |
1509 | // Float_t nPhotonHoughNorm = GetHoughPhotonsNorm(); | |
1510 | Float_t inRing = GetPortionOfRing(); | |
1511 | ||
1512 | Float_t massOfParticle = FindMassOfParticle(); | |
1513 | ||
1514 | Float_t houghArea = GetHoughArea(); | |
1515 | Float_t multiplicity = GetEventMultiplicity(); | |
1516 | ||
1517 | ||
1518 | Float_t var[20]; | |
1519 | ||
1520 | var[0] = 0; | |
1521 | var[1] = 0; | |
1522 | var[2] = vertZ; | |
1523 | var[3] = pmod; | |
1524 | var[4] = pt; | |
1525 | var[5] = trackEta; | |
1526 | var[6] = trackTheta; | |
1527 | var[7] = trackPhi; | |
1528 | var[8] = fittedTrackTheta; | |
1529 | var[9] = fittedTrackPhi; | |
1530 | var[10] = q; | |
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; | |
1540 | ||
1541 | fNtuple->Fill(var); | |
1542 | ||
1543 | ||
1544 | fittedTrackTheta = GetFittedTrackTheta(); | |
1545 | fittedTrackPhi = GetFittedTrackPhi(); | |
1546 | ||
1547 | ||
1548 | ||
1549 | if(thetaCerenkov > 0.505 && thetaCerenkov < 0.605) { | |
1550 | SetPhotonEnergy(6.85); | |
1551 | SetFreonRefractiveIndex(); | |
1552 | } | |
1553 | ||
1554 | Int_t nPhotons = GetPhotonsNumber(); | |
1555 | ||
1556 | for (Int_t j=0; j < nPhotons;j++) | |
1557 | SetPhotonIndex(j); | |
1558 | }//FillHistograms() | |
1559 | //__________________________________________________________________________________________________ | |
1560 | void AliRICHRecon::Minimization() | |
1561 | { | |
1562 | // minimization to find the best theta and phi of the track | |
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 | ||
1571 | Double_t trackThetaNew,trackPhiNew; | |
1572 | TString chname; | |
1573 | Double_t eps, b1, b2; | |
1574 | Int_t ierflg; | |
1575 | ||
1576 | gAliRICHminuit = new TMinuit(2); | |
1577 | gAliRICHminuit->SetObjectFit((TObject *)this); | |
1578 | gAliRICHminuit->SetFCN(fcnrecon); | |
1579 | gAliRICHminuit->mninit(5,10,7); | |
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 | ||
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); | |
1593 | ||
1594 | arglist = -1; | |
1595 | ||
1596 | // gAliRICHminuit->FixParameter(0); | |
1597 | ||
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); | |
1602 | arglist = 1; | |
1603 | gAliRICHminuit->mnexcm("SET ERR", &arglist, 1,ierflg); | |
1604 | arglist = -1; | |
1605 | ||
1606 | // gAliRICHminuit->mnscan(); | |
1607 | ||
1608 | // gAliRICHminuit->mnexcm("SIMPLEX",&arglist, 0, ierflag); | |
1609 | gAliRICHminuit->mnexcm("MIGRAD",&arglist, 0, ierflag); | |
1610 | gAliRICHminuit->mnexcm("EXIT" ,&arglist, 0, ierflag); | |
1611 | ||
1612 | gAliRICHminuit->mnpout(0,chname, trackThetaNew, eps , b1, b2, ierflg); | |
1613 | gAliRICHminuit->mnpout(1,chname, trackPhiNew, eps , b1, b2, ierflg); | |
1614 | ||
1615 | //values after the fit... | |
1616 | SetFittedTrackTheta((Float_t)trackThetaNew); | |
1617 | SetFittedTrackPhi((Float_t)trackPhiNew); | |
1618 | ||
1619 | delete gAliRICHminuit; | |
1620 | ||
1621 | } | |
1622 | ||
1623 | void AliRICHRecon::EstimationOfTheta() | |
1624 | { | |
1625 | // theta estimate | |
1626 | ||
1627 | Int_t nPhotons = 0; | |
1628 | ||
1629 | Float_t shiftX = GetShiftX(); | |
1630 | Float_t shiftY = GetShiftY(); | |
1631 | ||
1632 | Float_t *candidatePhotonX = GetCandidatePhotonX(); | |
1633 | Float_t *candidatePhotonY = GetCandidatePhotonY(); | |
1634 | ||
1635 | Int_t nPhotonsCandidates = GetCandidatePhotonsNumber(); | |
1636 | ||
1637 | // cout << "MINIM: Nphotons " << nPhotonsCandidates << endl; | |
1638 | ||
1639 | for (Int_t j=0; j < nPhotonsCandidates; j++) | |
1640 | { | |
1641 | ||
1642 | SetPhotonIndex(j); | |
1643 | ||
1644 | if(!GetPhotonFlag()) continue; | |
1645 | ||
1646 | Float_t xtoentr = candidatePhotonX[j] - shiftX; | |
1647 | Float_t ytoentr = candidatePhotonY[j] - shiftY; | |
1648 | ||
1649 | SetEntranceX(xtoentr); | |
1650 | SetEntranceY(ytoentr); | |
1651 | ||
1652 | FindPhiPoint(); | |
1653 | ||
1654 | FindThetaPhotonCerenkov(); | |
1655 | ||
1656 | Float_t thetaPhotonCerenkov = GetThetaPhotonCerenkov(); | |
1657 | ||
1658 | // cout << " ACCEPTED!!! " << thetaPhotonCerenkov << endl; | |
1659 | ||
1660 | SetPhotonEta(thetaPhotonCerenkov); | |
1661 | ||
1662 | nPhotons++; | |
1663 | ||
1664 | } | |
1665 | ||
1666 | Float_t xmean = 0.; | |
1667 | Float_t x2mean = 0.; | |
1668 | Int_t nev = 0; | |
1669 | ||
1670 | for (Int_t j=0; j < nPhotonsCandidates;j++) | |
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 | ||
1696 | Float_t vRMS = sqrt(x2mean - xmean*xmean); | |
1697 | ||
1698 | // cout << " RMS " << vRMS; | |
1699 | ||
1700 | SetEstimationOfTheta(xmean); | |
1701 | SetEstimationOfThetaRMS(vRMS); | |
1702 | } | |
1703 | ||
1704 | void fcnrecon(Int_t& /*npar*/, Double_t* /*gin*/, Double_t &f, Double_t *par, Int_t) | |
1705 | { | |
1706 | // function to be minimized | |
1707 | AliRICHRecon *gMyRecon = (AliRICHRecon*)gAliRICHminuit->GetObjectFit(); | |
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(); | |
1716 | Float_t vRMS = gMyRecon->GetEstimationOfThetaRMS(); | |
1717 | ||
1718 | Int_t houghPhotons = gMyRecon->GetHoughPhotons(); | |
1719 | ||
1720 | ||
1721 | f = (Double_t)(1000*vRMS/(Float_t)houghPhotons); | |
1722 | ||
1723 | // if(fDebug) cout << " f " << f | |
1724 | // << " theta " << par[0] << " phi " << par[1] | |
1725 | // << " HoughPhotons " << houghPhotons << endl; | |
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 | // } | |
1733 | } | |
1734 | ||
1735 | void AliRICHRecon::Waiting() | |
1736 | { | |
1737 | // wait, wait.... | |
1738 | if(!fIsDISPLAY) return; | |
1739 | cout << " Press any key to continue..."; | |
1740 | ||
1741 | // gSystem->ProcessEvents(); | |
1742 | getchar(); | |
1743 | ||
1744 | cout << endl; | |
1745 | ||
1746 | return; | |
1747 | } | |
1748 |