From f4bb06ee96c6d0cd25ce6deff4949f3445e416d3 Mon Sep 17 00:00:00 2001 From: mivanov Date: Tue, 10 Dec 2013 00:15:22 +0100 Subject: [PATCH] // // Adding Code to select identified V0 for the PID purposes // As an input, chain of filtered trees (output of AliAnalysisTaskFilteredTree task ) is used void filterPIDSelected( const char * chfinput="highptAll.list"){ // Parameter: // finput - name of the list file or the name of file itself // Output: // file - V0Selected.root --- PWGPP/TPC/macros/filterPIDSelected.C | 183 +++++++++++++++++++++++++++ 1 file changed, 183 insertions(+) create mode 100644 PWGPP/TPC/macros/filterPIDSelected.C diff --git a/PWGPP/TPC/macros/filterPIDSelected.C b/PWGPP/TPC/macros/filterPIDSelected.C new file mode 100644 index 00000000000..8c54362be43 --- /dev/null +++ b/PWGPP/TPC/macros/filterPIDSelected.C @@ -0,0 +1,183 @@ +/* + Macro to make additional filter of the V0 to select clean sample of identified particles. + As an input the + + .x $HOME/rootlogon.C + .L $ALICE_ROOT/PWGPP/TPC/macros/filterPIDSelected.C+ + + */ + + + + + +#include "TFile.h" +#include "TTree.h" +#include "TVectorD.h" +#include "TMatrixD.h" +#include "TH2.h" +#include "TF1.h" +#include "TTreeStream.h" +#include "AliMathBase.h" +#include "TSystem.h" +#include "TChain.h" +#include "TDatabasePDG.h" +#include "TRandom.h" +#include "AliTPCcalibBase.h" +#include "TCanvas.h" +#include "TLegend.h" +// +#include "AliESDv0.h" +#include "AliESDtrack.h" +#include "TMath.h" +#include "AliXRDPROOFtoolkit.h" +#include "TStatToolkit.h" +#include "TCut.h" + +TTree * tree = 0; +TTreeSRedirector *pcstream = 0; //new TTreeSRedirector("trend.root"); +// + + +void filterPIDSelected( const char * chfinput="highptAll.list"){ + // + // Code to select identified V0 for the PID + // As an input chain of filter trees is used + // Parameter: + // finput - name of the list file or the name of file itself + // Oputput: + // file - V0Selected.root + // + // + TTree * chain = 0; + if (TString(chfinput).Contains(".list")) { + chain = AliXRDPROOFtoolkit::MakeChainRandom(chfinput,"V0s",0,1000); + }else{ + TFile * finput= TFile::Open(chfinput); + if (!finput) finput= TFile::Open(TString::Format("%s#FilterEvents_Trees.root",finput)); + chain=(TTree*)finput->Get("V0s"); + } + chain->SetCacheSize(1000000000); + // + TDatabasePDG pdg; + Double_t massLambda = pdg.GetParticle("Lambda0")->Mass(); + Double_t massK0 = pdg.GetParticle("K0")->Mass(); + Double_t massPion = pdg.GetParticle("pi+")->Mass(); + Double_t massProton = pdg.GetParticle("proton")->Mass(); + // + // + chain->SetAlias("massPion",Form("(%f+0)",massPion)); + chain->SetAlias("massProton",Form("(%f+0)",massProton)); + chain->SetAlias("massK0",Form("(%f+0)",massK0)); + chain->SetAlias("massLambda",Form("(%f+0)",massLambda)); + // delta of mass + chain->SetAlias("K0Delta","(v0.GetEffMass(2,2)-massK0)"); + chain->SetAlias("LDelta","(v0.GetEffMass(4,2)-massLambda)"); + chain->SetAlias("ALDelta","(v0.GetEffMass(2,4)-massLambda)"); + chain->SetAlias("EDelta","(v0.GetEffMass(0,0))"); + // pull of the mass + chain->SetAlias("K0Pull","(v0.GetEffMass(2,2)-massK0)/v0.GetKFInfo(2,2,1)"); + chain->SetAlias("LPull","(v0.GetEffMass(4,2)-massLambda)/v0.GetKFInfo(4,2,1)"); + chain->SetAlias("ALPull","(v0.GetEffMass(2,4)-massLambda)/v0.GetKFInfo(2,4,1)"); + chain->SetAlias("EPull","EDelta/v0.GetKFInfo(0,0,1)"); + // effective pull of the mass - (empirical values form fits) + chain->SetAlias("K0PullEff","K0Delta/sqrt((3.63321e-03)**2+(5.68795e-04*v0.Pt())**2)"); + chain->SetAlias("LPullEff","LDelta/sqrt((1.5e-03)**2+(1.8e-04*v0.Pt())**2)"); + chain->SetAlias("ALPullEff","ALDelta/sqrt((1.5e-03)**2+(1.8e-04*v0.Pt())**2)"); + chain->SetAlias("EPullEff","v0.GetEffMass(0,0)/sqrt((5e-03)**2+(1.e-04*v0.Pt())**2)"); + // + // + chain->SetAlias("dEdx0DProton","AliMathBase::BetheBlochAleph(track0.fIp.P()/massProton)"); + chain->SetAlias("dEdx1DProton","AliMathBase::BetheBlochAleph(track1.fIp.P()/massProton)"); + chain->SetAlias("dEdx0DPion","AliMathBase::BetheBlochAleph(track0.fIp.P()/massPion)"); + chain->SetAlias("dEdx1DPion","AliMathBase::BetheBlochAleph(track1.fIp.P()/massPion)"); + // + // V0 - cuts -PID, + // + chain->SetAlias("cutDist","sqrt((track0.fIp.fP[0]-track1.fIp.fP[0])**2+(track0.fIp.fP[1]-track1.fIp.fP[1])**2)>3"); + chain->SetAlias("cutLong","track0.GetTPCClusterInfo(3,1,0)-5*abs(track0.fP[4])>130&&track1.GetTPCClusterInfo(3,1,0)>130-5*abs(track0.fP[4])"); + chain->SetAlias("cutPID","track0.fTPCsignal>0&&track1.fTPCsignal>0"); + chain->SetAlias("cutResol","sqrt(track0.fC[14]/track0.fP[4])<0.15&&sqrt(track1.fC[14]/track1.fP[4])<0.15"); + chain->SetAlias("cutV0","cutPID&&cutDist&&cutLong&&cutResol"); + // + // + chain->SetAlias("K0Selected", "abs(K0Pull)<3. &&abs(K0PullEff)<3. && abs(LPull)>3 && abs(ALPull)>3 &&v0.PtArmV0()>0.11"); + chain->SetAlias("LambdaSelected", "abs(LPull)<3. &&abs(LPullEff)<3. && abs(K0Pull)>3 && abs(EPull)>3 && abs(EDelta)>0.05"); + chain->SetAlias("ALambdaSelected", "abs(ALPull)<3. &&abs(ALPullEff)<3 && abs(K0Pull)>3 && abs(EPull)>3 &&abs(EDelta)>0.05"); + // + chain->SetAlias("GammaSelected", "abs(EPull)<3 && abs(K0Pull)>3 && abs(LPull)>3 && abs(ALPull)>3"); + // + // + TFile *fselected = TFile::Open("V0Selected.root","recreate"); + TTree * treeK0 = chain->CopyTree("type==8&&cutV0&&K0Selected"); + TTree * treeLambda = chain->CopyTree("type==4&&cutV0&&LambdaSelected"); + TTree * treeALambda = chain->CopyTree("type==2&&cutV0&&ALambdaSelected"); + TTree * treeGamma = chain->CopyTree("type==1&&cutV0&&GammaSelected"); + // + TTree * trees[4]={treeK0,treeLambda, treeGamma,treeALambda}; + TList * aliases = chain->GetListOfAliases(); + Int_t nalias= aliases->GetEntries(); + + for (Int_t i=0; i<4; i++){ + for (Int_t ialias=0; ialiasAt(ialias); + trees[i]->SetAlias(alias->GetName(),alias->GetTitle()); + } + } + treeK0->Write("treeK0"); + treeLambda->Write("treeLambda"); + treeALambda->Write("treeALambda"); + treeGamma->Write("treeGamma"); + fselected->Close(); + // +} + + +void FitPIDNCLSelected(){ + // + // fit the probability to find the cluster + // + Int_t kmarkers[5]={20,21,25,24,22}; + Int_t kcolors[5]={1,2,4,5,7}; + const char *chname[5]={"Proton","Pion","Electron"}; + + TFile f("V0Selected.root"); + TTree * treeLambda = (TTree*)f.Get("treeLambda"); + TTree * treeK0 = (TTree*)f.Get("treeK0"); + TTree * treeGamma = (TTree*)f.Get("treeGamma"); + // + // + TH2 *his3D=0; + TObjArray fitArray(3); + TH1D * hisNcldEdx[3]={0,0,0}; + TCut cutOut="abs(track0.fP[4])<2"; + 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"); + hisNcldEdx[0]=(TH1D*)treeLambda->GetHistogram()->Clone(); + 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"); + hisNcldEdx[1]=(TH1D*)treeK0->GetHistogram()->Clone(); + 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"); + hisNcldEdx[2]=(TH1D*)treeGamma->GetHistogram()->Clone(); + + TF1 fnclQ("fnclQ","[0]+[1]*exp(-[2]*abs(x))"); + fnclQ.SetParameters(0.02,5,3); + hisNcldEdx[0]->Fit(&fnclQ); + hisNcldEdx[1]->Fit(&fnclQ); + //hisNcldEdx[2]->Fit(&fnclQ); + + TCanvas * canvas = new TCanvas("canvasNCL0","canvasNCL0",600,500); + TLegend *legend = new TLegend(0.5,0.6,0.89,0.89,"Cluster finder eff."); + legend->SetBorderSize(0); + for (Int_t i=0; i<3; i++){ + hisNcldEdx[i]->SetMinimum(0); + hisNcldEdx[i]->SetMaximum(0.3); + hisNcldEdx[i]->SetMarkerColor(kcolors[i]); + hisNcldEdx[i]->SetMarkerStyle(kmarkers[i]); + hisNcldEdx[i]->GetXaxis()->SetTitle("Q ~ #sqrt{1+tan(#theta)^2}dE/dx"); + hisNcldEdx[i]->GetYaxis()->SetTitle("p_{cl}"); + if (i==0)hisNcldEdx[i]->Draw(""); + hisNcldEdx[i]->Draw("same"); + legend->AddEntry(hisNcldEdx[i],chname[i]); + } + legend->Draw(); + +} -- 2.43.0