From 84093f6f6bcade15879eaea0ab22cde192a407f2 Mon Sep 17 00:00:00 2001 From: kir Date: Fri, 19 May 2006 15:33:23 +0000 Subject: [PATCH] Prob vector is assigned in FillESD(), some clarification on calculation of heights --- RICH/AliRICHParam.cxx | 4 - RICH/AliRICHParam.h | 3 +- RICH/AliRICHRecon.cxx | 47 ++++++++--- RICH/AliRICHRecon.h | 6 +- RICH/AliRICHReconstructor.cxx | 44 +++++++++- RICH/AliRICHReconstructor.h | 14 +++- RICH/AliRICHTracker.cxx | 146 ++++++++++++++++------------------ RICH/AliRICHTracker.h | 4 +- RICH/RichConfig.C | 2 - 9 files changed, 166 insertions(+), 104 deletions(-) diff --git a/RICH/AliRICHParam.cxx b/RICH/AliRICHParam.cxx index c6f71e84572..e9f2a2cb957 100644 --- a/RICH/AliRICHParam.cxx +++ b/RICH/AliRICHParam.cxx @@ -36,10 +36,6 @@ ClassImp(AliRICHParam) AliRICHParam * AliRICHParam::fgInstance =0x0; //singleton pointer - -Double_t AliRICHParam::fgMass[5] ={0.00051,0.10566,0.13957,0.49360,0.93828}; - - //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ AliRICHParam::AliRICHParam():TNamed("RichParam","default version") { diff --git a/RICH/AliRICHParam.h b/RICH/AliRICHParam.h index 0edeada907f..33e50382cd4 100644 --- a/RICH/AliRICHParam.h +++ b/RICH/AliRICHParam.h @@ -6,6 +6,7 @@ #include //old style #include #include +#include //---------to be deleted------------- #include #include //Hit2SDigs() #include @@ -159,7 +160,6 @@ public: static void TestSeg(); //test the segmentation group of methodes static void TestTrans(); //test the transform group of methodes - static Double_t fgMass[5]; // mass array enum EPlaneId {kCenter,kPc,kRad,kAnod,kNch=7}; //4 planes in chamber and total number of chambers protected: AliRICHParam(); //default ctor is protected to enforce it to be singleton @@ -659,6 +659,7 @@ Double_t AliRICHParam::ErrCrom(Double_t thetaC, Double_t phiC, Double_t thetaM, Double_t refC6F14m = 1.29337; Double_t alpha =TMath::Cos(thetaM)-TMath::Tan(thetaC)*TMath::Cos(phiDelta)*TMath::Sin(thetaM); + //cout << "alpha : "<GetEntriesFast();iClu++){//clusters loop if(iClu == iMipId) continue; // do not consider MIP cluster as a photon candidate SetPhotonIndex(iClu); @@ -69,11 +73,10 @@ Double_t AliRICHRecon::ThetaCerenkov(AliRICHHelix *pHelix,TClonesArray *pCluster SetPhotonEta(-999.); SetPhotonWeight(0.); AliRICHCluster *pClu=(AliRICHCluster*)fClusters->UncheckedAt(iClu); //get pointer to current cluster -// if(pClu->Q()>AliRICHParam::QthMIP()) continue; //avoid MIP clusters from bkg + if(pClu->Q()>AliRICHParam::QthMIP()) continue; //avoid MIP clusters from bkg SetEntranceX(pClu->X() - GetShiftX()); SetEntranceY(pClu->Y() - GetShiftY()); //cluster position with respect to track intersection FindPhiPoint(); -// if(PhotonInBand()==0) continue; ???????????? - SetPhotonFlag(1); + SetPhotonFlag(1); FindThetaPhotonCerenkov(); Float_t thetaPhotonCerenkov = GetThetaPhotonCerenkov(); AliDebug(1,Form("Track Theta=%5.2f deg, Phi=%5.2f deg Photon clus=%2i ThetaCkov=%5.2f rad",GetTrackTheta()*TMath::RadToDeg(),GetTrackPhi()*TMath::RadToDeg() @@ -590,6 +593,7 @@ Float_t AliRICHRecon::FromEmissionToCathode() Float_t nfreon, nquartz, ngas; + //fParam->Print(); nfreon = fParam->IdxC6F14(fParam->EckovMean()); nquartz = fParam->IdxSiO2(fParam->EckovMean()); @@ -762,24 +766,41 @@ Double_t AliRICHRecon::HoughResponse() //__________________________________________________________________________________________________ void AliRICHRecon::FindThetaCerenkov() { - // manage with weight for photons +// Loops on all Ckov candidates and estimates the best Theta Ckov for a ring formed by those candidates. Also estimates an error for that Theat Ckov +// collecting errors for all single Ckov candidates thetas. (Assuming they are independent) +// Arguments: none +// Return: none Float_t wei = 0.; Float_t weightThetaCerenkov = 0.; Double_t etaMin=9999.,etaMax=0.; + Double_t sigma2 = 0; //to collect error squared for this ring + for(Int_t i=0;ietaMax) etaMax=photonEta; - Float_t photonWeight = GetPhotonWeight(); - weightThetaCerenkov += photonEta*photonWeight; - wei += photonWeight; - } + Float_t photonWeight = GetPhotonWeight(); + weightThetaCerenkov += photonEta*photonWeight; + wei += photonWeight; + //here comes sigma of the reconstructed ring + + //Double_t phiref=(GetPhiPoint()-GetTrackPhi()); + if(GetPhotonEta()<=0) continue;//?????????????????Flag photos = 2 may imply CkovEta = 0?????????????? + //??????????? Look at SetPhoton Flag method + Double_t phiref=GetTrackPhi(); + + Double_t beta = 1./(TMath::Cos(GetPhotonEta())*fParam->IdxC6F14(AliRICHParam::EckovMean())); + sigma2 += 1./AliRICHParam::SigmaSinglePhotonFormula(GetPhotonEta(),GetPhiPoint(),GetTrackTheta(),phiref,beta); + } } - + + if(sigma2>0) SetRingSigma2(1./sigma2); + else SetRingSigma2(1e10); + if(wei != 0.) weightThetaCerenkov /= wei; else weightThetaCerenkov = 0.; SetThetaCerenkov(weightThetaCerenkov); @@ -791,7 +812,6 @@ void AliRICHRecon::FindThetaCerenkov() Double_t nPhotBKG = (externalArea-internalArea)/effArea*fClusters->GetEntries(); if(nPhotBKG<0) nPhotBKG=0; //just protection from funny angles... SetPhotBKG(nPhotBKG); - // AliDebug(1,Form(" thetac weighted -> %f",weightThetaCerenkov)); } @@ -814,7 +834,10 @@ Int_t AliRICHRecon::FlagPhotons(Double_t thetaCkovHough) for(Int_t i=0;i= tmin && photonEta <= tmax) { SetPhotonFlag(2); iInsideCnt++;} + if(photonEta >= tmin && photonEta <= tmax) { + SetPhotonFlag(2); + iInsideCnt++; + } } return iInsideCnt; }//FlagPhotons diff --git a/RICH/AliRICHRecon.h b/RICH/AliRICHRecon.h index 50048cbfebe..484539bb036 100644 --- a/RICH/AliRICHRecon.h +++ b/RICH/AliRICHRecon.h @@ -84,7 +84,8 @@ public : Float_t GetFittedTrackTheta() const{ return fFittedTrackTheta;} // Float_t GetFittedTrackPhi() const{ return fFittedTrackPhi;} // Float_t GetFittedThetaCerenkov() const{ return fFittedThetaCerenkov;} // - void SetPhotonEnergy(Float_t e) { fEphot = e;} // + Float_t GetRingSigma2() const{ return fRingSigma2;} // + void SetPhotonEnergy(Float_t e) { fEphot = e;} // void SetEmissionPoint(Float_t LengthEmissionPoint) { fLengthEmissionPoint = LengthEmissionPoint;} // void SetEntranceX(Float_t Xtoentr) { fXtoentr = Xtoentr;} // void SetEntranceY(Float_t Ytoentr) { fYtoentr = Ytoentr;} // @@ -129,6 +130,7 @@ public : void SetFittedTrackPhi(Float_t FittedTrackPhi) { fFittedTrackPhi = FittedTrackPhi;} // void SetFittedThetaCerenkov(Float_t FittedThetaCerenkov) { fFittedThetaCerenkov = FittedThetaCerenkov;}// void SetFittedHoughPhotons(Int_t FittedHoughPhotons) { fFittedHoughPhotons = FittedHoughPhotons;} // + void SetRingSigma2(Float_t RingSigma2) { fRingSigma2 = RingSigma2;} // Float_t SnellAngle(Float_t n1, Float_t n2, Float_t theta1); // Float_t FromEmissionToCathode(); // @@ -201,7 +203,7 @@ protected: Float_t fThetaCerenkov; // Theta angle for Hough Float_t fWeightThetaCerenkov; // Theta Cerenkov angle weighted Float_t fThetaPeakPos; // Peak position - + Float_t fRingSigma2; // sigma2 of the reconstructed ring ClassDef(AliRICHRecon,0) }; diff --git a/RICH/AliRICHReconstructor.cxx b/RICH/AliRICHReconstructor.cxx index 0ae92b9680d..c7c70f5cca2 100644 --- a/RICH/AliRICHReconstructor.cxx +++ b/RICH/AliRICHReconstructor.cxx @@ -275,7 +275,7 @@ void AliRICHReconstructor::RichAna(Int_t iNevMin,Int_t iNevMax,Bool_t askPatRec) hnvec[12]=prob[0]+prob[1]+prob[2]; hnvec[13]=prob[3]; hnvec[14]=prob[4]; hnvec[15]=pTrRich->fErrPar[2]; hnvec[16]=pTrRich->fErrPar[3]; hnvec[17]=pTrRich->fErrPar[4]; for(Int_t i=0;i<3;i++) { - Double_t mass = AliRICHParam::fgMass[i+2]; + Double_t mass = AliPID::ParticleMass(i+2); Double_t refIndex=AliRICHParam::Instance()->IdxC6F14(AliRICHParam::EckovMean()); Double_t cosThetaTh = TMath::Sqrt(mass*mass+pTrack->GetP()*pTrack->GetP())/(refIndex*pTrack->GetP()); hnvec[18+i]=0; @@ -351,3 +351,45 @@ void AliRICHReconstructor::Test(TClonesArray *pDigLst,Bool_t isTryUnfold) delete pCluLst; }//Test() //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +void AliRICHReconstructor::FillESD(AliRunLoader *, AliESD *pESD) const +{ +// Calculates probability to be a electron-muon-pion-kaon-proton +// from the given Cerenkov angle and momentum assuming no initial particle composition +// (i.e. apriory probability to be the particle of the given sort is the same for all sorts) + + AliPID ppp; //needed + Double_t pid[AliPID::kSPECIES],h[AliPID::kSPECIES]; + Double_t refIndex=AliRICHParam::Instance()->IdxC6F14(AliRICHParam::EckovMean()); + + for(Int_t iTrk=0;iTrkGetNumberOfTracks();iTrk++){//ESD tracks loop + AliESDtrack *pTrack = pESD->GetTrack(iTrk);// get next reconstructed track + if(pTrack->GetRICHsignal()<=0){//RICH does not find anything reasonable for this track, assign 0.2 for all species + for(Int_t iPart=0;iPartSetRICHpid(pid); + continue; + } + Double_t pmod = pTrack->GetP(); + Double_t hTot=0; + for(Int_t iPart=0;iPartGetRICHsignal(),TMath::Sqrt(pTrack->GetRICHchi2()),kTRUE); + + else //beta < 1/ref. idx. => no light at all + h[iPart] =0 ; + hTot +=h[iPart]; //total height of all theoretical heights for normalization + }//species loop + + Double_t hMin=TMath::Gaus(5*TMath::Sqrt(pTrack->GetRICHchi2()),pTrack->GetRICHsignal(),TMath::Sqrt(pTrack->GetRICHchi2()),kTRUE);//5 sigma protection + + for(Int_t iPart=0;iParthMin) + pid[iPart]=h[iPart]/hTot; + else //all theoretical values are far away from experemental one + pid[iPart]=1.0/AliPID::kSPECIES; + pTrack->SetRICHpid(pid); + }//ESD tracks loop + //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 +}//FillESD diff --git a/RICH/AliRICHReconstructor.h b/RICH/AliRICHReconstructor.h index c790b4342e8..40029879b24 100644 --- a/RICH/AliRICHReconstructor.h +++ b/RICH/AliRICHReconstructor.h @@ -20,8 +20,15 @@ public: AliTracker* CreateTracker (AliRunLoader* )const{return new AliRICHTracker;} //from AliReconstructor for clusters->PID void Reconstruct (AliRunLoader* pAL )const; //from AliReconstruction for digits->clusters void Reconstruct (AliRunLoader* pAL,AliRawReader *pRR)const; //from AliReconstruction for raws->clusters - using AliReconstructor::Reconstruct; //to get rid of virtual hidden warning -//private part + virtual void FillESD (AliRunLoader* pAL,AliESD *pESD)const; //calculate pid for RICH + virtual void FillESD(AliRunLoader*, AliRawReader*, AliESD*) const { }; + virtual void FillESD(AliRawReader*, TTree*, AliESD*) const { }; + virtual void FillESD(TTree*, TTree*, AliESD*) const { }; + + + using AliReconstructor::Reconstruct; //to get rid of virtual hidden warning + + //private part static void Dig2Clu (TClonesArray*pDigLst,TClonesArray *pCluLst,Bool_t isTryUnfold=kTRUE );//digits list -> clusters list static void CluQA (AliRunLoader* pAL );//QA for clusters static void CheckPR ( );//utility-> run staff for stack @@ -30,7 +37,8 @@ public: static void RichAna (Int_t iNevMin=0,Int_t iNevMax=99999,Bool_t isPatRec=kFALSE );//utility-> create ntuples for analysis static void Test (Bool_t isTryUnfold=kTRUE );//test digits->clusters conversion static void Test (TClonesArray *pDigLst,Bool_t isTryUnfold=kTRUE );//test digits->clusters conversion -protected: + + protected: ClassDef(AliRICHReconstructor, 0) //class for the RICH reconstruction }; //__________________________________________________________________________________________________ diff --git a/RICH/AliRICHTracker.cxx b/RICH/AliRICHTracker.cxx index d675541d007..28526903d71 100644 --- a/RICH/AliRICHTracker.cxx +++ b/RICH/AliRICHTracker.cxx @@ -1,10 +1,13 @@ #include "AliRICHTracker.h" //class header #include "AliRICH.h" #include "AliRICHRecon.h" +#include "AliRICHParam.h" #include +#include //EsdQA() #include -#include //EsdPrint() -#include //EsdPrint() +#include //EsdQA() +#include //EsdQA() +#include //EsdQA() #include "AliRICHHelix.h" #include #include @@ -14,6 +17,9 @@ #include //RecWithStack(); #include //GetTrackPoint() #include //GetTrackPoint() +#include //EsdQA() +#include //EsdQA() +#include //EsdQA() ClassImp(AliRICHTracker) //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ AliRICHTracker::AliRICHTracker():AliTracker() @@ -92,41 +98,14 @@ Int_t AliRICHTracker::PropagateBack(AliESD *pESD) pTrack->SetRICHcluster(iMipId+1000000*iChamber); //set mip cluster index pTrack->SetRICHsignal(recon.ThetaCerenkov(&helix,pRich->Clus(iChamber),iMipId));//search for mean Cerenkov angle for this track + pTrack->SetRICHchi2(recon.GetRingSigma2()); pTrack->SetRICHnclusters(iMipId); //on return iMipId is number of photon clusters accepted in reconstruction AliDebug(1,Form("Ch=%i PC Intersection=(%5.2f,%5.2f) cm MIP cluster dist=(%5.2f,%5.2f)=%5.2f cm ThetaCkov=%f", iChamber,helix.PosPc().X(),helix.PosPc().Y(), mipDx,mipDy,mipDr, pTrack->GetRICHsignal())); //here comes PID calculations - if(pTrack->GetRICHsignal()>0) { - AliDebug(1,Form("Start to assign the probabilities")); - Double_t sigmaPID[AliPID::kSPECIES]; Double_t richPID[AliPID::kSPECIES]; - for (Int_t iPid=0;iPidClus(iChamber)->GetEntries();iCkov++){//Ckov candidates loop ????????????? somehting fomr AliRICHRecon must be - recon.SetPhotonIndex(iCkov); - if(recon.GetPhotonFlag() == 2) { - Double_t thetaCkov=0.6; //?????????????????? - Double_t phiCkov=0.6; //?????????????????? - Double_t thetaTrk=recon.GetTrackTheta(); - Double_t phiTrk=(recon.GetPhiPoint()-recon.GetTrackPhi()); - Double_t beta =pTrack->GetP()/TMath::Sqrt((pTrack->GetP()*pTrack->GetP()+pid.ParticleMass(iPid)*pid.ParticleMass(iPid))); - Double_t sigma2 = AliRICHParam::SigmaSinglePhotonFormula(thetaCkov,phiCkov,thetaTrk,phiTrk,beta); - if(sigma2>0) sigmaPID[iPid] += 1/sigma2; - } - }//Ckov candidates loop - - if (sigmaPID[iPid]>0) - sigmaPID[iPid] *= (Double_t)(iMipId-recon.GetPhotBKG())/(Double_t)(iMipId); // n total phots, m are background...the sigma are scaled.. - if(sigmaPID[iPid]>0) sigmaPID[iPid] = 1/TMath::Sqrt(sigmaPID[iPid])*0.001; // sigma from parametrization are in mrad... - else sigmaPID[iPid] = 0; - fErrPar[iPid]=sigmaPID[iPid]; - AliDebug(1,Form("sigma for %s is %f rad",AliPID::ParticleName(iPid),sigmaPID[iPid])); - } - CalcProb(pTrack->GetRICHsignal(),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])); - }//if(pTrack->GetRICHsignal()) +// CalcProb(pTrack->GetRICHsignal(),pTrack->GetP(),sigmaPID,richPID); }//ESD tracks loop AliDebug(1,"Stop pattern recognition"); return 0; // error code: 0=no error; @@ -219,42 +198,10 @@ void AliRICHTracker::RecWithStack(TNtupleD *hn) AliDebug(1,"Stop."); }//RecWithStack //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -void AliRICHTracker::CalcProb(Double_t thetaCer,Double_t pmod, Double_t *sigmaPID, Double_t *richPID) +void AliRICHTracker::EsdQA(Bool_t isPrint) { -// Calculates probability to be a electron-muon-pion-kaon-proton -// from the given Cerenkov angle and momentum assuming no initial particle composition -// (i.e. apriory probability to be the particle of the given sort is the same for all sorts) - Double_t height[AliPID::kSPECIES];Double_t totalHeight=0; - Double_t thetaTh[AliPID::kSPECIES]; - for(Int_t iPart=0;iPartIdxC6F14(AliRICHParam::EckovMean()); - Double_t cosThetaTh = TMath::Sqrt(mass*mass+pmod*pmod)/(refIndex*pmod); - thetaTh[iPart]=0; - if(cosThetaTh>=1) continue; - thetaTh[iPart] = TMath::ACos(cosThetaTh); -// 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); - if(sigmaPID[iPart]>0) height[iPart] = TMath::Gaus(thetaCer,thetaTh[iPart],sigmaPID[iPart],kTRUE); - else height[iPart] = 0; - totalHeight +=height[iPart]; - AliDebugClass(1,Form(" Particle %s with mass %f with height %f and thetaTH %f",AliPID::ParticleName(iPart),mass,height[iPart],thetaTh[iPart])); - AliDebugClass(1,Form(" partial height %15.14f total height %15.14f",height[iPart],totalHeight)); - } - if(totalHeight<1e-5) {for(Int_t iPart=0;iPart5.*sigmaPID[iPartNear]) for(Int_t iPart=0;iPartSetBranchAddress("ESD", &pEsd); Int_t iNevt=pTr->GetEntries(); Printf("This ESD contains %i events",iNevt); + //AliRICHParam *pPar; + + TH1D *pDx,*pDy,*pProbEl,*pProbMu,*pProbPi,*pProbKa,*pProbPr; TH2F *pThP; TProfile *pChiTh; + if(!isPrint){ + TH1::AddDirectory(kFALSE); + pDx =new TH1D("dX" ,"distance between Xmip and Xtrack;cm",300,-1.5,1.5); + pDy =new TH1D("dY" ,"distance between Ymip and Ytrack;cm",300,-1.5,1.5); + pThP =new TH2F("tvsp","#theta_{Ckov} radian;P GeV" ,65 ,-0.5,6.0,75,0,0.75); pThP->SetStats(0); + pProbPi=new TH1D("RichProbPion","HMPID PID probability for e #mu #pi" ,101,0.05,1.05); pProbPi->SetLineColor(kRed); + pProbEl=new TH1D("RichProbEle" ,"" ,101,0.05,1.05); pProbEl->SetLineColor(kGreen); + pProbMu=new TH1D("RichProbMuon" ,"" ,101,0.05,1.05); pProbMu->SetLineColor(kBlue); + pProbKa=new TH1D("RichProbKaon","HMPID PID probability for K" ,101,0.05,1.05); + pProbPr=new TH1D("RichProbProton","HMPID PID probability for p" ,101,0.05,1.05); + pChiTh =new TProfile("RichChiTh","#chi^{2};#theta_{C}" ,80 ,0,0.8 , -2,2); + + //if(!gGeoManager)TGeoManager::Import("geometry.root"); + //pPar = AliRICHParam::Instance(); + } + for(Int_t iEvt=0;iEvtGetEvent(iEvt); Int_t iNtracks=pEsd->GetNumberOfTracks(); Printf("ESD contains %i tracks created in Bz=%.2f Tesla",iNtracks,pEsd->GetMagneticField()/10.); @@ -270,20 +236,46 @@ void AliRICHTracker::EsdPrint() Float_t dx,dy; pTrack->GetRICHdxdy(dx,dy); Float_t theta,phi; pTrack->GetRICHthetaPhi(theta,phi); TString comment; - if(pTrack->GetRICHsignal()>0) - comment="OK"; - else if(pTrack->GetRICHsignal()==kMipQdcCut) - comment="no enough QDC"; - else if(pTrack->GetRICHsignal()==kMipDistCut) - comment="nearest cluster is too far"; - else if(pTrack->GetRICHsignal()==-1) - comment="no intersection"; - Printf("Track %2i Q=%4.1f P=%.3f GeV RICH: ChamClus %7i Track-Mip=(%7.2f,%7.2f)=%5.2f cm ThetaCer %7.1f rad %s", - iTrk,pTrack->GetSign(),pTrack->GetP(), - pTrack->GetRICHcluster(),dx,dy,TMath::Sqrt(dx*dx+dy*dy),pTrack->GetRICHsignal(), comment.Data()); + if(isPrint){ + if(pTrack->GetRICHsignal()>0) comment="OK"; + else if(pTrack->GetRICHsignal()==kMipQdcCut) comment="no enough QDC"; + else if(pTrack->GetRICHsignal()==kMipDistCut) comment="nearest cluster is too far"; + else if(pTrack->GetRICHsignal()==-1) comment="no intersection"; + Double_t pid[5]; pTrack->GetRICHpid(pid); + Printf("Tr %2i Q=%4.1f P=%.3f GeV Tr-Mip=%5.2f cm Th=%7.3f rad Prob : e=%.4f mu=%.4f pi=%.4f K=%.4f p=%.4f %s" , + iTrk,pTrack->GetSign(),pTrack->GetP(),TMath::Sqrt(dx*dx+dy*dy),pTrack->GetRICHsignal(), + pid[0],pid[1],pid[2],pid[3],pid[4], comment.Data()); + }else{//collect hists + + pDx->Fill(dx); pDy->Fill(dy); + pThP->Fill(pTrack->GetP(),pTrack->GetRICHsignal()); + pChiTh->Fill(pTrack->GetRICHsignal(),pTrack->GetRICHchi2()); + Double_t prob[5]; pTrack->GetRICHpid(prob); + pProbEl->Fill(prob[0]); pProbMu->Fill(prob[1]); pProbPi->Fill(prob[2]); pProbKa->Fill(prob[3]); pProbPr->Fill(prob[4]); + }//if plot }//ESD tracks loop }//ESD events loop delete pEsd; pFile->Close();//close AliESDs.root + if(!isPrint){ + TCanvas *pC=new TCanvas("c","Quality",1200,1500); pC->Divide(2,2); + TF1 *pPion = new TF1("RICHtheor","acos(sqrt(x*x+[0]*[0])/(x*[1]))",1.2,6); pPion->SetLineWidth(1);pPion->SetParameter(1,1.292); + AliPID ppp; pPion->SetParameter(0,AliPID::ParticleMass(AliPID::kPion)) ; pPion->SetLineColor(kRed); + TF1 *pKaon = (TF1*)pPion->Clone(); pKaon->SetParameter(0,AliPID::ParticleMass(AliPID::kKaon)) ; pKaon->SetLineColor(kGreen); + TF1 *pProt = (TF1*)pPion->Clone(); pProt->SetParameter(0,AliPID::ParticleMass(AliPID::kProton)); pProt->SetLineColor(kBlue); + + + pC->cd(1); pDx->Draw(); + pC->cd(2); pDy->Draw(); + pC->cd(3); pThP->Draw(); pPion->Draw("same"); pKaon->Draw("same"); pProt->Draw("same"); + pC->cd(4); pChiTh->Draw(); + + TCanvas *pC2=new TCanvas("c2","Quality 2",1200,1500); pC2->Divide(2,2); + pC2->cd(1); pProbPi->Draw(); pProbMu->Draw("same"); pProbEl->Draw("same"); + pC2->cd(2); pProbKa->Draw(); + pC2->cd(3); pProbPr->Draw(); + + + } } //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ void AliRICHTracker::MatrixPrint(Double_t probCut) diff --git a/RICH/AliRICHTracker.h b/RICH/AliRICHTracker.h index 9e4cf607bf0..efef96fb407 100644 --- a/RICH/AliRICHTracker.h +++ b/RICH/AliRICHTracker.h @@ -17,12 +17,12 @@ public: Int_t Clusters2Tracks(AliESD * ) {return 0;} //pure virtual from AliTracker Int_t LoadClusters (TTree *pCluTr ); //pure virtual from AliTracker Int_t PropagateBack (AliESD * ); //pure virtual from AliTracker invoked from AliReconstruction::RunTracking() + //void FillESD (AliESD *pESD ); //calculate pid for RICH Int_t RefitInward (AliESD * ) {return 0;} //pure virtual from AliTracker void UnloadClusters ( ) { } //pure virtual from AliTracker //private part void RecWithStack(TNtupleD *hn ); //recon from Stack in case ESD empty - static void CalcProb (Double_t thetaCer,Double_t pmod,Double_t *pidsigma, Double_t *pid); //calculate pid for RICH - static void EsdPrint ( ); //print ESD status + static void EsdQA (Bool_t isPrint=kFALSE ); //print QA info static void MatrixPrint (Double_t probCut=0.7 ); //print prob matrix with cut on probability Double_t fErrPar[5]; //Temporary stored for debug purpose enum ETrackingFlags {kMipDistCut=-990,kMipQdcCut=-999}; diff --git a/RICH/RichConfig.C b/RICH/RichConfig.C index 407bee42cc0..63c01af1c55 100644 --- a/RICH/RichConfig.C +++ b/RICH/RichConfig.C @@ -527,7 +527,6 @@ void RichConfig::WriteBatch() //reconstraction section - if(!fClusBG->GetButton(kNo)->GetState()){ fprintf(fp," AliReconstruction *pRec=new AliReconstruction;\n"); if (fClusBG->GetButton(kAll) ->GetState()) fprintf(fp," pRec->SetRunLocalReconstruction(\"ITS TPC TRD TOF RICH\"); //clusters created for these detectors\n"); else if(fClusBG->GetButton(kOnly) ->GetState()) fprintf(fp," pRec->SetRunLocalReconstruction(\"RICH\"); //clusters created for RICH only\n"); @@ -549,7 +548,6 @@ void RichConfig::WriteBatch() fprintf(fp," pRec->Run();delete pRec;\n\n"); - } //benchmarks fprintf(fp," cout<<\"!!!!!!!!!!!!Info in : Start time: \";time.Print();\n"); fprintf(fp," cout<<\"!!!!!!!!!!!!Info in : Stop time: \";time.Set(); time.Print();\n"); -- 2.43.0