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