PWGPP/TPC/macros/filterPIDSelected.C

   1 /*
   2   Macro to make additional filter of the V0 to select clean sample of identified particles.
   3   As an input the
   4
   5   .x $HOME/rootlogon.C
   6   .L $ALICE_ROOT/PWGPP/TPC/macros/filterPIDSelected.C+
   7
   8  */
   9
  10
  11
  12
  13
  14 #include "TFile.h"
  15 #include "TTree.h"
  16 #include "TVectorD.h"
  17 #include "TMatrixD.h"
  18 #include "TH2.h"
  19 #include "TF1.h"
  20 #include "TTreeStream.h"
  21 #include "AliMathBase.h"
  22 #include "TSystem.h"
  23 #include "TChain.h"
  24 #include "TDatabasePDG.h"
  25 #include "TRandom.h"
  26 #include "AliTPCcalibBase.h"
  27 #include "TCanvas.h"
  28 #include "TLegend.h"
  29 //
  30 #include "AliESDv0.h"
  31 #include "AliESDtrack.h"
  32 #include "TMath.h"
  33 #include "AliXRDPROOFtoolkit.h"
  34 #include "TStatToolkit.h"
  35 #include "TCut.h"
  36
  37 TTree * tree  = 0;
  38 TTreeSRedirector *pcstream = 0; //new TTreeSRedirector("trend.root");
  39 //
  40 void filterPIDSelectedTOF( const char * chfinput="highptAll.list");
  41 void filterPIDSelectedV0( const char * chfinput="highptAll.list");
  42
  43 void filterPIDSelected( const char * chfinput="highptAll.list"){
  44   filterPIDSelectedTOF(chfinput);
  45   filterPIDSelectedV0(chfinput);
  46 }
  47
  48 void filterPIDSelectedV0( const char * chfinput){
  49   //
  50   // Code to select identified V0 for the PID
  51   // As an input chain of filter trees is used
  52   // Parameter:
  53   //   finput - name of the list file or the name of file itself
  54   // Oputput:
  55   //   file - V0Selected.root
  56   //
  57   //
  58   TTree * chain  = 0;
  59   if (TString(chfinput).Contains(".list")) {
  60     chain = AliXRDPROOFtoolkit::MakeChainRandom(chfinput,"V0s",0,1000);
  61   }else{
  62     TFile * finput= TFile::Open(chfinput);
  63     if (!finput) finput= TFile::Open(TString::Format("%s#FilterEvents_Trees.root",finput));
  64     chain=(TTree*)finput->Get("V0s");
  65   }
  66   chain->SetCacheSize(1000000000);
  67   //
  68   TDatabasePDG pdg;
  69   Double_t massLambda = pdg.GetParticle("Lambda0")->Mass();
  70   Double_t massK0 = pdg.GetParticle("K0")->Mass();
  71   Double_t massPion = pdg.GetParticle("pi+")->Mass();
  72   Double_t massProton = pdg.GetParticle("proton")->Mass();
  73   //
  74   //
  75   chain->SetAlias("massPion",Form("(%f+0)",massPion));
  76   chain->SetAlias("massProton",Form("(%f+0)",massProton));
  77   chain->SetAlias("massK0",Form("(%f+0)",massK0));
  78   chain->SetAlias("massLambda",Form("(%f+0)",massLambda));
  79   // delta of mass
  80   chain->SetAlias("K0Delta","(v0.GetEffMass(2,2)-massK0)");
  81   chain->SetAlias("LDelta","(v0.GetEffMass(4,2)-massLambda)");
  82   chain->SetAlias("ALDelta","(v0.GetEffMass(2,4)-massLambda)");
  83   chain->SetAlias("EDelta","(v0.GetEffMass(0,0))");
  84   // pull of the mass
  85   chain->SetAlias("K0Pull","(v0.GetEffMass(2,2)-massK0)/v0.GetKFInfo(2,2,1)");
  86   chain->SetAlias("LPull","(v0.GetEffMass(4,2)-massLambda)/v0.GetKFInfo(4,2,1)");
  87   chain->SetAlias("ALPull","(v0.GetEffMass(2,4)-massLambda)/v0.GetKFInfo(2,4,1)");
  88   chain->SetAlias("EPull","EDelta/v0.GetKFInfo(0,0,1)");
  89   // effective pull of the mass - (empirical values form fits)
  90   chain->SetAlias("K0PullEff","K0Delta/sqrt((3.63321e-03)**2+(5.68795e-04*v0.Pt())**2)");
  91   chain->SetAlias("LPullEff","LDelta/sqrt((1.5e-03)**2+(1.8e-04*v0.Pt())**2)");
  92   chain->SetAlias("ALPullEff","ALDelta/sqrt((1.5e-03)**2+(1.8e-04*v0.Pt())**2)");
  93   chain->SetAlias("EPullEff","v0.GetEffMass(0,0)/sqrt((5e-03)**2+(1.e-04*v0.Pt())**2)");
  94   //
  95   //
  96   chain->SetAlias("dEdx0DProton","AliMathBase::BetheBlochAleph(track0.fIp.P()/massProton)");
  97   chain->SetAlias("dEdx1DProton","AliMathBase::BetheBlochAleph(track1.fIp.P()/massProton)");
  98   chain->SetAlias("dEdx0DPion","AliMathBase::BetheBlochAleph(track0.fIp.P()/massPion)");
  99   chain->SetAlias("dEdx1DPion","AliMathBase::BetheBlochAleph(track1.fIp.P()/massPion)");
 100   //
 101   // V0 - cuts -PID,
 102   //
 103   chain->SetAlias("cutDist","sqrt((track0.fIp.fP[0]-track1.fIp.fP[0])**2+(track0.fIp.fP[1]-track1.fIp.fP[1])**2)>3");
 104   chain->SetAlias("cutLong","track0.GetTPCClusterInfo(3,1,0)+5*abs(track0.fP[4])>130&&track1.GetTPCClusterInfo(3,1,0)>130-5*abs(track1.fP[4])");
 105   chain->SetAlias("cutPID","track0.fTPCsignal>0&&track1.fTPCsignal>0");
 106   chain->SetAlias("cutResol","sqrt(track0.fC[14]/track0.fP[4])<0.15&&sqrt(track1.fC[14]/track1.fP[4])<0.15");
 107   chain->SetAlias("cutV0","cutPID&&cutDist&&cutLong&&cutResol");
 108   //
 109   //
 110   chain->SetAlias("K0Selected",      "abs(K0Pull)<3. &&abs(K0PullEff)<3.  && abs(LPull)>3  && abs(ALPull)>3  &&v0.PtArmV0()>0.11");
 111   chain->SetAlias("LambdaSelected",  "abs(LPull)<3.  &&abs(LPullEff)<3.   && abs(K0Pull)>3 && abs(EPull)>3 && abs(EDelta)>0.05");
 112   chain->SetAlias("ALambdaSelected", "abs(ALPull)<3. &&abs(ALPullEff)<3   && abs(K0Pull)>3 && abs(EPull)>3 &&abs(EDelta)>0.05");
 113   //
 114   chain->SetAlias("GammaSelected", "abs(EPull)<3     && abs(K0Pull)>3 && abs(LPull)>3 && abs(ALPull)>3");
 115   //
 116   //
 117   TFile *fselected = TFile::Open("V0Selected.root","recreate");
 118   TTree * treeK0     =   chain->CopyTree("type==8&&cutV0&&K0Selected");
 119   TTree * treeLambda =   chain->CopyTree("type==4&&cutV0&&LambdaSelected");
 120   TTree * treeALambda =   chain->CopyTree("type==2&&cutV0&&ALambdaSelected");
 121   TTree * treeGamma =   chain->CopyTree("type==1&&cutV0&&GammaSelected");
 122   //
 123   TTree * trees[4]={treeK0,treeLambda, treeGamma,treeALambda};
 124   TList * aliases = chain->GetListOfAliases();
 125   Int_t nalias= aliases->GetEntries();
 126
 127   for (Int_t i=0; i<4; i++){
 128     for (Int_t ialias=0; ialias<nalias; ialias++){
 129       TNamed *alias = (TNamed*)aliases->At(ialias);
 130       trees[i]->SetAlias(alias->GetName(),alias->GetTitle());
 131     }
 132   }
 133   treeK0->Write("treeK0");
 134   treeLambda->Write("treeLambda");
 135   treeALambda->Write("treeALambda");
 136   treeGamma->Write("treeGamma");
 137   fselected->Close();
 138   //
 139 }
 140
 141 void filterPIDSelectedTOF( const char * chfinput){
 142   //
 143   // Code to select identified V0 for the PID
 144   // As an input chain of filter trees is used
 145    //
 146   TTree * chain  = 0;
 147   if (TString(chfinput).Contains(".list")) {
 148     chain = AliXRDPROOFtoolkit::MakeChainRandom(chfinput,"highPt",0,1000);
 149   }else{
 150     TFile * finput= TFile::Open(chfinput);
 151     if (!finput) finput= TFile::Open(TString::Format("%s#FilterEvents_Trees.root",finput));
 152     chain=(TTree*)finput->Get("highPt");
 153   }
 154   chain->SetCacheSize(1000000000);
 155   //
 156   TDatabasePDG pdg;
 157   Double_t massLambda = pdg.GetParticle("Lambda0")->Mass();
 158   Double_t massK0 = pdg.GetParticle("K0")->Mass();
 159   Double_t massPion = pdg.GetParticle("pi+")->Mass();
 160   Double_t massProton = pdg.GetParticle("proton")->Mass();
 161   chain->SetAlias("cutLong","esdTrack.GetTPCClusterInfo(3,1,0)+5*abs(esdTrack.fP[4])>130");
 162   TFile *fselected  = TFile::Open("TOFSelected.root","recreate");
 163
 164   TCut cutDeltaProton="abs((esdTrack.fTrackTime[4]-esdTrack.fTrackTime[3]))>400&&abs(esdTrack.fTOFsignalDz<3)";
 165   TCut cutDeltaKaon="abs((esdTrack.fTrackTime[3]-esdTrack.fTrackTime[2]))>400&&abs(esdTrack.fTOFsignalDz<3)";
 166   TCut cutDeltaPion="abs((esdTrack.fTrackTime[2]-esdTrack.fTrackTime[0]))>200&&abs(esdTrack.fTOFsignalDz<3)";
 167
 168   TTree * treeEl  = chain->CopyTree(cutDeltaPion+"cutLong&&esdTrack.fTOFr[0]>0.3+max(2*max(esdTrack.fTOFr[1],esdTrack.fTOFr[3]),esdTrack.fTOFr[4])");
 169   TTree * treePion  = chain->CopyTree(cutDeltaPion+"cutLong&&esdTrack.fTOFr[2]>0.3+max(max(esdTrack.fTOFr[0],esdTrack.fTOFr[3]),esdTrack.fTOFr[4])");
 170   TTree * treeKaon  = chain->CopyTree(cutDeltaKaon+"cutLong&&esdTrack.fTOFr[3]>0.3+max(max(esdTrack.fTOFr[0],2*esdTrack.fTOFr[2]),esdTrack.fTOFr[4])");
 171   TTree * treeProton  = chain->CopyTree(cutDeltaProton+"cutLong&&esdTrack.fTOFr[4]>0.3+max(max(esdTrack.fTOFr[0],2*esdTrack.fTOFr[2]),esdTrack.fTOFr[3])");
 172   treeEl->Write("treeEl");
 173   treePion->Write("treePion");
 174   treeKaon->Write("treeKaon");
 175   treeProton->Write("treeProton");
 176   //
 177   fselected->Close();
 178 }
 179
 180 void FitPIDNCLSelected(){
 181   //
 182   // fit the probability to find the cluster
 183   //
 184   Int_t kmarkers[5]={20,21,25,24,22};
 185   Int_t kcolors[5]={1,2,4,5,7};
 186   const char *chname[5]={"Proton","Pion","Electron"};
 187
 188   TFile f("V0Selected.root");
 189   TTree * treeLambda = (TTree*)f.Get("treeLambda");
 190   TTree * treeK0 = (TTree*)f.Get("treeK0");
 191   TTree * treeGamma = (TTree*)f.Get("treeGamma");
 192   //
 193   //
 194   TH2 *his3D=0;
 195   TObjArray  fitArray(3);
 196   TH1D * hisNcldEdx[3]={0,0,0};
 197   TCut cutOut="abs(track0.fP[4])<2";
 198   treeLambda->Draw("(1-track0.GetTPCClusterInfo(2,0)):sqrt(1+track0.fP[3]**2)*AliMathBase::BetheBlochAleph(track0.fIp.P()/massProton)>>hisNclProton(50,1,2.)",cutOut,"prof");
 199   hisNcldEdx[0]=(TH1D*)treeLambda->GetHistogram()->Clone();
 200   treeK0->Draw("(1-track0.GetTPCClusterInfo(2,0)):sqrt(1+track0.fP[3]**2)*AliMathBase::BetheBlochAleph(track0.fIp.P()/massPion)>>hisNclPion(50,1,2.)",cutOut,"prof");
 201   hisNcldEdx[1]=(TH1D*)treeK0->GetHistogram()->Clone();
 202   treeGamma->Draw("(1-track0.GetTPCClusterInfo(2,0)):sqrt(1+track0.fP[3]**2)*AliMathBase::BetheBlochAleph(track0.fIp.P()/0.0005)>>hisNclPion(50,1,2.5)",cutOut,"prof");
 203   hisNcldEdx[2]=(TH1D*)treeGamma->GetHistogram()->Clone();
 204
 205   TF1 fnclQ("fnclQ","[0]+[1]*exp(-[2]*abs(x))");
 206   fnclQ.SetParameters(0.02,5,3);
 207   hisNcldEdx[0]->Fit(&fnclQ);
 208   hisNcldEdx[1]->Fit(&fnclQ);
 209   //hisNcldEdx[2]->Fit(&fnclQ);
 210
 211   TCanvas * canvas = new TCanvas("canvasNCL0","canvasNCL0",600,500);
 212   TLegend *legend = new TLegend(0.5,0.6,0.89,0.89,"Cluster finder eff.");
 213   legend->SetBorderSize(0);
 214   for (Int_t i=0; i<3; i++){
 215     hisNcldEdx[i]->SetMinimum(0);
 216     hisNcldEdx[i]->SetMaximum(0.3);
 217     hisNcldEdx[i]->SetMarkerColor(kcolors[i]);
 218     hisNcldEdx[i]->SetMarkerStyle(kmarkers[i]);
 219     hisNcldEdx[i]->GetXaxis()->SetTitle("Q ~ #sqrt{1+tan(#theta)^2}dE/dx");
 220     hisNcldEdx[i]->GetYaxis()->SetTitle("p_{cl}");
 221     if (i==0)hisNcldEdx[i]->Draw("");
 222     hisNcldEdx[i]->Draw("same");
 223     legend->AddEntry(hisNcldEdx[i],chname[i]);
 224   }
 225   legend->Draw();
 226
 227 }