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 **************************************************************************/
18 Revision 1.7 2000/10/03 21:44:09 morsch
19 Use AliSegmentation and AliHit abstract base classes.
21 Revision 1.6 2000/10/02 21:28:12 fca
22 Removal of useless dependecies via forward declarations
24 Revision 1.5 2000/10/02 15:50:25 jbarbosa
25 Fixed forward declarations.
27 Revision 1.4 2000/06/30 16:33:43 dibari
28 Several changes (ring drawing, fiducial selection, etc.)
30 Revision 1.3 2000/06/15 15:47:12 jbarbosa
31 Corrected compilation errors on HP-UX (replaced pow with TMath::Power)
33 Revision 1.2 2000/06/12 15:26:09 jbarbosa
36 Revision 1.1 2000/06/09 14:53:01 jbarbosa
37 Bari's pattern recognition algorithm
41 #include "AliRICHHit.h"
42 #include "AliRICHCerenkov.h"
43 #include "AliRICHPadHit.h"
44 #include "AliRICHDigit.h"
45 #include "AliRICHRawCluster.h"
46 #include "AliRICHRecHit1D.h"
48 #include "AliDetector.h"
50 #include "AliRICHPoints.h"
51 #include "AliSegmentation.h"
52 #include "AliRICHPatRec.h"
54 #include "AliRICHConst.h"
55 #include "AliRICHPoints.h"
57 #include "AliHitMap.h"
59 #include <TParticle.h>
67 ClassImp(AliRICHPatRec)
68 //___________________________________________
69 AliRICHPatRec::AliRICHPatRec() : TObject()
71 // Default constructor
75 //___________________________________________
76 AliRICHPatRec::AliRICHPatRec(const char *name, const char *title)
79 //Constructor for Bari's pattern recogniton method object
82 void AliRICHPatRec::PatRec()
85 // Pattern recognition algorithm
87 AliRICHChamber* iChamber;
88 AliSegmentation* segmentation;
90 Int_t ntracks, ndigits[kNCH];
101 printf("PatRec started\n");
103 AliRICH *pRICH = (AliRICH*)gAlice->GetDetector("RICH");
104 TTree *treeH = gAlice->TreeH();
106 ntracks =(Int_t) treeH->GetEntries();
108 for (itr=0; itr<ntracks; itr++) {
110 status = TrackParam(itr,ich);
111 if(status==1) continue;
112 //printf(" theta %f phi %f track \n",fTrackTheta,fTrackPhi);
113 // ring->Fill(fTrackLoc[0],fTrackLoc[1],100.);
115 iChamber = &(pRICH->Chamber(ich));
116 segmentation=iChamber->GetSegmentationModel();
118 nent=(Int_t)gAlice->TreeD()->GetEntries();
119 gAlice->TreeD()->GetEvent(nent-1);
120 TClonesArray *pDigitss = pRICH->DigitsAddress(ich);
121 ndigits[ich] = pDigitss->GetEntriesFast();
122 printf("Digits in chamber %d: %d\n",ich,ndigits[ich]);
123 AliRICHDigit *padI = 0;
127 for (Int_t dig=0;dig<ndigits[ich];dig++) {
128 padI=(AliRICHDigit*) pDigitss->UncheckedAt(dig);
132 segmentation->GetPadC(x,y,rx,ry,rz);
134 //printf("Pad coordinates x:%d, Real coordinates x:%f\n",x,rx);
135 //printf("Pad coordinates y:%d, Real coordinates y:%f\n",y,ry);
141 fCerenkovAnglePad = PhotonCerenkovAngle();
142 if(fCerenkovAnglePad==-999) continue;
144 if(!PhotonInBand()) continue;
149 segmentation->GetPadI(fXpad,fYpad,0,xpad,ypad);
151 padsUsedX[goodPhotons]=xpad;
152 padsUsedY[goodPhotons]=ypad;
155 fEtaPhotons[goodPhotons-1] = fCerenkovAnglePad;
157 fNumEtaPhotons = goodPhotons;
159 BackgroundEstimation();
162 //CerenkovRingDrawing();
166 rechit[2] = fThetaCerenkov;
167 rechit[3] = fXshift + fTrackLoc[0];
168 rechit[4] = fYshift + fTrackLoc[1];
169 rechit[5] = fEmissPoint;
170 rechit[6] = goodPhotons;
172 //printf("Center coordinates:%f %f\n",rechit[3],rechit[4]);
174 pRICH->AddRecHit1D(ich,rechit,fEtaPhotons,padsUsedX,padsUsedY);
178 gAlice->TreeR()->Fill();
180 for (i=0;i<kNCH;i++) {
181 fRec=pRICH->RecHitsAddress1D(i);
182 int ndig=fRec->GetEntriesFast();
183 printf ("Chamber %d, rings %d\n",i,ndig);
185 pRICH->ResetRecHits1D();
190 Int_t AliRICHPatRec::TrackParam(Int_t itr, Int_t &ich)
192 // Get Local coordinates of track impact
194 AliRICHChamber* iChamber;
195 AliSegmentation* segmentation;
197 Float_t trackglob[3];
205 //printf("Calling TrackParam\n");
208 TTree *treeH = gAlice->TreeH();
209 treeH->GetEvent(itr);
211 AliRICH *pRICH = (AliRICH*)gAlice->GetDetector("RICH");
212 AliRICHHit* mHit=(AliRICHHit*)pRICH->FirstHit(-1);
213 if(mHit==0) return 1;
214 ich = mHit->fChamber-1;
215 trackglob[0] = mHit->X();
216 trackglob[1] = mHit->Y();
217 trackglob[2] = mHit->Z();
221 fTrackMom = sqrt(TMath::Power(pX,2)+TMath::Power(pY,2)+TMath::Power(pZ,2));
222 thetatr = (mHit->fTheta)*(Float_t)kDegrad;
223 phitr = mHit->fPhi*(Float_t)kDegrad;
225 part = mHit->fParticle;
227 iChamber = &(pRICH->Chamber(ich));
228 iChamber->GlobaltoLocal(trackglob,trackloc);
230 segmentation=iChamber->GetSegmentationModel();
232 // retrieve geometrical params
234 AliRICHGeometry* fGeometry=iChamber->GetGeometryModel();
236 fRw = fGeometry->GetFreonThickness();
237 fQw = fGeometry->GetQuartzThickness();
238 fTgap = fGeometry->GetGapThickness();
239 Float_t radiatorToPads= fGeometry->GetRadiatorToPads();
240 //+ fGeometry->GetProximityGapThickness();
242 //printf("Distance to pads. From geometry:%f, From calculations:%f\n",radiatorToPads,fRw + fQw + fTgap);
244 //Float_t apar = (fRw + fQw + fTgap)*tan(thetatr);
245 Float_t apar = radiatorToPads*tan(thetatr);
246 fTrackLoc[0] = apar*cos(phitr);
247 fTrackLoc[1] = apar*sin(phitr);
248 //fTrackLoc[2] = fRw + fQw + fTgap;
249 fTrackLoc[2] = radiatorToPads;
250 fTrackTheta = thetatr;
253 fXshift = trackloc[0] - fTrackLoc[0];
254 fYshift = trackloc[2] - fTrackLoc[1];
259 Float_t AliRICHPatRec::EstimationAtLimits(Float_t lim, Float_t radius,
263 // Estimation of emission point
265 Float_t nquartz = 1.585;
267 Float_t nfreon = 1.295;
270 // printf("Calling EstimationLimits\n");
272 Float_t apar = (fRw -fEmissPoint + fQw + fTgap)*tan(fTrackTheta);
273 Float_t b1 = (fRw-fEmissPoint)*tan(lim);
274 Float_t b2 = fQw / sqrt(TMath::Power(nquartz,2)-TMath::Power(nfreon*sin(lim),2));
275 Float_t b3 = fTgap / sqrt(TMath::Power(ngas,2)-TMath::Power(nfreon*sin(lim),2));
276 Float_t bpar = b1 + nfreon*sin(lim)*(b2+b3);
277 value = TMath::Power(radius,2)
278 -TMath::Power((apar*cos(fTrackPhi)-bpar*cos(phiphot)),2)
279 -TMath::Power((apar*sin(fTrackPhi)-bpar*sin(phiphot)),2);
284 Float_t AliRICHPatRec::PhotonCerenkovAngle()
286 // Cherenkov pad angle reconstruction
290 Float_t cherMax = 0.8;
292 Float_t eps = 0.0001;
293 Int_t niterEmiss = 0;
294 Int_t niterEmissMax = 0;
295 Float_t x1,x2,x3=0,p1,p2,p3;
299 // printf("Calling PhotonCerenkovAngle\n");
301 radius = sqrt(TMath::Power(fTrackLoc[0]-fXpad,2)+TMath::Power(fTrackLoc[1]-fYpad,2));
302 fEmissPoint = fRw/2.; //Start value of EmissionPoint
304 while(niterEmiss<=niterEmissMax) {
307 argY = fYpad - fEmissPoint*tan(fTrackTheta)*sin(fTrackPhi);
308 argX = fXpad - fEmissPoint*tan(fTrackTheta)*cos(fTrackPhi);
309 phiphot = atan2(argY,argX);
310 p1 = EstimationAtLimits(cherMin,radius,phiphot);
311 p2 = EstimationAtLimits(cherMax,radius,phiphot);
314 // printf("PhotonCerenkovAngle failed\n");
318 //start to find the Cherenkov pad angle
322 p3 = EstimationAtLimits(x3,radius,phiphot);
323 while(TMath::Abs(p3)>eps){
327 p1 = EstimationAtLimits(x1,radius,phiphot);
330 p3 = EstimationAtLimits(x3,radius,phiphot);
334 // printf(" max iterations in PhotonCerenkovAngle\n");
338 // printf("niterFun %i \n",niterFun);
340 if (niterEmiss != niterEmissMax+1) EmissionPoint();
343 printf(" phiphot %f fXpad %f fYpad %f fEmiss %f \n",
344 phiphot,fXpad,fYpad,fEmissPoint);
352 void AliRICHPatRec::EmissionPoint()
355 // Find emission point
357 Float_t absorbtionLength=7.83*fRw; //absorption length in the freon (cm)
358 // 7.83 = -1/ln(T0) where
359 // T0->Trasmission freon at 180nm = 0.88 (Eph=6.85eV)
360 Float_t photonLength, photonLengthMin, photonLengthMax;
362 photonLength=exp(-fRw/(absorbtionLength*cos(fCerenkovAnglePad)));
363 photonLengthMin=fRw*photonLength/(1.-photonLength);
364 photonLengthMax=absorbtionLength*cos(fCerenkovAnglePad);
365 fEmissPoint = fRw + photonLengthMin - photonLengthMax;
369 void AliRICHPatRec::PhotonSelection(Int_t track, Int_t &nphot, Float_t &thetamean)
372 // not implemented yet
374 printf("Calling PhotonSelection\n");
377 void AliRICHPatRec::BackgroundEstimation()
380 // estimate background noise
382 Float_t stepEta = 0.001;
383 Float_t etaMinBkg = 0.72;
384 Float_t etaMaxBkg = 0.75;
386 Float_t etaMax = 0.75;
388 Float_t nfreon = 1.295;
390 Float_t etaStepMin,etaStepMax,etaStepAvg;
392 Int_t numPhotBkg, numPhotonStep;
393 Float_t funBkg,areaBkg,normBkg;
394 Float_t densityBkg,storeBkg,numStore;
400 nstep = (int)((etaMaxBkg-etaMinBkg)/stepEta);
402 for (i=0;i<fNumEtaPhotons;i++) {
404 if(fEtaPhotons[i]>etaMinBkg && fEtaPhotons[i]<etaMaxBkg) {
408 if (numPhotBkg == 0) {
409 for (i=0;i<fNumEtaPhotons;i++) {
410 fWeightPhotons[i] = 1.;
415 // printf(" numPhotBkg %i ",numPhotBkg);
417 for (i=0;i<nstep;i++) {
418 etaStepMin = etaMinBkg + (Float_t)(i)*stepEta;
419 etaStepMax = etaMinBkg + (Float_t)(i+1)*stepEta;
420 etaStepAvg = 0.5*(etaStepMax + etaStepMin);
422 funBkg = tan(etaStepAvg)*TMath::Power((1.+TMath::Power(tan(etaStepAvg),2)),
423 5.52)-7.803 + 22.02*tan(etaStepAvg);
426 //printf("etaStepAvg: %f, etaStepMax: %f, etaStepMin: %f", etaStepAvg,etaStepMax,etaStepMin);
428 thetaSig = TMath::ASin(nfreon/ngas*TMath::Sin(etaStepAvg));
429 funBkg = tan(thetaSig)*(1.+TMath::Power(tan(thetaSig),2))*nfreon
430 /ngas*cos(etaStepAvg)/cos(thetaSig);
431 areaBkg += stepEta*funBkg;
434 densityBkg = 0.95*(Float_t)(numPhotBkg)/areaBkg;
435 // printf(" densityBkg %f \n",densityBkg);
437 nstep = (int)((etaMax-etaMin)/stepEta);
440 for (i=0;i<nstep;i++) {
441 etaStepMin = etaMinBkg + (Float_t)(i)*stepEta;
442 etaStepMax = etaMinBkg + (Float_t)(i+1)*stepEta;
443 etaStepAvg = 0.5*(etaStepMax + etaStepMin);
445 funBkg = tan(etaStepAvg)*TMath::Power((1.+TMath::Power(tan(etaStepAvg),2)),
446 5.52)-7.803 + 22.02*tan(etaStepAvg);
449 thetaSig = asin(nfreon/ngas*sin(etaStepAvg));
450 funBkg = tan(thetaSig)*(1.+TMath::Power(tan(thetaSig),2))*nfreon
451 /ngas*cos(etaStepAvg)/cos(thetaSig);
453 areaBkg = stepEta*funBkg;
454 normBkg = densityBkg*areaBkg;
456 for (ip=0;ip<fNumEtaPhotons;ip++) {
457 if(fEtaPhotons[ip]>etaStepMin && fEtaPhotons[ip]<etaStepMax) {
461 if (numPhotonStep == 0) {
469 if (numPhotonStep == 0) continue;
470 for (ip=0;ip<fNumEtaPhotons;ip++) {
471 if(fEtaPhotons[ip]>etaStepMin && fEtaPhotons[ip]<etaStepMax) {
475 fWeightPhotons[ip] = 1. - normBkg/(Float_t)(numPhotonStep);
477 printf(" normBkg %f numPhotonStep %i fW %f \n",
478 normBkg, numPhotonStep, fWeightPhotons[ip]);
480 if(fWeightPhotons[ip]<0) fWeightPhotons[ip] = 0.;
487 void AliRICHPatRec::FlagPhotons(Int_t track, Float_t theta)
490 // not implemented yet
492 printf("Calling FlagPhotons\n");
496 //////////////////////////////////////////
502 Int_t AliRICHPatRec::PhotonInBand()
504 //0=label for parameters giving internal band ellipse
505 //1=label for parameters giving external band ellipse
507 Float_t imp[2], mass[2], energy[2], beta[2];
508 Float_t emissPointLength[2];
509 Float_t e1, e2, f1, f2;
510 Float_t nfreon[2], nquartz[2];
512 Float_t pointsOnCathode[3];
514 Float_t phpad, thetacer[2];
515 Float_t bandradius[2], padradius;
517 imp[0] = 5.0; //threshold momentum for the proton Cherenkov emission
520 mass[0] = 0.938; //proton mass
521 mass[1] = 0.139; //pion mass
523 emissPointLength[0] = fRw-0.0001; //at the beginning of the radiator
524 emissPointLength[1] = 0.;//at the end of radiator
526 //parameters to calculate freon window refractive index vs. energy
530 //parameters to calculate quartz window refractive index vs. energy
544 for (times=0; times<=1; times++) {
546 nfreon[times] = a+b*energy[times];
549 nquartz[times] = sqrt(1+(f1/(TMath::Power(e1,2)-TMath::Power(energy[times],2)))+
550 (f2/(TMath::Power(e2,2)-TMath::Power(energy[times],2))));
552 beta[times] = imp[times]/sqrt(TMath::Power(imp[times],2)+TMath::Power(mass[times],2));
554 thetacer[times] = CherenkovAngle( nfreon[times], beta[times]);
556 bandradius[times] = DistanceFromMip( nfreon[times], nquartz[times],
557 emissPointLength[times],
558 thetacer[times], phpad, pointsOnCathode);
559 //printf(" ppp %f %f %f \n",pointsOnCathode);
562 bandradius[0] -= 1.6;
563 bandradius[1] += 1.6;
564 padradius = sqrt(TMath::Power(fXpad,2)+TMath::Power(fYpad,2));
565 // printf(" rmin %f r %f rmax %f \n",bandradius[0],padradius,bandradius[1]);
567 if(padradius>=bandradius[0] && padradius<=bandradius[1]) return 1;
571 Float_t AliRICHPatRec::DistanceFromMip(Float_t nfreon, Float_t nquartz,
572 Float_t emissPointLength, Float_t thetacer,
573 Float_t phpad, Float_t pointsOnCathode[3])
576 // Find the distance to MIP impact
578 Float_t distanceValue;
580 TVector3 radExitPhot(1,1,1);//photon impact at the radiator exit with respect
581 //to local reference sistem with the origin in the MIP entrance
583 TVector3 vectEmissPointLength(1,1,1);
584 Float_t magEmissPointLenght;
586 TVector3 radExitPhot2(1,1,1);//photon impact at the radiator exit with respect
587 Float_t magRadExitPhot2;
588 //to a reference sistem with origin in the photon emission point and
589 //axes parallel to the MIP reference sistem
591 TVector3 quarExitPhot(1,1,1);//photon impact at the quartz exit with respect
592 Float_t magQuarExitPhot;
594 TVector3 gapExitPhot(1,1,1) ;
595 Float_t magGapExitPhot;
597 TVector3 PhotocatExitPhot(1,1,1);
599 Double_t thetarad , phirad ;
600 Double_t thetaquar, phiquar;
601 Double_t thetagap , phigap ;
605 magEmissPointLenght = emissPointLength/cos(fTrackTheta);
607 vectEmissPointLength.SetMag(magEmissPointLenght);
608 vectEmissPointLength.SetTheta(fTrackTheta);
609 vectEmissPointLength.SetPhi(fTrackPhi);
612 radExitPhot2.SetTheta(thetacer);
613 radExitPhot2.SetPhi(phpad);
620 r1. RotateY(fTrackTheta);
621 r2. RotateZ(fTrackPhi);
625 r = r2 * r1;//rotation about the y axis by MIP theta incidence angle
626 //following by a rotation about the z axis by MIP phi incidence angle;
629 radExitPhot2 = r * radExitPhot2;
630 theta2 = radExitPhot2.Theta();
631 magRadExitPhot2 = (fRw - vectEmissPointLength(2))/cos(theta2);
632 radExitPhot2.SetMag(magRadExitPhot2);
635 radExitPhot = vectEmissPointLength + radExitPhot2;
636 thetarad = radExitPhot.Theta();
637 phirad = radExitPhot.Phi(); //check on the original file //
639 thetaquar = SnellAngle( nfreon, nquartz, theta2);
640 phiquar = radExitPhot2.Phi();
641 if(thetaquar == 999.) return thetaquar;
642 magQuarExitPhot = fQw/cos(thetaquar);
643 quarExitPhot.SetMag( magQuarExitPhot);
644 quarExitPhot.SetTheta(thetaquar);
645 quarExitPhot.SetPhi(phiquar);
647 thetagap = SnellAngle( nquartz, ngas, thetaquar);
649 if(thetagap == 999.) return thetagap;
650 magGapExitPhot = fTgap/cos(thetagap);
651 gapExitPhot.SetMag( magGapExitPhot);
652 gapExitPhot.SetTheta(thetagap);
653 gapExitPhot.SetPhi(phigap);
655 PhotocatExitPhot = radExitPhot + quarExitPhot + gapExitPhot;
657 distanceValue = sqrt(TMath::Power(PhotocatExitPhot(0),2)
658 +TMath::Power(PhotocatExitPhot(1),2));
659 pointsOnCathode[0] = (Float_t) PhotocatExitPhot(0) + fXshift - fTrackLoc[0];
660 pointsOnCathode[1] = (Float_t) PhotocatExitPhot(1) + fYshift - fTrackLoc[1];
661 pointsOnCathode[2] = (Float_t) PhotocatExitPhot(2);
663 //printf(" point in Distance.2. %f %f %f \n",pointsOnCathode[0],pointsOnCathode[1],pointsOnCathode[2]);
665 return distanceValue;
669 Float_t AliRICHPatRec::PhiPad()
675 Float_t thetapad, phipad;
676 Float_t thetarot, phirot;
678 zpad = fRw + fQw + fTgap;
680 TVector3 photonPad(fXpad, fYpad, zpad);
681 thetapad = photonPad.Theta();
682 phipad = photonPad.Phi();
688 thetarot = - fTrackTheta;
689 phirot = - fTrackPhi;
691 r2. RotateY(thetarot);
693 r = r2 * r1;//rotation about the z axis by MIP -phi incidence angle
694 //following by a rotation about the y axis by MIP -theta incidence angle;
696 photonPad = r * photonPad;
698 phipad = photonPad.Phi();
703 Float_t AliRICHPatRec:: SnellAngle(Float_t n1, Float_t n2, Float_t theta1)
706 // Compute the Snell angle
708 Float_t sinrefractangle;
709 Float_t refractangle;
711 sinrefractangle = (n1/n2)*sin(theta1);
713 if(sinrefractangle>1.) {
718 refractangle = asin(sinrefractangle);
722 Float_t AliRICHPatRec::CherenkovAngle(Float_t n, Float_t beta)
725 // Compute the cerenkov angle
734 thetacer = acos (1./(n*beta));
738 Float_t AliRICHPatRec::BetaCerenkov(Float_t n, Float_t theta)
745 beta = 1./(n*cos(theta));
752 void AliRICHPatRec::HoughResponse()
756 // Implement Hough response pat. rec. method
761 int i, j, k, nCorrBand;
764 float angle, thetaCerMean;
769 float stepEta = 0.001;
770 float windowEta = 0.040;
774 float etaPeakPos = -1;
775 Int_t etaPeakCount = -1;
780 nBin = (int)(0.5+etaMax/(stepEta));
781 nCorrBand = (int)(0.5+ windowEta/(2 * stepEta));
782 memset ((void *)hcs, 0, etaBin*sizeof(int));
784 for (k=0; k< fNumEtaPhotons; k++) {
786 angle = fEtaPhotons[k];
788 if (angle>=etaMin && angle<= etaMax) {
789 bin = (int)(0.5+angle/(stepEta));
793 if (bin2>nBin) bin2=nBin;
795 for (j=bin1; j<bin2; j++) {
796 hcs[j] += fWeightPhotons[k];
799 thetaCerMean += angle;
803 thetaCerMean /= fNumEtaPhotons;
807 for (bin=0; bin <nBin; bin++) {
808 angle = (bin+0.5) * (stepEta);
809 if (hcs[bin] && hcs[bin] > etaPeakPos) {
811 etaPeakPos = hcs[bin];
815 if (hcs[bin] == etaPeakPos) {
816 etaPeak[++etaPeakCount] = angle;
821 for (i=0; i<etaPeakCount+1; i++) {
822 fThetaCerenkov += etaPeak[i];
824 if (etaPeakCount>=0) {
825 fThetaCerenkov /= etaPeakCount+1;
826 fThetaPeakPos = etaPeakPos;
831 void AliRICHPatRec::HoughFiltering(float hcs[])
837 float k[5] = {0.05, 0.25, 0.4, 0.25, 0.05};
844 float stepEta = 0.001;
846 nBin = (int)(1+etaMax/stepEta);
847 sizeHCS = etaBin*sizeof(float);
849 memset ((void *)hcsFilt, 0, sizeHCS);
851 for (nx = 0; nx < nBin; nx++) {
852 for (i = 0; i < 5; i++) {
854 if (nxDx> -1 && nxDx<nBin)
855 hcsFilt[nx] += hcs[nxDx] * k[i];
859 for (nx = 0; nx < nBin; nx++) {
860 hcs[nx] = hcsFilt[nx];
864 /*void AliRICHPatRec::CerenkovRingDrawing()
867 //to draw Cherenkov ring by known Cherenkov angle
872 Float_t nfreonave, nquartzave;
875 Float_t e1, e2, f1, f2;
878 //parameters to calculate freon window refractive index vs. energy
883 //parameters to calculate quartz window refractive index vs. energy
895 for (Nphpad=0; Nphpad<nmaxdegrees;Nphpad++) {
897 phpad = (360./(Float_t)nmaxdegrees)*(Float_t)Nphpad;
899 aveEnerg = (energy[0]+energy[1])/2.;
901 nfreonave = a+b*aveEnerg;
902 nquartzave = sqrt(1+(f1/(TMath::Power(e1,2)-TMath::Power(aveEnerg,2)))+
903 (f2/(TMath::Power(e2,2)-TMath::Power(aveEnerg,2))));
905 bandradius = DistanceFromMip(nfreonave, nquartzave,
906 fEmissPoint,fThetaCerenkov, phpad);
908 fCoordEllipse[0][Nphpad] = fOnCathode[0];
909 fCoordEllipse[1][Nphpad] = fOnCathode[1];
910 printf(" values %f %f \n",fOnCathode[0],fOnCathode[1]);