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