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