]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/KINK/AliAnalysisKinkESDMC.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Update for Kink tasks from Evi Ganoti
[u/mrichter/AliRoot.git] / PWG2 / KINK / AliAnalysisKinkESDMC.cxx
CommitLineData
10eaad41 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 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16//-----------------------------------------------------------------
17// AliAnalysisKinkESDMC class
18// Example of an analysis task for kink topology study
19// Kaons from kink topology are 'identified' in this code
20//-----------------------------------------------------------------
21
22#include "TChain.h"
23#include "TTree.h"
24#include "TH1F.h"
25#include "TH2F.h"
26#include "TH1D.h"
27#include "TH2D.h"
28#include "TParticle.h"
29#include <TVector3.h>
30#include "TF1.h"
31
92adf4f6 32#include "AliAnalysisTask.h"
10eaad41 33#include "AliAnalysisManager.h"
34
92adf4f6 35#include "AliVEvent.h"
10eaad41 36#include "AliESDEvent.h"
10eaad41 37#include "AliMCEvent.h"
10eaad41 38#include "AliAnalysisKinkESDMC.h"
10eaad41 39#include "AliStack.h"
40#include "AliESDpid.h"
41#include "AliESDkink.h"
42
43ClassImp(AliAnalysisKinkESDMC)
44
10eaad41 45
46//________________________________________________________________________
47AliAnalysisKinkESDMC::AliAnalysisKinkESDMC(const char *name)
92adf4f6 48 : AliAnalysisTaskSE(name), fHistPtESD(0),fHistPt(0),fHistQtAll(0),fHistQt1(0),fHistQt2(0)
10eaad41 49 , fHistPtKaon(0),fHistPtKPDG(0),fHistEta(0),fHistEtaK(0),fptKMC(0),fMultiplMC(0),fESDMult(0),fgenpt(0),frad(0),
50 fKinkKaon(0), fKinkKaonBg(0), fM1kaon(0), fgenPtEtR(0),fPtKink(0), fptKink(0),
92adf4f6 51 fcodeH(0), fdcodeH(0), fAngMomK(0),fAngMomPi(0), fAngMomKC(0), fMultESDK(0), fMultMCK(0),
10eaad41 52 fRpr(0),fZpr(0),
10eaad41 53 fZvXv(0),fZvYv(0), fXvYv(0), fPtPrKink(0), f1(0), f2(0),
54 fListOfHistos(0)
10eaad41 55
56{
57 // Constructor
58
59 // Define input and output slots here
60 // Input slot #0 works with a TChain
61 DefineInput(0, TChain::Class());
62 DefineOutput(1, TList::Class());
63}
64
65//________________________________________________________________________
92adf4f6 66void AliAnalysisKinkESDMC::UserCreateOutputObjects()
10eaad41 67{
68 // Create histograms
69 // Called once
70
71 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);
72 f1->SetParameter(0,0.493677);
73 f1->SetParameter(1,0.9127037);
74 f1->SetParameter(2,TMath::Pi());
75
76
77 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);
78 f2->SetParameter(0,0.13957018);
79 f2->SetParameter(1,0.2731374);
80 f2->SetParameter(2,TMath::Pi());
81 //Open file 1= CAF
82 //OpenFile(1);
83
92adf4f6 84 fHistPtESD = new TH1F("fHistPtESD", "P_{T} distribution",100, 0.0,10.0);
10eaad41 85 fHistPtESD->GetXaxis()->SetTitle("P_{T} (GeV/c)");
86 fHistPtESD->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
87 fHistPtESD->SetMarkerStyle(kFullCircle);
92adf4f6 88 fHistPt = new TH1F("fHistPt", "P_{T} distribution",100, 0.0,10.0);
10eaad41 89 fHistQtAll = new TH1F("fHistQtAll", "Q_{T} distr All Kinks ",100, 0.0,.250);
90 fHistQt1= new TH1F("fHistQt1", "Q_{T} distribution",100, 0.0,.300);
91 fHistQt2= new TH1F("fHistQt2", "Q_{T} distribution",100, 0.0,.300);
92adf4f6 92 fHistPtKaon = new TH1F("fHistPtKaon", "P_{T}Kaon distribution",100, 0.0,10.0);
93 fHistPtKPDG = new TH1F("fHistPtKPDG", "P_{T}Kaon distribution",100, 0.0,10.0);
10eaad41 94 fHistEta= new TH1F("fHistEta", "Eta distribution", 26,-1.3, 1.3);
95 fHistEtaK= new TH1F("fHistEtaK", "EtaK distribution", 26,-1.3, 1.3);
92adf4f6 96 fptKMC= new TH1F("fptKMC", "P_{T}Kaon generated",100, 0.0,10.0);
97 fMultiplMC= new TH1F("fMultiplMC", "charge multiplicity MC",60, 0.5,120.5);
10eaad41 98 fESDMult= new TH1F("fESDMult", "charge multipliESD", 60, 0.5,120.5);
92adf4f6 99 fgenpt= new TH1F("fgenpt", "genpt K distribution",100, 0.0,10.0);
10eaad41 100 frad= new TH1F("frad", "radius K generated",100,0.,400.);
92adf4f6 101 fKinkKaon= new TH1F("fKinkKaon", "P_{T}Kaon kinks identi",100, 0.0,10.0);
102 fKinkKaonBg= new TH1F("fKinkKaonBg", "P_{T}Kaon kinks backgr",100, 0.0,10.0);
10eaad41 103 fM1kaon= new TH1F("fM1kaon","Invar m(kaon) from kink->mu+netrino decay",80,0.0, 0.8);
92adf4f6 104 fgenPtEtR= new TH1F("fgenPtEtR", "P_{T}Kaon distribution",100, 0.0,10.0);
105 fPtKink= new TH1F("fPtKink", "P_{T}Kaon Kink bution",100, 0.0,10.0);
106 fptKink= new TH1F("fptKink", "P_{T}Kaon Kink bution",100, 0.0,10.0);
10eaad41 107 fcodeH = new TH2F("fcodeH", "code vrs dcode dist. kinks,K",100,0.,2500.,100,0.,2500.);
108 fdcodeH = new TH2F("fdcodeH", "code vrs dcode dist. kinks,K",100,0.,2500.,100,0.,2500.);
109 fAngMomK= new TH2F("fAngMomK","Decay angle vrs Mother Mom,K",100,0.0,5.0,80,0.,80.);
110 fAngMomPi= new TH2F("fAngMomPi","Decay angle vrs Mother Mom,Pi",100,0.0,5.0,80,0.,80.);
92adf4f6 111 fAngMomKC= new TH2F("fAngMomKC","Decay angle vrs Mother Mom,K",100,0.0,5.0,80,0.,80.);
112 fMultESDK=new TH1F("fMultESDK", "charge multipliESD kaons", 60, 0.5,240.5);
113 fMultMCK=new TH1F("fMultMCK", "charge multipli MC kaons",100, 0.5,600.5);
10eaad41 114 fRpr = new TH1D("fRpr", "rad distribution PID pr",50,0.0, 2.5);
115 fZpr = new TH1D("fZpr", "z distribution PID pr ",60,-15.,15.);
10eaad41 116 fZvXv= new TH2F("fZvXv","Xv-Zv main vtx",60,-0.5,0.5,60, -15., 15.0);
117 fZvYv= new TH2F("fZvYv","Yv-Zv main vtx",60,-0.5,0.5, 60, -15., 15.);
118 fXvYv= new TH2F("fXvYv","Xv-Yv main vtx", 60,-1.5,1.5, 60, -1.5, 1.5);
92adf4f6 119 fPtPrKink=new TH1F("fPtPrKink","pt of ESD kaonKink tracks",100, 0.0,10.0);
10eaad41 120
121 fListOfHistos=new TList();
122
123 fListOfHistos->Add(fHistPtESD);
124 fListOfHistos->Add(fHistPt);
125 fListOfHistos->Add(fHistQtAll);
126 fListOfHistos->Add(fHistQt1);
127 fListOfHistos->Add(fHistQt2);
128 fListOfHistos->Add(fHistPtKaon);
129 fListOfHistos->Add(fHistPtKPDG);
130 fListOfHistos->Add(fHistEta);
131 fListOfHistos->Add(fHistEtaK);
132 fListOfHistos->Add(fptKMC);
133 fListOfHistos->Add(fMultiplMC);
134 fListOfHistos->Add(fESDMult);
135 fListOfHistos->Add(fgenpt);
136 fListOfHistos->Add(frad);
137 fListOfHistos->Add(fKinkKaon);
138 fListOfHistos->Add(fKinkKaonBg);
139 fListOfHistos->Add(fM1kaon);
140 fListOfHistos->Add(fgenPtEtR);
141 fListOfHistos->Add(fPtKink);
142 fListOfHistos->Add(fptKink);
143 fListOfHistos->Add(fcodeH);
144 fListOfHistos->Add(fdcodeH);
10eaad41 145 fListOfHistos->Add(fAngMomK);
146 fListOfHistos->Add(fAngMomPi);
92adf4f6 147 fListOfHistos->Add(fAngMomKC);
148 fListOfHistos->Add(fMultESDK);
149 fListOfHistos->Add(fMultMCK);
10eaad41 150 fListOfHistos->Add(fRpr);
151 fListOfHistos->Add(fZpr);
10eaad41 152 fListOfHistos->Add(fZvXv);
153 fListOfHistos->Add(fZvYv);
154 fListOfHistos->Add(fXvYv);
155 fListOfHistos->Add(fPtPrKink);
156}
157
158//________________________________________________________________________
92adf4f6 159void AliAnalysisKinkESDMC::UserExec(Option_t *)
10eaad41 160{
161 // Main loop
162 // Called for each event
163
164 // Process MC truth, therefore we receive the AliAnalysisManager and ask it for the AliMCEventHandler
165 // This handler can return the current MC event
92adf4f6 166
167 AliVEvent *event = InputEvent();
168 if (!event) {
169 Printf("ERROR: Could not retrieve event");
170 return;
171 }
172
173 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
174 if (!esd) {
175 Printf("ERROR: Could not retrieve esd");
176 return;
10eaad41 177 }
178
92adf4f6 179 AliMCEvent* mcEvent = MCEvent();
10eaad41 180 if (!mcEvent) {
181 Printf("ERROR: Could not retrieve MC event");
182 return;
183 }
184
185 Printf("MC particles: %d", mcEvent->GetNumberOfTracks());
186
10eaad41 187 AliStack* stack=mcEvent->Stack();
188
189//primary tracks in MC
190 Int_t nPrim = stack->GetNprimary();
191//
92adf4f6 192 const AliESDVertex *vertex=GetEventVertex(esd);
10eaad41 193 if(!vertex) return;
194//
195 Double_t vpos[3];
196 vertex->GetXYZ(vpos);
197 fZpr->Fill(vpos[2]);
92adf4f6 198// fZvXv->Fill( vpos[1], vpos[2]);
199// fZvYv->Fill( vpos[0], vpos[2]);
10eaad41 200
201 Double_t vtrack[3], ptrack[3];
202
203
92adf4f6 204 // Int_t nESDTracK = 0;
10eaad41 205
92adf4f6 206 Int_t nGoodTracks = esd->GetNumberOfTracks();
10eaad41 207 fESDMult->Fill(nGoodTracks);
92adf4f6 208
10eaad41 209//
210// track loop
92adf4f6 211 for (Int_t iTracks = 0; iTracks < esd->GetNumberOfTracks(); iTracks++) {
212// loop only on primary tracks
213
214 AliESDtrack* track = esd->GetTrack(iTracks);
10eaad41 215 if (!track) {
216 Printf("ERROR: Could not receive track %d", iTracks);
217 continue;
218 }
219
220
221 fHistPt->Fill(track->Pt());
222
92adf4f6 223
224 Int_t label = track->GetLabel();
225 label = TMath::Abs(label);
226 TParticle * part = stack->Particle(label);
227 if (!part) continue;
228// loop only on Primary tracks
229// gia oles tis troxies 5/7/2009 if (label > nPrim) continue;
230// if (label > nPrim) continue; /// primary tracks only
231 // ta nominal cauts tis apomakrynoun ols if (label < nPrim) continue; // gia tis secondary , 5/7/2009
232
10eaad41 233 // pt cut at 300 MeV
234 if ( (track->Pt())<.300)continue;
235
92adf4f6 236// UInt_t status=track->GetStatus();
10eaad41 237
238 // if((status&AliESDtrack::kITSrefit)==0) continue;
92adf4f6 239 // if((status&AliESDtrack::kTPCrefit)==0) continue;
240 // ayksanei to ratio BG/real if((track->GetTPCchi2()/track->GetTPCclusters(0))>4.5) continue;
241 if((track->GetTPCchi2()/track->GetTPCclusters(0))>4.5) continue;
10eaad41 242
243 Double_t extCovPos[15];
244 track->GetExternalCovariance(extCovPos);
92adf4f6 245 // if(extCovPos[0]>2) continue;
246 // if(extCovPos[2]>2) continue;
247 // if(extCovPos[5]>0.5) continue;
248 // if(extCovPos[9]>0.5) continue;
249 // if(extCovPos[14]>2) continue;
10eaad41 250
251
252 track->GetXYZ(vtrack);
253 fXvYv->Fill(vtrack[0],vtrack[1]);
92adf4f6 254 fZvYv->Fill(vtrack[0],vtrack[2]);
255 fZvXv->Fill(vtrack[1],vtrack[2]);
10eaad41 256
257// track momentum
258 track->GetPxPyPz(ptrack);
259
260 TVector3 trackMom(ptrack[0],ptrack[1],ptrack[2]);
261
262 Double_t trackEta=trackMom.Eta();
263
264
265
266 Float_t nSigmaToVertex = GetSigmaToVertex(track);
267
268 Float_t bpos[2];
269 Float_t bCovpos[3];
270 track->GetImpactParameters(bpos,bCovpos);
271
272 if (bCovpos[0]<=0 || bCovpos[2]<=0) {
273 Printf("Estimated b resolution lower or equal zero!");
274 bCovpos[0]=0; bCovpos[2]=0;
275 }
276
277 Float_t dcaToVertexXYpos = bpos[0];
278 Float_t dcaToVertexZpos = bpos[1];
279
280 fRpr->Fill(dcaToVertexZpos);
92adf4f6 281// test for secondary kinks , 5/7/2009
10eaad41 282 if(nSigmaToVertex>=4) continue;
92adf4f6 283 // if((dcaToVertexXYpos>3.0)||(dcaToVertexZpos>3.0)) continue; //arxiko
284 //if((dcaToVertexXYpos>0.4)||(dcaToVertexZpos>3.0)) continue;// me Zpos > 5 clen the secondaries
285 if((dcaToVertexXYpos>2.4)||(dcaToVertexZpos>5.0)) continue;// me Zpos > 5 clen the secondaries
286 // if((dcaToVertexXYpos>3.0)||(dcaToVertexZpos>3.0)) continue; // me ayto krataei 4 fakes apo ta 5
287 // if((dcaToVertexXYpos>4.0)||(dcaToVertexZpos>5.0)) continue;
10eaad41 288
10eaad41 289
290 // cut on eta
92adf4f6 291 if( (TMath::Abs(trackEta )) > 0.9 ) continue;
10eaad41 292 fHistPtESD->Fill(track->Pt());
293
294 // Add Kink analysis
295
296 Int_t indexKinkPos=track->GetKinkIndex(0);
297// loop on kinks
298 if(indexKinkPos<0){
299 fPtKink->Fill(track->Pt()); /// pt from track
300
301 // select kink class
302
92adf4f6 303 AliESDkink *kink=esd->GetKink(TMath::Abs(indexKinkPos)-1);
10eaad41 304//
305
306 Int_t eSDfLabel1=kink->GetLabel(0);
307 TParticle *particle1= stack->Particle(TMath::Abs(eSDfLabel1));
308 Int_t code1= particle1->GetPdgCode();
309
310 Int_t eSDfLabeld=kink->GetLabel(1);
311 TParticle *particled= stack->Particle(TMath::Abs(eSDfLabeld));
312 Int_t dcode1= particled->GetPdgCode();
313
314 const TVector3 motherMfromKink(kink->GetMotherP());
315 const TVector3 daughterMKink(kink->GetDaughterP());
316 Float_t qT=kink->GetQt();
317
92adf4f6 318 //fHistEta->Fill(trackEta) ; // Eta distr
10eaad41 319 fHistQtAll->Fill(qT) ; // Qt distr
320
321 fptKink->Fill(motherMfromKink.Pt()); /// pt from kink
322
92adf4f6 323 Float_t kinkAngle=TMath::RadToDeg()*kink->GetAngle(2);
324
10eaad41 325// fake kinks, small Qt and small kink angle
92adf4f6 326 if(( kinkAngle<1.)&&(TMath::Abs(code1)==321)&&(TMath::Abs(dcode1)==13)) fHistQt1 ->Fill(qT) ; // Qt distr
10eaad41 327 if ( qT<0.012) fcodeH->Fill( TMath::Abs(code1), TMath::Abs(dcode1) );
328 if(qT<0.012) continue; // remove fake kinks
92adf4f6 329 if( (kinkAngle<1.) ) continue;
10eaad41 330
10eaad41 331//remove the double taracks
332
92adf4f6 333 //if( (kinkAngle<1.) ) continue;
10eaad41 334//
10eaad41 335
92adf4f6 336 if((kink->GetR()>120.)&&(kink->GetR()<220.)&&(TMath::Abs(trackEta)<0.9)&&(label<nPrim)) {
337 if( ( (TMath::Abs(code1)==321)&&(TMath::Abs(dcode1)==13))||
338 ( (TMath::Abs(code1)==321)&&(TMath::Abs(dcode1)==11)) ||
339 ( (TMath::Abs(code1)==321)&&(TMath::Abs(dcode1)==211)) ) {
340 fHistPtKPDG->Fill(track->Pt()); // ALL KAONS (pdg) inside ESD kink sample
341 fHistEta->Fill(trackEta) ; // Eta distr of PDG kink ESD kaons
342 fMultESDK->Fill(nGoodTracks);
343 }
344 }
345 if( (qT>0.05)&& ( qT<0.25) ) fHistQt2->Fill(qT); // candidate kaon kinks
10eaad41 346
347// maximum decay angle at a given mother momentum
348 Double_t maxDecAngKmu=f1->Eval(motherMfromKink.Mag(),0.,0.,0.);
349 Double_t maxDecAngpimu=f2->Eval(motherMfromKink.Mag(),0.,0.,0.);
350// two dimensional plot
351 if(TMath::Abs(code1)==321) fAngMomK->Fill(motherMfromKink.Mag(), kinkAngle);
352 if(TMath::Abs(code1)==211) fAngMomPi->Fill(motherMfromKink.Mag(), kinkAngle);
353//
354// invariant mass of mother track decaying to mu
355 Float_t energyDaughterMu=TMath::Sqrt(daughterMKink.Mag()*daughterMKink.Mag()+0.105658*0.105658);
356 Float_t p1XM= motherMfromKink.Px();
357 Float_t p1YM= motherMfromKink.Py();
358 Float_t p1ZM= motherMfromKink.Pz();
359 Float_t p2XM= daughterMKink.Px();
360 Float_t p2YM= daughterMKink.Py();
361 Float_t p2ZM= daughterMKink.Pz();
362 Float_t p3Daughter=TMath::Sqrt(((p1XM-p2XM)*(p1XM-p2XM))+((p1YM-p2YM)*(p1YM-p2YM))+((p1ZM-p2ZM)*(p1ZM-p2ZM)));
363 Double_t invariantMassKmu= TMath::Sqrt((energyDaughterMu+p3Daughter)*(energyDaughterMu+p3Daughter)-motherMfromKink.Mag()*motherMfromKink.Mag());
364
365 fM1kaon->Fill(invariantMassKmu);
366
367// kaon selection from kinks
92adf4f6 368if((kinkAngle>maxDecAngpimu)&&(qT>0.05)&&(qT<0.25)&&((kink->GetR()>120.)&&(kink->GetR()<220.))&&(TMath::Abs(trackEta)<0.9)&&(invariantMassKmu<0.6)) {
369 // fHistEtaK->Fill(trackEta);
10eaad41 370
92adf4f6 371 // if(TMath::Abs(code1)==321) fAngMomKC->Fill(motherMfromKink.Mag(), kinkAngle); // real kaons
10eaad41 372
92adf4f6 373// fHistQt1 ->Fill(qT) ; // Qt distr
10eaad41 374 if ( (kinkAngle>maxDecAngKmu) && ( motherMfromKink.Mag()> 1.1 ) ) continue; // maximum angle selection
375
92adf4f6 376 if(TMath::Abs(code1)==321) fAngMomKC->Fill(motherMfromKink.Mag(), kinkAngle); // real kaons
10eaad41 377 fHistPtKaon->Fill(track->Pt()); //all PID kink-kaon
378
379// background inside the identified kaons, e.g KK , Kp
380 if((TMath::Abs(code1)==321)&&(( TMath::Abs(dcode1)) >=(TMath::Abs(code1)) )) fdcodeH->Fill( TMath::Abs(code1), TMath::Abs(dcode1));
92adf4f6 381// kaons from k to mun and k to pipi and to e decay
382 if( ( (TMath::Abs(code1)==321)&&(TMath::Abs(dcode1)==13))||
383 ( (TMath::Abs(code1)==321)&&(TMath::Abs(dcode1)==11)) ||
384 ( (TMath::Abs(code1)==321)&&(TMath::Abs(dcode1)==211)) ) {
385 // if((TMath::Abs(code1)==321) ) {
10eaad41 386
387 if ( label< nPrim ) fPtPrKink->Fill(track->Pt());
92adf4f6 388 fKinkKaon->Fill(track->Pt());
389 fHistEtaK->Fill(trackEta);
390 }
10eaad41 391 else {
392 fKinkKaonBg->Fill(track->Pt());
393 } // primary and secondary
394
395 } // kink selection
396
397
398 } //End Kink Information
399
400
401 } //track loop
402
403 // loop over mc particles
404
405 // variable to count tracks
406 Int_t nMCTracks = 0;
92adf4f6 407 Int_t nMCKinkKs = 0;
10eaad41 408
92adf4f6 409for (Int_t iMc = 0; iMc < nPrim; iMc++)
10eaad41 410 {
411 // AliDebug(AliLog::kDebug+1, Form("MC Loop: Processing particle %d.", iMc));
412
413 TParticle* particle = stack->Particle(iMc);
414
415 if (!particle)
416 {
417 // AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc));
418 continue;
419 }
420
421 //
422 Float_t eta = particle->Eta();
423
92adf4f6 424 if (TMath::Abs(eta) > 0.9)
10eaad41 425 continue;
426
427 Double_t ptK = particle->Pt();
428
429 if( ptK <0.300) continue;
430
431 Float_t code = particle->GetPdgCode();
92adf4f6 432 Int_t mcProcess=-1011;
433 if ((code==321)||(code==-321)){
10eaad41 434
92adf4f6 435
10eaad41 436 fptKMC->Fill(ptK);
92adf4f6 437
10eaad41 438
92adf4f6 439 // fMultiplMC->Fill(nPrim);
10eaad41 440
441 Int_t firstD=particle->GetFirstDaughter();
442 Int_t lastD=particle->GetLastDaughter();
443//loop on secondaries
444 for (Int_t k=firstD;k<lastD;k++) {
445 TParticle* daughter1=stack->Particle(k+1);
446 Float_t dcode = daughter1->GetPdgCode();
447
92adf4f6 448 mcProcess=daughter1->GetUniqueID();
449 if (mcProcess==4) {
10eaad41 450
451
92adf4f6 452 // frad->Fill(daughter1->R());
10eaad41 453
92adf4f6 454 if (((daughter1->R())>120)&&((daughter1->R())<220) ){
10eaad41 455 if (( ( code==321 )&& ( dcode ==-13 ))|| (( code == -321 )&& ( dcode == 13))) fgenPtEtR->Fill( ptK );//to muon
456 if (( ( code==321 )&& ( dcode ==211 ))|| (( code == -321 )&& ( dcode ==-211))) fgenPtEtR->Fill( ptK );//to pion
92adf4f6 457 if (( ( code==321 )&& ( dcode ==-11 ))|| (( code == -321 )&& ( dcode ==11))) fgenPtEtR->Fill( ptK );// to electr
10eaad41 458
92adf4f6 459 //fMultMCK->Fill(nPrim);
460 frad->Fill(daughter1->R());
461 nMCKinkKs++;
10eaad41 462 }
463// next only to mu decay
464 if (((code==321)&&(dcode==-13))||((code==-321)&&(dcode==13))){
465
466
92adf4f6 467 if (((daughter1->R())>120)&&((daughter1->R())<220)){
10eaad41 468 fgenpt->Fill(ptK);
469 }
470 }
92adf4f6 471 }// decay
472 }// daughters
10eaad41 473
92adf4f6 474 } /// kaons
10eaad41 475
476
477
478 nMCTracks++;
479 }// end of mc particle
480
481// printf("hello mc");
482 fMultiplMC->Fill(nMCTracks);
92adf4f6 483 if( nMCKinkKs>0) fMultMCK->Fill(nPrim);
10eaad41 484
485 PostData(1, fListOfHistos);
486
487}
488
489//________________________________________________________________________
490void AliAnalysisKinkESDMC::Terminate(Option_t *)
491{
492 // Draw result to the screen
493 // Called once at the end of the query
494
495}
496//____________________________________________________________________//
497
498Float_t AliAnalysisKinkESDMC::GetSigmaToVertex(AliESDtrack* esdTrack) const {
499 // Calculates the number of sigma to the vertex.
500
501 Float_t b[2];
502 Float_t bRes[2];
503 Float_t bCov[3];
504
505 esdTrack->GetImpactParameters(b,bCov);
506
507 if (bCov[0]<=0 || bCov[2]<=0) {
508 //AliDebug(1, "Estimated b resolution lower or equal zero!");
509 bCov[0]=0; bCov[2]=0;
510 }
511 bRes[0] = TMath::Sqrt(bCov[0]);
512 bRes[1] = TMath::Sqrt(bCov[2]);
513
514 if (bRes[0] == 0 || bRes[1] ==0) return -1;
515
516 Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
517
518 if (TMath::Exp(-d * d / 2) < 1e-10) return 1000;
519
520 d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
521
522 return d;
523 }
524//____________________________________________________________________//
525
526const AliESDVertex* AliAnalysisKinkESDMC::GetEventVertex(const AliESDEvent* esd) const
527 // Get the vertex from the ESD and returns it if the vertex is valid
528
529{
530 // Get the vertex
531
532 const AliESDVertex* vertex = esd->GetPrimaryVertex();
533
534 if((vertex->GetStatus()==kTRUE)&&(vertex->GetNContributors()>2)) return vertex;
535 else
536 {
537 vertex = esd->GetPrimaryVertexSPD();
538 if((vertex->GetStatus()==kTRUE)&&(vertex->GetNContributors()>2)) return vertex;
539 else
540 return 0;
541 }
f92b626a 542}