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