FILE *pChromErr, *pGeomErr, *pLocErr;
if(!count) {
- AliInfoClass("reading RICH error parameters...");
+ AliInfoGeneral("ReadErrFiles","reading RICH error parameters...");
pChromErr = fopen(Form("%s/RICH/RICHConfig/SigmaChromErr.txt",gSystem->Getenv("ALICE_ROOT")),"r");
pGeomErr = fopen(Form("%s/RICH/RICHConfig/SigmaGeomErr.txt",gSystem->Getenv("ALICE_ROOT")),"r");
pLocErr = fopen(Form("%s/RICH/RICHConfig/SigmaLocErr.txt",gSystem->Getenv("ALICE_ROOT")),"r");
fgErrLoc[2][i] = l2;
fgErrLoc[3][i] = l3;
}
- AliInfoClass("DONE successfully!");
+ AliInfoGeneral("ReadErrFiles","DONE successfully!");
fclose(pChromErr);
fclose(pGeomErr);
fclose(pLocErr);
// Pattern recognition method based on Hough transform
// Return theta Cerenkov for a given track and list of clusters which are set in ctor
- if(fpClusters->GetEntries()==0) return -1;//no clusters at all for a given track
+ if(fpClusters->GetEntries()==0) return -10;//no clusters at all for a given track
Bool_t kPatRec = kFALSE;
AliDebug(1,Form("---Track Parameters--- Theta: %f , Phi: %f ",GetTrackTheta()*TMath::RadToDeg(),GetTrackPhi()*TMath::RadToDeg()));
Int_t candidatePhotons = 0;
- SetThetaCerenkov(999.);
+ SetThetaCerenkov(-1);
SetHoughPhotons(0);
SetHoughPhotonsNorm(0);
if(nPhotonHough < 1)
{
- SetThetaCerenkov(999.);
- SetHoughPhotonsNorm(0.);
+ SetThetaCerenkov(-1);
+ SetHoughPhotonsNorm(0);
return -1;
}
#include <AliRawReader.h> //Reconstruct(...)
#include <AliRun.h> //ConvertDigits uses gAlice
#include <AliESD.h> //RichAna()
+#include <AliStack.h> //RichAna()
#include <TFile.h> //RichAna()
#include <TMinuit.h> //Dig2Clu()
+#include <TParticle.h> //TParticle()
ClassImp(AliRICHReconstructor)
void AliRICHReconstructor::Reconstruct(AliRunLoader *pAL,AliRawReader* pRR)const
if((pDig=UseDig(x[i],y[i],pDigList,pDigMap))) FormCluster(pClu,pDig,pDigList,pDigMap);
}//FormCluster()
//__________________________________________________________________________________________________
-void AliRICHReconstructor::RichAna(Int_t nNevMax,Bool_t askPatRec)
+void AliRICHReconstructor::CheckPR()
+{
+//Pattern recognition with stack particles
+ TFile *pFile = new TFile("$(HOME)/RPR.root","RECREATE","RICH Pattern Recognition");
+ TNtupleD *hn = new TNtupleD("hn","ntuple","Pmod:Charge:TrackTheta:TrackPhi:TrackX:TrackY:MinX:MinY:ChargeMIP:ThetaCerenkov:NPhotons:MipIndex:Chamber:Particle");
+// printf("\n\n");
+// printf("Pattern Recognition done for event %5i",0);
+ AliRICH *pRich=((AliRICH*)gAlice->GetDetector("RICH"));
+ AliMagF * magf = gAlice->Field();
+ AliTracker::SetFieldMap(magf,kTRUE);
+ for(Int_t iEvtN=0;iEvtN<pRich->GetLoader()->GetRunLoader()->GetNumberOfEvents();iEvtN++) {
+ pRich->GetLoader()->GetRunLoader()->GetEvent(iEvtN);
+ AliRICHTracker *tr = new AliRICHTracker();
+ tr->RecWithStack(hn);
+ AliInfoClass(Form("Pattern Recognition done for event %i \b",iEvtN));
+// printf("\b\b\b\b\b%5i",iEvtN+1);
+ }
+ printf("\n\n");
+ pFile->Write();pFile->Close();
+}
+//__________________________________________________________________________________________________
+void AliRICHReconstructor::RichAna(Int_t iNevMin,Int_t iNevMax,Bool_t askPatRec)
{
TFile *pFile=TFile::Open("AliESDs.root","read");
- if(!pFile || !pFile->IsOpen()) {return;} //open AliESDs.root
+ if(!pFile || !pFile->IsOpen()) {AliInfoClass("ESD file not open.");return;} //open AliESDs.root
TTree *pTree = (TTree*) pFile->Get("esdTree");
- if(!pTree){return;} //get ESD tree
-
+ if(!pTree){AliInfoClass("ESD not found.");return;} //get ESD tree
+ AliInfoClass("ESD found. Go ahead!");
+
AliRICH *pRich=((AliRICH*)gAlice->GetDetector("RICH"));
+
+ AliMagF * magf = gAlice->Field();
+ AliTracker::SetFieldMap(magf,kTRUE);
+ pRich->GetLoader()->GetRunLoader()->LoadHeader();
+ pRich->GetLoader()->GetRunLoader()->LoadKinematics();
+ TString var1 = "Pmod:Charge:TrackTheta:TrackPhi:MinX:MinY:ThetaCerenkov:NPhotons:";
+ TString var2 = "ChargeMIP:Chamber:TOF:LengthTOF:prob1:prob2:prob3:";
+ TString var3 = "ErrPar1:ErrPar2:ErrPar3:Th1:Th2:Th3:nPhotBKG:pdgCode";
+ TString varList = var1+var2+var3;
- AliMagF * magf = gAlice->Field(); AliTracker::SetFieldMap(magf,kTRUE); // AliTracker::SetFieldMap(magf);
+ Double_t dx,dy;
+ Double_t hnvec[30];
- 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:ErrPar1:ErrPar2:ErrPar3:Th1:Th2:Th3:nPhotBKG");
- if(nNevMax==0) nNevMax=999999;
- if(pRich->GetLoader()->GetRunLoader()->GetNumberOfEvents()<nNevMax) nNevMax = pRich->GetLoader()->GetRunLoader()->GetNumberOfEvents();
+// TFile *pFileRA = new TFile("$(HOME)/RichAna.root","RECREATE","RICH Pattern Recognition");
+ TFile *pFileRA = new TFile("./RichAna.root","RECREATE","RICH Pattern Recognition");
+ TNtupleD *hn = new TNtupleD("hn","ntuple",varList);
+ if(iNevMin<0) iNevMin=0;
+ if(iNevMin>iNevMax) {iNevMin=0;iNevMax=0;}
+ if(iNevMax==0) iNevMax=999999;
+ if(pRich->GetLoader()->GetRunLoader()->GetNumberOfEvents()<iNevMax) iNevMax = pRich->GetLoader()->GetRunLoader()->GetNumberOfEvents();
AliESD *pESD=new AliESD; pTree->SetBranchAddress("ESD", &pESD);
- for(Int_t iEvtN=0;iEvtN<nNevMax;iEvtN++) {
+ for(Int_t iEvtN=iNevMin;iEvtN<iNevMax;iEvtN++) {
pTree->GetEvent(iEvtN);
pRich->GetLoader()->GetRunLoader()->GetEvent(iEvtN);
pRich->GetLoader()->LoadRecPoints();
pRich->GetLoader()->TreeR()->GetEntry(0);
+ AliStack *pStack = pRich->GetLoader()->GetRunLoader()->Stack();
//Pattern recognition started
if(pESD->GetNumberOfTracks()) {
Int_t iNtracks=pESD->GetNumberOfTracks();
AliRICHTracker *pTrRich = new AliRICHTracker();
if(askPatRec==kTRUE) pTrRich->RecWithESD(pESD,pRich,iTrackN);
AliESDtrack *pTrack = pESD->GetTrack(iTrackN);// get next reconstructed track
- Double_t dx,dy;
- Double_t hnvec[30];
+
+ Int_t lab=TMath::Abs(pTrack->GetLabel());
+ TParticle *pPart=pStack->Particle(lab);
+ Int_t code=pPart->GetPdgCode();
+
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]);
if(cosThetaTh>=1) continue;
hnvec[18+i]= TMath::ACos(cosThetaTh);
}
- hnvec[21]=pTrRich->fnPhotBKG;
+ if(askPatRec==kTRUE) hnvec[21]=pTrRich->fnPhotBKG; else hnvec[21]=0;
+ hnvec[22]=code;
hn->Fill(hnvec);
}
}
}
pFileRA->Write();pFileRA->Close();// close RichAna.root
delete pESD; pFile->Close();//close AliESDs.root
-}//RichAna()
-//__________________________________________________________________________________________________
-void AliRICHReconstructor::CheckPR()
-{
-//Pattern recognition with stack particles
- TFile *pFile = new TFile("$(HOME)/RPR.root","RECREATE","RICH Pattern Recognition");
- TNtupleD *hn = new TNtupleD("hn","ntuple","Pmod:Charge:TrackTheta:TrackPhi:TrackX:TrackY:MinX:MinY:ChargeMIP:ThetaCerenkov:NPhotons:MipIndex:Chamber:Particle");
- AliRICH *pRich=((AliRICH*)gAlice->GetDetector("RICH"));
- AliMagF * magf = gAlice->Field();
- AliTracker::SetFieldMap(magf,kTRUE);
- for(Int_t iEvtN=0;iEvtN<pRich->GetLoader()->GetRunLoader()->GetNumberOfEvents();iEvtN++) {
- pRich->GetLoader()->GetRunLoader()->GetEvent(iEvtN);
- AliRICHTracker *tr = new AliRICHTracker();
- tr->RecWithStack(hn);
- AliInfoClass(Form("Pattern Recognition done for event %i \b",iEvtN));
- }
- pFile->Write();pFile->Close();
}
+//__________________________________________________________________________________________________
void FormCluster(AliRICHCluster *pClu,AliRICHDigit *pDig,TClonesArray *pDigList,TMatrixF *pDigMap)const;//form cluster recursive algorithm
inline AliRICHDigit *UseDig (Int_t padX,Int_t padY,TClonesArray *pDigList,TMatrixF *pDigMap )const;//use this pad's digit to form a cluster
static void CheckPR ( ); //utility-> run staff for stack
- static void RichAna (Int_t iNevMax=99999,Bool_t isPatRec=kFALSE ); //utility-> create ntuples for analysis
+ static void RichAna (Int_t iNevMin=0, Int_t iNevMax=99999,Bool_t isPatRec=kFALSE ); //utility-> create ntuples for analysis
protected:
ClassDef(AliRICHReconstructor, 0) //class for the RICH reconstruction
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) return;//intersection with no chamber found
-//find MIP cluster candidate (cluster which is closest to track intersection point)
+ if(!iChamber) {
+ pTrack->SetRICHsignal(-999); //to be improved by flags...
+ 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
Double_t chargeMip=0; //charge of the MIP
+ Bool_t kFound = kFALSE;
for(Int_t iClusN=0;iClusN<pRich->Clusters(iChamber)->GetEntries();iClusN++){//clusters loop for intersected chamber
AliRICHCluster *pClus=(AliRICHCluster*)pRich->Clusters(iChamber)->UncheckedAt(iClusN);//get pointer to current cluster
+ if(pClus->Q()<AliRICHParam::QthMIP()) continue;
Double_t distCurrent=pClus->DistTo(helix.PosPc());//distance between current cluster and helix intersection with PC
if(distCurrent<distMip){
+ kFound = kTRUE;
distMip=distCurrent;
iMipId=iClusN;
distX=pClus->DistX(helix.PosPc());
AliDebug(1,Form("Ploc (%f,%f,%f) dist= %f",helix.Ploc().Mag(),helix.Ploc().Theta()*TMath::RadToDeg(),
helix.Ploc().Phi()*TMath::RadToDeg(),pClus->DistTo(helix.PosPc())));
}//clusters loop for intersected chamber
-
+
+ if(!kFound) {
+ pTrack->SetRICHsignal(-999);
+ return;
+ }
AliDebug(1,Form("Min distance cluster: %i dist is %f",iMipId,distMip));
+ pTrack->SetRICHcluster(((Int_t)chargeMip)+1000000*iChamber);
+ pTrack->SetRICHdxdy(distX,distY);
+ pTrack->SetRICHthetaPhi(helix.Ploc().Theta(),helix.Ploc().Phi());
//
// HERE CUTS ON GOLD RINGS....
//
- if(distMip>AliRICHParam::DmatchMIP()||chargeMip<AliRICHParam::QthMIP()) {
+ if(distMip>AliRICHParam::DmatchMIP()) {
//track not accepted for pattern recognition
- pTrack->SetRICHsignal(-999.); //to be improved by flags...
+ pTrack->SetRICHsignal(-990); //to be improved by flags...
return;
}
//
Double_t thetaCerenkov=recon.ThetaCerenkov(); //search for mean Cerenkov angle for this track
- pTrack->SetRICHcluster(((Int_t)chargeMip)+1000000*iChamber);
- pTrack->SetRICHdxdy(distX,distY);
- pTrack->SetRICHthetaPhi(helix.Ploc().Theta(),helix.Ploc().Phi());
pTrack->SetRICHsignal(thetaCerenkov);
pTrack->SetRICHnclusters(recon.GetHoughPhotons());
if(recon.GetPhotonFlag() == 2) {
Double_t theta_g=recon.GetTrackTheta();
Double_t phi_g=(recon.GetPhiPoint()-recon.GetTrackPhi());
- Double_t sigma2 = AliRICHParam::SigmaSinglePhoton(iPart,pTrack->GetP(),theta_g,phi_g).Mag2();
- if (sigma2>0)
- sigmaPID[iPart] += 1/sigma2;
+ Double_t sigma2 = AliRICHParam::SigmaSinglePhoton(iPart,pTrack->GetP(),theta_g,phi_g).Mag2();
+ if(sigma2>0) sigmaPID[iPart] += 1/sigma2;
}
}
- sigmaPID[iPart] *= (Double_t)(recon.GetHoughPhotons()-fnPhotBKG)/(Double_t)(recon.GetHoughPhotons()); // n total phots, m are background...the sigma are scaled..
if (sigmaPID[iPart]>0)
- sigmaPID[iPart] = 1/TMath::Sqrt(sigmaPID[iPart])*0.001; // sigma from parametrization are in mrad...
- else
- sigmaPID[iPart] = 0;
- fErrPar[iPart]=sigmaPID[iPart];
+ sigmaPID[iPart] *= (Double_t)(recon.GetHoughPhotons()-fnPhotBKG)/(Double_t)(recon.GetHoughPhotons()); // n total phots, m are background...the sigma are scaled..
+ if(sigmaPID[iPart]>0) sigmaPID[iPart] = 1/TMath::Sqrt(sigmaPID[iPart])*0.001; // sigma from parametrization are in mrad...
+ else sigmaPID[iPart] = 0;
+ fErrPar[iPart]=sigmaPID[iPart];
AliDebug(1,Form("sigma for %s is %f rad",AliPID::ParticleName(iPart),sigmaPID[iPart]));
}
CalcProb(thetaCerenkov,pTrack->GetP(),sigmaPID,richPID);
// Double_t sinThetaThNorm = TMath::Sin(thetaTh)/TMath::Sqrt(1-1/(refIndex*refIndex));
// Double_t sigmaThetaTh = (0.014*(1/sinThetaThNorm-1) + 0.0043)*1.25;
// height[iPart] = TMath::Gaus(thetaCer,thetaTh,sigmaThetaTh);
- height[iPart] = TMath::Gaus(thetaCer,thetaTh[iPart],sigmaPID[iPart],kTRUE);
+ if(sigmaPID[iPart]>0) height[iPart] = TMath::Gaus(thetaCer,thetaTh[iPart],sigmaPID[iPart],kTRUE);
+ else height[iPart] = 0;
totalHeight +=height[iPart];
AliDebug(1,Form(" Particle %s with mass %f with height %f and thetaTH %f",AliPID::ParticleName(iPart),mass,height[iPart],thetaTh[iPart]));
AliDebug(1,Form(" partial height %15.14f total height %15.14f",height[iPart],totalHeight));
if(totalHeight<1e-5) {for(Int_t iPart=0;iPart<AliPID::kSPECIES;iPart++)richPID[iPart]=1.0/AliPID::kSPECIES;return;}
for(Int_t iPart=0;iPart<AliPID::kSPECIES;iPart++) richPID[iPart] = height[iPart]/totalHeight;
Int_t iPartNear = TMath::LocMax(AliPID::kSPECIES,richPID);
- if(TMath::Abs(thetaCer-thetaTh[iPartNear])
> 5*sigmaPID[iPartNear]) for(Int_t iPart=0;iPart<AliPID::kSPECIES;iPart++)richPID[iPart]=1.0/AliPID::kSPECIES;
+ if(TMath::Abs(thetaCer-thetaTh[iPartNear])
>5.*sigmaPID[iPartNear]) for(Int_t iPart=0;iPart<AliPID::kSPECIES;iPart++)richPID[iPart]=1.0/AliPID::kSPECIES;
//last line is to check if the nearest thetacerenkov to the teorethical one is within 5 sigma, otherwise no response (equal prob to every particle
}//CalcProb
if(AliceRead()){//it's from file, show some info
pMenu->AddButton("Display single chambers" ,"r->Display();" , "Display Fast");
pMenu->AddButton("Display ALL chambers" ,"r->DisplayEvent(0,0);" , "Display Fast");
- pMenu->AddButton("Hits QA" ,"hqa()" ,"????");
+ pMenu->AddButton("Hits QA" ,"hqa()" ,"QA plots for hits: hqa()");
pMenu->AddButton("Recon with stack" ,"AliRICHReconstructor::CheckPR( )","Create RSR.root with ntuple hn");
- pMenu->AddButton("RichAna no Recon" ,"AliRICHReconstructor::RichAna(0,kFALSE)","Create RichAna.root with ntuple hn without PatRec");
- pMenu->AddButton("RichAna with Recon" ,"AliRICHReconstructor::RichAna(0,kTRUE )","Create RichAna.root with ntuple hn with PatRec");
- pMenu->AddButton("Print hits" ,"h();" ,"????");
- pMenu->AddButton("Print sdigits" ,"s();" ,"????");
- pMenu->AddButton("Print digits" ,"d();" ,"????");
- pMenu->AddButton("Print clusters" ,"c();" ,"????");
- pMenu->AddButton("Print occupancy" ,"r->OccupancyPrint(-1);" ,"????");
- pMenu->AddButton("Print event summary " ,"r->SummaryOfEvent();" ,"????");
+ pMenu->AddButton("RichAna no Recon" ,"AliRICHReconstructor::RichAna(0,0,kFALSE)","Create RichAna.root with ntuple hn without PatRec");
+ pMenu->AddButton("RichAna with Recon" ,"AliRICHReconstructor::RichAna(0,0,kTRUE )","Create RichAna.root with ntuple hn with PatRec");
+ pMenu->AddButton("Print hits" ,"h();" ,"To print hits: h()");
+ pMenu->AddButton("Print sdigits" ,"s();" ,"To print sdigits: s()");
+ pMenu->AddButton("Print digits" ,"d();" ,"To print digits: d()");
+ pMenu->AddButton("Print clusters" ,"c();" ,"To print clusters: c()");
+ pMenu->AddButton("Print occupancy" ,"r->OccupancyPrint(-1);" ,"To print occupancy");
+ pMenu->AddButton("Print event summary " ,"r->SummaryOfEvent();" ,"To print a summary of the event");
}else{//it's aliroot, simulate
pMenu->AddButton("Debug ON", "DebugON();", "Switch debug on-off");
pMenu->AddButton("Debug OFF", "DebugOFF();", "Switch debug on-off");
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff