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