]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/SPECTRA/Kinks/AliAnalysisKinkESDat.cxx
kinks updates by Martha
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / Kinks / AliAnalysisKinkESDat.cxx
CommitLineData
828a6e0c 1/**************************************************************************
fcd98e58 2 * Authors: Martha Spyropoulou-Stassinaki and the members *
3 * of the Greek group at Physics Department of Athens University *
4 * Paraskevi Ganoti, Anastasia Belogianni and Filimon Roukoutakis. *
5 * The method is applied in pp and Pb-Pb real data. *
828a6e0c 6 * *
828a6e0c 7 **************************************************************************/
8
9//-----------------------------------------------------------------
10// AliAnalysisKinkESDat class
11// Example of an analysis task for kink topology study
12// Kaons from kink topology are 'identified' in this code
13//-----------------------------------------------------------------
14
15#include "TChain.h"
16#include "TTree.h"
17#include "TH1F.h"
18#include "TH2F.h"
19#include "TH3F.h"
20#include "TH1D.h"
21#include "TH2D.h"
22#include "TParticle.h"
23#include <TVector3.h>
24#include "TF1.h"
25
26#include "AliAnalysisTask.h"
27#include "AliAnalysisManager.h"
28
29#include "AliVEvent.h"
30#include "AliESDEvent.h"
31#include "AliMCEvent.h"
32#include "AliAnalysisKinkESDat.h"
33#include "AliStack.h"
34#include "AliESDpid.h"
35#include "AliPID.h"
36#include "AliESDkink.h"
37#include "AliESDtrack.h"
38#include "AliPhysicsSelectionTask.h"
39#include "AliInputEventHandler.h"
ad1036b2 40#include "AliESDInputHandler.h"
828a6e0c 41#include "AliAnalysisManager.h"
42#include "AliVEvent.h"
828a6e0c 43#include "AliESDtrackCuts.h"
fcd98e58 44#include "AliPIDResponse.h"
8e57b720 45/////////#include "AliTENDERSupplies.h"
828a6e0c 46ClassImp(AliAnalysisKinkESDat)
47
48
49//________________________________________________________________________
50AliAnalysisKinkESDat::AliAnalysisKinkESDat(const char *name)
51 : AliAnalysisTaskSE(name), fHistPtESD(0),fHistPt(0),fHistQtAll(0),fHistQt1(0),fHistQt2(0)
52 , fHistPtKaon(0),fHistPtKPDG(0),fHistEta(0),fHistEtaK(0),fptKMC(0),fMultiplMC(0),fESDMult(0),fgenpt(0),frad(0),
ad1036b2 53 fKinkKaon(0),fKinKRbn(0), fKinkKaonBg(0), fM1kaon(0), fPtKink(0), fptKink(0),
54 fAngMomK(0),fAngMomPi(0), fAngMomKC(0), fMultESDK(0), fMultMCK(0),
828a6e0c 55 fSignPtNcl(0), fSignPtEta(0), fEtaNcl(0), fSignPt(0), fChi2NclTPC(0), fRatChi2Ncl(0), fRadiusNcl(0), fTPCSgnlP(0),
56 fTPCSgnlPa(0), fRpr(0),fZpr(0), fdcatoVxXY(0), fnSigmToVx(0), fKinkMothDau(0),
ad1036b2 57 fZvXv(0),fZvYv(0), fXvYv(0), fHistPtKaoP(0), fHistPtKaoN(0),frapiKESD(0), flifetime(), fradLK(0),
828a6e0c 58 fradPtRpDt(0), fInvMuNuAll(0), fQtInvM(0),
ad1036b2 59 fDCAkink(0), fPosiKink(0), fPosiKinkK(0),fPosiKinKXZ(0), fPosiKinKYZ(0), fPosiKinKBg(0), fQtMothP(0), fTPCSgnlPtpc(0),
60 fTPCMomNSgnl(0), fMothKinkMomSgnl(0), fNSigmTPC(0), fTPCSgnlKinkDau(0), fPtKinkPos(0), fPtKinkNeg(0), fRadNclCln(0),
eb428cd9 61 fRatioCrossedRows(0), fRatioCrossedRowsKink(0),fRadiusPt(0), fRadiusPtcln(0), fInvMassMuNuPt(0), fInvMassMuNuPtAll(0),fPtCut1(0), fPtCut2(0),
c1e83dc8 62 fPtCut3(0), fAngMomKKinks(0),
828a6e0c 63 f1(0), f2(0),
5220a87d 64 fListOfHistos(0),fLowMulcut(-1),fUpMulcut(-1), fKinkRadUp(200), fKinkRadLow(130), fLowCluster(20), fLowQt(.12), fCutsMul(0), fMaxDCAtoVtxCut(0), fPIDResponse(0)
828a6e0c 65{
66 // Constructor
67
68 // Define input and output slots here
69 // Input slot #0 works with a TChain
70 // DefineInput(0, TChain::Class());
71//----------------------Marek multiplicity bins
72 fCutsMul=new AliESDtrackCuts("Mul","Mul");
73 fCutsMul->SetMinNClustersTPC(70);
74 fCutsMul->SetMaxChi2PerClusterTPC(4);
75 fCutsMul->SetAcceptKinkDaughters(kFALSE);
76 fCutsMul->SetRequireTPCRefit(kTRUE);
77 // ITS
78 fCutsMul->SetRequireITSRefit(kTRUE);
79 fCutsMul->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
80 AliESDtrackCuts::kAny);
81 fCutsMul->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
82
83 fCutsMul->SetMaxDCAToVertexZ(2);
84 fCutsMul->SetDCAToVertex2D(kFALSE);
85 fCutsMul->SetRequireSigmaToVertex(kFALSE);
86
87 fCutsMul->SetEtaRange(-0.8,+0.8);
88 fCutsMul->SetPtRange(0.15, 1e10);
89
8e57b720 90 fMaxDCAtoVtxCut=new AliESDtrackCuts("fMaxDCA", "fMaxDCA");
91 fMaxDCAtoVtxCut->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
92 fMaxDCAtoVtxCut->SetMaxChi2TPCConstrainedGlobal(36);
93// esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
94 // esdTrackCuts->SetMaxChi2TPCConstrainedGlobal(36);
95
96
828a6e0c 97
98 DefineOutput(1, TList::Class());
99}
100
101//________________________________________________________________________
102void AliAnalysisKinkESDat::UserCreateOutputObjects()
103{
104 // Create histograms
105 // Called once
106
107 f1=new TF1("f1","((atan([0]*[1]*(1.0/(sqrt((x^2)*(1.0-([1]^2))-([0]^2)*([1]^2))))))*180.)/[2]",1.1,10.0);
108 f1->SetParameter(0,0.493677);
109 f1->SetParameter(1,0.9127037);
110 f1->SetParameter(2,TMath::Pi());
111
112
113 f2=new TF1("f2","((atan([0]*[1]*(1.0/(sqrt((x^2)*(1.0-([1]^2))-([0]^2)*([1]^2))))))*180.)/[2]",0.1,10.0);
114 f2->SetParameter(0,0.13957018);
115 f2->SetParameter(1,0.2731374);
116 f2->SetParameter(2,TMath::Pi());
117 //Open file 1= CAF
118 //OpenFile(1);
ad1036b2 119 Double_t gPt7K0[45] = {0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8,0.9,1.0,
120 1.1, 1.2, 1.3, 1.4, 1.5, 1.6,1.7,1.8,1.9, 2.0,
121 2.2, 2.4, 2.6, 2.8, 3.0, 3.3, 3.6, 3.9,
fcd98e58 122 4.2, 4.6,5.0, 5.4, 5.9, 6.5, 7.0,7.5, 8.0,8.5, 9.2, 10., 11., 12., 13.5,15.0 }; // David K0
ad1036b2 123
124 Double_t gPt7TOF[47] = { 0.2,0.25, 0.3,0.35, 0.4,0.45, 0.5,0.55, 0.6,0.65, 0.7,0.75, 0.8, 0.85, 0.9, 0.95, 1.0,
828a6e0c 125 1.1, 1.2, 1.3, 1.4, 1.5, 1.6,1.7,1.8,1.9, 2.0,
126 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7,2.8, 2.9, 3.0,
ad1036b2 127 3.2, 3.4, 3.6, 3.8, 4.0, 4.2, 4.4, 4.6, 4.8,5.0 }; // Barbara TOF Kch
828a6e0c 128
ad1036b2 129 fHistPtESD = new TH1F("fHistPtESD", "P_{T} distribution",300, 0.0,15.0);
828a6e0c 130 fHistPtESD->GetXaxis()->SetTitle("P_{T} (GeV/c)");
131 fHistPtESD->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
132 fHistPtESD->SetMarkerStyle(kFullCircle);
ad1036b2 133 fHistPt = new TH1F("fHistPt", "P_{T} distribution",300, 0.0,15.0);
828a6e0c 134 fHistQtAll = new TH1F("fHistQtAll", "Q_{T} distr All Kinks ",100, 0.0,.300);
135 fHistQt1= new TH1F("fHistQt1", "Q_{T} distribution",100, 0.0,.300);
136 fHistQt2= new TH1F("fHistQt2", "Q_{T} distribution",100, 0.0,.300);
ad1036b2 137 fHistPtKaon = new TH1F("fHistPtKaon", "P_{T}Kaon distribution",300, 0.0,15.0);
138 fHistPtKPDG = new TH1F("fHistPtKPDG", "P_{T}Kaon distribution",300, 0.0,15.0);
828a6e0c 139 fHistEta= new TH1F("fHistEta", "Eta distribution", 26,-1.3, 1.3);
140 fHistEtaK= new TH1F("fHistEtaK", "EtaK distribution", 26,-1.3, 1.3);
ad1036b2 141 fptKMC= new TH1F("fptKMC", "P_{T}Kaon generated",300, 0.0,15.0);
828a6e0c 142 fMultiplMC= new TH1F("fMultiplMC", "charge multiplicity MC",100, 0.0,300.0);
143 fESDMult= new TH1F("fESDMult", "charge multipliESD",100, 0.0,300.0);
ad1036b2 144 fgenpt= new TH1F("fgenpt", "genpt K distribution",300, 0.0,15.0);
828a6e0c 145 //frad= new TH1F("frad", "radius K generated",100, 50., 250.0);
146 frad= new TH1F("frad", "radius K generated",100, 0.,1000.0);
ad1036b2 147 fKinkKaon= new TH1F("fKinkKaon", "P_{T}Kaon kinks identi",300, 0.0,15.0);
148 fKinKRbn= new TH1F("fKinKRbn", "p_{t}Kaon kinks identi[GeV/c],Entries",46,gPt7TOF);
149 fKinkKaonBg= new TH1F("fKinkKaonBg", "P_{T}Kaon kinks backgr",300, 0.0,15.0);
c1e83dc8 150 //fM1kaon= new TH1F("fM1kaon","Invar m(kaon) from kink->mu+netrino decay",180,0.10, 1.0);
151 fM1kaon= new TH1F("fM1kaon","Invar m(kaon) from kink->mu+netrino decay",600,0.10, 0.7); // 23/8/2013
ad1036b2 152 fPtKink= new TH1F("fPtKink", "P_{T}Kaon Kink distribution, counts",44, gPt7K0);
153 fptKink= new TH1F("fptKink", "P_{T}Kaon Kink bution",300, 0.0,15.0);
828a6e0c 154 fAngMomK= new TH2F("fAngMomK","Decay angle vrs Mother Mom,K",100,0.0,5.0,80,0.,80.);
155 fAngMomPi= new TH2F("fAngMomPi","Decay angle vrs Mother Mom,Pi",100,0.0,5.0,80,0.,80.);
156 fAngMomKC= new TH2F("fAngMomKC","Decay angle vrs Mother Mom,K",100,0.0,5.0,80,0.,80.);
157 fMultESDK=new TH1F("fMultESDK", "charge multipliESD kaons",100, 0.0,100.0);
158 fMultMCK=new TH1F("fMultMCK", "charge multipli MC kaons",100, 0.0,100.0);
159 fSignPtNcl= new TH2F("fSignPtNcl","SignPt vrs Ncl,K",80,-4.,4.0,70,20.,160.);
160 fSignPtEta= new TH2F("fSignPtEta","SignPt vrs Eta,K",80,-4.0,4.0,30,-1.5,1.5);
161 fEtaNcl= new TH2F("fEtaNcl","Eta vrs nclust,K",30,-1.5,1.5, 70,20, 160);
828a6e0c 162 fSignPt= new TH1F("fSignPt","SignPt ,K",80,-4.0,4.0);
163 fChi2NclTPC= new TH2F("fChi2NclTPC","Chi2vrs nclust,K",100,0.,500., 70,20, 160);
164 fRatChi2Ncl= new TH1F("fRatChi2Ncl","Ratio chi2/nclusters in TPC,K",50,0.0,5.0);
ad1036b2 165 fRadiusNcl = new TH2F("fRadiusNcl","kink radius vrs Nclust,K",75,100.,250., 80,0, 160);
8e57b720 166 fTPCSgnlP = new TH2F("fTPCSgnlP","TPC signal de/dx Mom,K",1000,0.0,20.0,150,0.,300.);
167 fTPCSgnlPa= new TH2F("fTPCSgnlPa","TPC signal de/dx Mom,K",1000,0.0,20.,150, 0.,300.);
828a6e0c 168 fRpr = new TH1D("fRpr", "rad distribution PID pr",100,-10.0, 10.0);
169 fZpr = new TH1D("fZpr", "z distribution PID pr ",80,-20.,20.);
170 fdcatoVxXY = new TH1D("fdcatoVxXY", "dca distribution PID ",20,-1.,1.);
171 fnSigmToVx = new TH1D("fnSigmToVx", "dca distribution PID ",80,0.,8.);
172 fKinkMothDau= new TH2F("fKinkMothDau","TPC kink Moth Daugh ,K",50,0.0,2.5,50, 0., 2.5);
173 fZvXv= new TH2F("fZvXv","Xv-Zv main vtx",60,-0.5,0.5,60, -15., 15.0);
174 fZvYv= new TH2F("fZvYv","Yv-Zv main vtx",60,-0.5,0.5, 60, -15., 15.);
175 fXvYv= new TH2F("fXvYv","Xv-Yv main vtx", 60,-1.5,1.5, 60, -1.5, 1.5);
ad1036b2 176 fHistPtKaoP = new TH1F("fHistPtKaoP", "P_{T}KaonP distribution",300, 0.0,15.0);
ad1036b2 177 fHistPtKaoN = new TH1F("fHistPtKaoN", "P_{T}KaonN distribution",300, 0.0,15.0);
828a6e0c 178 frapiKESD=new TH1F("frapiKESD","rapid Kdistribution", 26,-1.3, 1.3);
179 flifetime= new TH1F("flifetime", "ct study of K-kinks",100,0.,1000.);
180 fradLK= new TH1F("fradLK", "Length of K generated",100,0.,1000.);
181 fradPtRpDt=new TH3F("fradPtRpDt","rad pt rap dat",28,100.,240., 20, 0., 5., 20, -1., 1. );
c1e83dc8 182 //fInvMuNuAll= new TH1F("fInvMuNuAll", " Inv Mass MuNu all kink",180,0.1,1.0);
183 fInvMuNuAll= new TH1F("fInvMuNuAll", " Inv Mass MuNu all kink",600,0.1,0.7); // 23/8/2013
828a6e0c 184 fQtInvM= new TH2F("fQtInvM", "Q_{T} Versus Inv MuNu ",80, 0., 0.80 , 100 , 0., 0.300);
185 fDCAkink = new TH1F("fDCAkink ", "DCA kink vetrex ",50, 0.0,1.0);
828a6e0c 186 fPosiKink= new TH2F("fPosiKink", "Y vrx kink Vrex ",100, -300.0,300.0,100, -300, 300.);
187 fPosiKinkK= new TH2F("fPosiKinkK", "Y vrx kink VrexK ",100, -300.0,300.0,100, -300, 300.);
188 fPosiKinKXZ= new TH2F("fPosiKinKXZ", "Y vrx kink VrexK ",100, -300.0,300.0,100, -300, 300.);
189 fPosiKinKYZ= new TH2F("fPosiKinKYZ", "Y vrx kink VrexK ",100, -300.0,300.0,100, -300, 300.);
190 fPosiKinKBg= new TH2F("fPosiKinKBg", "z vrx kink rad ",100, -300.0,300.0,100,100., 300.);
191 fQtMothP = new TH2F("fQtMothP", " Qt vrs Mother P", 100, 0., 5.0,100, 0.,0.300);
ad1036b2 192 fTPCSgnlPtpc = new TH2F("fTPCSgnlPtpc","TPC signal de/dx Mom TPC,K ",300,0.0,15.0,100, 0., 250. );
193 fTPCMomNSgnl = new TH2F("fTPCMomNsgnl","TPC signal de/dx Mom TPC,K ",300,0.0,15.0,20 , -10., 10.);
8e57b720 194 fMothKinkMomSgnl = new TH2F("fMothKinkMomSgnl","TPC signal de/dx Mom TPC,Kink ",100,0.0,250.0,100, 0., 250. );
ad1036b2 195 fNSigmTPC = new TH1F("fNSigmTPC","TPC Nsigma de/dx TPC,K ", 30 , -7.5, 7.5);
8e57b720 196 fTPCSgnlKinkDau = new TH2F("fTPCSgnlKinkDau","TPC signal de/dx Mom,K",500,0.0,10.0,100,0.,250.);
ad1036b2 197 fPtKinkPos= new TH1F("fPtKinkPos", "Pos P_{T}Kaon Kink distribution, counts",44, gPt7K0);
198 fPtKinkNeg= new TH1F("fPtKinkNeg", "Neg P_{T}Kaon Kink distribution, counts",44, gPt7K0);
199 fRadNclCln = new TH2F("fRadNclCln","kink radius vrs Nclust,K Clean ",75,100.,250., 80,0, 160);
8e57b720 200 fRatioCrossedRows = new TH1F("fRatioCrossedRows","Ratio crossed rows in TPC",20,0.0,1.0);
201 fRatioCrossedRowsKink = new TH1F("fRatioCrossedRowsKink","Ratio crossed rows in TPC for kinks",20,0.0,1.0);
202 fRadiusPt =new TH2F("fRadiusPt","radius vs pt ",80, 90.,250.,100, 0.,10. );
203 fRadiusPtcln =new TH2F("fRadiusPtcln","radius vs pt clean ",80, 90.,250.,100, 0.,10. );
c1e83dc8 204 //fInvMassMuNuPt =new TH2F("fInvMassMuNuPt","Invariant mass-munu vs pt ",180, 0.10, 1.00, 100, 0.0, 10.0 );
205 fInvMassMuNuPt =new TH2F("fInvMassMuNuPt","Invariant mass-munu vs pt ",600, 0.10, 0.7, 100, 0.0, 10.0 );// 23/8/2013
206 fInvMassMuNuPtAll =new TH2F("fInvMassMuNuPtAll","Invariant mass-munu vs pt ",600, 0.10, 0.7, 100, 0.0, 10.0 );// 23/8/2013
8e57b720 207 fPtCut1 = new TH1F("fPtCut1", "P_{T}Kaon distribution",300, 0.0,15.0);
208 fPtCut2 = new TH1F("fPtCut2", "P_{T}Kaon distribution",300, 0.0,15.0);
209 fPtCut3 = new TH1F("fPtCut3", "P_{T}Kaon distribution",300, 0.0,15.0);
210 fAngMomKKinks = new TH2F("fAngMomKKinks","Decay angle vrs Mother Mom,Kinks",300,0.0,15.0,100,0.,100.);
828a6e0c 211
212 fListOfHistos=new TList();
213
214 fListOfHistos->Add(fHistPtESD);
215 fListOfHistos->Add(fHistPt);
216 fListOfHistos->Add(fHistQtAll);
217 fListOfHistos->Add(fHistQt1);
218 fListOfHistos->Add(fHistQt2);
219 fListOfHistos->Add(fHistPtKaon);
220 fListOfHistos->Add(fHistPtKPDG);
221 fListOfHistos->Add(fHistEta);
222 fListOfHistos->Add(fHistEtaK);
223 fListOfHistos->Add(fptKMC);
224 fListOfHistos->Add(fMultiplMC);
225 fListOfHistos->Add(fESDMult);
226 fListOfHistos->Add(fgenpt);
227 fListOfHistos->Add(frad);
228 fListOfHistos->Add(fKinkKaon);
229 fListOfHistos->Add(fKinKRbn);
230 fListOfHistos->Add(fKinkKaonBg);
231 fListOfHistos->Add(fM1kaon);
828a6e0c 232 fListOfHistos->Add(fPtKink);
233 fListOfHistos->Add(fptKink);
828a6e0c 234 fListOfHistos->Add(fAngMomK);
235 fListOfHistos->Add(fAngMomPi);
236 fListOfHistos->Add(fAngMomKC);
237 fListOfHistos->Add(fMultESDK);
238 fListOfHistos->Add(fMultMCK);
239 fListOfHistos->Add(fSignPtNcl);
240 fListOfHistos->Add(fSignPtEta);
241 fListOfHistos->Add(fEtaNcl);
242 fListOfHistos->Add(fSignPt);
243 fListOfHistos->Add(fChi2NclTPC);
244 fListOfHistos->Add(fRatChi2Ncl);
245 fListOfHistos->Add(fRadiusNcl);
246 fListOfHistos->Add(fTPCSgnlP);
247 fListOfHistos->Add(fTPCSgnlPa);
248 fListOfHistos->Add(fRpr);
249 fListOfHistos->Add(fZpr);
250 fListOfHistos->Add(fdcatoVxXY);
251 fListOfHistos->Add(fnSigmToVx);
252 fListOfHistos->Add(fKinkMothDau);
253 fListOfHistos->Add(fZvXv);
254 fListOfHistos->Add(fZvYv);
255 fListOfHistos->Add(fXvYv);
828a6e0c 256 fListOfHistos->Add(fHistPtKaoP);
257 fListOfHistos->Add(fHistPtKaoN);
258 fListOfHistos->Add(frapiKESD);
259 fListOfHistos->Add(flifetime);
260 fListOfHistos->Add(fradLK);
261 fListOfHistos->Add(fradPtRpDt);
262 fListOfHistos->Add(fInvMuNuAll);
263 fListOfHistos->Add(fQtInvM);
264 fListOfHistos->Add(fDCAkink);
265 fListOfHistos->Add(fPosiKink);
266 fListOfHistos->Add(fPosiKinkK);
267 fListOfHistos->Add(fPosiKinKXZ);
268 fListOfHistos->Add(fPosiKinKYZ);
269 fListOfHistos->Add(fPosiKinKBg);
270 fListOfHistos->Add(fQtMothP);
ad1036b2 271 fListOfHistos->Add(fTPCSgnlPtpc);
272 fListOfHistos->Add(fTPCMomNSgnl);
273 fListOfHistos->Add(fMothKinkMomSgnl);
274 fListOfHistos->Add(fNSigmTPC);
275 fListOfHistos->Add(fTPCSgnlKinkDau);
276 fListOfHistos->Add(fPtKinkPos);
277 fListOfHistos->Add(fPtKinkNeg);
278 fListOfHistos->Add(fRadNclCln);
8e57b720 279 fListOfHistos->Add(fRatioCrossedRows);
280 fListOfHistos->Add(fRatioCrossedRowsKink);
281 fListOfHistos->Add(fRadiusPt);
282 fListOfHistos->Add(fRadiusPtcln);
283 fListOfHistos->Add(fInvMassMuNuPt);
c1e83dc8 284 fListOfHistos->Add(fInvMassMuNuPtAll);
8e57b720 285 fListOfHistos->Add(fPtCut1);
286 fListOfHistos->Add(fPtCut2);
287 fListOfHistos->Add(fPtCut3);
288 fListOfHistos->Add(fAngMomKKinks);
828a6e0c 289
290}
828a6e0c 291//________________________________________________________________________
292void AliAnalysisKinkESDat::UserExec(Option_t *)
293{
294 // Main loop
295 // Called for each event
296
297 // Process MC truth, therefore we receive the AliAnalysisManager and ask it for the AliMCEventHandler
298 // This handler can return the current MC event
299
300 AliVEvent *event = InputEvent();
301 if (!event) {
302 Printf("ERROR: Could not retrieve event");
303 return;
304 }
305
306 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
307 if (!esd) {
308 Printf("ERROR: Could not retrieve esd");
309 return;
310 }
ad1036b2 311// Number ESD tracks
312 Int_t nESDTracks = esd->GetNumberOfTracks();
313 fMultiplMC->Fill(nESDTracks);
314//
828a6e0c 315//==================check of Physics selection?
316 Bool_t isSelected =
317((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()&AliVEvent::kMB;
318
319 if ( isSelected ==kFALSE) return; // 24/6/11 apo MF
320//
828a6e0c 321 fMultMCK->Fill(nESDTracks);
ad1036b2 322//
323//===============Marek's multiplicity
828a6e0c 324 Float_t refmultiplicity=fCutsMul->CountAcceptedTracks(esd);
325 if(fLowMulcut>-1)
326 {
327 if(refmultiplicity<fLowMulcut)
328 return;
329 }
330 if(fUpMulcut>-1)
331 {
332 if(refmultiplicity>fUpMulcut)
333 return;
334 }
335
336
337
338 fMultESDK->Fill(refmultiplicity);
339
340//
828a6e0c 341
342
343
344 const AliESDVertex *vertex=GetEventVertex(esd); // 22/8
ad1036b2 345 if(!vertex) return;
828a6e0c 346//
347 Double_t vpos[3];
348 vertex->GetXYZ(vpos);
349 fZpr->Fill(vpos[2]);
ad1036b2 350 if (TMath::Abs( vpos[2] ) > 10. ) return;
828a6e0c 351
352
353
354 Double_t vtrack[3], ptrack[3];
355
356
357 Int_t nESDTracK = 0;
5220a87d 358// Int_t nESDTrKink = 0;
828a6e0c 359
360 Int_t nGoodTracks = esd->GetNumberOfTracks();
361 fESDMult->Fill(nGoodTracks);
362
ad1036b2 363 Double_t nsigmall = 100.0;
828a6e0c 364 Double_t nsigma = 100.0;
ad1036b2 365 Double_t nsigmaPion =-100.0;
5220a87d 366 // Double_t nsigmaDau =-100.0;
8e57b720 367 Double_t dEdxKinkDau =0.0;
368 Double_t KinkDauCl =0.0;
fcd98e58 369// apo Eftihi
370 if(!fPIDResponse) {
371 AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
372 AliInputEventHandler* inputHandler =
373(AliInputEventHandler*)(man->GetInputEventHandler());
374 fPIDResponse = inputHandler->GetPIDResponse();
375 }
376
ad1036b2 377// loop on kink daughters
378 for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
379
380 AliESDtrack* trackD = esd->GetTrack(iTrack);
381 if (!trackD) {
382 Printf("ERROR: Could not receive track %d", iTrack);
383 continue;
384 }
8e57b720 385//
ad1036b2 386 Int_t indexKinkDau=trackD->GetKinkIndex(0);
387// daughter kink
388// AliESDkink *kink=esd->GetKink(TMath::Abs(indexKinkDau)-1);
8e57b720 389 if ( indexKinkDau > 0 ) {
390 Int_t labelD = trackD->GetLabel();
391 labelD = TMath::Abs(labelD);
392 nsigmaPion = (fPIDResponse->NumberOfSigmasTPC(trackD , AliPID::kPion));// 26/10 eftihis
393 dEdxKinkDau = (trackD->GetTPCsignal() ) ; // daughter kink dEdx
394 KinkDauCl=(trackD->GetTPCclusters(0) ) ;
395 }
fcd98e58 396//if((indexKinkDau >0)&& (nsigmaPion>1.2)) fTPCSgnlKinkDau->Fill(trackD->P(), (trackD->GetTPCsignal() ) ) ; // daughter kink
8e57b720 397// Ayto mexri 26/11/2012 if(indexKinkDau >0) fTPCSgnlKinkDau->Fill(trackD->P(), (trackD->GetTPCsignal() ) ) ; // daughter kink
ad1036b2 398 }
fcd98e58 399
828a6e0c 400// track loop
ad1036b2 401//
828a6e0c 402 for (Int_t iTracks = 0; iTracks < esd->GetNumberOfTracks(); iTracks++) {
403
404 AliESDtrack* track = esd->GetTrack(iTracks);
405 if (!track) {
406 Printf("ERROR: Could not receive track %d", iTracks);
407 continue;
408 }
409
410 fHistPt->Fill(track->Pt());
411
412
413 // sigmas
fcd98e58 414 nsigmall = (fPIDResponse->NumberOfSigmasTPC(track, AliPID::kKaon));
ad1036b2 415 // nsigmaPion= (fESDpid->NumberOfSigmasTPC(track,AliPID::kPion));
fcd98e58 416 nsigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kKaon));
828a6e0c 417
8e57b720 418//=======================new
419 Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
420 Float_t ratioCrossedRowsOverFindableClustersTPC = 1.0;
421 if (track->GetTPCNclsF()>0) {
422 ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC / track->GetTPCNclsF();
423 fRatioCrossedRows->Fill(ratioCrossedRowsOverFindableClustersTPC);
424 }
425//_______
828a6e0c 426
fcd98e58 427 Int_t indexKinkPos=track->GetKinkIndex(0); // kink index
828a6e0c 428
429 Int_t tpcNCl = track->GetTPCclusters(0);
430 Double_t tpcSign = track->GetSign();
431
432 Int_t label = track->GetLabel();
433 label = TMath::Abs(label);
434
435
436 UInt_t status=track->GetStatus();
437
438 if((status&AliESDtrack::kITSrefit)==0) continue;
439 if((status&AliESDtrack::kTPCrefit)==0) continue;
ad1036b2 440 if((track->GetTPCchi2()/track->GetTPCclusters(0))>4.0) continue;
828a6e0c 441
442 Double_t extCovPos[15];
443 track->GetExternalCovariance(extCovPos);
828a6e0c 444
445
446 track->GetXYZ(vtrack);
447 fXvYv->Fill(vtrack[0],vtrack[1]);
448 fZvYv->Fill(vtrack[0],vtrack[2]);
449 fZvXv->Fill(vtrack[1],vtrack[2]);
450
451// track momentum, rapidity calculation
452 track->GetPxPyPz(ptrack);
453
454 TVector3 trackMom(ptrack[0],ptrack[1],ptrack[2]);
455
ad1036b2 456// K-rapidity calcualtion
828a6e0c 457 Double_t etracK= TMath::Sqrt(trackMom.Mag()*trackMom.Mag() + 0.493677 *0.493677 );
458 Double_t rapiditK = 0.5 * (TMath::Log( (etracK + ptrack[2] ) / ( etracK - ptrack[2]) )) ;
459
460 Double_t trackEta=trackMom.Eta();
828a6e0c 461 Double_t trackPt = track->Pt();
462
463
464
465 Float_t bpos[2];
466 Float_t bCovpos[3];
467 track->GetImpactParameters(bpos,bCovpos);
468
469 if (bCovpos[0]<=0 || bCovpos[2]<=0) {
470 Printf("Estimated b resolution lower or equal zero!");
471 bCovpos[0]=0; bCovpos[2]=0;
472 }
473
474 Float_t dcaToVertexXYpos = bpos[0];
475 Float_t dcaToVertexZpos = bpos[1];
476
477 fRpr->Fill(dcaToVertexZpos);
478
0dc5d777 479// 14/2/13 /================/ if((TMath::Abs(dcaToVertexXYpos)>0.3)||(TMath::Abs(dcaToVertexZpos)>2.5)) nESDTrKink++; // count of second 23Jul11
8e57b720 480 //if((TMath::Abs(dcaToVertexXYpos)>0.3)||(TMath::Abs(dcaToVertexZpos)>2.5))
481// continue; // allagi 23Jul11
482
0dc5d777 483 if (!fMaxDCAtoVtxCut->AcceptTrack(track)) continue;
8e57b720 484
485 fdcatoVxXY->Fill(dcaToVertexXYpos);
828a6e0c 486//
828a6e0c 487
488// track Mult. after selection
489 nESDTracK++;
490 //
491//=========================================
492 fHistPtESD->Fill(track->Pt());
493
494 // Add Kink analysis =============================
495
ad1036b2 496// daughter kink
497//if(indexKinkPos >0)fTPCSgnlKinkDau->Fill(track->P(), (track->GetTPCsignal() ) ) ; // daughter kink
498
828a6e0c 499// loop on kinks
500 if(indexKinkPos<0){ ////mother kink
828a6e0c 501
fcd98e58 502 fptKMC ->Fill( track->Pt() ); // Pt from tracks , all kinks
503
8e57b720 504 fRatioCrossedRowsKink->Fill(ratioCrossedRowsOverFindableClustersTPC);
828a6e0c 505 // select kink class
506
507 AliESDkink *kink=esd->GetKink(TMath::Abs(indexKinkPos)-1);
508//
509 // DCA kink
510 Double_t Dist2 = kink->GetDistance();
ad1036b2 511 // fDCAkink->Fill( Dist2 );
512 // if (Dist2 > 0.2) continue; // systematics 11/8/11
828a6e0c 513//
514
ad1036b2 515// TPC mother momentum
828a6e0c 516
517 const TVector3 vposKink(kink->GetPosition());
518 fPosiKink ->Fill( vposKink[0], vposKink[1] );
828a6e0c 519 Double_t dzKink=vpos[2]-vposKink[2];
828a6e0c 520//
8e57b720 521// lifitime
828a6e0c 522 Double_t tanLamda = track->GetTgl(); // 25/6/2010
523
524 Double_t lifeKink= (TMath::Abs( dzKink ))*( TMath::Sqrt(1.+ tanLamda*tanLamda) ) / (TMath::Abs( tanLamda)) ;
8e57b720 525//
828a6e0c 526 const TVector3 motherMfromKink(kink->GetMotherP());
527 const TVector3 daughterMKink(kink->GetDaughterP());
528
529 Float_t qT=kink->GetQt();
530 Float_t motherPt=motherMfromKink.Pt();
ad1036b2 531// Kink mother momentum
5220a87d 532// Double_t trMomTPCKink=motherMfromKink.Mag();
ad1036b2 533// TPC mother momentun
534 Double_t trMomTPC=track->GetTPCmomentum();
8e57b720 535 // fTPCSgnlKinkDau->Fill( daughterMKink.Mag(), dEdxKinkDau ) ; // daughter kink
536 //
828a6e0c 537 fHistQtAll->Fill(qT) ; // Qt distr
538
539 fptKink->Fill(motherMfromKink.Pt()); /// pt from kink
540
541 Float_t kinkAngle=TMath::RadToDeg()*kink->GetAngle(2);
542
fcd98e58 543// rapiditya nd pt selection
828a6e0c 544 if( (TMath::Abs(rapiditK )) > 0.7 ) continue;
ad1036b2 545 if ( (track->Pt())<.200)continue; // new Feb 2012
828a6e0c 546
547 fQtMothP->Fill( track->P(), qT);
548
5220a87d 549 if ( qT> fLowQt ) fHistQt1 ->Fill(qT) ; // Qt distr
828a6e0c 550
551
552
fcd98e58 553 fHistEta->Fill(trackEta) ; // Eta distr
554 fHistQt2->Fill(qT); //
8e57b720 555 fKinkKaonBg->Fill(motherPt);
828a6e0c 556
557// maximum decay angle at a given mother momentum
558 //Double_t maxDecAngKmu=f1->Eval(motherMfromKink.Mag(),0.,0.,0.);
559 Double_t maxDecAngKmu=f1->Eval(track->P() ,0.,0.,0.);
560 Double_t maxDecAngpimu=f2->Eval( track->P(), 0.,0.,0.);
fcd98e58 561
562// fake kinks are removed
828a6e0c 563 if( (kinkAngle<2.) ) continue;
564
565
566 // BG ?????==============
567 if ( TMath::Abs(vposKink[2]) > 225. ) continue ;
568 if ( TMath::Abs(vposKink[2]) < 0.5 ) continue ;
828a6e0c 569//
8e57b720 570 fPtCut1 ->Fill(trackPt );
fcd98e58 571 fAngMomPi->Fill( track->P(), kinkAngle);
828a6e0c 572//
573// invariant mass of mother track decaying to mu
574 Float_t energyDaughterMu=TMath::Sqrt(daughterMKink.Mag()*daughterMKink.Mag()+0.105658*0.105658);
575 Float_t p1XM= motherMfromKink.Px();
576 Float_t p1YM= motherMfromKink.Py();
577 Float_t p1ZM= motherMfromKink.Pz();
578 Float_t p2XM= daughterMKink.Px();
579 Float_t p2YM= daughterMKink.Py();
580 Float_t p2ZM= daughterMKink.Pz();
581 Float_t p3Daughter=TMath::Sqrt(((p1XM-p2XM)*(p1XM-p2XM))+((p1YM-p2YM)*(p1YM-p2YM))+((p1ZM-p2ZM)*(p1ZM-p2ZM)));
582 Double_t invariantMassKmu= TMath::Sqrt((energyDaughterMu+p3Daughter)*(energyDaughterMu+p3Daughter)-motherMfromKink.Mag()*motherMfromKink.Mag());
583 fQtInvM -> Fill ( invariantMassKmu, qT);
584 fInvMuNuAll->Fill(invariantMassKmu);
c1e83dc8 585 fInvMassMuNuPtAll ->Fill(invariantMassKmu, trackPt);
fcd98e58 586//
8e57b720 587 fRadiusPt->Fill( kink->GetR(), trackPt); //
fcd98e58 588 // radius and Minv selection
b2fb60b8 589 //if( ( kink->GetR()> 120 ) && ( kink->GetR() < 210 ) ) {
590 if( ( kink->GetR()> fKinkRadLow ) && ( kink->GetR() <fKinkRadUp ) ) {
fcd98e58 591 // for systematics if( ( kink->GetR()> 130 ) && ( kink->GetR() < 200 ) ) {
5220a87d 592 if (qT>fLowQt ) fAngMomKC->Fill(track->P(), kinkAngle);
593 if ( qT> fLowQt ) fM1kaon->Fill(invariantMassKmu);
594 if ( qT > fLowQt)
828a6e0c 595 fRadiusNcl->Fill( (kink->GetR()) ,(track->GetTPCclusters(0) ) ) ;
596 }
fcd98e58 597// tails cleaning
5220a87d 598 if( ( tpcNCl<fLowCluster) ) continue; // test 27 feb 2012 ,, OK
599 // edw iatn !!! if( ( tpcNCl<50 ) ) continue; // test 15 March 13,, OK
fcd98e58 600// cleaning BG in tails
ad1036b2 601 Int_t tpcNClHigh = -51.67+ (11./12.) *( kink->GetR() ) ;
602 if ( tpcNCl > tpcNClHigh) continue;
603
ad1036b2 604 Int_t tpcNClMin = -85.5 + (65./95.) *( kink->GetR() ) ;
ad1036b2 605 if ( tpcNCl < tpcNClMin ) continue;
828a6e0c 606
8e57b720 607 // back, 20/1/2013 if (ratioCrossedRowsOverFindableClustersTPC< 0.5) continue;// test 14/1/2013
828a6e0c 608//
fcd98e58 609 fHistPtKPDG->Fill(track->Pt()); // ALL K-candidates until now
ad1036b2 610 // if((kinkAngle>maxDecAngpimu)&&(qT>0.12)&&(qT<0.30)&&((kink->GetR()>=120.)&&(kink->GetR()<=210.))&&(TMath::Abs(rapiditK)<0.7)&&(invariantMassKmu<0.6)){
b2fb60b8 611 //if((kinkAngle>maxDecAngpimu)&&(qT>0.12)&&(qT<0.30)&&((kink->GetR()>=120.)&&(kink->GetR()<=210.))&&(TMath::Abs(rapiditK)<0.7)&&(invariantMassKmu<0.8)){
5220a87d 612 //if((kinkAngle>maxDecAngpimu)&&(qT>0.04)&&(qT<0.30)&&((kink->GetR()>= fKinkRadLow )&&(kink->GetR()<= fKinkRadUp ))&&(TMath::Abs(rapiditK)<0.7)&&(invariantMassKmu<0.8)){
613 if((kinkAngle>maxDecAngpimu)&&(qT>fLowQt)&&(qT<0.30)&&((kink->GetR()>= fKinkRadLow )&&(kink->GetR()<= fKinkRadUp ))&&(TMath::Abs(rapiditK)<0.7)&&(invariantMassKmu<0.8)){
fcd98e58 614 // systematics if((kinkAngle>maxDecAngpimu)&&(qT>0.12)&&(qT<0.30)&&((kink->GetR()>=130.)&&(kink->GetR()<=200.))&&(TMath::Abs(rapiditK)<0.7)&&(invariantMassKmu<0.8)){
615//
8e57b720 616 fAngMomKKinks->Fill(track->P(), kinkAngle);
617 fPtCut2 ->Fill(trackPt );
fcd98e58 618// maximum angles selection with some error cut
828a6e0c 619 if( (kinkAngle<maxDecAngpimu*1.2) ) continue;
620 if ( (kinkAngle>maxDecAngKmu*.98) && ( track->P() >1.2 )) continue; ///5/5/2010
621
8e57b720 622 fPtCut3 ->Fill(trackPt );
fcd98e58 623// here the kaons selected by the decay features
ad1036b2 624 fTPCSgnlPa->Fill( track->GetInnerParam()->GetP() ,(track->GetTPCsignal() ) ) ;
fcd98e58 625//
8e57b720 626 // NO dEdx cut test 9/2/13 if ( nsigma > 3.5) continue;
5220a87d 627 if ( nsigma > 3.5) continue;
ad1036b2 628//
fcd98e58 629// next plots for the identified kaons by the kink analysis
630
631 fHistPtKaon->Fill(track->Pt()); //
632 if(tpcSign >0.) fHistPtKaoP->Fill( track->Pt() ) ; //
633 if ( tpcSign <0.) fHistPtKaoN->Fill( track->Pt() ) ; //
ad1036b2 634 fTPCSgnlP->Fill(track->GetInnerParam()->GetP(), (track->GetTPCsignal() ) ) ;
635 fRadNclCln->Fill( (kink->GetR()) ,(track->GetTPCclusters(0) ) ) ;
8e57b720 636 fRadiusPtcln->Fill( kink->GetR(), trackPt); //
637 fInvMassMuNuPt ->Fill(invariantMassKmu, trackPt);
ad1036b2 638
8e57b720 639 // fTPCSgnlPtpc->Fill(trMomTPC , (track->GetTPCsignal() ) ) ;
640 //fMothKinkMomSgnl ->Fill(trMomTPCKink , (track->GetTPCsignal() ) ) ;
641 fMothKinkMomSgnl ->Fill( dEdxKinkDau , (track->GetTPCsignal() ) ) ;
642//
643 fTPCSgnlKinkDau->Fill( daughterMKink.Mag(), dEdxKinkDau ) ; // daughter kink
ad1036b2 644//
645 fTPCMomNSgnl->Fill(trMomTPC ,nsigmall );
646 fNSigmTPC ->Fill(nsigmall );
647//
828a6e0c 648 frad->Fill(kink->GetR()); // kink
649 fradLK->Fill(lifeKink ); // kink
650 fHistEtaK->Fill(trackEta);
651 frapiKESD ->Fill(rapiditK); // rapidityof kaons
652 fPosiKinKBg->Fill( vposKink[2], kink->GetR() );
653
654 Float_t signPt= tpcSign*trackPt;
655 fSignPtNcl->Fill( signPt , tpcNCl ); /// 28/4/2010
656 fSignPtEta->Fill( signPt , rapiditK );
657 fEtaNcl->Fill( rapiditK, tpcNCl );
658 fSignPt->Fill( signPt );
659 fChi2NclTPC->Fill( (track->GetTPCchi2() ) , tpcNCl );
660 fRatChi2Ncl-> Fill ( (track->GetTPCchi2()/track->GetTPCclusters(0) )) ;
8e57b720 661 // fdcatoVxXY->Fill(dcaToVertexXYpos);
662 //if( dEdxKinkDau> 1.5* (track->GetTPCsignal() ) ) fKinkMothDau->Fill(track->P(),daughterMKink.Mag());
663 // if((dEdxKinkDau> 80. ) && (dEdxKinkDau > 4.*nsigmaPion) ) fKinkMothDau->Fill(track->P(),daughterMKink.Mag());
664 if (nsigmaPion> 3. ) fTPCSgnlPtpc->Fill( daughterMKink.Mag(), dEdxKinkDau ) ;
665 //if (TMath::Abs(dEdxKinkDau -(track->GetTPCsignal() )> 10. )) fTPCSgnlPtpc->Fill( daughterMKink.Mag(), dEdxKinkDau ) ;
828a6e0c 666 flifetime->Fill(( lifeKink*.493667 ) /track->P() ) ;
667 fKinkKaon->Fill(track->Pt());
ad1036b2 668 fDCAkink->Fill( Dist2 );
669
670 fPtKink->Fill(track->Pt()); /// K0 bins
671 if(tpcSign >0.) fPtKinkPos ->Fill( track->Pt() ) ; //K-plus bins K0
672 if(tpcSign <0.) fPtKinkNeg ->Fill( track->Pt() ) ; //K-minus bins K0
673 fKinKRbn->Fill(track->Pt()); // TOF
828a6e0c 674
828a6e0c 675 fradPtRpDt->Fill( kink->GetR(), 1./track->Pt(), rapiditK);
676 fAngMomK->Fill( track->P(), kinkAngle);
677 fPosiKinkK->Fill( vposKink[0], vposKink[1] );
678 fPosiKinKXZ->Fill( vposKink[2], vposKink[0] );
679 fPosiKinKYZ->Fill( vposKink[2], vposKink[1] );
680
681 } // kink selection
682
683
684 } //End Kink Information
685
686
687 } //track loop
ad1036b2 688// } //daughter loop
828a6e0c 689
690 // fMultiplMC->Fill(nESDTracK );
691
692 PostData(1, fListOfHistos);
693
694}
695
696//________________________________________________________________________
697void AliAnalysisKinkESDat::Terminate(Option_t *)
698{
699 // Draw result to the screen
700 // Called once at the end of the query
701
702}
703
704//____________________________________________________________________//
705
706
707const AliESDVertex* AliAnalysisKinkESDat::GetEventVertex(const AliESDEvent* esd) const
708 // Get the vertex from the ESD and returns it if the vertex is valid
709
710{
711 // Get the vertex
712
828a6e0c 713 const AliESDVertex* vertex = esd->GetPrimaryVertexTracks();
714
828a6e0c 715 if((vertex->GetStatus()==kTRUE)) return vertex;
716 else
717 {
718 vertex = esd->GetPrimaryVertexSPD();
719 if((vertex->GetStatus()==kTRUE)&&(vertex->GetNContributors()>0)) return vertex;
828a6e0c 720 else
721 return 0;
722 }
ad1036b2 723}