coverity fix
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / Kinks / AliAnalysisKinkESDat.cxx
CommitLineData
828a6e0c 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"
51ClassImp(AliAnalysisKinkESDat)
52
53
54//________________________________________________________________________
55AliAnalysisKinkESDat::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//________________________________________________________________________
98void 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//________________________________________________________________________
269void 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//________________________________________________________________________
591void 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
601const 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}