Int_t totPadsPerChamber = AliRICHParam::NpadsX()*AliRICHParam::NpadsY();
- Int_t nDigCh[kNchambers]={0,0,0,0,0,0,0};
- Int_t iChHits[kNchambers]={0,0,0,0,0,0,0};
- Int_t nPrim[kNchambers]={0,0,0,0,0,0,0};
- Int_t nSec[kNchambers]={0,0,0,0,0,0,0};
- for(Int_t iEvtN=iEvtNmin;iEvtN<iEvtNmax;iEvtN++){
- if(iEvtN%10==0) AliInfo(Form("events processed %i",iEvtN));
+ for(Int_t iEvtN=iEvtNmin;iEvtN<iEvtNmax;iEvtN++){
+ Int_t nDigCh[kNchambers]={0,0,0,0,0,0,0};
+ Int_t iChHits[kNchambers]={0,0,0,0,0,0,0};
+ Int_t nPrim[kNchambers]={0,0,0,0,0,0,0};
+ Int_t nSec[kNchambers]={0,0,0,0,0,0,0};
+ AliInfo(Form("events processed %i",iEvtN));
if(GetLoader()->GetRunLoader()->GetEvent(iEvtN)) return;
AliStack *pStack = GetLoader()->GetRunLoader()->Stack();
for(Int_t iPrimN=0;iPrimN<GetLoader()->TreeH()->GetEntries();iPrimN++){//prims loop
}
}
GetLoader()->TreeD()->GetEntry(0);
- for(Int_t iChamber=1;iChamber<=kNchambers;iChamber++) nDigCh[iChamber-1]= Digits(iChamber)->GetEntries();
- }
- for(Int_t iChamber=1;iChamber<=kNchambers;iChamber++){
- Double_t occupancy = (Double_t)nDigCh[iChamber-1]/(Double_t)totPadsPerChamber;
- Info("Occupancy","for chamber %i = %4.2f %% and charged prim tracks %i and sec. tracks %i with total %i",
+ for(Int_t iChamber=1;iChamber<=kNchambers;iChamber++) {
+ nDigCh[iChamber-1]= Digits(iChamber)->GetEntries();
+ Double_t occupancy = (Double_t)nDigCh[iChamber-1]/(Double_t)totPadsPerChamber;
+ Info("Occupancy","for chamber %i = %4.2f %% and charged prim tracks %i and sec. tracks %i with total %i",
iChamber,occupancy*100.,nPrim[iChamber-1],nSec[iChamber-1],iChHits[iChamber-1]);
- }
+ }
+ }
GetLoader()->UnloadHits();
GetLoader()->UnloadDigits();
GetLoader()->GetRunLoader()->UnloadHeader();
pFile->Write();pFile->Close();
}
//__________________________________________________________________________________________________
+void AliRICH::RichAna()
+{
+ TFile *pFile=TFile::Open("AliESDs.root","read");
+ if(!pFile || !pFile->IsOpen()) {AliInfo("ESD file not open.");return;} //open AliESDs.root
+ TTree *pTree = (TTree*) pFile->Get("esdTree");
+ if(!pTree){AliInfo("ESD not found.");return;} //get ESD tree
+ AliInfo("ESD found. Go ahead!");
+
+ AliMagF * magf = gAlice->Field();
+ AliTracker::SetFieldMap(magf,kTRUE);
+// AliTracker::SetFieldMap(magf);
+
+ TFile *pFileRA = new TFile("$(HOME)/RichAna.root","RECREATE","RICH Pattern Recognition");
+ TNtupleD *hn = new TNtupleD("hn","ntuple","Pmod:Charge:TrackTheta:TrackPhi:MinX:MinY:ThetaCerenkov:NPhotons:ChargeMIP:Chamber:TOF:LengthTOF:prob1:prob2:prob3");
+
+ AliESD *pESD=new AliESD; pTree->SetBranchAddress("ESD", &pESD);
+// for(Int_t iEvtN=0;iEvtN<GetLoader()->GetRunLoader()->GetNumberOfEvents();iEvtN++) {
+ for(Int_t iEvtN=0;iEvtN<2;iEvtN++) {
+ pTree->GetEvent(iEvtN);
+ AliRICH *pRich=((AliRICH*)gAlice->GetDetector("RICH"));
+ pRich->GetLoader()->GetRunLoader()->GetEvent(iEvtN);
+ pRich->GetLoader()->LoadRecPoints();
+ pRich->GetLoader()->TreeR()->GetEntry(0);
+//Pattern recognition started
+ if(pESD->GetNumberOfTracks()) {
+ Int_t iNtracks=pESD->GetNumberOfTracks();
+ Info("RichAna",Form("Start with %i tracks",iNtracks));
+ for(Int_t iTrackN=0;iTrackN<iNtracks;iTrackN++){//ESD tracks loop
+ if(iTrackN%100==0)Info("RichAna",Form("Track %i to be processed",iTrackN));
+ AliRICHTracker *pTrRich = new AliRICHTracker();
+ pTrRich->RecWithESD(pESD,pRich,iTrackN);
+ AliESDtrack *pTrack = pESD->GetTrack(iTrackN);// get next reconstructed track
+ Double_t dx,dy;
+ Double_t hnvec[20];
+ pTrack->GetRICHdxdy(dx,dy);
+ hnvec[0]=pTrack->GetP();
+ hnvec[1]=pTrack->GetSign();
+// cout << " Track momentum " << pTrack->GetP() << " charge " << pTrack->GetSign() << endl;
+
+ pTrack->GetRICHthetaPhi(hnvec[2],hnvec[3]);
+ pTrack->GetRICHdxdy(hnvec[4],hnvec[5]);
+ hnvec[6]=pTrack->GetRICHsignal();
+ hnvec[7]=pTrack->GetRICHnclusters();
+ hnvec[9]=pTrack->GetRICHcluster()/1000000;
+ hnvec[8]=pTrack->GetRICHcluster()-hnvec[9]*1000000;
+ hnvec[10]=pTrack->GetTOFsignal();
+ hnvec[11]=pTrack->GetIntegratedLength();
+ Double_t prob[5];
+ pTrack->GetRICHpid(prob);
+ hnvec[12]=prob[0]+prob[1]+prob[2];
+ hnvec[13]=prob[3];
+ hnvec[14]=prob[4];
+ hn->Fill(hnvec);
+ }
+ }
+ Info("RichAna","Pattern Recognition done for event %i",iEvtN);
+ }
+ pFileRA->Write();pFileRA->Close();// close RichAna.root
+ delete pESD; pFile->Close();//close AliESDs.root
+}
+//__________________________________________________________________________________________________
void AliRICH::DisplayEvent(Int_t iEvtNmin,Int_t iEvtNmax)const
{
TH2F *pDigitsH2[8];
pHitsH2->SetMarkerColor(kRed); pHitsH2->SetMarkerStyle(29); pHitsH2->SetMarkerSize(0.4);
ReadESD(iEventN,iChamber);
pHitsH2->Draw();
- ReadESD(iEventN,iChamber);
+// ReadESD(iEventN,iChamber);
AliRICHParam::DrawSectors();
TLatex l; l.SetNDC(); l.SetTextSize(0.02);
if(!isHits) {l.SetTextColor(kRed) ;l.DrawLatex(0.1,0.01,"No Hits" );}
void PrintTracks (Int_t iEvent=0); //prints a list of tracks for a given event
void CheckPR ()const; //utility-> run staff for stack ??????
void ReadESD(Int_t iEventN, Int_t iChamber)const;
+ void RichAna (); //utility-> reconstruct ESD rings outside ESD
void DrawRing(TVector3 entrance,TVector3 vectorTrack,Double_t thetaCer)const;
protected:
#include <TMath.h>
#include <TRotation.h>
#include <TVector3.h>
+#include <TH1F.h>
#include "AliRICH.h"
#include "AliRICHParam.h"
{
// main ctor
SetFreonScaleFactor(1);
- fIsWEIGHT = kFALSE;
+// fIsWEIGHT = kFALSE;
+ fIsWEIGHT = kTRUE;
fThetaBin=750; fThetaMin = 0.0; fThetaMax = 0.75;
fDTheta = 0.001; fWindowWidth = 0.045;
+ fMinNumPhots = 3;
fRadiatorWidth = AliRICHParam::Zfreon();
fQuartzWidth = AliRICHParam::Zwin();
fGapWidth = AliRICHParam::Freon2Pc() - fRadiatorWidth - fQuartzWidth;
- fXmin = -AliRICHParam::PcSizeX()/2.;
- fXmax = AliRICHParam::PcSizeX()/2.;
- fYmin = -AliRICHParam::PcSizeY()/2.;
- fYmax = AliRICHParam::PcSizeY()/2.;
+ fXmin = 0.;
+ fXmax = AliRICHParam::PcSizeX();
+ fYmin = 0.;
+ fYmax = AliRICHParam::PcSizeY();
SetTrackTheta(pHelix->Ploc().Theta());
SetTrackPhi(pHelix->Ploc().Phi());
SetMipIndex(iMipId);
Int_t zone = CheckDetectorAcceptance();
+ AliDebug(1,Form("acceptance to detector zone -> %d",zone));
if (zone != 0)
{
//__________________________________________________________________________________________________
void AliRICHRecon::HoughResponse()
{
- //Hough response
-
-// Implement Hough response pat. rec. method
-
- Float_t *hCSspace;
-
- int bin=0;
- int bin1=0;
- int bin2=0;
- int i, j, k, nCorrBand;
- float hcs[750],hcsw[750];
- float angle, weight;
- float lowerlimit,upperlimit;
-
- float etaPeak[100];
-
- int nBin;
-
- float etaPeakPos = -1;
-
- Int_t etaPeakCount = -1;
-
- Float_t thetaCerenkov = 0.;
-
- nBin = (int)(0.5+fThetaMax/(fDTheta));
- nCorrBand = (int)(0.5+ fWindowWidth/(2 * fDTheta));
-
- memset ((void *)hcs, 0, fThetaBin*sizeof(float));
- memset ((void *)hcsw, 0, fThetaBin*sizeof(float));
-
+ TH1F *phots = new TH1F("phots","phots",750,0.,0.75);
+ TH1F *photsw = new TH1F("photsw","photsw",750,0.,0.75);
+ TH1F *resultw = new TH1F("resultw","resultw",750,0.,0.75);
Int_t nPhotons = GetPhotonsNumber();
-
- Int_t weightFlag = 0;
-
- for (k=0; k< nPhotons; k++) {
-
- SetPhotonIndex(k);
-
- angle = GetPhotonEta();
-
+ Int_t nBin = (Int_t)(fThetaMax/fDTheta);
+ Int_t nCorrBand = (Int_t)(fWindowWidth/(2*fDTheta));
+ AliDebug(1,Form("Ring reconstruction for track with theta %f",GetTrackTheta()*TMath::RadToDeg()));
+ for (Int_t kPhot=0; kPhot< nPhotons; kPhot++){
+ SetPhotonIndex(kPhot);
+ Double_t angle = GetPhotonEta();
if(angle == -999.) continue;
-
- if (angle>=fThetaMin && angle<= fThetaMax)
-
- {
-
- bin = (int)(0.5+angle/(fDTheta));
-
- bin1= bin-nCorrBand;
- bin2= bin+nCorrBand;
-
- // calculate weights
-
- if(fIsWEIGHT)
- {
- lowerlimit = ((Float_t)bin1)*fDTheta + 0.5*fDTheta;
- SetThetaOfRing(lowerlimit);
- FindAreaAndPortionOfRing();
- Float_t area1 = GetAreaOfRing();
-
- upperlimit = ((Float_t)bin2)*fDTheta + 0.5*fDTheta;
- SetThetaOfRing(upperlimit);
- FindAreaAndPortionOfRing();
- Float_t area2 = GetAreaOfRing();
-
- // cout << "lowerlimit" << lowerlimit << "upperlimit " << upperlimit << endl;
- Float_t diffarea = area2 - area1;
-
- if(diffarea>0)
- {
- weight = 1./(area2-area1);
- }
- else
- {
- weightFlag = 1;
- weight = 1.;
- }
-
- // cout <<" low "<< lowerlimit << " up " << upperlimit <<
- // " area1 " << area1 << " area2 " << area2 << " weight " << weight << endl;
-
- }
- else
- {
- weight = 1.;
- }
-
- SetPhotonWeight(weight);
-
- // cout << "weight..." << weight << endl;
-
-
- if (bin1<0) bin1=0;
- if (bin2>nBin) bin2=nBin;
-
- for (j=bin1; j<bin2; j++)
- {
- hcs[j] += 1;
- hcsw[j] += weight;
- }
- }
- }
-
-
- if(weightFlag == 0)
- {
- hCSspace = hcsw;
- }
- else
- {
- hCSspace = hcs;
- // cout << " probems with weight...normal procedure adopted " << endl;
- }
-
- HoughFiltering(hCSspace);
-
- for (bin=0; bin <nBin; bin++) {
- angle = (bin+0.5) * (fDTheta);
- if (hCSspace[bin] && hCSspace[bin] > etaPeakPos) {
- etaPeakCount = 0;
- etaPeakPos = hCSspace[bin];
- etaPeak[0]=angle;
- }
- else {
- if (hCSspace[bin] == etaPeakPos) {
- etaPeak[++etaPeakCount] = angle;
- }
- }
- }
-
- for (i=0; i<etaPeakCount+1; i++) {
- thetaCerenkov += etaPeak[i];
- }
- if (etaPeakCount>=0) {
- thetaCerenkov /= etaPeakCount+1;
- fThetaPeakPos = etaPeakPos;
- }
-
- SetThetaCerenkov(thetaCerenkov);
-}
-//__________________________________________________________________________________________________
-void AliRICHRecon::HoughFiltering(float hcs[])
-{
- // filter for Hough
-
-// hough filtering
-
- float hcsFilt[750];
- float k[5] = {0.05, 0.25, 0.4, 0.25, 0.05};
- int nx, i, nxDx;
- int sizeHCS;
- int nBin;
-
- nBin = (int)(1+fThetaMax/fDTheta);
- sizeHCS = fThetaBin*sizeof(float);
-
- memset ((void *)hcsFilt, 0, sizeHCS);
-
- for (nx = 0; nx < nBin; nx++) {
- for (i = 0; i < 5; i++) {
- nxDx = nx + (i-2);
- if (nxDx> -1 && nxDx<nBin)
- hcsFilt[nx] += hcs[nxDx] * k[i];
- }
- }
-
- for (nx = 0; nx < nBin; nx++) {
- hcs[nx] = hcsFilt[nx];
- }
-}
+ phots->Fill(angle);
+ Int_t bin = (Int_t)(0.5+angle/(fDTheta));
+ Double_t weight=1.;
+ if(fIsWEIGHT){
+ Double_t lowerlimit = ((Float_t)bin)*fDTheta - 0.5*fDTheta;
+ SetThetaOfRing(lowerlimit);
+ FindAreaAndPortionOfRing();
+ Float_t area1 = GetAreaOfRing();
+ Double_t upperlimit = ((Float_t)bin)*fDTheta + 0.5*fDTheta;
+ SetThetaOfRing(upperlimit);
+ FindAreaAndPortionOfRing();
+ Float_t area2 = GetAreaOfRing();
+ AliDebug(1,Form("lowerlimit %f area %f ; upperlimit %f area %f",lowerlimit,area1,upperlimit,area2));
+ Float_t diffarea = area2 - area1;
+ if(diffarea>0){weight = 1./(area2-area1);}else{weight = 1.;}
+ }
+ photsw->Fill(angle,weight);
+ SetPhotonWeight(weight);
+ }
+ for (Int_t i=1; i<=nBin;i++){
+ Int_t bin1= i-nCorrBand;
+ Int_t bin2= i+nCorrBand;
+ if(bin1<1) bin1=1;
+ if(bin2>nBin)bin2=nBin;
+ Double_t sumPhots=phots->Integral(bin1,bin2);
+ if(sumPhots<fMinNumPhots) continue; // cut on minimum n. of photons per ring
+ Double_t sumPhotsw=photsw->Integral(bin1,bin2);
+ resultw->Fill((Float_t)(i*fDTheta),sumPhotsw);
+}
+// evaluate the "BEST" thetacerenkov....
+ Float_t *pVec = resultw->GetArray();
+ Int_t locMax = TMath::LocMax(nBin,pVec);
+ SetThetaCerenkov((Double_t)(locMax*fDTheta+0.5*fDTheta));
+
+// Reset and delete objects
+ phots->Delete();photsw->Delete();resultw->Delete();
+}//HoughResponse
//__________________________________________________________________________________________________
void AliRICHRecon::FindThetaCerenkov()
{
}
SetHoughPhotons(nPhotonHough);
}
-
Int_t fCandidatePhotonsNumber; // number of candidate photons
Int_t fHoughPhotons; // n. photons after Hough
Int_t fFittedHoughPhotons; // n. photons after Hough and after minimization
+ Int_t fMinNumPhots; // minimum number of photons for a given ring
Float_t fTrackTheta; // Theta of track at RICH
Float_t fTrackPhi; // Phi of track at RICH
// 1. AliESD - reconstructed tracks are used
// 2. RICH private ntuple for debug- stack particles used instead of reconstructed tracks
AliDebug(1,"Start pattern recognition");
- if(pESD->GetNumberOfTracks())
- RecWithESD(pESD);
- else
- RecWithStack(0);
+ if(pESD->GetNumberOfTracks()) {
+ Int_t iNtracks=pESD->GetNumberOfTracks();
+ AliDebug(1,Form("Start with %i tracks",iNtracks));
+ AliRICH *pRich=((AliRICH*)gAlice->GetDetector("RICH"));
+ for(Int_t iTrackN=0;iTrackN<iNtracks;iTrackN++){//ESD tracks loop
+ RecWithESD(pESD,pRich,iTrackN);
+ }
+ }
+ else RecWithStack(0);
AliDebug(1,"Stop pattern recognition");
return 0; // error code: 0=no error;
}//PropagateBack()
//__________________________________________________________________________________________________
-void AliRICHTracker::RecWithESD(AliESD *pESD)
+void AliRICHTracker::RecWithESD(AliESD *pESD,AliRICH *pRich,Int_t iTrackN)
{
//recontruction from ESD- primary way to reconstruct particle ID signal from tracks provided by core detectors
- Int_t iNtracks=pESD->GetNumberOfTracks();
- Double_t b=GetFieldMap()->SolenoidField()/10;// magnetic field in Tesla
- AliDebug(1,Form("Start with %i tracks in %f Tesla field",iNtracks,b));
-
- AliRICH *pRich=((AliRICH*)gAlice->GetDetector("RICH"));
-
- for(Int_t iTrackN=0;iTrackN<iNtracks;iTrackN++){//ESD tracks loop
+ Double_t fField=GetFieldMap()->SolenoidField()/10;// magnetic field in Tesla
AliESDtrack *pTrack = pESD->GetTrack(iTrackN);// get next reconstructed track
// if((pTrack->GetStatus()&AliESDtrack::kTOFout)==0) continue; //ignore tracks not recontructed by TOF
// pTrack->GetXYZ(xb);
// pTrack->GetPxPyPz(pb);
Int_t status=pTrack->GetStatus()&AliESDtrack::kTOFout;//get running track parameters
- Int_t charge = (Int_t)(-TMath::Sign(1.,pTrack->GetSign()*b));
+ Int_t charge = (Int_t)(-TMath::Sign(1.,pTrack->GetSign()*fField));
AliDebug(1,Form("Track %i pmod=%f charge=%i stat=%i",iTrackN,pTrack->GetP(),charge,status));
- AliRICHHelix helix(pTrack->X3(),pTrack->P3(),charge,b);
+ AliRICHHelix helix(pTrack->X3(),pTrack->P3(),charge,fField);
Int_t iChamber=helix.RichIntersect(pRich->P());
AliDebug(1,Form("intersection with %i chamber found",iChamber));
- if(!iChamber) continue;//intersection with no chamber found
+ if(!iChamber) return;//intersection with no chamber found
//find MIP cluster candidate (cluster which is closest to track intersection point)
Double_t distMip=9999,distX=0,distY=0; //min distance between clusters and track position on PC
Int_t iMipId=0; //index of that min distance cluster
if(distMip>1||chargeMip<100) {
//track not accepted for pattern recognition
pTrack->SetRICHsignal(-999.); //to be improved by flags...
- continue;
+ return;
}
//
AliRICHRecon recon(&helix,pRich->Clusters(iChamber),iMipId); //actual job is done there
Double_t thetaCerenkov=recon.ThetaCerenkov(); //search for mean Cerenkov angle for this track
- pTrack->SetRICHcluster(iMipId+100000*iChamber);
+ pTrack->SetRICHcluster(((Int_t)chargeMip)+1000000*iChamber);
pTrack->SetRICHdxdy(distX,distY);
pTrack->SetRICHthetaPhi(helix.Ploc().Theta(),helix.Ploc().Phi());
pTrack->SetRICHsignal(thetaCerenkov);
Double_t richPID[AliPID::kSPECIES];
for (Int_t iPart=0;iPart<AliPID::kSPECIES;iPart++) {
sigmaPID[iPart] = 0;
- richPID[iPart] = 0;
for(Int_t iphot=0;iphot<pRich->Clusters(iChamber)->GetEntries();iphot++) {
recon.SetPhotonIndex(iphot);
if(recon.GetPhotonFlag() == 2) {
- Double_t sigma2 = AliRICHParam::SigmaSinglePhoton(iPart,pTrack->GetP(),recon.GetTrackTheta(),recon.GetPhiPoint()-recon.GetTrackPhi()).Mag2();
- if (sigma2 >0)
- sigmaPID[iPart] += 1/sigma2;
+ Double_t sigma = AliRICHParam::SigmaSinglePhoton(iPart,pTrack->GetP(),recon.GetTrackTheta(),recon.GetPhiPoint()-recon.GetTrackPhi()).Mag();
+ sigmaPID[iPart] += 1/(sigma*sigma);
}
}
if (sigmaPID[iPart]>0)
CalcProb(thetaCerenkov,pTrack->GetP(),sigmaPID,richPID);
pTrack->SetRICHpid(richPID);
AliDebug(1,Form("PROBABILITIES ---> %f - %f - %f - %f - %f",richPID[0],richPID[1],richPID[2],richPID[3],richPID[4]));
- }
-
-
- }//ESD tracks loop
+ }
AliDebug(1,"Stop.");
} //RecWithESD
//__________________________________________________________________________________________________
class AliCluster;
class AliESD;
+class AliRICH;
class TTree;
class AliRICHTracker : public AliTracker
void UnloadClusters() {AliDebug(1,"Start.");} //pure virtual from AliTracker
AliCluster *GetCluster(Int_t )const {AliDebug(1,"Start.");return 0;} //pure virtual from AliTracker
Int_t PropagateBack(AliESD *); //pure virtual from AliTracker
- void RecWithESD(AliESD *); //recon with ESD
+ void RecWithESD(AliESD *,AliRICH *,Int_t iTrackN); //recon with ESD
void RecWithStack(TNtupleD *hn); //recon from Stack in case ESD empty
void CalcProb(Double_t thetaCer,Double_t pmod,Double_t *sigmaPID, Double_t *richPID); // calculate pid for RICH
Int_t LoadClusters(TTree *); //pure virtual from AliTracker
protected:
+ Double_t fField; // magnetic field stored
+
ClassDef(AliRICHTracker,0)
};//class AliRICHTracker
pMenu->AddButton("Print occupancy" ,"r->OccupancyPrint(-1);" ,"????");
pMenu->AddButton("Event Summary " ,"r->SummaryOfEvent();" ,"????");
pMenu->AddButton("Hits plots" ,"r->ControlPlots()" ,"????");
- pMenu->AddButton("Recon with stack" ,"r->CheckPR()" , "Create RSR.root with ntuple hn");
+ pMenu->AddButton("Recon with stack" ,"r->CheckPR()", "Create RSR.root with ntuple hn");
+ pMenu->AddButton("RichAna" ,"r->RichAna()", "Create RichAna.root with ntuple hn");
pMenu->Show();
}//TestMenu()
//__________________________________________________________________________________________________
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff