/*
$Log$
+ Revision 1.10 2001/02/27 15:21:06 jbarbosa
+ Transition to SDigits.
+
+ Revision 1.9 2001/02/13 20:38:48 jbarbosa
+ Changes to make it work with new IO.
+
+ Revision 1.8 2000/11/01 15:37:18 jbarbosa
+ Updated to use its own rec. point object.
+
+ Revision 1.7 2000/10/03 21:44:09 morsch
+ Use AliSegmentation and AliHit abstract base classes.
+
+ Revision 1.6 2000/10/02 21:28:12 fca
+ Removal of useless dependecies via forward declarations
+
+ Revision 1.5 2000/10/02 15:50:25 jbarbosa
+ Fixed forward declarations.
+
+ Revision 1.4 2000/06/30 16:33:43 dibari
+ Several changes (ring drawing, fiducial selection, etc.)
+
+ Revision 1.3 2000/06/15 15:47:12 jbarbosa
+ Corrected compilation errors on HP-UX (replaced pow with TMath::Power)
+
+ Revision 1.2 2000/06/12 15:26:09 jbarbosa
+ Cleaned up version.
+
Revision 1.1 2000/06/09 14:53:01 jbarbosa
Bari's pattern recognition algorithm
#include "AliRICHHit.h"
#include "AliRICHCerenkov.h"
-#include "AliRICHPadHit.h"
+#include "AliRICHSDigit.h"
#include "AliRICHDigit.h"
#include "AliRICHRawCluster.h"
-#include "AliRICHRecHit.h"
+#include "AliRICHRecHit1D.h"
#include "AliRun.h"
#include "AliDetector.h"
#include "AliRICH.h"
#include "AliRICHPoints.h"
-#include "AliRICHSegmentation.h"
+#include "AliSegmentation.h"
#include "AliRICHPatRec.h"
#include "AliRICH.h"
#include "AliRICHConst.h"
#include "AliRICHPoints.h"
#include "AliConst.h"
+#include "AliHitMap.h"
#include <TParticle.h>
#include <TMath.h>
#include <TRandom.h>
#include <TCanvas.h>
#include <TH2.h>
+#include <TTree.h>
ClassImp(AliRICHPatRec)
// Pattern recognition algorithm
AliRICHChamber* iChamber;
- AliRICHSegmentation* segmentation;
+ AliSegmentation* segmentation;
Int_t ntracks, ndigits[kNCH];
Int_t itr, ich, i;
Int_t goodPhotons;
Int_t x,y,q;
- Float_t rx,ry;
+ Float_t rx,ry,rz;
Int_t nent,status;
+ Int_t padsUsedX[100];
+ Int_t padsUsedY[100];
- Float_t gamma,massCer,betaCer;
-
- Float_t rechit[5];
-
- printf("PatRec started\n");
+ Float_t rechit[7];
- TCanvas *c1 = new TCanvas("c1","Alice RICH pad hits",50,10,700,700);
+ //printf("PatRec started\n");
- TH2F *ring = new TH2F("ring","map",90,-30.,30.,90,-30.,30.);
- TH2F *ringband = new TH2F("ringband","map",90,-65.,65.,90,-65.,65.);
- TH1F *cerangle = new TH1F("cerangle","phot",300,0.45,0.75);
- TH1F *ceranglew= new TH1F("ceranglew","phot",300,0.45,0.75);
- TH1F *hough = new TH1F("hough","hough",75,0.45,0.75);
- TH1F *mass = new TH1F("mass","mass",100,50.,1050.);
-
AliRICH *pRICH = (AliRICH*)gAlice->GetDetector("RICH");
TTree *treeH = gAlice->TreeH();
status = TrackParam(itr,ich);
if(status==1) continue;
- printf(" theta %f phi %f track \n",fTrackTheta,fTrackPhi);
+ //printf(" theta %f phi %f track \n",fTrackTheta,fTrackPhi);
// ring->Fill(fTrackLoc[0],fTrackLoc[1],100.);
iChamber = &(pRICH->Chamber(ich));
segmentation=iChamber->GetSegmentationModel();
nent=(Int_t)gAlice->TreeD()->GetEntries();
- gAlice->TreeD()->GetEvent(nent-1);
+ gAlice->TreeD()->GetEvent(0);
TClonesArray *pDigitss = pRICH->DigitsAddress(ich);
ndigits[ich] = pDigitss->GetEntriesFast();
- printf("ndigits %d in chamber %d\n",ndigits[ich],ich);
+ printf("Digits in chamber %d: %d\n",ich,ndigits[ich]);
AliRICHDigit *padI = 0;
goodPhotons = 0;
x=padI->fPadX;
y=padI->fPadY;
q=padI->fSignal;
- segmentation->GetPadCxy(x,y,rx,ry);
+ segmentation->GetPadC(x,y,rx,ry,rz);
+
+ //printf("Pad coordinates x:%d, Real coordinates x:%f\n",x,rx);
+ //printf("Pad coordinates y:%d, Real coordinates y:%f\n",y,ry);
fXpad = rx-fXshift;
fYpad = ry-fYshift;
fQpad = q;
- ringband->Fill(x,y,1.);
fCerenkovAnglePad = PhotonCerenkovAngle();
if(fCerenkovAnglePad==-999) continue;
if(!PhotonInBand()) continue;
- ring->Fill(fXpad,fYpad,1.);
- cerangle->Fill(fCerenkovAnglePad,1.);
+ Int_t xpad;
+ Int_t ypad;
+
+ segmentation->GetPadI(fXpad,fYpad,0,xpad,ypad);
+ padsUsedX[goodPhotons]=xpad;
+ padsUsedY[goodPhotons]=ypad;
+
goodPhotons++;
- fEtaPhotons[goodPhotons] = fCerenkovAnglePad;
+ fEtaPhotons[goodPhotons-1] = fCerenkovAnglePad;
}
fNumEtaPhotons = goodPhotons;
BackgroundEstimation();
- for(i=0;i<goodPhotons;i++) {
- ceranglew->Fill(fEtaPhotons[i],fWeightPhotons[i]);
- // printf(" Eta %f weight %f \n",fEtaPhotons[i],fWeightPhotons[i]);
- }
-
HoughResponse();
+ //CerenkovRingDrawing();
rechit[0] = 0;
rechit[1] = 0;
rechit[2] = fThetaCerenkov;
- rechit[3] = 0;
- rechit[4] = 0;
+ rechit[3] = fXshift + fTrackLoc[0];
+ rechit[4] = fYshift + fTrackLoc[1];
+ rechit[5] = fEmissPoint;
+ rechit[6] = goodPhotons;
-
- hough->Fill(fThetaCerenkov,1.);
+ //printf("Center coordinates:%f %f\n",rechit[3],rechit[4]);
- pRICH->AddRecHit(ich,rechit);
+ pRICH->AddRecHit1D(ich,rechit,fEtaPhotons,padsUsedX,padsUsedY);
- betaCer = BetaCerenkov(1.29,fThetaCerenkov);
- gamma = 1./sqrt(1.-pow(betaCer,2));
- massCer = fTrackMom/(betaCer*gamma);
- // printf(" mass %f \n",massCer);
- mass->Fill(massCer*1000,1.);
}
gAlice->TreeR()->Fill();
TClonesArray *fRec;
for (i=0;i<kNCH;i++) {
- fRec=pRICH->RecHitsAddress(i);
+ fRec=pRICH->RecHitsAddress1D(i);
int ndig=fRec->GetEntriesFast();
printf ("Chamber %d, rings %d\n",i,ndig);
}
- pRICH->ResetRecHits();
+ pRICH->ResetRecHits1D();
-
- c1->Divide(2,2);
- c1->cd(1);
- ring->Draw("box");
- c1->cd(2);
- ringband->Draw("box");
- c1->cd(3);
- cerangle->Draw();
- c1->cd(4);
- hough->Draw();
-
}
// Get Local coordinates of track impact
AliRICHChamber* iChamber;
- AliRICHSegmentation* segmentation;
+ AliSegmentation* segmentation;
Float_t trackglob[3];
Float_t trackloc[3];
Float_t part;
Float_t pX, pY, pZ;
- printf("Calling TrackParam\n");
+ //printf("Calling TrackParam\n");
gAlice->ResetHits();
TTree *treeH = gAlice->TreeH();
AliRICHHit* mHit=(AliRICHHit*)pRICH->FirstHit(-1);
if(mHit==0) return 1;
ich = mHit->fChamber-1;
- trackglob[0] = mHit->fX;
- trackglob[1] = mHit->fY;
- trackglob[2] = mHit->fZ;
+ trackglob[0] = mHit->X();
+ trackglob[1] = mHit->Y();
+ trackglob[2] = mHit->Z();
pX = mHit->fMomX;
pY = mHit->fMomY;
pZ = mHit->fMomZ;
- fTrackMom = sqrt(pow(pX,2)+pow(pY,2)+pow(pZ,2));
- thetatr = (180 - mHit->fTheta)*(Float_t)kDegrad;
+ fTrackMom = sqrt(TMath::Power(pX,2)+TMath::Power(pY,2)+TMath::Power(pZ,2));
+ thetatr = (mHit->fTheta)*(Float_t)kDegrad;
phitr = mHit->fPhi*(Float_t)kDegrad;
iloss = mHit->fLoss;
part = mHit->fParticle;
fRw = fGeometry->GetFreonThickness();
fQw = fGeometry->GetQuartzThickness();
- fTgap = fGeometry->GetGapThickness()
- + fGeometry->GetProximityGapThickness();
+ fTgap = fGeometry->GetGapThickness();
+ Float_t radiatorToPads= fGeometry->GetRadiatorToPads();
+ //+ fGeometry->GetProximityGapThickness();
- Float_t apar = (fRw + fQw + fTgap)*tan(thetatr);
+ //printf("Distance to pads. From geometry:%f, From calculations:%f\n",radiatorToPads,fRw + fQw + fTgap);
+
+ //Float_t apar = (fRw + fQw + fTgap)*tan(thetatr);
+ Float_t apar = radiatorToPads*tan(thetatr);
fTrackLoc[0] = apar*cos(phitr);
fTrackLoc[1] = apar*sin(phitr);
- fTrackLoc[2] = fRw + fQw + fTgap;
+ //fTrackLoc[2] = fRw + fQw + fTgap;
+ fTrackLoc[2] = radiatorToPads;
fTrackTheta = thetatr;
fTrackPhi = phitr;
Float_t apar = (fRw -fEmissPoint + fQw + fTgap)*tan(fTrackTheta);
Float_t b1 = (fRw-fEmissPoint)*tan(lim);
- Float_t b2 = fQw / sqrt(pow(nquartz,2)-pow(nfreon*sin(lim),2));
- Float_t b3 = fTgap / sqrt(pow(ngas,2)-pow(nfreon*sin(lim),2));
+ Float_t b2 = fQw / sqrt(TMath::Power(nquartz,2)-TMath::Power(nfreon*sin(lim),2));
+ Float_t b3 = fTgap / sqrt(TMath::Power(ngas,2)-TMath::Power(nfreon*sin(lim),2));
Float_t bpar = b1 + nfreon*sin(lim)*(b2+b3);
- value = pow(radius,2)
- -pow((apar*cos(fTrackPhi)-bpar*cos(phiphot)),2)
- -pow((apar*sin(fTrackPhi)-bpar*sin(phiphot)),2);
+ value = TMath::Power(radius,2)
+ -TMath::Power((apar*cos(fTrackPhi)-bpar*cos(phiphot)),2)
+ -TMath::Power((apar*sin(fTrackPhi)-bpar*sin(phiphot)),2);
return value;
}
// printf("Calling PhotonCerenkovAngle\n");
- radius = sqrt(pow(fTrackLoc[0]-fXpad,2)+pow(fTrackLoc[1]-fYpad,2));
+ radius = sqrt(TMath::Power(fTrackLoc[0]-fXpad,2)+TMath::Power(fTrackLoc[1]-fYpad,2));
fEmissPoint = fRw/2.; //Start value of EmissionPoint
while(niterEmiss<=niterEmissMax) {
etaStepMax = etaMinBkg + (Float_t)(i+1)*stepEta;
etaStepAvg = 0.5*(etaStepMax + etaStepMin);
/*
- funBkg = tan(etaStepAvg)*pow((1.+pow(tan(etaStepAvg),2)),
+ funBkg = tan(etaStepAvg)*TMath::Power((1.+TMath::Power(tan(etaStepAvg),2)),
5.52)-7.803 + 22.02*tan(etaStepAvg);
*/
- thetaSig = asin(nfreon/ngas*sin(etaStepAvg));
- funBkg = tan(thetaSig)*(1.+pow(tan(thetaSig),2))*nfreon
+
+ //printf("etaStepAvg: %f, etaStepMax: %f, etaStepMin: %f", etaStepAvg,etaStepMax,etaStepMin);
+
+ thetaSig = TMath::ASin(nfreon/ngas*TMath::Sin(etaStepAvg));
+ funBkg = tan(thetaSig)*(1.+TMath::Power(tan(thetaSig),2))*nfreon
/ngas*cos(etaStepAvg)/cos(thetaSig);
areaBkg += stepEta*funBkg;
}
etaStepMax = etaMinBkg + (Float_t)(i+1)*stepEta;
etaStepAvg = 0.5*(etaStepMax + etaStepMin);
/*
- funBkg = tan(etaStepAvg)*pow((1.+pow(tan(etaStepAvg),2)),
+ funBkg = tan(etaStepAvg)*TMath::Power((1.+TMath::Power(tan(etaStepAvg),2)),
5.52)-7.803 + 22.02*tan(etaStepAvg);
*/
thetaSig = asin(nfreon/ngas*sin(etaStepAvg));
- funBkg = tan(thetaSig)*(1.+pow(tan(thetaSig),2))*nfreon
+ funBkg = tan(thetaSig)*(1.+TMath::Power(tan(thetaSig),2))*nfreon
/ngas*cos(etaStepAvg)/cos(thetaSig);
areaBkg = stepEta*funBkg;
Float_t e1, e2, f1, f2;
Float_t nfreon[2], nquartz[2];
Int_t times;
-
+ Float_t pointsOnCathode[3];
Float_t phpad, thetacer[2];
Float_t bandradius[2], padradius;
f1 = 46.411;
f2 = 228.71;
-
phpad = PhiPad();
for (times=0; times<=1; times++) {
nfreon[times] = a+b*energy[times];
+ //nfreon[times] = 1;
- nquartz[times] = sqrt(1+(f1/(pow(e1,2)-pow(energy[times],2)))+
- (f2/(pow(e2,2)-pow(energy[times],2))));
+ nquartz[times] = sqrt(1+(f1/(TMath::Power(e1,2)-TMath::Power(energy[times],2)))+
+ (f2/(TMath::Power(e2,2)-TMath::Power(energy[times],2))));
- beta[times] = imp[times]/sqrt(pow(imp[times],2)+pow(mass[times],2));
+ beta[times] = imp[times]/sqrt(TMath::Power(imp[times],2)+TMath::Power(mass[times],2));
thetacer[times] = CherenkovAngle( nfreon[times], beta[times]);
bandradius[times] = DistanceFromMip( nfreon[times], nquartz[times],
- emissPointLength[times],
- thetacer[times], phpad);
- }
+ emissPointLength[times],
+ thetacer[times], phpad, pointsOnCathode);
+ //printf(" ppp %f %f %f \n",pointsOnCathode);
+}
bandradius[0] -= 1.6;
bandradius[1] += 1.6;
- padradius = sqrt(pow(fXpad,2)+pow(fYpad,2));
+ padradius = sqrt(TMath::Power(fXpad,2)+TMath::Power(fYpad,2));
// printf(" rmin %f r %f rmax %f \n",bandradius[0],padradius,bandradius[1]);
if(padradius>=bandradius[0] && padradius<=bandradius[1]) return 1;
Float_t AliRICHPatRec::DistanceFromMip(Float_t nfreon, Float_t nquartz,
Float_t emissPointLength, Float_t thetacer,
- Float_t phpad)
+ Float_t phpad, Float_t pointsOnCathode[3])
{
// Find the distance to MIP impact
TVector3 gapExitPhot(1,1,1) ;
Float_t magGapExitPhot;
//
- TVector3 fPhotocatExitPhot(1,1,1);
+ TVector3 PhotocatExitPhot(1,1,1);
Double_t theta2;
Double_t thetarad , phirad ;
Double_t thetaquar, phiquar;
radExitPhot = vectEmissPointLength + radExitPhot2;
thetarad = radExitPhot.Theta();
-
phirad = radExitPhot.Phi(); //check on the original file //
thetaquar = SnellAngle( nfreon, nquartz, theta2);
gapExitPhot.SetTheta(thetagap);
gapExitPhot.SetPhi(phigap);
- fPhotocatExitPhot = radExitPhot + quarExitPhot + gapExitPhot;
+ PhotocatExitPhot = radExitPhot + quarExitPhot + gapExitPhot;
- distanceValue = sqrt(pow(fPhotocatExitPhot(0),2)
- +pow(fPhotocatExitPhot(1),2));
- return distanceValue ;
-}
+ distanceValue = sqrt(TMath::Power(PhotocatExitPhot(0),2)
+ +TMath::Power(PhotocatExitPhot(1),2));
+ pointsOnCathode[0] = (Float_t) PhotocatExitPhot(0) + fXshift - fTrackLoc[0];
+ pointsOnCathode[1] = (Float_t) PhotocatExitPhot(1) + fYshift - fTrackLoc[1];
+ pointsOnCathode[2] = (Float_t) PhotocatExitPhot(2);
+
+ //printf(" point in Distance.2. %f %f %f \n",pointsOnCathode[0],pointsOnCathode[1],pointsOnCathode[2]);
+
+ return distanceValue;
+
+ }
Float_t AliRICHPatRec::PhiPad()
{
}
}
-Float_t AliRICHPatRec::CherenkovRingDrawing(Float_t fixedthetacer)
+/*void AliRICHPatRec::CerenkovRingDrawing()
{
//to draw Cherenkov ring by known Cherenkov angle
- Int_t nmaxdegrees, nstepdegrees;
- Float_t phpad, thetacer;
+ Int_t nmaxdegrees;
+ Int_t Nphpad;
+ Float_t phpad;
Float_t nfreonave, nquartzave;
Float_t aveEnerg;
Float_t energy[2];
Float_t e1, e2, f1, f2;
Float_t bandradius;
- Float_t coordPadRing;
-
+
//parameters to calculate freon window refractive index vs. energy
+
Float_t a = 1.177;
Float_t b = 0.0172;
//parameters to calculate quartz window refractive index vs. energy
-/*
- Energ[0] = 5.6;
- Energ[1] = 7.7;
-*/
+
energy[0] = 5.0;
energy[1] = 8.0;
e1 = 10.666;
f2 = 228.71;
- nmaxdegrees = 360;
-
- nstepdegrees = 36;
+ nmaxdegrees = 36;
- for (phpad=0; phpad<nmaxdegrees;phpad++) {
+ for (Nphpad=0; Nphpad<nmaxdegrees;Nphpad++) {
+
+ phpad = (360./(Float_t)nmaxdegrees)*(Float_t)Nphpad;
aveEnerg = (energy[0]+energy[1])/2.;
nfreonave = a+b*aveEnerg;
- nquartzave = sqrt(1+(f1/(pow(e1,2)-pow(aveEnerg,2)))+
- (f2/(pow(e2,2)-pow(aveEnerg,2))));
-
- thetacer = fixedthetacer;
+ nquartzave = sqrt(1+(f1/(TMath::Power(e1,2)-TMath::Power(aveEnerg,2)))+
+ (f2/(TMath::Power(e2,2)-TMath::Power(aveEnerg,2))));
bandradius = DistanceFromMip(nfreonave, nquartzave,
- fEmissPoint,thetacer, phpad);
+ fEmissPoint,fThetaCerenkov, phpad);
- coordPadRing=fPhotocatExitPhot;
+ fCoordEllipse[0][Nphpad] = fOnCathode[0];
+ fCoordEllipse[1][Nphpad] = fOnCathode[1];
+ printf(" values %f %f \n",fOnCathode[0],fOnCathode[1]);
- phpad = (nmaxdegrees/nstepdegrees)*phpad;
-
- return coordPadRing;
}
- return coordPadRing;
-}
+}*/
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff