Fixing Coding violations
[u/mrichter/AliRoot.git] / RICH / AliRICHRecon.cxx
CommitLineData
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 47TMinuit *gAliRICHminuit ;
5cb4dfc3 48
8cc78f7f 49void fcnrecon(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
b068561d 50//__________________________________________________________________________________________________
51AliRICHRecon::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 75void 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 295void 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 303void 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
487void 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
507Int_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 615void 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
682void 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
797void 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
890void 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 973Int_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 1004Float_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
1029void 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
1062Float_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
1132void 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
1152Float_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
1170Float_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
1192void 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
1336void 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
1366void 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
1402void 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 1447void 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
1454Float_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
1478void 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 1560void 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
1623void 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 1704void 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 1735void 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