task for ESD kink analysis (Martha)
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / Kinks / AliAnalysisKinkESDat.cxx
1 /**************************************************************************
2  * of the Greek group at Physics Department of Athens University
3  * Paraskevi Ganoti, Anastasia Belogianni and Filimon Roukoutakis 
4  * Contributors are mentioned in the code where appropriate.              *
5  *                                                                        *
6  * Permission to use, copy, modify and distribute this software and its   *
7  * documentation strictly for non-commercial purposes is hereby granted   *
8  * without fee, provided that the above copyright notice appears in all   *
9  * copies and that both the copyright notice and this permission notice   *
10  * appear in the supporting documentation. The authors make no claims     *
11  * about the suitability of this software for any purpose. It is          *
12  * provided "as is" without express or implied warranty.                  *
13  **************************************************************************/
14
15 //-----------------------------------------------------------------
16 //                 AliAnalysisKinkESDat class
17 //       Example of an analysis task for kink topology study
18 //      Kaons from kink topology are 'identified' in this code
19 //-----------------------------------------------------------------
20
21 #include "TChain.h"
22 #include "TTree.h"
23 #include "TH1F.h"
24 #include "TH2F.h"
25 #include "TH3F.h"
26 #include "TH1D.h"
27 #include "TH2D.h"
28 #include "TParticle.h"
29 #include <TVector3.h>
30 #include "TF1.h"
31
32 #include "AliAnalysisTask.h"
33 #include "AliAnalysisManager.h"
34
35 #include "AliVEvent.h"
36 #include "AliESDEvent.h"
37 #include "AliMCEvent.h"
38 #include "AliAnalysisKinkESDat.h"
39 #include "AliStack.h"
40 #include "AliESDpid.h"
41 #include "AliPID.h"
42 #include "AliESDkink.h"
43 #include "AliESDtrack.h"
44 #include "AliPhysicsSelectionTask.h"
45 #include "AliInputEventHandler.h"
46 #include "AliAnalysisManager.h"
47 #include "AliVEvent.h"
48 //#include "AddTaskTender.h"
49 //  #include "AliTPCpidESD.h"
50 #include "AliESDtrackCuts.h"
51 ClassImp(AliAnalysisKinkESDat)
52
53
54 //________________________________________________________________________
55 AliAnalysisKinkESDat::AliAnalysisKinkESDat(const char *name) 
56   : AliAnalysisTaskSE(name), fHistPtESD(0),fHistPt(0),fHistQtAll(0),fHistQt1(0),fHistQt2(0)
57   , fHistPtKaon(0),fHistPtKPDG(0),fHistEta(0),fHistEtaK(0),fptKMC(0),fMultiplMC(0),fESDMult(0),fgenpt(0),frad(0),
58   fKinkKaon(0),fKinKRbn(0), fKinkKaonBg(0), fM1kaon(0),  fgenPtEtR(0),fPtKink(0),  fptKink(0),
59    fcodeH(0), fdcodeH(0), fAngMomK(0),fAngMomPi(0), fAngMomKC(0),  fMultESDK(0), fMultMCK(0),
60  fSignPtNcl(0), fSignPtEta(0), fEtaNcl(0), fSignPt(0), fChi2NclTPC(0), fRatChi2Ncl(0), fRadiusNcl(0), fTPCSgnlP(0),
61    fTPCSgnlPa(0), fRpr(0),fZpr(0), fdcatoVxXY(0), fnSigmToVx(0), fKinkMothDau(0),
62  fZvXv(0),fZvYv(0), fXvYv(0), fPtPrKink(0), fHistPtKaoP(0), fHistPtKaoN(0),frapiKESD(0), flifetime(), fradLK(0),
63     fradPtRpDt(0), fInvMuNuAll(0), fQtInvM(0), 
64          fDCAkink(0), fPosiKink(0),  fPosiKinkK(0),fPosiKinKXZ(0), fPosiKinKYZ(0), fPosiKinKBg(0), fQtMothP(0),
65  f1(0), f2(0),
66       fListOfHistos(0),fLowMulcut(-1),fUpMulcut(-1),fCutsMul(0)
67
68 {
69   // Constructor
70
71   // Define input and output slots here
72   // Input slot #0 works with a TChain
73  // DefineInput(0, TChain::Class());
74 //----------------------Marek multiplicity bins 
75  fCutsMul=new AliESDtrackCuts("Mul","Mul");
76         fCutsMul->SetMinNClustersTPC(70);
77         fCutsMul->SetMaxChi2PerClusterTPC(4);
78         fCutsMul->SetAcceptKinkDaughters(kFALSE);
79         fCutsMul->SetRequireTPCRefit(kTRUE);
80         // ITS
81         fCutsMul->SetRequireITSRefit(kTRUE);
82         fCutsMul->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
83                                                 AliESDtrackCuts::kAny);
84         fCutsMul->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
85
86         fCutsMul->SetMaxDCAToVertexZ(2);
87         fCutsMul->SetDCAToVertex2D(kFALSE);
88         fCutsMul->SetRequireSigmaToVertex(kFALSE);
89
90         fCutsMul->SetEtaRange(-0.8,+0.8);
91         fCutsMul->SetPtRange(0.15, 1e10);
92
93
94   DefineOutput(1, TList::Class());
95 }
96
97 //________________________________________________________________________
98 void AliAnalysisKinkESDat::UserCreateOutputObjects() 
99 {
100   // Create histograms
101   // Called once
102   
103    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);
104    f1->SetParameter(0,0.493677);
105    f1->SetParameter(1,0.9127037);
106    f1->SetParameter(2,TMath::Pi());
107
108
109    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);
110    f2->SetParameter(0,0.13957018);
111    f2->SetParameter(1,0.2731374);
112    f2->SetParameter(2,TMath::Pi());
113    //Open file  1= CAF 
114     //OpenFile(1); 
115  //  Double_t gPt[31] = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8,0.9,1.0,
116    //                     1.1, 1.2, 1.3, 1.4, 1.5, 1.6,1.7,1.8,1.9,  2.0,
117      //                    2.2, 2.4, 2.6, 2.8, 3.0, 3.3, 3.6,3.9, 4.2, 4.5, 4.8};
118    Double_t gPt7[43] = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8,0.9,1.0,
119                         1.1, 1.2, 1.3, 1.4, 1.5, 1.6,1.7,1.8,1.9,  2.0,
120                          2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7,2.8, 2.9, 3.0, 
121                          3.2, 3.4, 3.6, 3.8, 4.0, 4.4, 4.8,5.2, 5.6, 6.0,  7.0, 8.0,10.0 };
122
123   fHistPtESD = new TH1F("fHistPtESD", "P_{T} distribution",50, 0.0,5.0);
124   fHistPtESD->GetXaxis()->SetTitle("P_{T} (GeV/c)");
125   fHistPtESD->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
126   fHistPtESD->SetMarkerStyle(kFullCircle);
127   fHistPt = new TH1F("fHistPt", "P_{T} distribution",50, 0.0,5.0); 
128   fHistQtAll = new TH1F("fHistQtAll", "Q_{T} distr All Kinks ",100, 0.0,.300); 
129   fHistQt1= new TH1F("fHistQt1", "Q_{T} distribution",100, 0.0,.300); 
130   fHistQt2= new TH1F("fHistQt2", "Q_{T} distribution",100, 0.0,.300); 
131   //fHistPtKaon = new TH1F("fHistPtKaon", "P_{T}Kaon distribution",50, 0.0,5.0); 
132   fHistPtKaon = new TH1F("fHistPtKaon", "P_{T}Kaon distribution",100, 0.0,10.0); 
133   fHistPtKPDG = new TH1F("fHistPtKPDG", "P_{T}Kaon distribution",50, 0.0,5.0); 
134   fHistEta= new TH1F("fHistEta", "Eta distribution", 26,-1.3, 1.3); 
135   fHistEtaK= new TH1F("fHistEtaK", "EtaK distribution", 26,-1.3, 1.3); 
136   fptKMC= new TH1F("fptKMC", "P_{T}Kaon generated",100, 0.0,10.0); 
137   fMultiplMC= new TH1F("fMultiplMC", "charge multiplicity MC",100, 0.0,300.0); 
138   fESDMult= new TH1F("fESDMult", "charge multipliESD",100, 0.0,300.0); 
139   fgenpt= new TH1F("fgenpt", "genpt   K distribution",50, 0.0,5.0); 
140    //frad= new TH1F("frad", "radius  K generated",100, 50., 250.0);
141    frad= new TH1F("frad", "radius  K generated",100, 0.,1000.0);
142   // fKinkKaon= new TH1F("fKinkKaon", "P_{T}Kaon kinks identi",50, 0.0,5.0); 
143   fKinkKaon= new TH1F("fKinkKaon", "P_{T}Kaon kinks identi",100, 0.0,10.0); 
144   fKinKRbn= new TH1F("fKinKRbn", "p_{t}Kaon kinks identi[GeV/c],Entries",42,gPt7); 
145   fKinkKaonBg= new TH1F("fKinkKaonBg", "P_{T}Kaon kinks backgr",50, 0.0,5.0); 
146   fM1kaon= new TH1F("fM1kaon","Invar m(kaon) from kink->mu+netrino decay",80,0.0, 0.8); 
147   fgenPtEtR= new TH1F("fgenPtEtR", "P_{T}Kaon distribution",50, 0.0,5.0); 
148   fPtKink= new TH1F("fPtKink", "P_{T}Kaon Kink  bution",50, 0.0,5.0); 
149   fptKink= new TH1F("fptKink", "P_{T}Kaon Kink  bution",50, 0.0,5.0); 
150   fcodeH   = new TH2F("fcodeH", "code vrs dcode dist. kinks,K",100,0.,2500.,100,0.,2500.);
151   fdcodeH = new TH2F("fdcodeH", "code vrs dcode dist. kinks,K",100,0.,2500.,100,0.,2500.);
152   fAngMomK= new TH2F("fAngMomK","Decay angle vrs Mother Mom,K",100,0.0,5.0,80,0.,80.);
153   fAngMomPi= new TH2F("fAngMomPi","Decay angle vrs Mother Mom,Pi",100,0.0,5.0,80,0.,80.);
154   fAngMomKC= new TH2F("fAngMomKC","Decay angle vrs Mother Mom,K",100,0.0,5.0,80,0.,80.);
155   fMultESDK=new TH1F("fMultESDK", "charge multipliESD kaons",100, 0.0,100.0); 
156   fMultMCK=new TH1F("fMultMCK", "charge multipli MC kaons",100, 0.0,100.0); 
157   fSignPtNcl= new TH2F("fSignPtNcl","SignPt vrs Ncl,K",80,-4.,4.0,70,20.,160.);
158   fSignPtEta= new TH2F("fSignPtEta","SignPt vrs Eta,K",80,-4.0,4.0,30,-1.5,1.5);
159   fEtaNcl= new TH2F("fEtaNcl","Eta vrs nclust,K",30,-1.5,1.5, 70,20, 160);
160   //fSignPt= new TH1F("fSignPt","SignPt ,K",40,-4.0,4.0);
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., 70,20, 160);
165     fTPCSgnlP = new TH2F("fTPCSgnlP","TPC signal de/dx Mom,K",100,0.0,4.0,100,0.,250.);
166   fTPCSgnlPa= new TH2F("fTPCSgnlPa","TPC signal de/dx Mom,K",100,0.0,4.,100, 0.,250.);
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   fPtPrKink=new TH1F("fPtPrKink","pt of ESD  kaonKink tracks",100, 0.0,10.0);
176   // fHistPtKaoP = new TH1F("fHistPtKaoP", "P_{T}KaonP  distribution",50, 0.0,5.0); 
177   fHistPtKaoP = new TH1F("fHistPtKaoP", "P_{T}KaonP  distribution",100, 0.0,10.0); 
178   // fHistPtKaoN = new TH1F("fHistPtKaoN", "P_{T}KaonN  distribution",50, 0.0,5.0); 
179   fHistPtKaoN = new TH1F("fHistPtKaoN", "P_{T}KaonN  distribution",100, 0.0,10.0); 
180
181   frapiKESD=new TH1F("frapiKESD","rapid Kdistribution", 26,-1.3, 1.3); 
182   flifetime= new TH1F("flifetime", "ct study of K-kinks",100,0.,1000.); 
183   fradLK= new TH1F("fradLK", "Length of   K generated",100,0.,1000.); 
184   fradPtRpDt=new TH3F("fradPtRpDt","rad pt rap dat",28,100.,240., 20, 0., 5., 20, -1., 1. );
185   fInvMuNuAll= new TH1F("fInvMuNuAll", " Inv Mass MuNu all kink",80,0.,0.8); 
186   fQtInvM= new TH2F("fQtInvM", "Q_{T} Versus Inv MuNu ",80, 0., 0.80 , 100 , 0., 0.300); 
187     fDCAkink = new TH1F("fDCAkink ", "DCA kink vetrex ",50, 0.0,1.0);
188
189   fPosiKink= new TH2F("fPosiKink", "Y vrx kink Vrex ",100, -300.0,300.0,100, -300, 300.);
190   fPosiKinkK= new TH2F("fPosiKinkK", "Y vrx kink VrexK ",100, -300.0,300.0,100, -300, 300.);
191   fPosiKinKXZ= new TH2F("fPosiKinKXZ", "Y vrx kink VrexK ",100, -300.0,300.0,100, -300, 300.);
192   fPosiKinKYZ= new TH2F("fPosiKinKYZ", "Y vrx kink VrexK ",100, -300.0,300.0,100, -300, 300.);
193   fPosiKinKBg= new TH2F("fPosiKinKBg", "z vrx kink rad    ",100, -300.0,300.0,100,100., 300.);
194   fQtMothP = new TH2F("fQtMothP", " Qt vrs Mother P", 100, 0., 5.0,100, 0.,0.300);
195
196    fListOfHistos=new TList();
197
198    fListOfHistos->Add(fHistPtESD);
199    fListOfHistos->Add(fHistPt);
200    fListOfHistos->Add(fHistQtAll);
201    fListOfHistos->Add(fHistQt1);
202    fListOfHistos->Add(fHistQt2);
203    fListOfHistos->Add(fHistPtKaon);
204    fListOfHistos->Add(fHistPtKPDG);
205    fListOfHistos->Add(fHistEta);
206    fListOfHistos->Add(fHistEtaK);
207    fListOfHistos->Add(fptKMC);
208    fListOfHistos->Add(fMultiplMC);
209    fListOfHistos->Add(fESDMult);
210    fListOfHistos->Add(fgenpt);
211    fListOfHistos->Add(frad);
212    fListOfHistos->Add(fKinkKaon);
213    fListOfHistos->Add(fKinKRbn);
214    fListOfHistos->Add(fKinkKaonBg);
215    fListOfHistos->Add(fM1kaon);
216    fListOfHistos->Add(fgenPtEtR);
217    fListOfHistos->Add(fPtKink);
218    fListOfHistos->Add(fptKink);
219    fListOfHistos->Add(fcodeH);
220    fListOfHistos->Add(fdcodeH);
221    fListOfHistos->Add(fAngMomK);
222    fListOfHistos->Add(fAngMomPi);
223    fListOfHistos->Add(fAngMomKC);
224    fListOfHistos->Add(fMultESDK);
225    fListOfHistos->Add(fMultMCK);
226    fListOfHistos->Add(fSignPtNcl);
227    fListOfHistos->Add(fSignPtEta);
228    fListOfHistos->Add(fEtaNcl);
229    fListOfHistos->Add(fSignPt);
230    fListOfHistos->Add(fChi2NclTPC);
231    fListOfHistos->Add(fRatChi2Ncl);
232    fListOfHistos->Add(fRadiusNcl);
233    fListOfHistos->Add(fTPCSgnlP);
234    fListOfHistos->Add(fTPCSgnlPa);
235    fListOfHistos->Add(fRpr);
236    fListOfHistos->Add(fZpr);
237    fListOfHistos->Add(fdcatoVxXY);
238    fListOfHistos->Add(fnSigmToVx);
239    fListOfHistos->Add(fKinkMothDau);
240    fListOfHistos->Add(fZvXv);
241    fListOfHistos->Add(fZvYv);
242    fListOfHistos->Add(fXvYv);
243    fListOfHistos->Add(fPtPrKink);
244    fListOfHistos->Add(fHistPtKaoP);
245    fListOfHistos->Add(fHistPtKaoN);
246    fListOfHistos->Add(frapiKESD);
247    fListOfHistos->Add(flifetime);
248    fListOfHistos->Add(fradLK);
249    fListOfHistos->Add(fradPtRpDt);
250    fListOfHistos->Add(fInvMuNuAll);
251    fListOfHistos->Add(fQtInvM);
252    fListOfHistos->Add(fDCAkink);
253    fListOfHistos->Add(fPosiKink);
254    fListOfHistos->Add(fPosiKinkK);
255    fListOfHistos->Add(fPosiKinKXZ);
256    fListOfHistos->Add(fPosiKinKYZ);
257    fListOfHistos->Add(fPosiKinKBg);
258    fListOfHistos->Add(fQtMothP);
259
260
261 }
262 //=======================new thing
263 //     Float_t nCrossedRowsTPC = esdTrack->GetTPCClusterInfo(2,1);
264 //  Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
265 //  if (esdTrack->GetTPCNclsF()>0) {
266 //    ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC / esdTrack->GetTPCNclsF();
267 //   }
268 //________________________________________________________________________
269 void AliAnalysisKinkESDat::UserExec(Option_t *) 
270 {
271   // Main loop
272   // Called for each event
273
274   // Process MC truth, therefore we receive the AliAnalysisManager and ask it for the AliMCEventHandler
275   // This handler can return the current MC event
276   
277    AliVEvent *event = InputEvent();
278   if (!event) {
279      Printf("ERROR: Could not retrieve event");
280      return;
281   }
282
283   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
284   if (!esd) {
285      Printf("ERROR: Could not retrieve esd");
286      return;
287   }
288 //==================check of Physics selection?
289        Bool_t isSelected =
290 ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()&AliVEvent::kMB;
291
292          if ( isSelected ==kFALSE) return;   //  24/6/11 apo MF
293 //
294 //
295    Int_t nESDTracks =  esd->GetNumberOfTracks();
296       fMultMCK->Fill(nESDTracks);
297 //===============Marek multiplicity
298    Float_t refmultiplicity=fCutsMul->CountAcceptedTracks(esd);
299         if(fLowMulcut>-1)
300         {
301                 if(refmultiplicity<fLowMulcut)
302                         return;
303         }
304         if(fUpMulcut>-1)
305         {
306                 if(refmultiplicity>fUpMulcut)
307                         return;
308         }
309
310
311
312        fMultESDK->Fill(refmultiplicity);
313
314 //
315 //   Int_t nESDTracks =  esd->GetNumberOfTracks();
316   //   if ( nESDTracks>0 ) fMultMCK->Fill(nESDTracks);
317
318
319
320   const AliESDVertex *vertex=GetEventVertex(esd);    // 22/8
321     if(!vertex) return;
322       fMultiplMC->Fill(nESDTracks);
323 //
324   Double_t vpos[3];
325   vertex->GetXYZ(vpos);
326     fZpr->Fill(vpos[2]);         
327      if (TMath::Abs( vpos[2] ) > 10. ) return;   
328
329     
330
331   Double_t vtrack[3], ptrack[3];
332   
333      
334  Int_t nESDTracK = 0;
335  Int_t nESDTrKink = 0;
336
337    Int_t nGoodTracks =  esd->GetNumberOfTracks();
338     fESDMult->Fill(nGoodTracks);
339       
340   Double_t fAlephParameters[5] = {0.0283086,
341                                   2.63394e+01,
342                                   5.04114e-11,
343                                   2.12543e+00,
344                                   4.88663e+00};
345        Double_t nsigma = 100.0;
346        AliESDpid *fESDpid = new AliESDpid();                  
347     fESDpid->GetTPCResponse().SetBetheBlochParameters(fAlephParameters[0],
348                                                     fAlephParameters[1],
349                                                    fAlephParameters[2],
350                                                  fAlephParameters[3],
351                                                fAlephParameters[4]);
352 //
353 // track loop
354    for (Int_t iTracks = 0; iTracks < esd->GetNumberOfTracks(); iTracks++) {
355
356     AliESDtrack* track = esd->GetTrack(iTracks);
357     if (!track) {
358       Printf("ERROR: Could not receive track %d", iTracks);
359       continue;
360     }
361     
362     fHistPt->Fill(track->Pt());
363
364
365      //    sigmas
366     nsigma = TMath::Abs(fESDpid->NumberOfSigmasTPC(track,AliPID::kKaon));
367
368
369
370       Int_t tpcNCl = track->GetTPCclusters(0);  
371       Double_t tpcSign = track->GetSign();  
372     
373     Int_t label = track->GetLabel();
374     label = TMath::Abs(label);
375
376
377     UInt_t status=track->GetStatus();
378
379     if((status&AliESDtrack::kITSrefit)==0) continue;   
380     if((status&AliESDtrack::kTPCrefit)==0) continue;
381       if((track->GetTPCchi2()/track->GetTPCclusters(0))>3.8) continue;  
382
383       Double_t extCovPos[15];
384       track->GetExternalCovariance(extCovPos);    
385       if(extCovPos[0]>2) continue;
386      if(extCovPos[2]>2) continue;    
387      if(extCovPos[5]>0.5) continue;  
388      if(extCovPos[9]>0.5) continue;
389      if(extCovPos[14]>2) continue;
390
391
392     track->GetXYZ(vtrack);
393  fXvYv->Fill(vtrack[0],vtrack[1]);  
394  fZvYv->Fill(vtrack[0],vtrack[2]);  
395  fZvXv->Fill(vtrack[1],vtrack[2]);  
396
397 // track momentum, rapidity calculation
398      track->GetPxPyPz(ptrack);
399     
400     TVector3 trackMom(ptrack[0],ptrack[1],ptrack[2]);
401     
402           Double_t   etracK= TMath::Sqrt(trackMom.Mag()*trackMom.Mag() + 0.493677 *0.493677  );
403          Double_t rapiditK = 0.5 * (TMath::Log(  (etracK + ptrack[2]  ) / ( etracK - ptrack[2])  ))  ;
404     
405     Double_t trackEta=trackMom.Eta();
406      Double_t trMoment=trackMom.Mag();       
407     Double_t trackPt = track->Pt();
408     
409     
410       
411     Float_t bpos[2];
412     Float_t bCovpos[3];
413     track->GetImpactParameters(bpos,bCovpos);
414     
415     if (bCovpos[0]<=0 || bCovpos[2]<=0) {
416      Printf("Estimated b resolution lower or equal zero!");
417      bCovpos[0]=0; bCovpos[2]=0;
418     }
419
420     Float_t dcaToVertexXYpos = bpos[0];
421     Float_t dcaToVertexZpos = bpos[1];
422     
423     fRpr->Fill(dcaToVertexZpos);
424  
425    if((TMath::Abs(dcaToVertexXYpos)>0.3)||(TMath::Abs(dcaToVertexZpos)>3.0))  nESDTrKink++;  //  count of second  23Jul11    
426 //
427      //   if((dcaToVertexXYpos>0.3)||(dcaToVertexZpos>0.3)) continue;   //    allagi-dokini    3/6                 
428      // if((TMath::Abs(dcaToVertexXYpos)>0.4)||(dcaToVertexZpos>2.5)) continue;   //    allagi  23Jul11               
429      if((TMath::Abs(dcaToVertexXYpos)>0.3)||(TMath::Abs(dcaToVertexZpos)>2.5)) continue;   //    allagi  23Jul11               
430     
431 //  track Mult. after selection 
432     nESDTracK++;        
433   //    
434 //=========================================
435     fHistPtESD->Fill(track->Pt());
436
437    // Add Kink analysis           =============================
438    
439                 Int_t indexKinkPos=track->GetKinkIndex(0);
440 //  loop on kinks
441                 if(indexKinkPos<0){     ////mother kink
442                fPtKink->Fill(track->Pt()); /// pt from track
443
444         // select kink class    
445
446           AliESDkink *kink=esd->GetKink(TMath::Abs(indexKinkPos)-1);
447 //
448          // DCA kink
449           Double_t  Dist2 = kink->GetDistance();
450           fDCAkink->Fill( Dist2   );
451 //
452         
453       
454            const TVector3 vposKink(kink->GetPosition());
455  fPosiKink ->Fill( vposKink[0], vposKink[1]  );
456 //   Double_t  lengthK = TMath::Sqrt( vposKink[0]*vposKink[0] + vposKink[1]*vposKink[1] + vposKink[2]*vposKink[2] ) ;
457    // Double_t dxKink = vpos[0]-vposKink[0], dyKink=vpos[1]-vposKink[1], dzKink=vpos[2]-vposKink[2]; 
458    Double_t  dzKink=vpos[2]-vposKink[2]; 
459    // Double_t lifeKink= TMath::Sqrt( dxKink*dxKink + dyKink*dyKink + dzKink*dzKink ) ;
460 //
461             Double_t tanLamda = track->GetTgl();  // 25/6/2010
462
463    Double_t lifeKink= (TMath::Abs( dzKink ))*( TMath::Sqrt(1.+ tanLamda*tanLamda) ) / (TMath::Abs( tanLamda)) ;
464
465            const TVector3 motherMfromKink(kink->GetMotherP());
466            const TVector3 daughterMKink(kink->GetDaughterP());
467
468            Float_t qT=kink->GetQt();
469             Float_t motherPt=motherMfromKink.Pt();
470        //     Float_t etaMother=motherMfromKink.Eta();
471
472            fHistQtAll->Fill(qT) ;  //  Qt   distr
473                   
474            fptKink->Fill(motherMfromKink.Pt()); /// pt from kink
475
476            Float_t kinkAngle=TMath::RadToDeg()*kink->GetAngle(2);
477
478           if(  (TMath::Abs(rapiditK )) > 0.7 ) continue;
479         if ( (track->Pt())<.250)continue;
480
481                 fQtMothP->Fill( track->P(), qT);
482
483     if ( qT> 0.04)  fHistQt1  ->Fill(qT) ;  //  Qt   distr
484
485
486
487             fHistEta->Fill(trackEta) ;  //   Eta distr of PDG kink ESD  kaons
488       fHistQt2->Fill(qT);  // PDG ESD kaons            
489
490 //          maximum decay angle at a given mother momentum
491            //Double_t maxDecAngKmu=f1->Eval(motherMfromKink.Mag(),0.,0.,0.);
492            Double_t maxDecAngKmu=f1->Eval(track->P()          ,0.,0.,0.);
493            Double_t maxDecAngpimu=f2->Eval(    track->P(),       0.,0.,0.);
494          if( (kinkAngle<2.)  ) continue;
495
496            
497       //  BG  ?????==============
498               if ( TMath::Abs(vposKink[2]) >  225. ) continue ;
499               if ( TMath::Abs(vposKink[2]) <  0.5 ) continue ;
500 //                if(( vposKink[2] >0. )&& (vposKink[2]< 5.)  ) continue; 
501 //
502             fKinkKaonBg->Fill(motherPt);     
503  fAngMomPi->Fill( track->P(),           kinkAngle); 
504 //
505 // invariant mass of mother track decaying to mu
506          Float_t energyDaughterMu=TMath::Sqrt(daughterMKink.Mag()*daughterMKink.Mag()+0.105658*0.105658);
507          Float_t p1XM= motherMfromKink.Px();
508          Float_t p1YM= motherMfromKink.Py();
509          Float_t p1ZM= motherMfromKink.Pz();
510          Float_t p2XM= daughterMKink.Px();
511          Float_t p2YM= daughterMKink.Py();
512          Float_t p2ZM= daughterMKink.Pz();
513          Float_t p3Daughter=TMath::Sqrt(((p1XM-p2XM)*(p1XM-p2XM))+((p1YM-p2YM)*(p1YM-p2YM))+((p1ZM-p2ZM)*(p1ZM-p2ZM)));
514          Double_t invariantMassKmu= TMath::Sqrt((energyDaughterMu+p3Daughter)*(energyDaughterMu+p3Daughter)-motherMfromKink.Mag()*motherMfromKink.Mag());
515            fQtInvM -> Fill ( invariantMassKmu,  qT);
516            fInvMuNuAll->Fill(invariantMassKmu);
517    //     Float_t ptKink=TMath::Sqrt(p1XM*p1XM + p1YM*p1YM);
518   
519       if( ( kink->GetR()> 120 ) && ( kink->GetR() < 210 )  )  {
520       if (qT>0.12)  fAngMomKC->Fill(track->P(), kinkAngle); 
521           if ( qT>0.12) fM1kaon->Fill(invariantMassKmu);
522              if ( qT > 0.12) 
523          fRadiusNcl->Fill( (kink->GetR()) ,(track->GetTPCclusters(0)  ) ) ;
524   }    
525         //  fPosiKinKBg->Fill( vposKink[2], kink->GetR() );
526           if(  ( tpcNCl<30) ) continue;
527          if( ( ( track->GetTPCclusters(0) ) / (kink->GetR() ) ) > 0.63 ) continue;
528             if( ( ( track->GetTPCclusters(0) ) / (kink->GetR() ) ) < 0.20 ) continue;
529
530 //
531                fHistPtKPDG->Fill(track->Pt());  // ALL KAONS (pdg) inside ESD  kink sample
532 //if((kinkAngle>maxDecAngpimu)&&(qT>0.04)&&(qT<0.30)&&((kink->GetR()>120.)&&(kink->GetR()<200.))&&(TMath::Abs(trackEta)<0.9)&&(invariantMassKmu<0.6)) {
533  // if((kinkAngle>maxDecAngpimu)&&(qT>0.04)&&(qT<0.30)&&((kink->GetR()>=120.)&&(kink->GetR()<=200.))&&(TMath::Abs(etaMother)<0.9)&&(invariantMassKmu<0.6)){
534     if((kinkAngle>maxDecAngpimu)&&(qT>0.12)&&(qT<0.30)&&((kink->GetR()>=120.)&&(kink->GetR()<=210.))&&(TMath::Abs(rapiditK)<0.7)&&(invariantMassKmu<0.6)){
535  // 16/10  if((kinkAngle>maxDecAngpimu)&&(qT>0.04)&&(qT<0.30)&&((kink->GetR()>=120.)&&(kink->GetR()<=210.))&&(TMath::Abs(rapiditK)<0.7)&&(invariantMassKmu<0.6)){
536 //  if((kinkAngle>maxDecAngpimu)&&(qT>0.04)&&(qT<0.30)&&((kink->GetR()>=133.)&&(kink->GetR()<=179.))&&(TMath::Abs(rapiditK)<0.5)&&(invariantMassKmu<0.6)){      // STAR 
537
538         if( (kinkAngle<maxDecAngpimu*1.2) ) continue; 
539                  if ( (kinkAngle>maxDecAngKmu*.98) && ( track->P() >1.2 )) continue;  ///5/5/2010
540
541 /*
542 */
543
544          fTPCSgnlPa->Fill( trMoment ,(track->GetTPCsignal()  ) ) ;
545                           if ( nsigma               > 3.5) continue;
546                                      fHistPtKaon->Fill(track->Pt());   //all PID kink-kaon
547               if(tpcSign >0.)        fHistPtKaoP->Fill( track->Pt()         ) ;   //all PID kink-kaon
548                if ( tpcSign <0.)    fHistPtKaoN->Fill(  track->Pt()        ) ;   //all PID kink-kaon
549                 frad->Fill(kink->GetR());  // kink 
550                fradLK->Fill(lifeKink    );  // kink 
551              fHistEtaK->Fill(trackEta);
552             frapiKESD ->Fill(rapiditK);  //  rapidityof kaons 
553           fPosiKinKBg->Fill( vposKink[2], kink->GetR() );
554
555                      Float_t signPt= tpcSign*trackPt;
556                   fSignPtNcl->Fill( signPt  ,   tpcNCl   );   ///  28/4/2010
557                   fSignPtEta->Fill( signPt  , rapiditK  );
558                   fEtaNcl->Fill( rapiditK, tpcNCl    );
559                   fSignPt->Fill( signPt );
560                   fChi2NclTPC->Fill( (track->GetTPCchi2() ) ,  tpcNCl );
561          fRatChi2Ncl-> Fill (  (track->GetTPCchi2()/track->GetTPCclusters(0)  )) ;
562     fdcatoVxXY->Fill(dcaToVertexXYpos);
563                   fKinkMothDau->Fill(track->P(),daughterMKink.Mag());
564          fTPCSgnlP->Fill(track->P(), (track->GetTPCsignal()  ) ) ;
565                         flifetime->Fill(( lifeKink*.493667   )  /track->P()   ) ;
566              fKinkKaon->Fill(track->Pt());        
567
568              fKinKRbn->Fill(track->Pt());        
569              fptKMC   ->Fill(  track->Pt()    );        
570               fradPtRpDt->Fill( kink->GetR(), 1./track->Pt(), rapiditK); 
571                fAngMomK->Fill(    track->P(),        kinkAngle); 
572        fPosiKinkK->Fill( vposKink[0], vposKink[1]  );
573           fPosiKinKXZ->Fill( vposKink[2], vposKink[0]  );
574           fPosiKinKYZ->Fill( vposKink[2], vposKink[1]  );
575
576         }  //  kink selection 
577                   
578
579         }  //End Kink Information    
580   
581
582   } //track loop 
583
584  //     fMultiplMC->Fill(nESDTracK );
585
586   PostData(1, fListOfHistos);
587
588 }      
589
590 //________________________________________________________________________
591 void AliAnalysisKinkESDat::Terminate(Option_t *) 
592 {
593   // Draw result to the screen 
594   // Called once at the end of the query
595
596 }
597
598 //____________________________________________________________________//
599
600
601 const AliESDVertex* AliAnalysisKinkESDat::GetEventVertex(const AliESDEvent* esd) const
602   // Get the vertex from the ESD and returns it if the vertex is valid
603   
604 {
605   // Get the vertex 
606   
607 // 24/3  const AliESDVertex* vertex = esd->GetPrimaryVertex();
608    const AliESDVertex* vertex = esd->GetPrimaryVertexTracks();
609
610   // if((vertex->GetStatus()==kTRUE)&&(vertex->GetNContributors()>2)) return vertex;
611   if((vertex->GetStatus()==kTRUE)) return vertex;
612   else
613   { 
614      vertex = esd->GetPrimaryVertexSPD();
615       if((vertex->GetStatus()==kTRUE)&&(vertex->GetNContributors()>0)) return vertex;
616 //     if((vertex->GetStatus()==kTRUE)) return vertex;
617      else
618      return 0;
619   }
620 }