1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //////////////////////////////////////////////////////////////////////////
20 // RICH class to perfom pattern recognition based on Hough transfrom //
22 //////////////////////////////////////////////////////////////////////////
26 #include "AliRICHRecon.h"
27 #include "AliRICHParam.h"
28 #include <AliLoader.h>
31 #include <Riostream.h>
32 #include <TParticle.h>
39 #include <TRotation.h>
43 #define NPointsOfRing 201
45 TMinuit *gAliRICHminuit ;
47 void fcnrecon(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
48 //__________________________________________________________________________________________________
49 AliRICHRecon::AliRICHRecon(const char*name, const char*title)
53 fThetaBin=750; fThetaMin = 0.0; fThetaMax = 0.75;
54 fDTheta = 0.001; fWindowWidth = 0.060;
55 fNpadX = AliRICHParam::NpadsY();
56 fNpadY = AliRICHParam::NpadsX();
57 fPadSizeX = AliRICHParam::PadSizeY();
58 fPadSizeY = AliRICHParam::PadSizeX();
59 fRadiatorWidth = AliRICHParam::FreonThickness();
60 fQuartzWidth = AliRICHParam::QuartzThickness();
61 fGapWidth = AliRICHParam::RadiatorToPads();
62 fXmin = -AliRICHParam::PcSizeY()/2.;
63 fXmax = AliRICHParam::PcSizeY()/2.;
64 fYmin = -AliRICHParam::PcSizeX()/2.;
65 fYmax = AliRICHParam::PcSizeX()/2.;
66 fRich = (AliRICH*)gAlice->GetDetector("RICH");
67 fOutFile=new TFile("Anal.root","RECREATE","My Analysis histos");
68 if(fIsDISPLAY) fDisplay = new TCanvas("Display","RICH Display",0,0,1200,750);
69 fNtuple=new TNtuple("hn","ntuple",
70 "Run:Trig:VertZ:Pmod:Pt:Eta:TrackTheta:TrackPhi:TrackThetaFit:TrackPhiFit:Charge::NPhotons:NPhotonsFit:InRing:MassOfParticle:HoughArea:Multiplicity:TPCLastZ");
72 //__________________________________________________________________________________________________
73 void AliRICHRecon::StartProcessEvent()
75 //start to process for pattern recognition
77 Float_t trackThetaStored = 0;
78 Float_t trackPhiStored = 0;
79 Float_t thetaCerenkovStored = 0;
80 Int_t houghPhotonsStored = 0;
82 SetFreonScaleFactor(0.994);
90 Rich()->GetLoader()->LoadHits();
91 Rich()->GetLoader()->LoadRecPoints();
92 Rich()->GetLoader()->LoadDigits();
93 gAlice->GetRunLoader()->LoadHeader();
94 gAlice->GetRunLoader()->LoadKinematics();
96 Rich()->GetLoader()->TreeR()->GetEntry(0);
98 Float_t clusX[7][500],clusY[7][500];
99 Int_t clusQ[7][500],clusMul[7][500];
102 for (Int_t ich=0;ich<7;ich++) {
103 nClusters[ich] = Rich()->ClustersOld(ich+1)->GetEntries();
104 for(Int_t k=0;k<nClusters[ich];k++) {
105 AliRICHRawCluster *pCluster = (AliRICHRawCluster *)Rich()->ClustersOld(ich+1)->At(k);
106 clusX[ich][k] = pCluster->fX;
107 clusY[ich][k] = pCluster->fY;
108 clusQ[ich][k] = pCluster->fQ;
109 clusMul[ich][k] = pCluster->fMultiplicity;
114 Int_t nPrimaries = (Int_t)Rich()->GetLoader()->TreeH()->GetEntries();
116 cout << " N. primaries " << nPrimaries << endl;
118 for(Int_t i=0;i<nPrimaries;i++){
120 Rich()->GetLoader()->TreeH()->GetEntry(i);
122 Rich()->Hits()->Print();
127 for(Int_t j=0;j<Rich()->Hits()->GetEntries();j++) {
129 pHit = (AliRICHhit*)Rich()->Hits()->At(j);
130 if(pHit->GetTrack() < nPrimaries) break;
134 cout << " iPrim " << iPrim << " pHit " << pHit << endl;
140 TParticle *pParticle = gAlice->GetRunLoader()->Stack()->Particle(pHit->GetTrack());
141 Float_t pmod = pParticle->P();
142 Float_t pt = pParticle->Pt();
143 Float_t trackEta = pParticle->Eta();
144 Int_t q = (Int_t)TMath::Sign(1.,pParticle->GetPDG()->Charge());
148 cout << " pmod " << pmod << " pt " << pt << " Eta " << trackEta << " charge " << q << endl;
150 SetTrackMomentum(pmod);
152 SetTrackEta(trackEta);
155 TVector3 pGlob(pHit->MomFreoX(),pHit->MomFreoY(),pHit->MomFreoZ());
156 TVector3 pLocal = Rich()->C(pHit->Chamber())->Glob2Loc(pGlob,1);
158 Float_t primGlobalX = pHit->X();
159 Float_t primGlobalY = pHit->Y();
160 Float_t primGlobalZ = pHit->Z();
161 TVector3 primGlobal(primGlobalX,primGlobalY,primGlobalZ);
162 TVector3 primLocal = Rich()->C(pHit->Chamber())->Glob2Loc(primGlobal);
164 // Float_t pmodFreo = pLocal.Mag();
165 Float_t trackTheta = pLocal.Theta();
166 Float_t trackPhi = pLocal.Phi();
168 // cout << " trackTheta " << trackTheta << " trackPhi " << trackPhi << endl;
170 SetTrackTheta(trackTheta);
171 SetTrackPhi(trackPhi);
174 Float_t minDist = 999.;
176 // cout << " n Clusters " << nClusters[pHit->Chamber()-1] << " for chamber n. " << pHit->Chamber() << endl;
178 for(Int_t j=0;j<nClusters[pHit->Chamber()-1];j++)
180 Float_t diffx = primLocal.X() - clusX[pHit->Chamber()-1][j];
181 Float_t diffy = primLocal.Y() - clusY[pHit->Chamber()-1][j];
184 Float_t diff = sqrt(diffx*diffx + diffy*diffy);
194 Float_t diffx = primLocal.X() - clusX[pHit->Chamber()-1][maxInd];
195 Float_t diffy = primLocal.Y() - clusY[pHit->Chamber()-1][maxInd];
197 cout << " diffx " << diffx << " diffy " << diffy << endl;
203 Float_t shiftX = primLocal.X()/primLocal.Z()*(fRadiatorWidth+fQuartzWidth+fGapWidth) + primLocal.X();
204 Float_t shiftY = primLocal.Y()/primLocal.Z()*(fRadiatorWidth+fQuartzWidth+fGapWidth) + primLocal.Y();
209 Float_t *pclusX = &clusX[pHit->Chamber()-1][0];
210 Float_t *pclusY = &clusY[pHit->Chamber()-1][0];
212 SetCandidatePhotonX(pclusX);
213 SetCandidatePhotonY(pclusY);
214 SetCandidatePhotonsNumber(nClusters[pHit->Chamber()-1]);
216 Int_t qch = clusQ[pHit->Chamber()-1][maxInd];
219 if(minDist < 3.0 && qch > 120 && maxInd !=0)
225 Float_t xrndm = fXmin + (fXmax-fXmin)*gRandom->Rndm(280964);
226 Float_t yrndm = fYmin + (fYmax-fYmin)*gRandom->Rndm(280964);
234 trackThetaStored = GetTrackTheta();
235 trackPhiStored = GetTrackPhi();
236 thetaCerenkovStored = GetThetaCerenkov();
237 houghPhotonsStored = GetHoughPhotons();
239 Int_t diffNPhotons = 999;
241 Float_t diffTrackTheta = 999.;
242 Float_t diffTrackPhi = 999.;
244 while(fIsMINIMIZER && GetHoughPhotons() > 2
246 && diffTrackTheta > 0.0001
250 Int_t houghPhotonsBefore = GetHoughPhotons();
252 Float_t trackThetaBefore = GetTrackTheta();
253 Float_t trackPhiBefore = GetTrackPhi();
259 diffNPhotons = TMath::Abs(houghPhotonsBefore - GetHoughPhotons());
261 Float_t trackThetaAfter = GetTrackTheta();
262 Float_t trackPhiAfter = GetTrackPhi();
264 diffTrackTheta = TMath::Abs(trackThetaAfter - trackThetaBefore);
265 diffTrackPhi = TMath::Abs(trackPhiAfter - trackPhiBefore);
268 cout << " houghPhotonsBefore " << houghPhotonsBefore
269 << " GetHoughPhotons() " << GetHoughPhotons();
274 SetFittedThetaCerenkov(GetThetaCerenkov());
275 SetFittedHoughPhotons(GetHoughPhotons());
277 SetTrackTheta(trackThetaStored);
278 SetTrackPhi(trackPhiStored);
279 SetThetaCerenkov(thetaCerenkovStored);
280 SetHoughPhotons(houghPhotonsStored);
286 if(fIsDISPLAY) DrawEvent(1);
292 if(fIsDISPLAY) fDisplay->Print("display.ps");
293 }//StartProcessEvent()
294 //__________________________________________________________________________________________________
295 void AliRICHRecon::EndProcessEvent()
297 // function called at the end of the event loop
302 //__________________________________________________________________________________________________
303 void AliRICHRecon::PatRec()
305 //pattern recognition method based on Hough transform
308 Float_t trackTheta = GetTrackTheta();
309 Float_t trackPhi = GetTrackPhi();
310 Float_t pmod = GetTrackMomentum();
311 Int_t iMipIndex = GetMipIndex();
313 Bool_t kPatRec = kFALSE;
315 Int_t candidatePhotons = 0;
317 Float_t shiftX = GetShiftX();
318 Float_t shiftY = GetShiftY();
320 Float_t* candidatePhotonX = GetCandidatePhotonX();
321 Float_t* candidatePhotonY = GetCandidatePhotonY();
323 Int_t candidatePhotonsNumber = GetCandidatePhotonsNumber();
325 if(fDebug) cout << " n " << candidatePhotonsNumber << endl;
327 SetThetaCerenkov(999.);
329 SetHoughPhotonsNorm(0);
332 for (Int_t j=0; j < candidatePhotonsNumber; j++)
341 if (j == iMipIndex) continue;
344 if(candidatePhotonX[j] < -64.) continue; /* avoid artificial clusters from edge uesd by Yale.... */
346 Float_t xtoentr = candidatePhotonX[j] - shiftX;
347 Float_t ytoentr = candidatePhotonY[j] - shiftY;
349 // Float_t chargehit = fHits_charge[j];
350 // if(chargehit > 150) continue;
352 SetEntranceX(xtoentr);
353 SetEntranceY(ytoentr);
357 Int_t photonStatus = PhotonInBand();
361 cout << " Photon n. " << j << " Status " << photonStatus << " accepted " << endl;
362 cout << " CandidatePhotonX[j] " << candidatePhotonX[j] << " CandidatePhotonY[j] " << candidatePhotonY[j] << endl;
365 if(photonStatus == 0) continue;
369 FindThetaPhotonCerenkov();
371 Float_t thetaPhotonCerenkov = GetThetaPhotonCerenkov();
373 if(fDebug) cout << " theta photon " << thetaPhotonCerenkov << endl;
375 SetPhotonEta(thetaPhotonCerenkov);
382 if(candidatePhotons >= 1) kPatRec = kTRUE;
386 SetThetaCerenkov(999.);
389 SetPhotonsNumber(candidatePhotonsNumber);
396 Int_t nPhotonHough = GetHoughPhotons();
400 SetThetaCerenkov(999.);
401 SetHoughPhotonsNorm(0.);
405 if(fIsWEIGHT) FindWeightThetaCerenkov();
407 Float_t thetaCerenkov = GetThetaCerenkov();
409 SetThetaOfRing(thetaCerenkov);
410 FindAreaAndPortionOfRing();
412 Float_t nPhotonHoughNorm = ((Float_t)nPhotonHough)/GetPortionOfRing();
413 SetHoughPhotonsNorm(nPhotonHoughNorm);
415 // Calculate the area where the photon are accepted...
417 Float_t thetaInternal = thetaCerenkov - 0.5*fWindowWidth;
418 SetThetaOfRing(thetaInternal);
419 FindAreaAndPortionOfRing();
420 Float_t internalArea = GetAreaOfRing();
422 Float_t thetaExternal = thetaCerenkov + 0.5*fWindowWidth;
423 SetThetaOfRing(thetaExternal);
424 FindAreaAndPortionOfRing();
425 Float_t externalArea = GetAreaOfRing();
427 Float_t houghArea = externalArea - internalArea;
429 SetHoughArea(houghArea);
433 cout << " ----- SUMMARY OF RECONSTRUCTION ----- " << endl;
434 cout << " Rings found " << fNrings << " with thetac " << thetaCerenkov << endl;
437 cout << " Nphotons " << GetPhotonsNumber()
438 << " Hough " << nPhotonHough
439 << " norm " << nPhotonHoughNorm << endl;
441 cout << " In PatRec:p " << pmod << " theta " << trackTheta << " phi " << trackPhi << endl;
442 cout << " ------------------------------------- " << endl;
445 Int_t nPhotons = GetPhotonsNumber();
451 for (Int_t j=0; j < nPhotons;j++)
455 Float_t eta = GetPhotonEta();
459 if(GetPhotonFlag() == 2)
472 xmean /=(Float_t)nev;
473 x2mean /=(Float_t)nev;
479 Float_t vRMS = sqrt(x2mean - xmean*xmean);
483 if(fDebug) cout << " RMS " << vRMS << endl;
487 void AliRICHRecon::FindEmissionPoint()
489 //estimate the emission point in radiator
491 // Find emission point
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;
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;
503 SetEmissionPoint(emissionPoint);
507 Int_t AliRICHRecon::PhotonInBand()
509 //search band fro photon candidates
511 // Float_t massOfParticle;
517 Float_t xtoentr = GetEntranceX();
518 Float_t ytoentr = GetEntranceY();
523 Float_t phpad = GetPhiPoint();
525 // Float_t pmod = GetTrackMomentum();
526 // Float_t trackTheta = GetTrackTheta();
527 // Float_t trackPhi = GetTrackPhi();
530 SetPhotonEnergy(5.6);
531 SetEmissionPoint(fRadiatorWidth -0.0001);
532 SetMassHypotesis(0.93828);
535 SetFreonRefractiveIndex();
537 beta = GetBetaOfParticle();
538 nfreon = GetFreonRefractiveIndex();
540 thetacer = Cerenkovangle(nfreon,beta);
544 if(fDebug) cout << " thetacer in photoninband min " << thetacer << endl;
546 FindThetaAtQuartz(thetacer);
548 if(thetacer == 999. || GetThetaAtQuartz() == 999.)
551 SetXInnerRing(-999.);
552 SetYInnerRing(-999.);
553 SetRadiusInnerRing(-999.);
557 SetThetaPhotonInDRS(GetThetaAtQuartz());
558 SetPhiPhotonInDRS(phpad);
560 innerRadius = FromEmissionToCathode();
561 if(innerRadius == 999.) innerRadius = -999.;
563 SetXInnerRing(GetXPointOnCathode());
564 SetYInnerRing(GetYPointOnCathode());
565 SetRadiusInnerRing(innerRadius);
569 SetPhotonEnergy(7.7);
570 SetEmissionPoint(0.);
571 // SetMassHypotesis(0.139567);
572 SetMassHypotesis(0.);
575 SetFreonRefractiveIndex();
577 beta = GetBetaOfParticle();
578 nfreon = GetFreonRefractiveIndex();
580 thetacer = Cerenkovangle(nfreon,beta);
584 if(fDebug) cout << " thetacer in photoninband max " << thetacer << endl;
586 FindThetaAtQuartz(thetacer);
588 if(thetacer == 999. || GetThetaAtQuartz() == 999.)
593 SetRadiusOuterRing(999.);
597 SetThetaPhotonInDRS(GetThetaAtQuartz());
598 SetPhiPhotonInDRS(phpad);
600 outerRadius = FromEmissionToCathode();
601 // cout << " outerRadius " << outerRadius << endl;
602 SetXOuterRing(GetXPointOnCathode());
603 SetYOuterRing(GetYPointOnCathode());
604 SetRadiusOuterRing(outerRadius);
607 Float_t padradius = sqrt(TMath::Power(xtoentr,2)+TMath::Power(ytoentr,2));
609 if(fDebug) printf(" rmin %f r %f rmax %f \n",innerRadius,padradius,outerRadius);
611 if(padradius>=innerRadius && padradius<=outerRadius) return 1;
615 void AliRICHRecon::FindThetaAtQuartz(Float_t thetaCerenkov)
617 //find the theta at the quartz plate
619 if(thetaCerenkov == 999.)
621 SetThetaAtQuartz(999.);
625 Float_t thetaAtQuartz = 999.;
627 Float_t trackTheta = GetTrackTheta();
629 if(trackTheta == 0) {
631 if(fDebug) cout << " Theta sol unique " << thetaCerenkov << endl;
633 thetaAtQuartz = thetaCerenkov;
634 SetThetaAtQuartz(thetaAtQuartz);
638 Float_t trackPhi = GetTrackPhi();
639 Float_t phiPoint = GetPhiPoint();
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;
650 Double_t underSqrt = 1 + b*b - c*c;
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;
662 if(fDebug) cout << " sqrt negative !!!!" << underSqrt << endl;
663 SetThetaAtQuartz(999.);
667 Double_t sol1 = (1+TMath::Sqrt(underSqrt))/(b-c);
668 Double_t sol2 = (1-TMath::Sqrt(underSqrt))/(b-c);
670 Double_t thetaSol1 = 2*TMath::ATan(sol1);
671 Double_t thetaSol2 = 2*TMath::ATan(sol2);
673 if(fDebug) cout << " Theta sol 1 " << thetaSol1
674 << " Theta sol 2 " << thetaSol2 << endl;
676 if(thetaSol1>0 && thetaSol1 < TMath::Pi()) thetaAtQuartz = (Float_t)thetaSol1;
677 if(thetaSol2>0 && thetaSol2 < TMath::Pi()) thetaAtQuartz = (Float_t)thetaSol2;
679 SetThetaAtQuartz(thetaAtQuartz);
682 void AliRICHRecon::FindThetaPhotonCerenkov()
684 //find theta cerenkov of ring
686 Float_t thetaCerMin = 0.;
687 Float_t thetaCerMax = 0.75;
688 Float_t thetaCerMean;
690 Float_t radiusMin, radiusMax, radiusMean;
691 Int_t nIteration = 0;
693 const Float_t kTollerance = 0.05;
695 // Float_t pmod = GetTrackMomentum();
696 // Float_t trackTheta = GetTrackTheta();
697 // Float_t trackPhi = GetTrackPhi();
699 Float_t phiPoint = GetPhiPoint();
701 SetPhotonEnergy(6.85);
702 SetEmissionPoint(fRadiatorWidth/2);
704 Float_t xPoint = GetEntranceX();
705 Float_t yPoint = GetEntranceY();
706 Float_t distPoint = sqrt(xPoint*xPoint + yPoint*yPoint);
708 if(fDebug) cout << " DistPoint " << distPoint << endl;
710 // Star minimization...
714 FindThetaAtQuartz(thetaCerMin);
716 if(GetThetaAtQuartz() == 999.)
722 SetThetaPhotonInDRS(GetThetaAtQuartz());
723 SetPhiPhotonInDRS(phiPoint);
725 radiusMin = FromEmissionToCathode();
730 FindThetaAtQuartz(thetaCerMax);
731 if(GetThetaAtQuartz() == 999.)
737 SetThetaPhotonInDRS(GetThetaAtQuartz());
738 SetPhiPhotonInDRS(phiPoint);
740 radiusMax = FromEmissionToCathode();
744 thetaCerMean = (thetaCerMax + thetaCerMin)/2;
746 FindThetaAtQuartz(thetaCerMean);
747 if(GetThetaAtQuartz() == 999.)
753 SetThetaPhotonInDRS(GetThetaAtQuartz());
754 SetPhiPhotonInDRS(phiPoint);
756 radiusMean = FromEmissionToCathode();
759 if(fDebug) cout << " r1 " << radiusMin << " rmean "
760 << radiusMean << " r2 " << radiusMax << endl;
762 while (TMath::Abs(radiusMean-distPoint) > kTollerance)
765 if((radiusMin-distPoint)*(radiusMean-distPoint) < 0) thetaCerMax = thetaCerMean;
766 if((radiusMin-distPoint)*(radiusMean-distPoint) > 0) {
768 thetaCerMin = thetaCerMean;
770 FindThetaAtQuartz(thetaCerMin);
771 SetThetaPhotonInDRS(GetThetaAtQuartz());
772 SetPhiPhotonInDRS(phiPoint);
774 radiusMin =FromEmissionToCathode();
777 thetaCerMean = (thetaCerMax + thetaCerMin)/2;
779 FindThetaAtQuartz(thetaCerMean);
780 SetThetaPhotonInDRS(GetThetaAtQuartz());
781 SetPhiPhotonInDRS(phiPoint);
783 radiusMean = FromEmissionToCathode();
787 if(fDebug) printf(" max iterations in FindPhotonCerenkov\n");
788 SetThetaPhotonCerenkov(999.);
793 SetThetaPhotonCerenkov(thetaCerMean);
797 void AliRICHRecon::FindAreaAndPortionOfRing()
799 //find fraction of the ring accepted by the RICH
801 Float_t xPoint[NPointsOfRing], yPoint[NPointsOfRing];
803 // Float_t xtoentr = GetEntranceX();
804 // Float_t ytoentr = GetEntranceY();
805 Float_t shiftX = GetShiftX();
806 Float_t shiftY = GetShiftY();
808 Float_t xemiss = GetXCoordOfEmission();
809 Float_t yemiss = GetYCoordOfEmission();
811 Float_t x0 = xemiss + shiftX;
812 Float_t y0 = yemiss + shiftY;
814 // Float_t pmod = GetTrackMomentum();
815 // Float_t trackTheta = GetTrackTheta();
816 // Float_t trackPhi = GetTrackPhi();
818 SetPhotonEnergy(6.85);
819 SetFreonRefractiveIndex();
821 SetEmissionPoint(fRadiatorWidth/2.);
823 Float_t theta = GetThetaOfRing();
826 Int_t nPsiAccepted = 0;
829 for(Int_t i=0;i<NPointsOfRing-1;i++)
832 Float_t psi = 2*TMath::Pi()*i/NPointsOfRing;
834 SetThetaPhotonInTRS(theta);
835 SetPhiPhotonInTRS(psi);
836 FindPhotonAnglesInDRS();
838 Float_t radius = FromEmissionToCathode();
839 if (radius == 999.) continue;
843 Float_t xPointRing = GetXPointOnCathode() + shiftX;
844 Float_t yPointRing = GetYPointOnCathode() + shiftY;
846 SetDetectorWhereX(xPointRing);
847 SetDetectorWhereY(yPointRing);
849 Int_t zone = CheckDetectorAcceptance();
854 FindIntersectionWithDetector();
855 xPoint[nPoints] = GetIntersectionX();
856 yPoint[nPoints] = GetIntersectionY();
860 xPoint[nPoints] = xPointRing;
861 yPoint[nPoints] = yPointRing;
869 xPoint[nPoints] = xPoint[0];
870 yPoint[nPoints] = yPoint[0];
876 for (Int_t i = 0; i < nPoints; i++)
878 area += TMath::Abs((xPoint[i]-x0)*(yPoint[i+1]-y0) - (xPoint[i+1]-x0)*(yPoint[i]-y0));
883 Float_t portionOfRing = ((Float_t)nPsiAccepted)/((Float_t)(nPsiTotal));
887 SetPortionOfRing(portionOfRing);
890 void AliRICHRecon::FindIntersectionWithDetector()
892 // find ring intersection with CsI edges
894 Float_t xIntersect, yIntersect;
895 Float_t x1, x2, y1, y2;
897 Float_t shiftX = GetShiftX();
898 Float_t shiftY = GetShiftY();
900 Float_t xPoint = GetXPointOnCathode() + shiftX;
901 Float_t yPoint = GetYPointOnCathode() + shiftY;
903 Float_t xemiss = GetXCoordOfEmission();
904 Float_t yemiss = GetYCoordOfEmission();
906 Float_t phi = GetPhiPhotonInDRS();
907 Float_t m = tan(phi);
909 Float_t x0 = xemiss + shiftX;
910 Float_t y0 = yemiss + shiftY;
934 yIntersect = m*(xIntersect - x0) + y0;
935 if (yIntersect >= fYmin && yIntersect <= fYmax && xIntersect >= x1 && xIntersect <= x2)
937 SetIntersectionX(xIntersect);
938 SetIntersectionY(yIntersect);
943 yIntersect = m*(xIntersect - x0) + y0;
944 if (yIntersect >= fYmin && yIntersect <= fYmax && xIntersect >= x1 && xIntersect <= x2)
946 SetIntersectionX(xIntersect);
947 SetIntersectionY(yIntersect);
952 xIntersect = (yIntersect - y0)/m + x0;
953 if (xIntersect >= fXmin && xIntersect <= fXmax && yIntersect >= y1 && yIntersect <= y2)
955 SetIntersectionX(xIntersect);
956 SetIntersectionY(yIntersect);
961 xIntersect = (yIntersect - y0)/m + x0;
962 if (xIntersect >= fXmin && xIntersect <= fXmax && yIntersect >= y1 && yIntersect <= y2)
964 SetIntersectionX(xIntersect);
965 SetIntersectionY(yIntersect);
969 cout << " sono fuori!!!!!!" << endl;
972 //__________________________________________________________________________________________________
973 Int_t AliRICHRecon::CheckDetectorAcceptance() const
975 // check for the acceptance
977 // crosses X -2.6 2.6 cm
980 Float_t xcoord = GetDetectorWhereX();
981 Float_t ycoord = GetDetectorWhereY();
985 if(ycoord > fYmax) return 2;
986 if(ycoord > fYmin && ycoord < fYmax) return 3;
987 if(ycoord < fYmin) return 4;
991 if(ycoord > fYmax) return 8;
992 if(ycoord > fYmin && ycoord < fYmax) return 7;
993 if(ycoord < fYmin) return 6;
995 if(xcoord > fXmin && xcoord < fXmax)
997 if(ycoord > fYmax) return 1;
998 if(ycoord > fYmin && ycoord < fYmax) return 0;
999 if(ycoord < fYmin) return 5;
1003 //__________________________________________________________________________________________________
1004 Float_t AliRICHRecon::PhotonPositionOnCathode()
1006 // find the photon position on the CsI
1007 // Float_t massOfParticle;
1012 SetPhotonEnergy(6.85);
1013 SetEmissionPoint(fRadiatorWidth/2.);
1014 SetMassHypotesis(0.139567);
1016 SetBetaOfParticle();
1017 SetFreonRefractiveIndex();
1019 beta = GetBetaOfParticle();
1020 nfreon = GetFreonRefractiveIndex();
1023 Float_t radius = FromEmissionToCathode();
1024 if (radius == 999.) return 999.;
1029 void AliRICHRecon::FindPhotonAnglesInDRS()
1031 // Setup the rotation matrix of the track...
1038 Float_t trackTheta = GetTrackTheta();
1039 Float_t trackPhi = GetTrackPhi();
1041 mtheta.RotateY(trackTheta);
1042 mphi.RotateZ(trackPhi);
1044 mrot = mphi * mtheta;
1045 // minv = mrot.Inverse();
1047 TVector3 photonInRadiator(1,1,1);
1049 Float_t thetaCerenkov = GetThetaPhotonInTRS();
1050 Float_t phiCerenkov = GetPhiPhotonInTRS();
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);
1062 Float_t AliRICHRecon::FromEmissionToCathode()
1064 // trace from emission point to cathode
1066 Float_t nfreon, nquartz, ngas;
1068 SetFreonRefractiveIndex();
1069 SetQuartzRefractiveIndex();
1070 SetGasRefractiveIndex();
1072 nfreon = GetFreonRefractiveIndex();
1073 nquartz = GetQuartzRefractiveIndex();
1074 ngas = GetGasRefractiveIndex();
1076 Float_t trackTheta = GetTrackTheta();
1077 Float_t trackPhi = GetTrackPhi();
1078 Float_t lengthOfEmissionPoint = GetEmissionPoint();
1080 Float_t theta = GetThetaPhotonInDRS();
1081 Float_t phi = GetPhiPhotonInDRS();
1083 // cout << " Theta " << Theta << " Phi " << Phi << endl;
1085 Float_t xemiss = lengthOfEmissionPoint*tan(trackTheta)*cos(trackPhi);
1086 Float_t yemiss = lengthOfEmissionPoint*tan(trackTheta)*sin(trackPhi);
1088 SetXCoordOfEmission(xemiss);
1089 SetYCoordOfEmission(yemiss);
1091 Float_t thetaquar = SnellAngle(nfreon, nquartz, theta);
1093 if(thetaquar == 999.)
1095 SetXPointOnCathode(999.);
1096 SetYPointOnCathode(999.);
1100 Float_t thetagap = SnellAngle( nquartz, ngas, thetaquar);
1102 if(thetagap == 999.)
1104 SetXPointOnCathode(999.);
1105 SetYPointOnCathode(999.);
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);
1117 Float_t xtot = xemiss + xw + xq + xg;
1118 Float_t ytot = yemiss + yw + yq + yg;
1120 SetXPointOnCathode(xtot);
1121 SetYPointOnCathode(ytot);
1124 Float_t distanceFromEntrance = sqrt(TMath::Power(fPhotonLimitX,2)
1125 +TMath::Power(fPhotonLimitY,2));
1127 return distanceFromEntrance;
1132 void AliRICHRecon::FindPhiPoint()
1134 //find phi of generated point
1136 Float_t xtoentr = GetEntranceX();
1137 Float_t ytoentr = GetEntranceY();
1139 Float_t trackTheta = GetTrackTheta();
1140 Float_t trackPhi = GetTrackPhi();
1142 Float_t emissionPoint = GetEmissionPoint();
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);
1148 SetPhiPoint(phipad);
1152 Float_t AliRICHRecon::Cerenkovangle(Float_t n, Float_t beta)
1154 // cerenkov angle from n and beta
1156 // Compute the cerenkov angle
1162 // cout << " warning in Cerenkoangle !!!!!! " << endl;
1166 thetacer = acos (1./(n*beta));
1170 Float_t AliRICHRecon::SnellAngle(Float_t n1, Float_t n2, Float_t theta1)
1174 // Compute the Snell angle
1176 Float_t sinrefractangle;
1177 Float_t refractangle;
1179 sinrefractangle = (n1/n2)*sin(theta1);
1181 if(sinrefractangle>1.) {
1182 // cout << " PROBLEMS IN SNELL ANGLE !!!!! " << endl;
1183 refractangle = 999.;
1184 return refractangle;
1187 refractangle = asin(sinrefractangle);
1188 return refractangle;
1192 void AliRICHRecon::HoughResponse()
1196 // Implement Hough response pat. rec. method
1203 int i, j, k, nCorrBand;
1204 float hcs[750],hcsw[750];
1205 float angle, weight;
1206 float lowerlimit,upperlimit;
1212 float etaPeakPos = -1;
1214 Int_t etaPeakCount = -1;
1216 Float_t thetaCerenkov = 0.;
1218 nBin = (int)(0.5+fThetaMax/(fDTheta));
1219 nCorrBand = (int)(0.5+ fWindowWidth/(2 * fDTheta));
1221 memset ((void *)hcs, 0, fThetaBin*sizeof(float));
1222 memset ((void *)hcsw, 0, fThetaBin*sizeof(float));
1224 Int_t nPhotons = GetPhotonsNumber();
1226 Int_t weightFlag = 0;
1228 for (k=0; k< nPhotons; k++) {
1232 angle = GetPhotonEta();
1234 if(angle == -999.) continue;
1236 if (angle>=fThetaMin && angle<= fThetaMax)
1240 bin = (int)(0.5+angle/(fDTheta));
1242 bin1= bin-nCorrBand;
1243 bin2= bin+nCorrBand;
1245 // calculate weights
1249 lowerlimit = ((Float_t)bin1)*fDTheta + 0.5*fDTheta;
1250 SetThetaOfRing(lowerlimit);
1251 FindAreaAndPortionOfRing();
1252 Float_t area1 = GetAreaOfRing();
1254 upperlimit = ((Float_t)bin2)*fDTheta + 0.5*fDTheta;
1255 SetThetaOfRing(upperlimit);
1256 FindAreaAndPortionOfRing();
1257 Float_t area2 = GetAreaOfRing();
1259 // cout << "lowerlimit" << lowerlimit << "upperlimit " << upperlimit << endl;
1260 Float_t diffarea = area2 - area1;
1264 weight = 1./(area2-area1);
1272 // cout <<" low "<< lowerlimit << " up " << upperlimit <<
1273 // " area1 " << area1 << " area2 " << area2 << " weight " << weight << endl;
1281 SetPhotonWeight(weight);
1283 // cout << "weight..." << weight << endl;
1287 if (bin2>nBin) bin2=nBin;
1289 for (j=bin1; j<bin2; j++)
1305 // cout << " probems with weight...normal procedure adopted " << endl;
1308 HoughFiltering(hCSspace);
1310 for (bin=0; bin <nBin; bin++) {
1311 angle = (bin+0.5) * (fDTheta);
1312 if (hCSspace[bin] && hCSspace[bin] > etaPeakPos) {
1314 etaPeakPos = hCSspace[bin];
1318 if (hCSspace[bin] == etaPeakPos) {
1319 etaPeak[++etaPeakCount] = angle;
1324 for (i=0; i<etaPeakCount+1; i++) {
1325 thetaCerenkov += etaPeak[i];
1327 if (etaPeakCount>=0) {
1328 thetaCerenkov /= etaPeakCount+1;
1329 fThetaPeakPos = etaPeakPos;
1332 SetThetaCerenkov(thetaCerenkov);
1336 void AliRICHRecon::HoughFiltering(float hcs[])
1343 float k[5] = {0.05, 0.25, 0.4, 0.25, 0.05};
1348 nBin = (int)(1+fThetaMax/fDTheta);
1349 sizeHCS = fThetaBin*sizeof(float);
1351 memset ((void *)hcsFilt, 0, sizeHCS);
1353 for (nx = 0; nx < nBin; nx++) {
1354 for (i = 0; i < 5; i++) {
1356 if (nxDx> -1 && nxDx<nBin)
1357 hcsFilt[nx] += hcs[nxDx] * k[i];
1361 for (nx = 0; nx < nBin; nx++) {
1362 hcs[nx] = hcsFilt[nx];
1366 void AliRICHRecon::FindWeightThetaCerenkov()
1368 // manage with weight for photons
1371 Float_t weightThetaCerenkov = 0.;
1373 Int_t nPhotons = GetPhotonsNumber();
1374 for(Int_t i=0;i<nPhotons;i++)
1378 if(GetPhotonFlag() == 2)
1380 Float_t photonEta = GetPhotonEta();
1381 Float_t photonWeight = GetPhotonWeight();
1382 weightThetaCerenkov += photonEta*photonWeight;
1383 wei += photonWeight;
1389 weightThetaCerenkov /= wei;
1393 weightThetaCerenkov = 0.;
1396 SetThetaCerenkov(weightThetaCerenkov);
1398 cout << " thetac weighted -> " << weightThetaCerenkov << endl;
1402 void AliRICHRecon::FlagPhotons()
1406 Int_t nPhotonHough = 0;
1408 Float_t thetaCerenkov = GetThetaCerenkov();
1409 if(fDebug) cout << " fThetaCerenkov " << thetaCerenkov << endl;
1411 Float_t thetaDist= thetaCerenkov - fThetaMin;
1412 Int_t steps = (Int_t)(thetaDist / fDTheta);
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);
1418 tmin = tavg - 0.5*fWindowWidth;
1419 tmax = tavg + 0.5*fWindowWidth;
1421 if(fDebug) cout << " tmin " << tmin << " tmax " << tmax << endl;
1422 if(fDebug) cout << " thetac " << thetaCerenkov << endl;
1424 // Int_t candidatePhotonsNumber = GetCandidatePhotonsNumber();
1426 Int_t nPhotons = GetPhotonsNumber();
1428 // for(Int_t i=0;i<candidatePhotonsNumber;i++)
1430 for(Int_t i=0;i<nPhotons;i++)
1434 Float_t photonEta = GetPhotonEta();
1436 if(photonEta == -999.) continue;
1438 if(photonEta >= tmin && photonEta <= tmax)
1444 SetHoughPhotons(nPhotonHough);
1447 void AliRICHRecon::DrawEvent(Int_t flag) const
1449 // draw event with rings
1451 flag=1; // dummy to be removed...
1454 Float_t AliRICHRecon::FindMassOfParticle()
1456 // find mass of the particle from theta cerenkov
1458 Float_t pmod = GetTrackMomentum();
1460 SetPhotonEnergy(6.85);
1461 SetFreonRefractiveIndex();
1463 Float_t thetaCerenkov = GetThetaCerenkov();
1464 FindBetaFromTheta(thetaCerenkov);
1466 Double_t beta = (Double_t)(GetBetaOfParticle());
1467 Double_t den = 1. - beta*beta;
1468 if(den<=0.) return 999.;
1470 Double_t gamma = 1./TMath::Sqrt(den);
1472 Float_t mass = pmod/(beta*(Float_t)gamma);
1478 void AliRICHRecon::FillHistograms()
1480 // fill histograms..
1482 Float_t fittedTrackTheta, fittedTrackPhi;
1484 Float_t thetaCerenkov = GetThetaCerenkov();
1485 if(thetaCerenkov == 999.) return;
1487 Float_t vertZ = GetEventVertexZ();
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();
1498 fittedTrackTheta = GetFittedTrackTheta();
1499 fittedTrackPhi = GetFittedTrackPhi();
1500 Int_t fittednPhotonHough = GetFittedHoughPhotons();
1504 cout << " p " << pmod << " ThetaC " << thetaCerenkov
1505 << " rings " << fNrings << endl;
1508 Int_t nPhotonHough = GetHoughPhotons();
1509 // Float_t nPhotonHoughNorm = GetHoughPhotonsNorm();
1510 Float_t inRing = GetPortionOfRing();
1512 Float_t massOfParticle = FindMassOfParticle();
1514 Float_t houghArea = GetHoughArea();
1515 Float_t multiplicity = GetEventMultiplicity();
1526 var[6] = trackTheta;
1528 var[8] = fittedTrackTheta;
1529 var[9] = fittedTrackPhi;
1531 var[11] = thetaCerenkov;
1532 var[12] = (Float_t)nPhotonHough;
1533 var[13] = (Float_t)fittednPhotonHough;
1535 var[15] = massOfParticle;
1536 var[16] = houghArea;
1537 var[17] = multiplicity;
1544 fittedTrackTheta = GetFittedTrackTheta();
1545 fittedTrackPhi = GetFittedTrackPhi();
1549 if(thetaCerenkov > 0.505 && thetaCerenkov < 0.605) {
1550 SetPhotonEnergy(6.85);
1551 SetFreonRefractiveIndex();
1554 Int_t nPhotons = GetPhotonsNumber();
1556 for (Int_t j=0; j < nPhotons;j++)
1559 //__________________________________________________________________________________________________
1560 void AliRICHRecon::Minimization()
1562 // minimization to find the best theta and phi of the track
1567 static Double_t vstart[2];
1568 static Double_t lower[2], upper[2];
1569 static Double_t step[2]={0.001,0.001};
1571 Double_t trackThetaNew,trackPhiNew;
1573 Double_t eps, b1, b2;
1576 gAliRICHminuit = new TMinuit(2);
1577 gAliRICHminuit->SetObjectFit((TObject *)this);
1578 gAliRICHminuit->SetFCN(fcnrecon);
1579 gAliRICHminuit->mninit(5,10,7);
1581 vstart[0] = (Double_t)GetTrackTheta();
1582 vstart[1] = (Double_t)GetTrackPhi();
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;
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);
1596 // gAliRICHminuit->FixParameter(0);
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);
1603 gAliRICHminuit->mnexcm("SET ERR", &arglist, 1,ierflg);
1606 // gAliRICHminuit->mnscan();
1608 // gAliRICHminuit->mnexcm("SIMPLEX",&arglist, 0, ierflag);
1609 gAliRICHminuit->mnexcm("MIGRAD",&arglist, 0, ierflag);
1610 gAliRICHminuit->mnexcm("EXIT" ,&arglist, 0, ierflag);
1612 gAliRICHminuit->mnpout(0,chname, trackThetaNew, eps , b1, b2, ierflg);
1613 gAliRICHminuit->mnpout(1,chname, trackPhiNew, eps , b1, b2, ierflg);
1615 //values after the fit...
1616 SetFittedTrackTheta((Float_t)trackThetaNew);
1617 SetFittedTrackPhi((Float_t)trackPhiNew);
1619 delete gAliRICHminuit;
1623 void AliRICHRecon::EstimationOfTheta()
1629 Float_t shiftX = GetShiftX();
1630 Float_t shiftY = GetShiftY();
1632 Float_t *candidatePhotonX = GetCandidatePhotonX();
1633 Float_t *candidatePhotonY = GetCandidatePhotonY();
1635 Int_t nPhotonsCandidates = GetCandidatePhotonsNumber();
1637 // cout << "MINIM: Nphotons " << nPhotonsCandidates << endl;
1639 for (Int_t j=0; j < nPhotonsCandidates; j++)
1644 if(!GetPhotonFlag()) continue;
1646 Float_t xtoentr = candidatePhotonX[j] - shiftX;
1647 Float_t ytoentr = candidatePhotonY[j] - shiftY;
1649 SetEntranceX(xtoentr);
1650 SetEntranceY(ytoentr);
1654 FindThetaPhotonCerenkov();
1656 Float_t thetaPhotonCerenkov = GetThetaPhotonCerenkov();
1658 // cout << " ACCEPTED!!! " << thetaPhotonCerenkov << endl;
1660 SetPhotonEta(thetaPhotonCerenkov);
1667 Float_t x2mean = 0.;
1670 for (Int_t j=0; j < nPhotonsCandidates;j++)
1674 Float_t eta = GetPhotonEta();
1678 if(GetPhotonFlag() == 2)
1689 xmean /=(Float_t)nev;
1690 x2mean /=(Float_t)nev;
1696 Float_t vRMS = sqrt(x2mean - xmean*xmean);
1698 // cout << " RMS " << vRMS;
1700 SetEstimationOfTheta(xmean);
1701 SetEstimationOfThetaRMS(vRMS);
1704 void fcnrecon(Int_t& /*npar*/, Double_t* /*gin*/, Double_t &f, Double_t *par, Int_t)
1706 // function to be minimized
1707 AliRICHRecon *gMyRecon = (AliRICHRecon*)gAliRICHminuit->GetObjectFit();
1709 Float_t p0 = (Float_t)par[0];
1710 Float_t p1 = (Float_t)par[1];
1712 gMyRecon->SetTrackTheta(p0);
1713 gMyRecon->SetTrackPhi(p1);
1715 gMyRecon->EstimationOfTheta();
1716 Float_t vRMS = gMyRecon->GetEstimationOfThetaRMS();
1718 Int_t houghPhotons = gMyRecon->GetHoughPhotons();
1721 f = (Double_t)(1000*vRMS/(Float_t)houghPhotons);
1723 // if(fDebug) cout << " f " << f
1724 // << " theta " << par[0] << " phi " << par[1]
1725 // << " HoughPhotons " << houghPhotons << endl;
1727 // if(fDebug&&iflag == 3)
1729 // cout << " --- end convergence...summary --- " << endl;
1730 // cout << " theta " << par[0] << endl;
1731 // cout << " phi " << par[1] << endl;
1735 void AliRICHRecon::Waiting()
1738 if(!fIsDISPLAY) return;
1739 cout << " Press any key to continue...";
1741 // gSystem->ProcessEvents();