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 **************************************************************************/
20 #include "DataStructures.h"
22 #include "AliDetector.h"
24 #include "AliRICHPoints.h"
25 #include "AliRICHSegResV0.h"
26 #include "AliRICHPatRec.h"
28 #include "AliRICHConst.h"
29 #include "AliRICHPoints.h"
31 #include "TParticle.h"
38 ClassImp(AliRICHPatRec)
39 //___________________________________________
40 AliRICHPatRec::AliRICHPatRec() : TObject()
44 //___________________________________________
45 AliRICHPatRec::AliRICHPatRec(const char *name, const char *title)
51 void AliRICHPatRec::PatRec()
54 AliRICHChamber* iChamber;
55 AliRICHSegmentation* segmentation;
57 Int_t ntracks, ndigits[7];
64 Float_t gamma,MassCer,BetaCer;
68 printf("PatRec started\n");
70 TCanvas *c1 = new TCanvas("c1","Alice RICH pad hits",50,10,700,700);
72 TH2F *ring = new TH2F("ring","map",90,-30.,30.,90,-30.,30.);
73 TH2F *ringband = new TH2F("ringband","map",90,-65.,65.,90,-65.,65.);
74 TH1F *cerangle = new TH1F("cerangle","phot",300,0.45,0.75);
75 TH1F *ceranglew= new TH1F("ceranglew","phot",300,0.45,0.75);
76 TH1F *hough = new TH1F("hough","hough",75,0.45,0.75);
77 TH1F *mass = new TH1F("mass","mass",100,50.,1050.);
79 AliRICH *RICH = (AliRICH*)gAlice->GetDetector("RICH");
80 TTree *TH = gAlice->TreeH();
82 ntracks =(Int_t) TH->GetEntries();
84 for (itr=0; itr<ntracks; itr++) {
86 status = TrackParam(itr,ich);
87 if(status==1) continue;
88 printf(" theta %f phi %f track \n",fTrackTheta,fTrackPhi);
89 // ring->Fill(fTrackLoc[0],fTrackLoc[1],100.);
91 iChamber = &(RICH->Chamber(ich));
92 segmentation=iChamber->GetSegmentationModel();
94 nent=(Int_t)gAlice->TreeD()->GetEntries();
95 gAlice->TreeD()->GetEvent(nent-1);
96 TClonesArray *Digits = RICH->DigitsAddress(ich);
97 ndigits[ich] = Digits->GetEntriesFast();
98 printf("ndigits %d in chamber %d\n",ndigits[ich],ich);
99 AliRICHDigit *padI = 0;
103 for (Int_t dig=0;dig<ndigits[ich];dig++) {
104 padI=(AliRICHDigit*) Digits->UncheckedAt(dig);
108 segmentation->GetPadCxy(x,y,rx,ry);
114 ringband->Fill(x,y,1.);
115 fCerenkovAnglePad = PhotonCerenkovAngle();
116 if(fCerenkovAnglePad==-999) continue;
118 if(!PhotonInBand()) continue;
120 ring->Fill(fXpad,fYpad,1.);
121 cerangle->Fill(fCerenkovAnglePad,1.);
124 fEtaPhotons[GoodPhotons] = fCerenkovAnglePad;
126 fNumEtaPhotons = GoodPhotons;
128 BackgroundEstimation();
130 for(i=0;i<GoodPhotons;i++) {
131 ceranglew->Fill(fEtaPhotons[i],fWeightPhotons[i]);
132 // printf(" Eta %f weight %f \n",fEtaPhotons[i],fWeightPhotons[i]);
139 rechit[2] = fThetaCerenkov;
144 hough->Fill(fThetaCerenkov,1.);
146 RICH->AddRecHit(ich,rechit);
148 BetaCer = BetaCerenkov(1.29,fThetaCerenkov);
149 gamma = 1./sqrt(1.-pow(BetaCer,2));
150 MassCer = fTrackMom/(BetaCer*gamma);
151 // printf(" mass %f \n",MassCer);
152 mass->Fill(MassCer*1000,1.);
155 gAlice->TreeR()->Fill();
158 fRec=RICH->RecHitsAddress(i);
159 int ndig=fRec->GetEntriesFast();
160 printf ("Chamber %d, rings %d\n",i,ndig);
162 RICH->ResetRecHits();
169 ringband->Draw("box");
178 Int_t AliRICHPatRec::TrackParam(Int_t itr, Int_t &ich)
180 // Get Local coordinates of track impact
182 AliRICHChamber* iChamber;
183 AliRICHSegmentation* segmentation;
185 Float_t trackglob[3];
193 printf("Calling TrackParam\n");
196 TTree *TH = gAlice->TreeH();
199 AliRICH *RICH = (AliRICH*)gAlice->GetDetector("RICH");
200 AliRICHHit* mHit=(AliRICHHit*)RICH->FirstHit(-1);
201 if(mHit==0) return 1;
202 ich = mHit->fChamber-1;
203 trackglob[0] = mHit->fX;
204 trackglob[1] = mHit->fY;
205 trackglob[2] = mHit->fZ;
209 fTrackMom = sqrt(pow(pX,2)+pow(pY,2)+pow(pZ,2));
210 thetatr = (180 - mHit->fTheta)*(Float_t)kDegrad;
211 phitr = mHit->fPhi*(Float_t)kDegrad;
213 part = mHit->fParticle;
215 iChamber = &(RICH->Chamber(ich));
216 iChamber->GlobaltoLocal(trackglob,trackloc);
218 segmentation=iChamber->GetSegmentationModel();
220 // retrieve geometrical params
222 AliRICHGeometry* fGeometry=iChamber->GetGeometryModel();
224 fRw = fGeometry->GetFreonThickness();
225 fQw = fGeometry->GetQuartzThickness();
226 fTgap = fGeometry->GetGapThickness()
227 + fGeometry->GetProximityGapThickness();
229 Float_t apar = (fRw + fQw + fTgap)*tan(thetatr);
230 fTrackLoc[0] = apar*cos(phitr);
231 fTrackLoc[1] = apar*sin(phitr);
232 fTrackLoc[2] = fRw + fQw + fTgap;
233 fTrackTheta = thetatr;
236 fXshift = trackloc[0] - fTrackLoc[0];
237 fYshift = trackloc[2] - fTrackLoc[1];
242 Float_t AliRICHPatRec::EstimationAtLimits(Float_t lim, Float_t radius,
245 Float_t nquartz = 1.585;
247 Float_t nfreon = 1.295;
250 // printf("Calling EstimationLimits\n");
252 Float_t apar = (fRw -fEmissPoint + fQw + fTgap)*tan(fTrackTheta);
253 Float_t b1 = (fRw-fEmissPoint)*tan(lim);
254 Float_t b2 = fQw / sqrt(pow(nquartz,2)-pow(nfreon*sin(lim),2));
255 Float_t b3 = fTgap / sqrt(pow(ngas,2)-pow(nfreon*sin(lim),2));
256 Float_t bpar = b1 + nfreon*sin(lim)*(b2+b3);
257 value = pow(radius,2)
258 -pow((apar*cos(fTrackPhi)-bpar*cos(phiphot)),2)
259 -pow((apar*sin(fTrackPhi)-bpar*sin(phiphot)),2);
264 Float_t AliRICHPatRec::PhotonCerenkovAngle()
266 // Cherenkov pad angle reconstruction
270 Float_t cherMax = 0.8;
272 Float_t eps = 0.0001;
273 Int_t niterEmiss = 0;
274 Int_t niterEmissMax = 0;
275 Float_t x1,x2,x3,p1,p2,p3;
279 // printf("Calling PhotonCerenkovAngle\n");
281 radius = sqrt(pow(fTrackLoc[0]-fXpad,2)+pow(fTrackLoc[1]-fYpad,2));
282 fEmissPoint = fRw/2.; //Start value of EmissionPoint
284 while(niterEmiss<=niterEmissMax) {
287 argY = fYpad - fEmissPoint*tan(fTrackTheta)*sin(fTrackPhi);
288 argX = fXpad - fEmissPoint*tan(fTrackTheta)*cos(fTrackPhi);
289 phiphot = atan2(argY,argX);
290 p1 = EstimationAtLimits(cherMin,radius,phiphot);
291 p2 = EstimationAtLimits(cherMax,radius,phiphot);
294 // printf("PhotonCerenkovAngle failed\n");
298 //start to find the Cherenkov pad angle
302 p3 = EstimationAtLimits(x3,radius,phiphot);
303 while(TMath::Abs(p3)>eps){
307 p1 = EstimationAtLimits(x1,radius,phiphot);
310 p3 = EstimationAtLimits(x3,radius,phiphot);
314 // printf(" max iterations in PhotonCerenkovAngle\n");
318 // printf("niterFun %i \n",niterFun);
320 if (niterEmiss != niterEmissMax+1) EmissionPoint();
323 printf(" phiphot %f fXpad %f fYpad %f fEmiss %f \n",
324 phiphot,fXpad,fYpad,fEmissPoint);
332 void AliRICHPatRec::EmissionPoint()
334 Float_t AbsorLength=7.83*fRw; //absorption length in the freon (cm)
335 // 7.83 = -1/ln(T0) where
336 // T0->Trasmission freon at 180nm = 0.88 (Eph=6.85eV)
337 Float_t PhotonLength, PhotonLengthMin, PhotonLengthMax;
339 PhotonLength=exp(-fRw/(AbsorLength*cos(fCerenkovAnglePad)));
340 PhotonLengthMin=fRw*PhotonLength/(1.-PhotonLength);
341 PhotonLengthMax=AbsorLength*cos(fCerenkovAnglePad);
342 fEmissPoint = fRw + PhotonLengthMin - PhotonLengthMax;
346 void AliRICHPatRec::PhotonSelection(Int_t track, Int_t &nphot, Float_t &thetamean)
348 printf("Calling PhotonSelection\n");
351 void AliRICHPatRec::BackgroundEstimation()
353 Float_t StepEta = 0.001;
354 Float_t EtaMinBkg = 0.72;
355 Float_t EtaMaxBkg = 0.75;
357 Float_t EtaMax = 0.75;
359 Float_t nfreon = 1.295;
361 Float_t EtaStepMin,EtaStepMax,EtaStepAvg;
363 Int_t NumPhotBkg, NumPhotonStep;
364 Float_t FunBkg,AreaBkg,NormBkg;
365 Float_t DensityBkg,StoreBkg,NumStore;
371 nstep = (int)((EtaMaxBkg-EtaMinBkg)/StepEta);
373 for (i=0;i<fNumEtaPhotons;i++) {
375 if(fEtaPhotons[i]>EtaMinBkg && fEtaPhotons[i]<EtaMaxBkg) {
379 if (NumPhotBkg == 0) {
380 for (i=0;i<fNumEtaPhotons;i++) {
381 fWeightPhotons[i] = 1.;
386 // printf(" NumPhotBkg %i ",NumPhotBkg);
388 for (i=0;i<nstep;i++) {
389 EtaStepMin = EtaMinBkg + (Float_t)(i)*StepEta;
390 EtaStepMax = EtaMinBkg + (Float_t)(i+1)*StepEta;
391 EtaStepAvg = 0.5*(EtaStepMax + EtaStepMin);
393 FunBkg = tan(EtaStepAvg)*pow((1.+pow(tan(EtaStepAvg),2)),
394 5.52)-7.803 + 22.02*tan(EtaStepAvg);
396 ThetaSig = asin(nfreon/ngas*sin(EtaStepAvg));
397 FunBkg = tan(ThetaSig)*(1.+pow(tan(ThetaSig),2))*nfreon
398 /ngas*cos(EtaStepAvg)/cos(ThetaSig);
399 AreaBkg += StepEta*FunBkg;
402 DensityBkg = 0.95*(Float_t)(NumPhotBkg)/AreaBkg;
403 // printf(" DensityBkg %f \n",DensityBkg);
405 nstep = (int)((EtaMax-EtaMin)/StepEta);
408 for (i=0;i<nstep;i++) {
409 EtaStepMin = EtaMinBkg + (Float_t)(i)*StepEta;
410 EtaStepMax = EtaMinBkg + (Float_t)(i+1)*StepEta;
411 EtaStepAvg = 0.5*(EtaStepMax + EtaStepMin);
413 FunBkg = tan(EtaStepAvg)*pow((1.+pow(tan(EtaStepAvg),2)),
414 5.52)-7.803 + 22.02*tan(EtaStepAvg);
417 ThetaSig = asin(nfreon/ngas*sin(EtaStepAvg));
418 FunBkg = tan(ThetaSig)*(1.+pow(tan(ThetaSig),2))*nfreon
419 /ngas*cos(EtaStepAvg)/cos(ThetaSig);
421 AreaBkg = StepEta*FunBkg;
422 NormBkg = DensityBkg*AreaBkg;
424 for (ip=0;ip<fNumEtaPhotons;ip++) {
425 if(fEtaPhotons[ip]>EtaStepMin && fEtaPhotons[ip]<EtaStepMax) {
429 if (NumPhotonStep == 0) {
437 if (NumPhotonStep == 0) continue;
438 for (ip=0;ip<fNumEtaPhotons;ip++) {
439 if(fEtaPhotons[ip]>EtaStepMin && fEtaPhotons[ip]<EtaStepMax) {
443 fWeightPhotons[ip] = 1. - NormBkg/(Float_t)(NumPhotonStep);
445 printf(" NormBkg %f NumPhotonStep %i fW %f \n",
446 NormBkg, NumPhotonStep, fWeightPhotons[ip]);
448 if(fWeightPhotons[ip]<0) fWeightPhotons[ip] = 0.;
455 void AliRICHPatRec::FlagPhotons(Int_t track, Float_t theta)
457 printf("Calling FlagPhotons\n");
461 //////////////////////////////////////////
467 Int_t AliRICHPatRec::PhotonInBand()
469 //0=label for parameters giving internal band ellipse
470 //1=label for parameters giving external band ellipse
472 Float_t imp[2], mass[2], Energ[2], beta[2];
473 Float_t EmissPointLength[2];
474 Float_t E1, E2, F1, F2;
475 Float_t nfreon[2], nquartz[2];
479 Float_t phpad, thetacer[2];
480 Float_t bandradius[2], padradius;
482 imp[0] = 5.0; //threshold momentum for the proton Cherenkov emission
485 mass[0] = 0.938; //proton mass
486 mass[1] = 0.139; //pion mass
488 EmissPointLength[0] = fRw-0.0001; //at the beginning of the radiator
489 EmissPointLength[1] = 0.;//at the end of radiator
491 //parameters to calculate freon window refractive index vs. energy
495 //parameters to calculate quartz window refractive index vs. energy
510 for (times=0; times<=1; times++) {
512 nfreon[times] = a+b*Energ[times];
514 nquartz[times] = sqrt(1+(F1/(pow(E1,2)-pow(Energ[times],2)))+
515 (F2/(pow(E2,2)-pow(Energ[times],2))));
517 beta[times] = imp[times]/sqrt(pow(imp[times],2)+pow(mass[times],2));
519 thetacer[times] = CherenkovAngle( nfreon[times], beta[times]);
521 bandradius[times] = DistanceFromMip( nfreon[times], nquartz[times],
522 EmissPointLength[times],
523 thetacer[times], phpad);
526 bandradius[0] -= 1.6;
527 bandradius[1] += 1.6;
528 padradius = sqrt(pow(fXpad,2)+pow(fYpad,2));
529 // printf(" rmin %f r %f rmax %f \n",bandradius[0],padradius,bandradius[1]);
531 if(padradius>=bandradius[0] && padradius<=bandradius[1]) return 1;
535 Float_t AliRICHPatRec::DistanceFromMip(Float_t nfreon, Float_t nquartz,
536 Float_t EmissPointLength, Float_t thetacer,
539 Float_t DistanceValue;
541 TVector3 RadExitPhot(1,1,1);//photon impact at the radiator exit with respect
542 //to local reference sistem with the origin in the MIP entrance
544 TVector3 VectEmissPointLength(1,1,1);
545 Float_t MagEmissPointLenght;
547 TVector3 RadExitPhot2(1,1,1);//photon impact at the radiator exit with respect
548 Float_t MagRadExitPhot2;
549 //to a reference sistem with origin in the photon emission point and
550 //axes parallel to the MIP reference sistem
552 TVector3 QuarExitPhot(1,1,1);//photon impact at the quartz exit with respect
553 Float_t MagQuarExitPhot;
555 TVector3 GapExitPhot(1,1,1) ;
556 Float_t MagGapExitPhot;
558 TVector3 fPhotocatExitPhot(1,1,1);
560 Double_t thetarad , phirad ;
561 Double_t thetaquar, phiquar;
562 Double_t thetagap , phigap ;
566 MagEmissPointLenght = EmissPointLength/cos(fTrackTheta);
568 VectEmissPointLength.SetMag(MagEmissPointLenght);
569 VectEmissPointLength.SetTheta(fTrackTheta);
570 VectEmissPointLength.SetPhi(fTrackPhi);
573 RadExitPhot2.SetTheta(thetacer);
574 RadExitPhot2.SetPhi(phpad);
581 r1. RotateY(fTrackTheta);
582 r2. RotateZ(fTrackPhi);
586 r = r2 * r1;//rotation about the y axis by MIP theta incidence angle
587 //following by a rotation about the z axis by MIP phi incidence angle;
590 RadExitPhot2 = r * RadExitPhot2;
591 theta2 = RadExitPhot2.Theta();
592 MagRadExitPhot2 = (fRw - VectEmissPointLength(2))/cos(theta2);
593 RadExitPhot2.SetMag(MagRadExitPhot2);
596 RadExitPhot = VectEmissPointLength + RadExitPhot2;
597 thetarad = RadExitPhot.Theta();
599 phirad = RadExitPhot.Phi(); //check on the original file //
601 thetaquar = SnellAngle( nfreon, nquartz, theta2);
602 phiquar = RadExitPhot2.Phi();
603 if(thetaquar == 999.) return thetaquar;
604 MagQuarExitPhot = fQw/cos(thetaquar);
605 QuarExitPhot.SetMag( MagQuarExitPhot);
606 QuarExitPhot.SetTheta(thetaquar);
607 QuarExitPhot.SetPhi(phiquar);
609 thetagap = SnellAngle( nquartz, ngas, thetaquar);
611 if(thetagap == 999.) return thetagap;
612 MagGapExitPhot = fTgap/cos(thetagap);
613 GapExitPhot.SetMag( MagGapExitPhot);
614 GapExitPhot.SetTheta(thetagap);
615 GapExitPhot.SetPhi(phigap);
617 fPhotocatExitPhot = RadExitPhot + QuarExitPhot + GapExitPhot;
619 DistanceValue = sqrt(pow(fPhotocatExitPhot(0),2)
620 +pow(fPhotocatExitPhot(1),2));
621 return DistanceValue ;
624 Float_t AliRICHPatRec::PhiPad()
627 Float_t thetapad, phipad;
628 Float_t thetarot, phirot;
630 zpad = fRw + fQw + fTgap;
632 TVector3 PhotonPad(fXpad, fYpad, zpad);
633 thetapad = PhotonPad.Theta();
634 phipad = PhotonPad.Phi();
640 thetarot = - fTrackTheta;
641 phirot = - fTrackPhi;
643 r2. RotateY(thetarot);
645 r = r2 * r1;//rotation about the z axis by MIP -phi incidence angle
646 //following by a rotation about the y axis by MIP -theta incidence angle;
648 PhotonPad = r * PhotonPad;
650 phipad = PhotonPad.Phi();
655 Float_t AliRICHPatRec:: SnellAngle(Float_t n1, Float_t n2, Float_t theta1)
657 Float_t sinrefractangle;
658 Float_t refractangle;
660 sinrefractangle = (n1/n2)*sin(theta1);
662 if(sinrefractangle>1.) {
667 refractangle = asin(sinrefractangle);
671 Float_t AliRICHPatRec::CherenkovAngle(Float_t n, Float_t beta)
680 thetacer = acos (1./(n*beta));
684 Float_t AliRICHPatRec::BetaCerenkov(Float_t n, Float_t theta)
688 beta = 1./(n*cos(theta));
695 void AliRICHPatRec::HoughResponse()
701 int i, j, k, NcorrBand;
704 float angle, ThetaCerMean;
709 float StepEta = 0.001;
710 float WindowEta = 0.040;
714 float EtaPeakPos = -1;
715 Int_t EtaPeakCount = -1;
720 Nbin = (int)(0.5+EtaMax/(StepEta));
721 NcorrBand = (int)(0.5+ WindowEta/(2 * StepEta));
722 memset ((void *)HCS, 0, EtaBin*sizeof(int));
724 for (k=0; k< fNumEtaPhotons; k++) {
726 angle = fEtaPhotons[k];
728 if (angle>=EtaMin && angle<= EtaMax) {
729 bin = (int)(0.5+angle/(StepEta));
733 if (bin2>Nbin) bin2=Nbin;
735 for (j=bin1; j<bin2; j++) {
736 HCS[j] += fWeightPhotons[k];
739 ThetaCerMean += angle;
743 ThetaCerMean /= fNumEtaPhotons;
747 for (bin=0; bin <Nbin; bin++) {
748 angle = (bin+0.5) * (StepEta);
749 if (HCS[bin] && HCS[bin] > EtaPeakPos) {
751 EtaPeakPos = HCS[bin];
755 if (HCS[bin] == EtaPeakPos) {
756 EtaPeak[++EtaPeakCount] = angle;
761 for (i=0; i<EtaPeakCount+1; i++) {
762 fThetaCerenkov += EtaPeak[i];
764 if (EtaPeakCount>=0) {
765 fThetaCerenkov /= EtaPeakCount+1;
766 fThetaPeakPos = EtaPeakPos;
771 void AliRICHPatRec::HoughFiltering(float HCS[])
774 float K[5] = {0.05, 0.25, 0.4, 0.25, 0.05};
781 float StepEta = 0.001;
783 Nbin = (int)(1+EtaMax/StepEta);
784 sizeHCS = EtaBin*sizeof(float);
786 memset ((void *)HCS_filt, 0, sizeHCS);
788 for (nx = 0; nx < Nbin; nx++) {
789 for (i = 0; i < 5; i++) {
791 if (nx_dx> -1 && nx_dx<Nbin)
792 HCS_filt[nx] += HCS[nx_dx] * K[i];
796 for (nx = 0; nx < Nbin; nx++) {
797 HCS[nx] = HCS_filt[nx];
801 Float_t AliRICHPatRec::CherenkovRingDrawing(Float_t fixedthetacer)
805 //to draw Cherenkov ring by known Cherenkov angle
807 Int_t nmaxdegrees, nstepdegrees;
808 Float_t phpad, thetacer;
809 Float_t nfreonave, nquartzave;
812 Float_t E1, E2, F1, F2;
814 Float_t CoordPadRing;
816 //parameters to calculate freon window refractive index vs. energy
820 //parameters to calculate quartz window refractive index vs. energy
837 for (phpad=0; phpad<nmaxdegrees;phpad++) {
839 AveEnerg = (Energ[0]+Energ[1])/2.;
841 nfreonave = a+b*AveEnerg;
842 nquartzave = sqrt(1+(F1/(pow(E1,2)-pow(AveEnerg,2)))+
843 (F2/(pow(E2,2)-pow(AveEnerg,2))));
845 thetacer = fixedthetacer;
847 bandradius = DistanceFromMip(nfreonave, nquartzave,
848 fEmissPoint,thetacer, phpad);
850 CoordPadRing=fPhotocatExitPhot;
852 phpad = (nmaxdegrees/nstepdegrees)*phpad;