RICH/AliRICHTracker.cxx

   1 #include "AliRICHTracker.h" //class header
   2 #include "AliRICH.h"
   3 #include "AliRICHRecon.h"
   4 #include "AliRICHParam.h"
   5 #include <AliESD.h>
   6 #include <TGeoManager.h>     //EsdQA()
   7 #include <TVector3.h>
   8 #include <TTree.h>          //EsdQA()
   9 #include <TFile.h>          //EsdQA()
  10 #include <TProfile.h>   //EsdQA()
  11 #include "AliRICHHelix.h"
  12 #include <AliMagF.h>
  13 #include <AliStack.h>
  14 #include <TParticle.h>
  15 #include <TMath.h>
  16 #include <AliRun.h>
  17 #include <TNtupleD.h>            //RecWithStack();
  18 #include <AliTrackPointArray.h>  //GetTrackPoint()
  19 #include <AliAlignObj.h>         //GetTrackPoint()
  20 #include <TH1F.h>                //EsdQA()
  21 #include <TH2F.h>                //EsdQA()
  22 #include <TCanvas.h>             //EsdQA()
  23 ClassImp(AliRICHTracker)
  24 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  25 AliRICHTracker::AliRICHTracker():AliTracker()
  26 {
  27 // AliRICHTracker is created from AliReconstraction::Run() which invokes AliReconstraction::CreateTrackers()
  28 // which in turn invokes AliRICHReconstructor::CreateTracker().
  29 // Note that this is done just once per session before AliReconstruction::Run() goes to events loop.
  30   AliRICHParam::Instance()->CdbRead(0,0);
  31   for(Int_t i=0;i<5;i++)fErrPar[i]=0;
  32 }
  33 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  34 Bool_t AliRICHTracker::GetTrackPoint(Int_t idx, AliTrackPoint& point) const
  35 {
  36 // Interface callback methode invoked from AliReconstruction::WriteAlignmentData() to get position of MIP cluster in MARS associated to a current track.
  37 // MIP cluster is reffered by index which is stored in AliESDtrack  ???????
  38 // Arguments: idx- cluster index which is stored by RICH in AliESDtrack
  39 //            point- reference to the object where to store the point
  40 //   Returns: status of operation  if FALSE then AliReconstruction::WriteAlignmentData() do not store this point to array of points for current track.
  41   if(idx<0) return kFALSE; //no MIP cluster assigned to this track in PropagateBack()
  42   Int_t iCham=idx/1000000;
  43   Int_t iClu=idx%1000000;
  44   point.SetVolumeID(AliAlignObj::LayerToVolUID(AliAlignObj::kRICH,iCham-1));//layer and chamber number
  45   AliRICH *pRich=((AliRICH*)gAlice->GetDetector("RICH"));
  46   AliRICHCluster *pClu=(AliRICHCluster*)pRich->Clus(iCham)->UncheckedAt(iClu);//get pointer to cluster
  47   TVector3 mars=AliRICHParam::Instance()->Lors2Mars(iCham,pClu->X(),pClu->Y());
  48   point.SetXYZ(mars.X(),mars.Y(),mars.Z());
  49   return kTRUE;
  50 }
  51 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  52 Int_t AliRICHTracker::LoadClusters(TTree *pCluTree)
  53 {
  54 // Interface callback methode invoked from AliReconstruction::RunTracking() to load RICH clusters for RICH
  55 // Arguments: pCluTree- pointer to clusters tree got by AliRICHLoader::LoadRecPoints("read") then AliRICHLoader::TreeR()
  56 //   Returns: error code (currently ignored in AliReconstruction::RunTraking())
  57   AliDebug(1,"Start.");  pCluTree->GetEntry(0);  AliDebug(1,"Stop."); return 0;
  58 }
  59 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  60 Int_t AliRICHTracker::PropagateBack(AliESD *pESD)
  61 {
  62 // Interface callback methode invoked by AliRecontruction::RunTracking() during tracking after TOF. It's done just once per event
  63 // Arguments: pESD - pointer to Event Summary Data class instance which contains a list of tracks
  64 //   Returns: error code, 0 if no errors
  65   Int_t iNtracks=pESD->GetNumberOfTracks();  AliDebug(1,Form("Start with %i tracks",iNtracks));
  66   AliRICH *pRich=((AliRICH*)gAlice->GetDetector("RICH"));
  67   AliRICHRecon recon;
  68   Int_t nphots =0;
  69
  70   for(Int_t iTrk=0;iTrk<iNtracks;iTrk++){//ESD tracks loop
  71     AliESDtrack *pTrk = pESD->GetTrack(iTrk);// get next reconstructed track
  72     Double_t mom[3], pos[3];
  73     pTrk->GetPxPyPz(mom); TVector3 mom3(mom[0],mom[1],mom[2]);
  74     pTrk->GetXYZ(pos);    TVector3 pos3(pos[0],pos[1],pos[2]);
  75     AliRICHHelix helix(pos3,mom3,(Int_t)pTrk->GetSign(),-0.1*GetBz()); //construct helix out of track running parameters
  76      //Printf(" magnetic field %f charged %f\n",GetBz(),pTrack->GetSign()); helix.Print("Track");
  77     Int_t iChamber=helix.RichIntersect(AliRICHParam::Instance());
  78     if(!iChamber) continue;                                         //no intersection with chambers, ignore this track go after the next one
  79
  80     //find MIP cluster candidate (closest to track intersection point cluster with large enough QDC)
  81     Double_t    dR=9999,   dX=9999,   dY=9999;                                             //distance between track-PC intersection point and current cluster
  82     Double_t mipDr=9999,mipDx=9999,mipDy=9999,mipX=9999,mipY=9999; Int_t mipQ=0;           //nearest cluster parameters
  83     Int_t   iMipId=-1;                                                                     //index of this nearest cluster
  84     for(Int_t iClu=0;iClu<pRich->Clus(iChamber)->GetEntries();iClu++){                     //clusters loop for intersected chamber
  85       AliRICHCluster *pClu=(AliRICHCluster*)pRich->Clus(iChamber)->UncheckedAt(iClu);      //get pointer to current cluster
  86       if(pClu->Q()<AliRICHParam::QthMIP()) continue;                                       //to low QDC, go after another one
  87       pClu->DistXY(helix.PosPc(),dX,dY); dR=TMath::Sqrt(dX*dX+dY*dY);                      //get distance for current cluster
  88       if(dR<mipDr){iMipId=iClu; mipDr=dR; mipDx=dX; mipDy=dY; mipX=pClu->X(); mipY=pClu->Y(); mipQ=pClu->Q();} //current cluster is closer, overwrite data for min cluster
  89     }//clusters loop for intersected chamber
  90
  91     pTrk->SetRICHthetaPhi(helix.Ploc().Theta(),helix.Ploc().Phi()); //store track impact angles with respect to RICH planes
  92     pTrk->SetRICHdxdy(mipDx,mipDy);                                 //distance between track-PC intersection and closest cluster with Qdc>100
  93     pTrk->SetRICHmipXY(mipX,mipY);                                  //position of that closest cluster with Qdc>100
  94     pTrk->SetRICHnclusters(1000000*mipQ);                           //charge of that closest cluster with Qdc>100
  95
  96     if(iMipId==-1)                        {pTrk->SetRICHsignal(kMipQdcCut);  continue;} //no cluster with enough QDC found
  97     if(mipDr>AliRICHParam::DmatchMIP())   {pTrk->SetRICHsignal(kMipDistCut); continue;} //closest cluster with enough carge is still too far
  98
  99     pTrk->SetRICHcluster(iMipId+1000000*iChamber);                                //set mip cluster index
 100     pTrk->SetRICHsignal(recon.ThetaCerenkov(&helix,pRich->Clus(iChamber),nphots));//search for mean Cerenkov angle for this track
 101     pTrk->SetRICHnclusters(1000000*mipQ+nphots);                                  //on return nphots is number of photon clusters accepted in reconstruction
 102     pTrk->SetRICHchi2(recon.GetRingSigma2());
 103
 104     AliDebug(1,Form("Ch=%i PC Intersection=(%5.2f,%5.2f) cm MIP cluster dist=(%5.2f,%5.2f)=%5.2f cm ThetaCkov=%f",
 105                      iChamber,helix.PosPc().X(),helix.PosPc().Y(),            mipDx,mipDy,mipDr,     pTrk->GetRICHsignal()));
 106
 107 //here comes PID calculations
 108 //    CalcProb(pTrack->GetRICHsignal(),pTrack->GetP(),sigmaPID,richPID);
 109   }//ESD tracks loop
 110   AliDebug(1,"Stop pattern recognition");
 111   return 0; // error code: 0=no error;
 112 }//PropagateBack()
 113 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 114 void AliRICHTracker::RecWithStack(TNtupleD *hn)
 115 {
 116 // Reconstruction for particles from STACK. This methode is to be used for RICH standalone when no other detectors are switched on, so normal tracking is not available.
 117 // Arguments: hn- output ntuple where to store all variables
 118 //   Returns: none
 119   AliDebug(1,"Start.");
 120   AliRICH *pRich=((AliRICH*)gAlice->GetDetector("RICH"));
 121
 122 //  pRich->GetLoader()->GetRunLoader()->LoadHeader();
 123   if(!pRich->GetLoader()->GetRunLoader()->TreeK()) pRich->GetLoader()->GetRunLoader()->LoadKinematics();
 124   AliStack *pStack =   pRich->GetLoader()->GetRunLoader()->Stack();
 125   if(!pStack) {AliDebug(1,Form("No STACK found in AliRoot"));return;}
 126   Int_t iNtracks=pStack->GetNtrack();
 127   AliDebug(1,Form(" Start reconstruction with %i track(s) from Stack",iNtracks));
 128
 129   Double_t hnvec[20];
 130
 131   Double_t b=GetFieldMap()->SolenoidField()/10;// magnetic field in Tesla
 132   AliDebug(1,Form("Start with simulated %i tracks in %f Tesla field",iNtracks,b));
 133   TVector3 x0(0,0,0); TVector3 p0(0,0,0);//tmp storage for AliRICHHelix
 134
 135
 136   if(pRich->GetLoader()->LoadRecPoints()) {AliDebug(1,Form("No clusters found in RICH"));return;}
 137   pRich->GetLoader()->TreeR()->GetEntry(0);
 138
 139   AliRICHRecon recon;
 140   for(Int_t iTrackN=0;iTrackN<iNtracks;iTrackN++){//stack particles loop
 141     TParticle *pParticle = pStack->Particle(iTrackN);
 142     if(!pParticle) {AliDebug(1,Form("Not a valid TParticle pointer. Track skipped"));continue;}
 143     AliDebug(1,Form(" PDG code : %i",pParticle->GetPdgCode()));
 144 //
 145 // problem of PDG code of some extra particles to be solved!!!!!!!!!
 146 //
 147 // found problem! Look in TRD directory : codes from Fluka are :
 148 //
 149 //    if ((pdg_code == 10010020) ||
 150 //        (pdg_code == 10010030) ||
 151 //        (pdg_code == 50000050) ||
 152 //        (pdg_code == 50000051) ||
 153 //        (pdg_code == 10020040)) {
 154 //
 155     if(pParticle->GetPdgCode()>=50000050||pParticle->GetPdgCode()==0||pParticle->GetPdgCode()>10000) {AliDebug(1,Form("A photon as track... Track skipped"));continue;}
 156     if(!pParticle->GetPDG()) continue;
 157 //
 158 // to be updated for us!!
 159 //
 160     AliDebug(1,Form("Track %i is a %s with charge %i and momentum %f",
 161             iTrackN,pParticle->GetPDG()->GetName(),Int_t(pParticle->GetPDG()->Charge()),pParticle->P()));
 162 //    if(pParticle->GetMother(0)!=-1) continue; //consider only primaries
 163     if(pParticle->GetPDG()->Charge()==0||TMath::Abs(Int_t(pParticle->GetPDG()->Charge()))!=3) continue; //to avoid photons from stack...
 164     hnvec[0]=pParticle->P();
 165     hnvec[1]=pParticle->GetPDG()->Charge();
 166
 167     p0.SetMagThetaPhi(pParticle->P(),pParticle->Theta(),pParticle->Phi());   x0.SetXYZ(pParticle->Vx(),pParticle->Vy(),pParticle->Vz());
 168     AliRICHHelix helix(x0,p0,TMath::Sign(1,(Int_t)pParticle->GetPDG()->Charge()),b);
 169     Int_t iChamber=helix.RichIntersect(AliRICHParam::Instance());
 170     if(!iChamber) continue;// no intersection with RICH found
 171
 172     hnvec[2]=helix.Ploc().Theta();
 173     hnvec[3]=helix.Ploc().Phi();
 174     hnvec[4]=helix.PosPc().X();
 175     hnvec[5]=helix.PosPc().Y();
 176
 177     Double_t dX,dY,dR,dRmip=9999;      //min distance between clusters and track position on PC
 178     Int_t iMipId=-1;       //index of that min distance cluster
 179     for(Int_t iClu=0;iClu<pRich->Clus(iChamber)->GetEntries();iClu++){//clusters loop for intersected chamber
 180       AliRICHCluster *pClu=(AliRICHCluster*)pRich->Clus(iChamber)->UncheckedAt(iClu);//get pointer to current cluster
 181       pClu->DistXY(helix.PosPc(),dX,dY); dR=TMath::Sqrt(dX*dX+dY*dY);//ditance between current cluster and helix intersection with PC
 182       if(dR<dRmip){dRmip=dR; hnvec[6]=pClu->X();hnvec[7]=pClu->Y();hnvec[8]=pClu->Q();
 183                                                   iMipId=1000000*iChamber+iClu;}//find cluster nearest to the track
 184
 185     }//clusters loop for intersected chamber
 186
 187
 188     hnvec[9]=recon.ThetaCerenkov(&helix,pRich->Clus(iChamber),iMipId); //search for mean Cerenkov angle for this track
 189     hnvec[10]=iMipId;//on return from ThetaCerenkov() contains number of photon candidates accepted
 190     hnvec[11]=(Double_t)iMipId;
 191     hnvec[12]=(Double_t)iChamber;
 192     hnvec[13]=(Double_t)pParticle->GetPdgCode();
 193     if(hn) hn->Fill(hnvec);
 194     AliDebug(1,Form("FINAL Theta Cerenkov=%f",hnvec[9]));
 195   }//stack particles loop
 196
 197   pRich->GetLoader()->UnloadRecPoints();
 198   AliDebug(1,"Stop.");
 199 }//RecWithStack
 200 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 201 void AliRICHTracker::EsdQA(Bool_t isPrint)
 202 {
 203 // Reads ESD file and print out or plot some information for QA
 204 // Arguments: isPrint is a flag to choose between printing (isPrint = kTRUE) and plotting (isPrint = kFALSE)
 205 //   Returns: none
 206
 207   TFile *pFile=TFile::Open("AliESDs.root","read");              if(!pFile) {Printf("ERROR: AliESDs.root does not exist!");return;}
 208   TTree *pTr=(TTree*)pFile->Get("esdTree");                     if(!pTr)   {Printf("ERROR: AliESDs.root, no ESD tree inside!");return;}
 209   AliESD *pEsd=new AliESD;  pTr->SetBranchAddress("ESD", &pEsd);
 210
 211   TH1D     *pProbEl=0,*pProbMu=0,*pProbPi=0,*pProbKa=0,*pProbPr=0,*pMom=0,*pMipQ=0;
 212   TH2F     *pThP=0,*pDxDy=0;
 213   TProfile *pChiTh=0;
 214   if(!isPrint){
 215     TH1::AddDirectory(kFALSE);
 216     pProbEl=new TH1D("RiProbE" ,"Prob e"       ,101  ,0   ,1.05); pProbEl->SetLineColor(kGreen);
 217     pProbPi=new TH1D("RiProbPi","Prob #pi"     ,101  ,0   ,1.05); pProbPi->SetLineColor(kRed);
 218     pProbMu=new TH1D("RiProbMu","Prob #mu"     ,101  ,0   ,1.05); pProbMu->SetLineColor(kBlue);
 219     pProbKa=new TH1D("RiProbK" ,"Prob K"       ,101  ,0   ,1.05);
 220     pProbPr=new TH1D("RiProbP" ,"Prob p"       ,101  ,0   ,1.05);
 221     pMom   =new TH1D("pMom"    ,"Track P, GeV" ,200  ,0   ,20  );
 222     pMipQ  =new TH1D("RiMipQ"  ,"Mip Q, ADC"   ,2000 ,0   ,4000);
 223     pThP   =new TH2F("RiThP"   ,"#theta_{Ckov} radian;P GeV"              ,65 ,-0.5,6.0,75,0,0.75); pThP->SetStats(0);
 224     pDxDy  =new TH2F("RiDxDy"  ,"distance between mip and track;cm",300,-2.5,2.5, 300,-2.5,2.5);
 225     pChiTh =new TProfile("RiChiTh","#chi^{2};#theta_{C}"            ,80 ,0,0.8 , -2,2);
 226   }
 227
 228   Int_t iEvtCnt=0,iTrkCnt=0,iGoodCnt=0;   Float_t bz=0;
 229   for(Int_t iEvt=0;iEvt<pTr->GetEntries();iEvt++){//ESD events loop
 230     pTr->GetEvent(iEvt); iEvtCnt++; if(isPrint) Printf("");
 231     bz=pEsd->GetMagneticField()/10.;
 232     for(Int_t iTrk=0;iTrk<pEsd->GetNumberOfTracks();iTrk++){//ESD tracks loop
 233       AliESDtrack *pTrk=pEsd->GetTrack(iTrk); iTrkCnt++;         //get next reconstructed track and increment total tracks counter
 234
 235       Float_t mom           =pTrk->GetP();                       //track momentum
 236       Double_t sign         =pTrk->GetSign();                    //track sign
 237       Float_t ckov          =pTrk->GetRICHsignal();              //Theta ckov for this track, rad
 238       Float_t chi2          =pTrk->GetRICHchi2();                //Theta ckov error for this track, rad^2
 239       Int_t   qdc           =pTrk->GetRICHnclusters()/1000000;   //Mip candidate charge, qdc
 240       Int_t   nphot         =pTrk->GetRICHnclusters()%1000000;   //number of photon candidates
 241       Float_t dx,dy;         pTrk->GetRICHdxdy(dx,dy);           //distance between mip position and track instersection
 242       Float_t theta,phi;     pTrk->GetRICHthetaPhi(theta,phi);   //track inclination angles in LORS
 243       Double_t pid[5];       pTrk->GetRICHpid(pid);              //pid vector
 244
 245       if(ckov>0) iGoodCnt++;
 246       if(isPrint){
 247         TString comment;
 248         if(ckov>0)                         comment="OK";
 249         else if(ckov==kMipQdcCut)          comment="small QDC";
 250         else if(ckov==kMipDistCut)         comment="mip too far";
 251         else if(ckov==-1)                  comment="no intersection";
 252         Printf("Tr=%2i Q=%4.1f P=%.3f R=%4.2f Th=%6.3f MipQ= %4i Nph=%2i" " rad Prob : e=%.4f mu=%.4f pi=%.4f K=%.4f p=%.4f %s" ,
 253                      iTrk,sign,mom,TMath::Sqrt(dx*dx+dy*dy),ckov,qdc,nphot,           pid[0],pid[1],pid[2],pid[3],pid[4], comment.Data());
 254       }else{//collect hists
 255                                  pMom->Fill(mom);
 256         pMipQ->Fill(qdc);
 257                                  pDxDy->Fill(dx,dy);
 258         pThP->Fill(mom,ckov);
 259                                  pChiTh->Fill(ckov,chi2);
 260         pProbEl->Fill(pid[0]);
 261                                  pProbMu->Fill(pid[1]);
 262         pProbPi->Fill(pid[2]);
 263                                  pProbKa->Fill(pid[3]);
 264         pProbPr->Fill(pid[4]);
 265       }//if plot
 266     }//ESD tracks loop
 267   }//ESD events loop
 268   delete pEsd;  pFile->Close();//close AliESDs.root
 269
 270   TString summary=Form("Events: %i Tracks %i Good RICH: %i Mag Fld %.2f Tesla",iEvtCnt,iTrkCnt,iGoodCnt,bz);
 271   if(isPrint){
 272     Printf(summary.Data());
 273   }else{
 274     TCanvas *pC=new TCanvas("c",summary.Data()); pC->Divide(2,2);
 275     TF1 *pPion = new TF1("RITheor","acos(sqrt(x*x+[0]*[0])/(x*[1]))",1.2,6); pPion->SetLineWidth(1);
 276                                                                     pPion->SetParameter(1,1.292);                                  //ref idx
 277     AliPID ppp;                        pPion->SetLineColor(kRed);   pPion->SetParameter(0,AliPID::ParticleMass(AliPID::kPion));    //mass
 278     TF1 *pKaon = (TF1*)pPion->Clone(); pKaon->SetLineColor(kGreen); pKaon->SetParameter(0,AliPID::ParticleMass(AliPID::kKaon));
 279     TF1 *pProt = (TF1*)pPion->Clone(); pProt->SetLineColor(kBlue);  pProt->SetParameter(0,AliPID::ParticleMass(AliPID::kProton));
 280
 281     pC->cd(1);        pDxDy->Draw();                                                                     //distance between mip and track intersection
 282     pC->cd(2);        pMipQ->Draw();
 283     pC->cd(3);        pThP->Draw();       pPion->Draw("same"); pKaon->Draw("same"); pProt->Draw("same"); //Theta Ckov versus p + theoretical curves
 284     pC->cd(4);        pChiTh->Draw();                                                                    //Theta Ckov error versus theta Ckov
 285
 286     TCanvas *pC2=new TCanvas("c2",summary.Data()); pC2->Divide(2,2);
 287     pC2->cd(1);        pProbPi->Draw(); pProbMu->Draw("same"); pProbEl->Draw("same");
 288     pC2->cd(2);        pProbKa->Draw();
 289     pC2->cd(3);        pProbPr->Draw();
 290     pC2->cd(4);        pMom->Draw();
 291   }
 292 }//EsdQA()
 293 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 294 void AliRICHTracker::MatrixPrint(Double_t probCut)
 295 {
 296 // Reads a set of 3  ESD files from current directory and prints out the matrix of probabilities to pion kaon or proton completely blindly withou nay assumption on the contents of files.
 297 // Normally it implies that those 3 ESDs contain only particles of the same sort namly pions, kaons and protons in that order.
 298 // Arguments: probCut - cut on probability
 299 //   Returns: none
 300   for(Int_t iFile=0;iFile<3;iFile++){
 301     TFile *pFile=TFile::Open(Form("Esd%1i.root",iFile+1),"read"); if(!pFile) {Printf("ERROR: Esd%1i.root does not exist!",iFile+1);return;}
 302     TTree *pTr=(TTree*)pFile->Get("esdTree");                     if(!pTr)   {Printf("ERROR: Esd%1i.root, no ESD tree inside!",iFile+1);return;}
 303     AliESD *pEsd=new AliESD;  pTr->SetBranchAddress("ESD", &pEsd);
 304     Int_t iProtCnt=0,iKaonCnt=0,iPionCnt=0,iUnreconCnt=0,iTrkCnt=0; //counters
 305
 306     for(Int_t iEvt=0;iEvt<pTr->GetEntries();iEvt++){//ESD events loop
 307       pTr->GetEvent(iEvt);
 308       iTrkCnt+=pEsd->GetNumberOfTracks();
 309       for(Int_t iTrk=0;iTrk<pEsd->GetNumberOfTracks();iTrk++){//ESD tracks loop
 310         AliESDtrack *pTrack = pEsd->GetTrack(iTrk);// get next reconstructed track
 311         Float_t dx,dy;         pTrack->GetRICHdxdy(dx,dy);
 312         Float_t theta,phi;     pTrack->GetRICHthetaPhi(theta,phi);
 313         Double_t prob[5];       pTrack->GetRICHpid(prob);
 314         if(pTrack->GetRICHsignal()>0){
 315           if(prob[4]>probCut)                         iProtCnt++;
 316           if(prob[3]>probCut)                         iKaonCnt++;
 317           if((prob[0]+prob[1]+prob[2])>probCut)       iPionCnt++;
 318         } else
 319           iUnreconCnt++;
 320       }//ESD tracks loop
 321
 322     }//ESD events loop
 323     Printf("Bz=%5.2f Events=%i Total tracks=%i No recognized tracks=%i Pion=%i Kaon=%i Proton=%i ProbCut=%.2f",
 324         0.1*pEsd->GetMagneticField(),pTr->GetEntries(),iTrkCnt,iUnreconCnt,iPionCnt,iKaonCnt,iProtCnt,probCut);
 325     delete pEsd;  pFile->Close();//close AliESDs.root
 326   }//files loop
 327 }
 328 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++