]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/SPECTRA/IdentifiedHighPt/grid/AliAnalysisTaskHighPtDeDxV0.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
TENDER becomes Tender
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / IdentifiedHighPt / grid / AliAnalysisTaskHighPtDeDxV0.cxx
CommitLineData
4ebdd20e 1#include "AliAnalysisTaskHighPtDeDxV0.h"
2
3// ROOT includes
4#include <TList.h>
5#include <TTree.h>
6#include <TMath.h>
7#include <TH1.h>
8#include <TParticle.h>
9#include <TFile.h>
10
11// AliRoot includes
12#include <AliAnalysisManager.h>
13#include <AliAnalysisFilter.h>
14#include <AliESDInputHandler.h>
15#include <AliESDEvent.h>
16#include <AliESDVertex.h>
17#include <AliLog.h>
18#include <AliExternalTrackParam.h>
19
20#include <AliMCEventHandler.h>
21#include <AliMCEvent.h>
22#include <AliStack.h>
23
24#include <AliHeader.h>
25#include <AliGenPythiaEventHeader.h>
26#include <AliGenDPMjetEventHeader.h>
27
28#include <AliAODVertex.h>
29#include <AliESDv0.h>
30#include <AliKFVertex.h>
31
32#include "AliCentrality.h"
33#include "AliESDtrackCuts.h"
34
35#include <AliAODTrack.h>
36#include <AliAODPid.h>
37#include <AliAODMCHeader.h>
38
39// STL includes
40#include <iostream>
41using namespace std;
42
43//_____________________________________________________________________________
44//
45// Responsible:
46// Peter Christiansen (Lund)
47//
48//_____________________________________________________________________________
49
50
51ClassImp(AliAnalysisTaskHighPtDeDxV0)
52
53const Double_t AliAnalysisTaskHighPtDeDxV0::fgkClight = 2.99792458e-2;
54
55//_____________________________________________________________________________
56AliAnalysisTaskHighPtDeDxV0::AliAnalysisTaskHighPtDeDxV0():
57 AliAnalysisTaskSE(),
58 fESD(0x0),
59 fAOD(0x0),
60 fMC(0x0),
61 fMCStack(0x0),
62 fMCArray(0x0),
63 fTrackFilter(0x0),
64 fTrackFilterGolden(0x0),
65 fTrackFilterTPC(0x0),
66 fAnalysisType("ESD"),
67 fAnalysisMC(kFALSE),
68 fAnalysisPbPb(kFALSE),
69 ftrigBit1(0x0),
70 ftrigBit2(0x0),
71 fEvent(0x0),
72 fV0ArrayGlobalPar(0x0),
73 fV0ArrayTPCPar(0x0),
74 fTrackArrayMC(0x0),
75 fVtxCut(10.0),
76 fEtaCut(0.9),
77 fMinPt(0.1),
78 fMinCent(0.0),
79 fMaxCent(100.0),
80 fMassCut(0.1),
81 fTreeOption(0),
82 fRequireRecV0(kTRUE),
83 fStoreMcIn(kFALSE),
84 fMcProcessType(-999),
85 fTriggeredEventMB(-999),
86 fVtxStatus(-999),
87 fZvtx(-999),
88 fZvtxMC(-999),
89 fRun(-999),
90 fEventId(-999),
91 fListOfObjects(0),
92 fEvents(0x0), fVtx(0x0), fVtxMC(0x0), fVtxBeforeCuts(0x0), fVtxAfterCuts(0x0),
93 fTree(0x0),
94 fnv0(0)
95{
96 // Default constructor (should not be used)
97
98}
99
100//______________________________________________________________________________
101AliAnalysisTaskHighPtDeDxV0::AliAnalysisTaskHighPtDeDxV0(const char *name):
102 AliAnalysisTaskSE(name),
103 fESD(0x0),
104 fAOD(0x0),
105 fMC(0x0),
106 fMCStack(0x0),
107 fMCArray(0x0),
108 fTrackFilter(0x0),
109 fTrackFilterGolden(0x0),
110 fTrackFilterTPC(0x0),
111 fAnalysisType("ESD"),
112 fAnalysisMC(kFALSE),
113 fAnalysisPbPb(kFALSE),
114 ftrigBit1(0x0),
115 ftrigBit2(0x0),
116 fEvent(0x0),
117 fV0ArrayGlobalPar(0x0),
118 fV0ArrayTPCPar(0x0),
119 fTrackArrayMC(0x0),
120 fVtxCut(10.0),
121 fEtaCut(0.9),
122 fMinPt(0.1),
123 fMinCent(0.0),
124 fMaxCent(100.0),
125 fMassCut(0.1),
126 fTreeOption(0),
127 fRequireRecV0(kTRUE),
128 fStoreMcIn(kFALSE),
129 fMcProcessType(-999),
130 fTriggeredEventMB(-999),
131 fVtxStatus(-999),
132 fZvtx(-999),
133 fZvtxMC(-999),
134 fRun(-999),
135 fEventId(-999),
136 fListOfObjects(0),
137 fEvents(0x0), fVtx(0x0), fVtxMC(0x0), fVtxBeforeCuts(0x0), fVtxAfterCuts(0x0),
138 fTree(0x0),
139 fnv0(0)
140{
141 // Output slot #1 writes into a TList
142 DefineOutput(1, TList::Class());
143}
144
145//_____________________________________________________________________________
146AliAnalysisTaskHighPtDeDxV0::~AliAnalysisTaskHighPtDeDxV0()
147{
148 // Destructor
149 // histograms are in the output list and deleted when the output
150 // list is deleted by the TSelector dtor
151 if (fListOfObjects && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
152 delete fListOfObjects;
153 fListOfObjects = 0;
154 }
155
156 // //for proof running; I cannot create tree do to memory limitations -> work with THnSparse
157 // if (fListOfObjects && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fOutputList;
158
159
160
161}
162
163//______________________________________________________________________________
164void AliAnalysisTaskHighPtDeDxV0::UserCreateOutputObjects()
165{
166 // This method is called once per worker node
167 // Here we define the output: histograms and debug tree if requested
168
169 OpenFile(1);
170 fListOfObjects = new TList();
171 fListOfObjects->SetOwner();
172
173 //
174 // Histograms
175 //
176 fEvents = new TH1I("fEvents","Number of analyzed events; Events; Counts", 1, 0, 1);
177 fListOfObjects->Add(fEvents);
178
179 fVtx = new TH1I("fVtx","Vtx info (0=no, 1=yes); Vtx; Counts", 2, -0.5, 1.5);
180 fListOfObjects->Add(fVtx);
181
182 fnv0 = new TH1F("fnv0","# V0's", 10000,0,10000);
183 fListOfObjects->Add(fnv0);
184
185 fVtxBeforeCuts = new TH1F("fVtxBeforeCuts", "Vtx distribution (before cuts); Vtx z [cm]; Counts", 120, -30, 30);
186 fListOfObjects->Add(fVtxBeforeCuts);
187
188 fVtxAfterCuts = new TH1F("fVtxAfterCuts", "Vtx distribution (before cuts); Vtx z [cm]; Counts", 120, -30, 30);
189 fListOfObjects->Add(fVtxAfterCuts);
190
191 if (fAnalysisMC) {
192 fVtxMC = new TH1I("fVtxMC","Vtx info - no trigger cut (0=no, 1=yes); Vtx; Counts", 2, -0.5, 1.5);
193 fListOfObjects->Add(fVtxMC);
194
195 }
196
197 if (fTreeOption) {
198
199 fTree = new TTree("tree","Event data");
200 fEvent = new DeDxEvent();
201 fTree->Branch("event", &fEvent);
202
203 fV0ArrayGlobalPar = new TClonesArray("DeDxV0", 1000);
204 fTree->Bronch("v0GlobalPar", "TClonesArray", &fV0ArrayGlobalPar);
205
206 fV0ArrayTPCPar = new TClonesArray("DeDxV0", 1000);
207 fTree->Bronch("v0TPCPar", "TClonesArray", &fV0ArrayTPCPar);
208
209 if (fAnalysisMC && fStoreMcIn) {
210 fTrackArrayMC = new TClonesArray("DeDxTrackMC", 1000);
211 fTree->Bronch("trackMC", "TClonesArray", &fTrackArrayMC);
212 }
213
214 fTree->SetDirectory(0);
215
216 fListOfObjects->Add(fTree);
217
218 }
219
220 // Post output data.
221 PostData(1, fListOfObjects);
222}
223
224//______________________________________________________________________________
225void AliAnalysisTaskHighPtDeDxV0::UserExec(Option_t *)
226{
227 // Main loop
228
229 //
230 // Comment: This method matches completely the same method for the high pT
231 // tracks
232 //
233
234
235 //
236 // First we make sure that we have valid input(s)!
237 //
238 if (fAnalysisType == "ESD"){
239 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
240 if(!fESD){
241 Printf("%s:%d ESDEvent not found in Input Manager",(char*)__FILE__,__LINE__);
242 this->Dump();
243 return;
244 }
245 } else {
246 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
247 if(!fAOD){
248 Printf("%s:%d AODEvent not found in Input Manager",(char*)__FILE__,__LINE__);
249 this->Dump();
250 return;
251 }
252 }
253
254 if (fAnalysisMC) {
255
256 if (fAnalysisType == "ESD"){
257 fMC = dynamic_cast<AliMCEvent*>(MCEvent());
258 if(!fMC){
259 Printf("%s:%d MCEvent not found in Input Manager",(char*)__FILE__,__LINE__);
260 this->Dump();
261 return;
262 }
263
264 fMCStack = fMC->Stack();
265
266 if(!fMCStack){
267 Printf("%s:%d MCStack not found in Input Manager",(char*)__FILE__,__LINE__);
268 this->Dump();
269 return;
270 }
271 } else { // AOD
272
273 fMC = dynamic_cast<AliMCEvent*>(MCEvent());
274 if(fMC)
275 fMC->Dump();
276
277 fMCArray = (TClonesArray*)fAOD->FindListObject("mcparticles");
278 if(!fMCArray){
279 Printf("%s:%d AOD MC array not found in Input Manager",(char*)__FILE__,__LINE__);
280 this->Dump();
281 return;
282 }
283 }
284 }
285
286 // Get trigger decision
287 fTriggeredEventMB = 0; //init
288
289
290
291 if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
292 ->IsEventSelected() & ftrigBit1 ){
293 fTriggeredEventMB = 1; //event triggered as minimum bias
294 }
295 if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
296 ->IsEventSelected() & ftrigBit2 ){
297 // From AliVEvent:
298 // kINT7 = BIT(1), // V0AND trigger, offline V0 selection
299 fTriggeredEventMB += 2;
300 }
301
302
303 // Get process type for MC
304 fMcProcessType = 0; // -1=invalid, 0=data, 1=ND, 2=SD, 3=DD
305
306 // real data that are not triggered we skip
307 if(!fAnalysisMC && !fTriggeredEventMB)
308 return;
309
310 if (fAnalysisMC) {
311
312 if (fAnalysisType == "ESD"){
313
314 AliHeader* headerMC = fMC->Header();
315 if (headerMC) {
316
317 AliGenEventHeader* genHeader = headerMC->GenEventHeader();
318 TArrayF vtxMC(3); // primary vertex MC
319 vtxMC[0]=9999; vtxMC[1]=9999; vtxMC[2]=9999; //initialize with dummy
320 if (genHeader) {
321 genHeader->PrimaryVertex(vtxMC);
322 }
323 fZvtxMC = vtxMC[2];
324
325 // PYTHIA:
326 AliGenPythiaEventHeader* pythiaGenHeader =
327 dynamic_cast<AliGenPythiaEventHeader*>(headerMC->GenEventHeader());
328 if (pythiaGenHeader) { //works only for pythia
329 fMcProcessType = GetPythiaEventProcessType(pythiaGenHeader->ProcessType());
330 }
331 // PHOJET:
332 AliGenDPMjetEventHeader* dpmJetGenHeader =
333 dynamic_cast<AliGenDPMjetEventHeader*>(headerMC->GenEventHeader());
334 if (dpmJetGenHeader) {
335 fMcProcessType = GetDPMjetEventProcessType(dpmJetGenHeader->ProcessType());
336 }
337 }
338 } else { // AOD
339
340 AliAODMCHeader* mcHeader = dynamic_cast<AliAODMCHeader*>(fAOD->FindListObject("mcHeader"));
341 if(mcHeader) {
342 fZvtxMC = mcHeader->GetVtxZ();
343
344 if(strstr(mcHeader->GetGeneratorName(), "Pythia")) {
345 fMcProcessType = GetPythiaEventProcessType(mcHeader->GetEventType());
346 } else {
347 fMcProcessType = GetDPMjetEventProcessType(mcHeader->GetEventType());
348 }
349 }
350 }
351 }
352
353 // There are 3 cases
354 // Vertex: NO - status -1
355 // Vertex: YES : outside cut - status 0
356 // : inside cut - status 1
357 // We have to be careful how we normalize because we probably want to
358 // normalize to:
359 // Nevents=(No vertex + outside + inside)/(out + in)*in
360
361 if (fAnalysisType == "ESD")
362 fZvtx = GetVertex(fESD);
363 else // AOD
364 fZvtx = GetVertex(fAOD);
365
366 fVtxStatus = -999;
367
368 if(fZvtx<-990) {
369
370 fVtxStatus = -1;
371 if(fTriggeredEventMB)
372 fVtx->Fill(0);
373 if(fAnalysisMC)
374 fVtxMC->Fill(0);
375 } else {
376
377 if(fTriggeredEventMB)
378 fVtx->Fill(1);
379 if(fAnalysisMC)
380 fVtxMC->Fill(1);
381 fVtxBeforeCuts->Fill(fZvtx);
382 fVtxStatus = 0;
383 if (TMath::Abs(fZvtx) < fVtxCut) {
384 fVtxAfterCuts->Fill(fZvtx);
385 fVtxStatus = 1;
386 }
387 }
388
389 // store MC event data no matter what
390 if(fAnalysisMC && fStoreMcIn) {
391
392 if (fAnalysisType == "ESD") {
393 ProcessMCTruthESD();
394 } else { // AOD
395 ProcessMCTruthAOD();
396 }
397 }
398 Float_t centrality = -10;
399 // only analyze triggered events
400 if(fTriggeredEventMB) {
401
402 if (fAnalysisType == "ESD"){
403 if(fAnalysisPbPb){
404 AliCentrality *centObject = fESD->GetCentrality();
405 centrality = centObject->GetCentralityPercentile("V0M");
406 if((centrality>fMaxCent)||(centrality<fMinCent))return;
407 }
408 Int_t nv0test=fESD->GetNumberOfV0s();
409 cout<<"&&&&&&&&&&&&&&&& hello world"<<endl;
410 fnv0->Fill(nv0test);
411 AnalyzeESD(fESD);
412
413 } else { // AOD
414 if(fAnalysisPbPb){
415 AliCentrality *centObject = fAOD->GetCentrality();
416 centrality = centObject->GetCentralityPercentile("V0M");
417 if((centrality>fMaxCent)||(centrality<fMinCent))return;
418 }
419 Int_t nv0test=fAOD->GetNumberOfV0s();
420 fnv0->Fill(nv0test);
421 AnalyzeAOD(fAOD);
422
423 }
424 }
425
426
427
428 if( fTreeOption) {
429
430 fEvent->process = fMcProcessType;
431 fEvent->trig = fTriggeredEventMB;
432 fEvent->zvtxMC = fZvtxMC;
433 fEvent->cent = centrality;
434
435 if(!fRequireRecV0 || fEvent->n > 0) // only fill if there are accepted V0s
436 fTree->Fill();
437 fV0ArrayGlobalPar->Clear();
438 fV0ArrayTPCPar->Clear();
439 if (fAnalysisMC && fStoreMcIn)
440 fTrackArrayMC->Clear();
441 }
442
443 // Post output data.
444 PostData(1, fListOfObjects);
445}
446
447//________________________________________________________________________
448void AliAnalysisTaskHighPtDeDxV0::AnalyzeESD(AliESDEvent* esdEvent)
449{
450 fRun = esdEvent->GetRunNumber();
451 fEventId = 0;
452 if(esdEvent->GetHeader())
453 fEventId = GetEventIdAsLong(esdEvent->GetHeader());
454
455 // Int_t event = esdEvent->GetEventNumberInFile();
456 UInt_t time = esdEvent->GetTimeStamp();
457 // ULong64_t trigger = esdEvent->GetTriggerMask();
458 Float_t magf = esdEvent->GetMagneticField();
459
460
461
462 Short_t isPileup = esdEvent->IsPileupFromSPD();
463
464 // centrality
465 Float_t centrality = 120;
466 AliCentrality *centObject = esdEvent->GetCentrality();
467 if(centObject)
468 centrality = centObject->GetCentralityPercentile("V0M");
469
470 if(fTriggeredEventMB) {// Only MC case can we have not triggered events
471
472 if(fVtxStatus==1) { // accepted vertex
473
474 // accepted event
475 fEvents->Fill(0);
476
477 ProduceArrayTrksESD( esdEvent, kGlobalTrk );//produce array with global track parameters
478 ProduceArrayTrksESD( esdEvent, kTPCTrk );//array with tpc track parametes
479
480
481 } // end if vtx status
482 } // end if triggered
483
484 if(fTreeOption) {
485
486 fEvent->run = fRun;
487 fEvent->eventid = fEventId;
488 fEvent->time = time;
489 //fEvent->cent = centrality;
490 fEvent->mag = magf;
491 fEvent->zvtx = fZvtx;
492 fEvent->vtxstatus = fVtxStatus;
493 //fEvent->trackmult = trackmult;
494 //fEvent->n = nadded;
495 fEvent->pileup = isPileup;
496 }
497}
498
499//________________________________________________________________________
500void AliAnalysisTaskHighPtDeDxV0::AnalyzeAOD(AliAODEvent* aodEvent)
501{
502 fRun = aodEvent->GetRunNumber();
503 fEventId = 0;
504 if(aodEvent->GetHeader())
505 fEventId = GetEventIdAsLong(aodEvent->GetHeader());
506
507 // Int_t event = aodEvent->GetEventNumberInFile();
508 UInt_t time = 0; // Missing AOD info? aodEvent->GetTimeStamp();
509 // ULong64_t trigger = aodEvent->GetTriggerMask();
510 Float_t magf = aodEvent->GetMagneticField();
511
512 Int_t trackmult = 0; // no pt cuts
513 Int_t nadded = 0;
514
515 Short_t isPileup = aodEvent->IsPileupFromSPD();
516
517 // centrality
518 Float_t centrality = 120;
519 AliCentrality *centObject = aodEvent->GetCentrality();
520 if(centObject)
521 centrality = centObject->GetCentralityPercentile("V0M");
522
523 if(fTriggeredEventMB) {// Only MC case can we have not triggered events
524
525 if(fVtxStatus==1) { // accepted vertex
526
527 // accepted event
528 fEvents->Fill(0);
529
530 ProduceArrayTrksAOD( aodEvent, kGlobalTrk );//produce array with global track parameters
531 ProduceArrayTrksAOD( aodEvent, kTPCTrk );//array with tpc track parametes
532
533
534 } // end if vtx status
535 } // end if triggered
536
537 if(fTreeOption) {
538
539 fEvent->run = fRun;
540 fEvent->eventid = fEventId;
541 fEvent->time = time;
542 //fEvent->cent = centrality;
543 fEvent->mag = magf;
544 fEvent->zvtx = fZvtx;
545 fEvent->vtxstatus = fVtxStatus;
546 //fEvent->trackmult = trackmult;
547 //fEvent->n = nadded;
548 fEvent->pileup = isPileup;
549 }
550
551
552
553
554
555
556
557}
558
559//_____________________________________________________________________________
560Float_t AliAnalysisTaskHighPtDeDxV0::GetVertex(const AliVEvent* event) const
561{
562 Float_t zvtx = -999;
563
564 const AliVVertex* primaryVertex = event->GetPrimaryVertex();
565
566 if(primaryVertex->GetNContributors()>0)
567 zvtx = primaryVertex->GetZ();
568
569 return zvtx;
570}
571
572//_____________________________________________________________________________
573Short_t AliAnalysisTaskHighPtDeDxV0::GetPidCode(Int_t pdgCode) const
574{
575 // return our internal code for pions, kaons, and protons
576
577 Short_t pidCode = 6;
578
579 switch (TMath::Abs(pdgCode)) {
580 case 211:
581 pidCode = 1; // pion
582 break;
583 case 321:
584 pidCode = 2; // kaon
585 break;
586 case 2212:
587 pidCode = 3; // proton
588 break;
589 case 11:
590 pidCode = 4; // electron
591 break;
592 case 13:
593 pidCode = 5; // proton
594 break;
595 default:
596 pidCode = 6; // something else?
597 };
598
599 return pidCode;
600}
601
602
603//_____________________________________________________________________________
604void AliAnalysisTaskHighPtDeDxV0::ProcessMCTruthESD()
605{
606 // Fill the special MC histogram with the MC truth info
607
608 Short_t trackmult = 0;
609 Short_t nadded = 0;
610 const Int_t nTracksMC = fMCStack->GetNtrack();
611
612 for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
613
614 //Cuts
615 if(!(fMCStack->IsPhysicalPrimary(iTracks)))
616 continue;
617
618 TParticle* trackMC = fMCStack->Particle(iTracks);
619
620 Double_t chargeMC = trackMC->GetPDG()->Charge();
621 if (chargeMC != 0) // select charge 0 particles!
622 continue;
623
624 if (TMath::Abs(trackMC->Eta()) > fEtaCut )
625 continue;
626
627 trackmult++;
628
629 // "charge:pt:p:eta:phi:pidCode"
630 Float_t ptMC = trackMC->Pt();
631 Float_t pMC = trackMC->P();
632 Float_t etaMC = trackMC->Eta();
633 Float_t phiMC = trackMC->Phi();
634
635 Int_t pdgCode = trackMC->GetPdgCode();
636 if(!(pdgCode==310 || pdgCode == 3122 || pdgCode == -3122))
637 continue;
638 Short_t pidCodeMC = 0;
639 pidCodeMC = GetPidCode(pdgCode);
640
641 // Warning: In the data code we cut on the daughters.....
642 if (trackMC->Pt() < fMinPt)
643 continue;
644
645 if(fTreeOption) {
646
647 DeDxTrackMC* track = new((*fTrackArrayMC)[nadded]) DeDxTrackMC();
648 nadded++;
649
650 track->pMC = pMC;
651 track->ptMC = ptMC;
652 track->etaMC = etaMC;
653 track->phiMC = phiMC;
654 track->qMC = Short_t(chargeMC);
655 track->pidMC = pidCodeMC;
656 track->pdgMC = pdgCode;
657 }
658
659 }//MC track loop
660
661 if(fTreeOption) {
662
663 Sort(fTrackArrayMC, kTRUE);
664
665 fEvent->trackmultMC = trackmult;
666 fEvent->nMC = nadded;
667 }
668
669}
670
671//_____________________________________________________________________________
672void AliAnalysisTaskHighPtDeDxV0::ProcessMCTruthAOD()
673{
674 // Fill the special MC histogram with the MC truth info
675
676 Short_t trackmult = 0;
677 Short_t nadded = 0;
678 const Int_t nTracksMC = fMCArray->GetEntriesFast();
679
680 for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
681
682 AliAODMCParticle* trackMC = dynamic_cast<AliAODMCParticle*>(fMCArray->At(iTracks));
683
684 //Cuts
685 if(!(trackMC->IsPhysicalPrimary()))
686 continue;
687
688
689 Double_t chargeMC = trackMC->Charge();
690 if (chargeMC != 0) // Keep the neutral particles
691 continue;
692
693 if (TMath::Abs(trackMC->Eta()) > fEtaCut )
694 continue;
695
696 trackmult++;
697
698 // "charge:pt:p:eta:phi:pidCode"
699 Float_t ptMC = trackMC->Pt();
700 Float_t pMC = trackMC->P();
701 Float_t etaMC = trackMC->Eta();
702 Float_t phiMC = trackMC->Phi();
703
704 Int_t pdgCode = trackMC->PdgCode();
705 if(!(pdgCode==310 || pdgCode == 3122 || pdgCode == -3122))
706 continue;
707 Short_t pidCodeMC = 0;
708 pidCodeMC = GetPidCode(pdgCode);
709
710 if (trackMC->Pt() < fMinPt)
711 continue;
712
713 if(fTreeOption) {
714
715 DeDxTrackMC* track = new((*fTrackArrayMC)[nadded]) DeDxTrackMC();
716 nadded++;
717
718 track->pMC = pMC;
719 track->ptMC = ptMC;
720 track->etaMC = etaMC;
721 track->phiMC = phiMC;
722 track->qMC = Short_t(chargeMC);
723 track->pidMC = pidCodeMC;
724 track->pdgMC = pdgCode;
725 }
726
727 }//MC track loop
728
729 if(fTreeOption) {
730
731 Sort(fTrackArrayMC, kTRUE);
732
733 fEvent->trackmultMC = trackmult;
734 fEvent->nMC = nadded;
735 }
736
737}
738
739//_____________________________________________________________________________
740Short_t AliAnalysisTaskHighPtDeDxV0::GetPythiaEventProcessType(Int_t pythiaType) {
741 //
742 // Get the process type of the event. PYTHIA
743 //
744 // source PWG0 dNdpt
745
746 Short_t globalType = -1; //init
747
748 if(pythiaType==92||pythiaType==93){
749 globalType = 2; //single diffractive
750 }
751 else if(pythiaType==94){
752 globalType = 3; //double diffractive
753 }
754 //else if(pythiaType != 91){ // also exclude elastic to be sure... CKB??
755 else {
756 globalType = 1; //non diffractive
757 }
758 return globalType;
759}
760
761//_____________________________________________________________________________
762Short_t AliAnalysisTaskHighPtDeDxV0::GetDPMjetEventProcessType(Int_t dpmJetType) {
763 //
764 // get the process type of the event. PHOJET
765 //
766 //source PWG0 dNdpt
767 // can only read pythia headers, either directly or from cocktalil header
768 Short_t globalType = -1;
769
770 if (dpmJetType == 1 || dpmJetType == 4) { // explicitly inelastic plus central diffraction
771 globalType = 1;
772 }
773 else if (dpmJetType==5 || dpmJetType==6) {
774 globalType = 2;
775 }
776 else if (dpmJetType==7) {
777 globalType = 3;
778 }
779 return globalType;
780}
781
782//_____________________________________________________________________________
783ULong64_t AliAnalysisTaskHighPtDeDxV0::GetEventIdAsLong(AliVHeader* header) const
784{
785 // To have a unique id for each event in a run!
786 // Modified from AliRawReader.h
787 return ((ULong64_t)header->GetBunchCrossNumber()+
788 (ULong64_t)header->GetOrbitNumber()*3564+
789 (ULong64_t)header->GetPeriodNumber()*16777215*3564);
790}
791
792
793//____________________________________________________________________
794TParticle* AliAnalysisTaskHighPtDeDxV0::FindPrimaryMother(AliStack* stack, Int_t label)
795{
796 //
797 // Finds the first mother among the primary particles of the particle identified by <label>,
798 // i.e. the primary that "caused" this particle
799 //
800 // Taken from AliPWG0Helper class
801 //
802
803 Int_t nSteps = 0;
804
805 Int_t motherLabel = FindPrimaryMotherLabel(stack, label, nSteps);
806 if (motherLabel < 0)
807 return 0;
808
809 return stack->Particle(motherLabel);
810}
811
812//____________________________________________________________________
813Int_t AliAnalysisTaskHighPtDeDxV0::FindPrimaryMotherLabel(AliStack* stack, Int_t label, Int_t& nSteps)
814{
815 //
816 // Finds the first mother among the primary particles of the particle identified by <label>,
817 // i.e. the primary that "caused" this particle
818 //
819 // returns its label
820 //
821 // Taken from AliPWG0Helper class
822 //
823 nSteps = 0;
824 const Int_t nPrim = stack->GetNprimary();
825
826 while (label >= nPrim) {
827
828 //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0));
829
830 nSteps++; // 1 level down
831
832 TParticle* particle = stack->Particle(label);
833 if (!particle) {
834
835 AliDebugGeneral("FindPrimaryMotherLabel", AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack.", label));
836 return -1;
837 }
838
839 // find mother
840 if (particle->GetMother(0) < 0) {
841
842 AliDebugGeneral("FindPrimaryMotherLabel", AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label));
843 return -1;
844 }
845
846 label = particle->GetMother(0);
847 }
848
849 return label;
850}
851
852//____________________________________________________________________
853AliAODMCParticle* AliAnalysisTaskHighPtDeDxV0::FindPrimaryMotherAOD(AliAODMCParticle* startParticle, Int_t& nSteps)
854{
855 //
856 // Finds the first mother among the primary particles of the particle identified by <label>,
857 // i.e. the primary that "caused" this particle
858 //
859 // Taken from AliPWG0Helper class
860 //
861
862 nSteps = 0;
863
864 AliAODMCParticle* mcPart = startParticle;
865
866 while (mcPart)
867 {
868
869 if(mcPart->IsPrimary())
870 return mcPart;
871
872 Int_t mother = mcPart->GetMother();
873
874 mcPart = dynamic_cast<AliAODMCParticle*>(fMCArray->At(mother));
875 nSteps++; // 1 level down
876 }
877
878 return 0;
879}
880
881//_____________________________________________________________________________
882void AliAnalysisTaskHighPtDeDxV0::Sort(TClonesArray* array, Bool_t isMC)
883{
884 const Int_t n = array->GetEntriesFast();
885 if(n==0) {
886
887 if(isMC)
888 fEvent->ptmaxMC = 0;
889 else
890 fEvent->ptmax = 0;
891
892 return;
893 }
894
895 Float_t ptArray[n];
896 Int_t indexArray[n];
897
898 for(Int_t i = 0; i < n; i++) {
899
900 Float_t pt = 0;
901 if(isMC) {
902 DeDxTrackMC* track = (DeDxTrackMC*)array->At(i);
903 pt = track->ptMC;
904 } else {
905 DeDxTrack* track = (DeDxTrack*)array->At(i);
906 pt = track->pt;
907 }
908 ptArray[i] = pt;
909 indexArray[i] = i;
910 }
911
912 TMath::Sort(n, ptArray, indexArray);
913
914 // set max pt
915 if(isMC) {
916 DeDxTrackMC* track = (DeDxTrackMC*)array->At(indexArray[0]);
917 fEvent->ptmaxMC = track->ptMC;
918 } else {
919 DeDxTrack* track = (DeDxTrack*)array->At(indexArray[0]);
920 fEvent->ptmax = track->pt;
921 }
922
923 // set order of each track
924 for(Int_t i = 0; i < n; i++) {
925
926 if(isMC) {
927 DeDxTrackMC* track = (DeDxTrackMC*)array->At(indexArray[i]);
928 track->orderMC = i;
929 } else {
930 DeDxTrack* track = (DeDxTrack*)array->At(indexArray[i]);
931 track->order = i;
932 }
933 }
934}
935//__________________________________________________________________
936void AliAnalysisTaskHighPtDeDxV0::ProduceArrayTrksESD( AliESDEvent *ESDevent, AnalysisMode analysisMode ){
937 Int_t nv0s = ESDevent->GetNumberOfV0s();
938 if(nv0s<1)return;
939 Int_t trackmult = 0; // no pt cuts
940 Int_t nadded = 0;
941 const AliESDVertex *myBestPrimaryVertex = ESDevent->GetPrimaryVertex();
942 if (!myBestPrimaryVertex) return;
943 if (!(myBestPrimaryVertex->GetStatus())) return;
944 Double_t lPrimaryVtxPosition[3];
945 myBestPrimaryVertex->GetXYZ(lPrimaryVtxPosition);
946
947 Double_t lPrimaryVtxCov[6];
948 myBestPrimaryVertex->GetCovMatrix(lPrimaryVtxCov);
949 Double_t lPrimaryVtxChi2 = myBestPrimaryVertex->GetChi2toNDF();
950
951 AliAODVertex* myPrimaryVertex = new AliAODVertex(lPrimaryVtxPosition, lPrimaryVtxCov, lPrimaryVtxChi2, NULL, -1, AliAODVertex::kPrimary);
952
953
954 if( analysisMode == kGlobalTrk ){
955 if(fV0ArrayGlobalPar)
956 fV0ArrayGlobalPar->Clear();
957 } else if( analysisMode == kTPCTrk ){
958 if(fV0ArrayTPCPar)
959 fV0ArrayTPCPar->Clear();
960 }
961
962 if( analysisMode==kGlobalTrk ){
963
964
965
966
967 // ################################
968 // #### BEGINNING OF V0 CODE ######
969 // ################################
970
971
972 for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {
973
974 // This is the begining of the V0 loop
975 AliESDv0 *esdV0 = ESDevent->GetV0(iV0);
976 if (!esdV0) continue;
977
978 // AliESDTrack (V0 Daughters)
979 UInt_t lKeyPos = (UInt_t)TMath::Abs(esdV0->GetPindex());
980 UInt_t lKeyNeg = (UInt_t)TMath::Abs(esdV0->GetNindex());
981
982 AliESDtrack *pTrack = ESDevent->GetTrack(lKeyPos);
983 AliESDtrack *nTrack = ESDevent->GetTrack(lKeyNeg);
984 if (!pTrack || !nTrack) {
985 Printf("ERROR: Could not retreive one of the daughter track");
986 continue;
987 }
988
989 // Remove like-sign
990 if (pTrack->GetSign() == nTrack->GetSign()) {
991 //cout<< "like sign, continue"<< endl;
992 continue;
993 }
994
995 // Eta cut on decay products
996 if(TMath::Abs(pTrack->Eta()) > fEtaCut || TMath::Abs(nTrack->Eta()) > fEtaCut)
997 continue;
998
999 // Pt cut on decay products
1000 if (esdV0->Pt() < fMinPt)
1001 // if (pTrack->Pt() < fMinPt && nTrack->Pt() < fMinPt)
1002 continue;
1003
1004 //filter for positive track
1005 UShort_t filterFlag_p = 0;
1006 Bool_t filterCut_Set1_p = kFALSE;//parameters from global tracks, with TPC cuts (filter bit =1 in AOD)
1007 Bool_t filterCut_Set2_p = kFALSE;//parameters from global tracks, cuts tpc+its 2010 W/O golden cuts
1008 Bool_t filterCut_Set3_p = kFALSE;//parameters from global tracks, cuts its+tpc 2010 WITH golden cuts
1009
1010 UInt_t selectDebug_p = 0;
1011 if (fTrackFilterGolden) {
1012 selectDebug_p = fTrackFilterGolden->IsSelected(pTrack);
1013 if (selectDebug_p) {
1014 filterFlag_p +=1;
1015 filterCut_Set3_p=kTRUE;
1016 }
1017 }
1018
1019 if (fTrackFilterTPC) {
1020
1021 selectDebug_p = fTrackFilterTPC->IsSelected(pTrack);
1022 if (selectDebug_p){//only tracks which pass the TPC-only track cuts
1023 filterFlag_p +=2;
1024 filterCut_Set1_p=kTRUE;
1025
1026 }
1027
1028 }
1029
1030 if (fTrackFilter) {
1031 selectDebug_p = fTrackFilter->IsSelected(pTrack);
1032 if (selectDebug_p) {
1033 filterCut_Set2_p=kTRUE;
1034 }
1035 }
1036
1037 if(filterFlag_p ==0 )
1038 continue;
1039
1040
1041 //filter for negative track
1042 UShort_t filterFlag_n = 0;
1043 Bool_t filterCut_Set1_n = kFALSE;//parameters from global tracks, with TPC cuts (filter bit =1 in AOD)
1044 Bool_t filterCut_Set2_n = kFALSE;//parameters from global tracks, cuts tpc+its 2010 W/O golden cuts
1045 Bool_t filterCut_Set3_n = kFALSE;//parameters from global tracks, cuts its+tpc 2010 WITH golden cuts
1046
1047 UInt_t selectDebug_n = 0;
1048 if (fTrackFilterGolden) {
1049 selectDebug_n = fTrackFilterGolden->IsSelected(nTrack);
1050 if (selectDebug_n) {
1051 filterFlag_n +=1;
1052 filterCut_Set3_n=kTRUE;
1053 }
1054 }
1055
1056 if (fTrackFilterTPC) {
1057
1058 selectDebug_n = fTrackFilterTPC->IsSelected(nTrack);
1059 if (selectDebug_n){//only tracks which pass the TPC-only track cuts
1060 filterFlag_n +=2;
1061 filterCut_Set1_n=kTRUE;
1062
1063 }
1064
1065 }
1066
1067 if (fTrackFilter) {
1068 selectDebug_n = fTrackFilter->IsSelected(nTrack);
1069 if (selectDebug_n) {
1070 filterCut_Set2_n=kTRUE;
1071 }
1072 }
1073
1074
1075 // Check if switch does anything!
1076 Bool_t isSwitched = kFALSE;
1077 if (pTrack->GetSign() < 0) { // switch
1078
1079 isSwitched = kTRUE;
1080 AliESDtrack* helpTrack = nTrack;
1081 nTrack = pTrack;
1082 pTrack = helpTrack;
1083 }
1084
1085 Float_t alpha = esdV0->AlphaV0();
1086 Float_t ptarm = esdV0->PtArmV0();
1087 // Double_t pVtxPos= v0->PrimaryVtxPosition();
1088
1089 Double_t lV0Position[3];
1090 esdV0->GetXYZ(lV0Position[0], lV0Position[1], lV0Position[2]);
1091
1092 Double_t lV0Radius = TMath::Sqrt(lV0Position[0]*lV0Position[0]+lV0Position[1]*lV0Position[1]);
1093 Double_t lV0DecayLength = TMath::Sqrt(TMath::Power(lV0Position[0] - lPrimaryVtxPosition[0],2) +
1094 TMath::Power(lV0Position[1] - lPrimaryVtxPosition[1],2) +
1095 TMath::Power(lV0Position[2] - lPrimaryVtxPosition[2],2 ));
1096 AliKFVertex primaryVtxKF( *myPrimaryVertex );
1097 AliKFParticle::SetField(ESDevent->GetMagneticField());
1098
1099 // Also implement switch here!!!!!!
1100 AliKFParticle* negEKF = 0; // e-
1101 AliKFParticle* posEKF = 0; // e+
1102 AliKFParticle* negPiKF = 0; // pi -
1103 AliKFParticle* posPiKF = 0; // pi +
1104 AliKFParticle* posPKF = 0; // p
1105 AliKFParticle* negAPKF = 0; // p-bar
1106
1107 if(!isSwitched) {
1108 negEKF = new AliKFParticle( *(esdV0->GetParamN()) , 11);
1109 posEKF = new AliKFParticle( *(esdV0->GetParamP()) ,-11);
1110 negPiKF = new AliKFParticle( *(esdV0->GetParamN()) ,-211);
1111 posPiKF = new AliKFParticle( *(esdV0->GetParamP()) , 211);
1112 posPKF = new AliKFParticle( *(esdV0->GetParamP()) , 2212);
1113 negAPKF = new AliKFParticle( *(esdV0->GetParamN()) ,-2212);
1114 } else { // switch + and -
1115 negEKF = new AliKFParticle( *(esdV0->GetParamP()) , 11);
1116 posEKF = new AliKFParticle( *(esdV0->GetParamN()) ,-11);
1117 negPiKF = new AliKFParticle( *(esdV0->GetParamP()) ,-211);
1118 posPiKF = new AliKFParticle( *(esdV0->GetParamN()) , 211);
1119 posPKF = new AliKFParticle( *(esdV0->GetParamN()) , 2212);
1120 negAPKF = new AliKFParticle( *(esdV0->GetParamP()) ,-2212);
1121 }
1122
1123 AliKFParticle v0GKF; // Gamma e.g. from pi0
1124 v0GKF+=(*negEKF);
1125 v0GKF+=(*posEKF);
1126 v0GKF.SetProductionVertex(primaryVtxKF);
1127
1128 AliKFParticle v0K0sKF; // K0 short
1129 v0K0sKF+=(*negPiKF);
1130 v0K0sKF+=(*posPiKF);
1131 v0K0sKF.SetProductionVertex(primaryVtxKF);
1132
1133 AliKFParticle v0LambdaKF; // Lambda
1134 v0LambdaKF+=(*negPiKF);
1135 v0LambdaKF+=(*posPKF);
1136 v0LambdaKF.SetProductionVertex(primaryVtxKF);
1137
1138 AliKFParticle v0AntiLambdaKF; // Lambda-bar
1139 v0AntiLambdaKF+=(*posPiKF);
1140 v0AntiLambdaKF+=(*negAPKF);
1141 v0AntiLambdaKF.SetProductionVertex(primaryVtxKF);
1142
1143 Double_t deltaInvMassG = v0GKF.GetMass();
1144 Double_t deltaInvMassK0s = v0K0sKF.GetMass()-0.498;
1145 Double_t deltaInvMassL = v0LambdaKF.GetMass()-1.116;
1146 Double_t deltaInvMassAntiL = v0AntiLambdaKF.GetMass()-1.116;
1147
1148 if(TMath::Abs(deltaInvMassK0s) > fMassCut &&
1149 TMath::Abs(deltaInvMassL) > fMassCut &&
1150 TMath::Abs(deltaInvMassAntiL) > fMassCut)
1151 continue;
1152
1153 // Extract track information
1154
1155 Short_t pcharge = pTrack->Charge();
1156 Float_t ppt = pTrack->Pt();
1157 Float_t pp = pTrack->P();
1158 Float_t peta = pTrack->Eta();
1159 Float_t pphi = pTrack->Phi();
1160 Short_t pncl = pTrack->GetTPCsignalN();
1161 Short_t pneff = Short_t(pTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
1162 Float_t pdedx = pTrack->GetTPCsignal();
1163
1164 Float_t ptpcchi = 0;
1165 if(pTrack->GetTPCNcls() > 0)
1166 ptpcchi = pTrack->GetTPCchi2()/Float_t(pTrack->GetTPCNcls());
1167
1168 Short_t ncharge = nTrack->Charge();
1169 Float_t npt = nTrack->Pt();
1170 Float_t np = nTrack->P();
1171 Float_t neta = nTrack->Eta();
1172 Float_t nphi = nTrack->Phi();
1173 Short_t nncl = nTrack->GetTPCsignalN();
1174 Short_t nneff = Short_t(nTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
1175 Float_t ndedx = nTrack->GetTPCsignal();
1176 Float_t ntpcchi = 0;
1177 if(nTrack->GetTPCNcls() > 0)
1178 ntpcchi = nTrack->GetTPCchi2()/Float_t(nTrack->GetTPCNcls());
1179
1180 Float_t b[2];
1181 Float_t bCov[3];
1182 pTrack->GetImpactParameters(b,bCov);
1183 Float_t pdcaxy = b[0];
1184 Float_t pdcaz = b[1];
1185 nTrack->GetImpactParameters(b,bCov);
1186 Float_t ndcaxy = b[0];
1187 Float_t ndcaz = b[1];
1188
1189 Int_t primaryV0 = 0; // 0 means that the tracks are not both daughters of a primary particle (1 means they are)
1190 Int_t pdgV0 = 0; // 0 means that they don't have same origin for MC (1 means they have the same original mother)
1191 Float_t p_ptMC = 0;
1192 Short_t p_pidCode = 0; // 0 = real data / no mc track!
1193 Short_t p_primaryFlag = 0; // 0 = real data / not primary mc track
1194 Int_t p_pdgMother = 0;
1195 Float_t n_ptMC = 0;
1196 Short_t n_pidCode = 0; // 0 = real data / no mc track!
1197 Short_t n_primaryFlag = 0; // 0 = real data / not primary mc track
1198 Int_t n_pdgMother = 0;
1199 if(fAnalysisMC) {
1200
1201 Int_t p_mother_label = 0;
1202 Int_t p_mother_steps = 0;
1203 Int_t n_mother_label = 0;
1204 Int_t n_mother_steps = 0;
1205
1206 // positive track
1207 const Int_t p_label = TMath::Abs(pTrack->GetLabel());
1208 TParticle* p_mcTrack = fMCStack->Particle(p_label);
1209 if (p_mcTrack){
1210
1211 if(fMCStack->IsPhysicalPrimary(p_label))
1212 p_primaryFlag = 1;
1213
1214 Int_t p_pdgCode = p_mcTrack->GetPdgCode();
1215 p_pidCode = GetPidCode(p_pdgCode);
1216
1217 p_ptMC = p_mcTrack->Pt();
1218
1219 p_mother_label = FindPrimaryMotherLabel(fMCStack, p_label,
1220 p_mother_steps);
1221 if(p_mother_label>0) {
1222 TParticle* p_mother = fMCStack->Particle(p_mother_label);
1223 p_pdgMother = p_mother->GetPdgCode();
1224 }
1225 }
1226
1227 // negative track
1228 const Int_t n_label = TMath::Abs(nTrack->GetLabel());
1229 TParticle* n_mcTrack = fMCStack->Particle(n_label);
1230 if (n_mcTrack){
1231
1232 if(fMCStack->IsPhysicalPrimary(n_label))
1233 n_primaryFlag = 1;
1234
1235 Int_t n_pdgCode = n_mcTrack->GetPdgCode();
1236 n_pidCode = GetPidCode(n_pdgCode);
1237
1238 n_ptMC = n_mcTrack->Pt();
1239
1240 n_mother_label = FindPrimaryMotherLabel(fMCStack, n_label,
1241 n_mother_steps);
1242 if(n_mother_label>0) {
1243 TParticle* n_mother = fMCStack->Particle(n_mother_label);
1244 n_pdgMother = n_mother->GetPdgCode();
1245 }
1246 }
1247
1248 // Check if V0 is primary = first and the same mother of both partciles
1249 if(p_mother_label>0 && n_mother_label>0 && p_mother_label == n_mother_label) {
1250 pdgV0 = p_pdgMother;
1251 if(p_mother_steps == 1 && n_mother_steps == 1) {
1252 primaryV0 = 1;
1253 }
1254 }
1255 }
1256
1257
1258
1259 if(fTreeOption) {
1260
1261
1262
1263 DeDxV0* v0data = new((*fV0ArrayGlobalPar)[nadded]) DeDxV0();
1264 nadded++;
1265
1266 // v0 data
1267 v0data->p = esdV0->P();
1268 v0data->pt = esdV0->Pt();
1269 v0data->eta = esdV0->Eta();
1270 v0data->phi = esdV0->Phi();
1271 v0data->pdca = TMath::Sqrt(pdcaxy*pdcaxy + pdcaz*pdcaz);
1272 v0data->ndca = TMath::Sqrt(ndcaxy*ndcaxy + ndcaz*ndcaz);
1273 v0data->dmassG = deltaInvMassG;
1274 v0data->dmassK0 = deltaInvMassK0s;
1275 v0data->dmassL = deltaInvMassL;
1276 v0data->dmassAL = deltaInvMassAntiL;
1277 v0data->alpha = alpha;
1278 v0data->ptarm = ptarm;
1279 v0data->decayr = lV0Radius;
1280 v0data->decayl = lV0DecayLength;
1281
1282 // New parameters
1283 v0data->status = esdV0->GetOnFlyStatus();
1284 v0data->chi2 = esdV0->GetChi2V0();
1285 v0data->cospt = esdV0->GetV0CosineOfPointingAngle();
1286 // cospt: as I understand this means that the pointing to the vertex
1287 // is fine so I remove the dcaxy and dcaz for the V= class
1288 v0data->dcadaughters = esdV0->GetDcaV0Daughters();
1289 v0data->primary = primaryV0;
1290 v0data->pdg = pdgV0;
1291
1292 // positive track
1293 v0data->ptrack.p = pp;
1294 v0data->ptrack.pt = ppt;
1295 // v0data->ptrack.ptcon = ppt_con;
1296 // v0data->ptrack.tpcchi = ptpcchi;
1297 v0data->ptrack.eta = peta;
1298 v0data->ptrack.phi = pphi;
1299 v0data->ptrack.q = pcharge;
1300 v0data->ptrack.ncl = pncl;
1301 v0data->ptrack.neff = pneff;
1302 v0data->ptrack.dedx = pdedx;
1303 v0data->ptrack.dcaxy = pdcaxy;
1304 v0data->ptrack.dcaz = pdcaz;
1305 v0data->ptrack.pid = p_pidCode;
1306 v0data->ptrack.primary = p_primaryFlag;
1307 v0data->ptrack.pttrue = p_ptMC;
1308 v0data->ptrack.mother = p_pdgMother;
1309 v0data->ptrack.filter = filterFlag_p;
1310 v0data->ptrack.filterset1 = filterCut_Set1_p;
1311 v0data->ptrack.filterset2 = filterCut_Set2_p;
1312 v0data->ptrack.filterset3 = filterCut_Set3_p;
1313
1314 // negative track
1315 v0data->ntrack.p = np;
1316 v0data->ntrack.pt = npt;
1317 // v0data->ntrack.ptcon = npt_con;
1318 // v0data->ntrack.tpcchi = ntpcchi;
1319 v0data->ntrack.eta = neta;
1320 v0data->ntrack.phi = nphi;
1321 v0data->ntrack.q = ncharge;
1322 v0data->ntrack.ncl = nncl;
1323 v0data->ntrack.neff = nneff;
1324 v0data->ntrack.dedx = ndedx;
1325 v0data->ntrack.dcaxy = ndcaxy;
1326 v0data->ntrack.dcaz = ndcaz;
1327 v0data->ntrack.pid = n_pidCode;
1328 v0data->ntrack.primary = n_primaryFlag;
1329 v0data->ntrack.pttrue = n_ptMC;
1330 v0data->ntrack.mother = n_pdgMother;
1331 v0data->ntrack.filter = filterFlag_n;
1332 v0data->ntrack.filterset1 = filterCut_Set1_n;
1333 v0data->ntrack.filterset2 = filterCut_Set2_n;
1334 v0data->ntrack.filterset3 = filterCut_Set3_n;
1335
1336
1337 }
1338
1339 // clean up loop over v0
1340
1341 delete negPiKF;
1342 delete posPiKF;
1343 delete posPKF;
1344 delete negAPKF;
1345
1346
1347
1348 }
1349
1350 // clean up event
1351 //delete myPrimaryVertex;
1352 }
1353 else if( analysisMode==kTPCTrk ){
1354 cout<<"&&&&&&&&&&&&&&&&&&&&& Hello world"<<endl;
1355 const AliESDVertex *vtxSPD = ESDevent->GetPrimaryVertexSPD();
1356 if( vtxSPD->GetNContributors() < 1 || TMath::Abs(vtxSPD->GetZ()) > 10.0 ) return;
1357
1358
1359 // ################################
1360 // #### BEGINNING OF V0 CODE ######
1361 // ################################
1362
1363
1364 for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {
1365
1366 // This is the begining of the V0 loop
1367 AliESDv0 *esdV0 = ESDevent->GetV0(iV0);
1368 if (!esdV0) continue;
1369
1370 // AliESDTrack (V0 Daughters)
1371 UInt_t lKeyPos = (UInt_t)TMath::Abs(esdV0->GetPindex());
1372 UInt_t lKeyNeg = (UInt_t)TMath::Abs(esdV0->GetNindex());
1373
1374 AliESDtrack *pTrack = ESDevent->GetTrack(lKeyPos);
1375 AliESDtrack *nTrack = ESDevent->GetTrack(lKeyNeg);
1376 if (!pTrack || !nTrack) {
1377 Printf("ERROR: Could not retreive one of the daughter track");
1378 continue;
1379 }
1380
1381 // Remove like-sign
1382 if (pTrack->GetSign() == nTrack->GetSign()) {
1383 //cout<< "like sign, continue"<< endl;
1384 continue;
1385 }
1386
1387 AliESDtrack *pTrackTPC = AliESDtrackCuts::GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(ESDevent),pTrack->GetID());
1388 AliESDtrack *nTrackTPC = AliESDtrackCuts::GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(ESDevent),nTrack->GetID());
1389
1390 if (!pTrackTPC || !nTrackTPC) {
1391 Printf("ERROR: Could not retreive one of the daughter TPC track");
1392 continue;
1393 }
1394
1395 //filter for positive track
1396 UInt_t selectDebug_p = 0;
1397 UShort_t filterFlag_p = 0;
1398 selectDebug_p = fTrackFilterTPC->IsSelected(pTrackTPC);
1399 //if(selectDebug_p==0) continue;
1400
1401 //filter for negative track
1402 UInt_t selectDebug_n = 0;
1403 UShort_t filterFlag_n = 0;
1404 selectDebug_n = fTrackFilterTPC->IsSelected(nTrackTPC);
1405 if(selectDebug_n==0 || selectDebug_p==0) continue;
1406
1407 if(selectDebug_p)
1408 filterFlag_p += 1;
1409
1410 if(pTrackTPC->Pt()>0.){
1411 // only constrain tracks above threshold
1412 AliExternalTrackParam exParamp;
1413 // take the B-field from the ESD, no 3D fieldMap available at this point
1414 Bool_t relate_p = false;
1415 relate_p = pTrackTPC->RelateToVertexTPC(vtxSPD,ESDevent->GetMagneticField(),
1416 kVeryBig,&exParamp);
1417 if(!relate_p){
1418 delete pTrackTPC;
1419 continue;
1420 }
1421 pTrackTPC->Set(exParamp.GetX(),exParamp.GetAlpha(),exParamp.GetParameter(),
1422 exParamp.GetCovariance());
1423 }
1424 else continue;
1425
1426
1427 //filter for negative track
1428
1429 if(selectDebug_n)
1430 filterFlag_n += 1;
1431
1432
1433 if(nTrackTPC->Pt()>0.){
1434 // only constrain tracks above threshold
1435 AliExternalTrackParam exParamn;
1436 // take the B-field from the ESD, no 3D fieldMap available at this point
1437 Bool_t relate_n = false;
1438 relate_n = nTrackTPC->RelateToVertexTPC(vtxSPD,ESDevent->GetMagneticField(),
1439 kVeryBig,&exParamn);
1440 if(!relate_n){
1441 delete nTrackTPC;
1442 continue;
1443 }
1444 nTrackTPC->Set(exParamn.GetX(),exParamn.GetAlpha(),exParamn.GetParameter(),
1445 exParamn.GetCovariance());
1446 }
1447 else continue;
1448
1449
1450 // Eta cut on decay products
1451 if(TMath::Abs(pTrackTPC->Eta()) > fEtaCut || TMath::Abs(nTrackTPC->Eta()) > fEtaCut)
1452 continue;
1453
1454 // Pt cut on decay products
1455 if (esdV0->Pt() < fMinPt)
1456 // if (pTrack->Pt() < fMinPt && nTrack->Pt() < fMinPt)
1457 continue;
1458
1459
1460 // Check if switch does anything!
1461 Bool_t isSwitched = kFALSE;
1462 if (pTrackTPC->GetSign() < 0) { // switch
1463
1464 isSwitched = kTRUE;
1465 AliESDtrack* helpTrack = nTrack;
1466 nTrackTPC = pTrackTPC;
1467 pTrackTPC = helpTrack;
1468 }
1469
1470
1471 // Extract track information
1472
1473 Short_t pcharge = pTrackTPC->Charge();
1474 Float_t ppt = pTrackTPC->Pt();
1475 Float_t pp = pTrackTPC->P();
1476 Float_t peta = pTrackTPC->Eta();
1477 Float_t pphi = pTrackTPC->Phi();
1478 Short_t pncl = pTrackTPC->GetTPCsignalN();
1479 Short_t pneff = Short_t(pTrackTPC->GetTPCClusterInfo(2, 1)); // effective track length for pT res
1480 Float_t pdedx = pTrackTPC->GetTPCsignal();
1481
1482 Float_t ptpcchi = 0;
1483 if(pTrackTPC->GetTPCNcls() > 0)
1484 ptpcchi = pTrackTPC->GetTPCchi2()/Float_t(pTrackTPC->GetTPCNcls());
1485
1486 Short_t ncharge = nTrackTPC->Charge();
1487 Float_t npt = nTrackTPC->Pt();
1488 Float_t np = nTrackTPC->P();
1489 Float_t neta = nTrackTPC->Eta();
1490 Float_t nphi = nTrackTPC->Phi();
1491 Short_t nncl = nTrackTPC->GetTPCsignalN();
1492 Short_t nneff = Short_t(nTrackTPC->GetTPCClusterInfo(2, 1)); // effective track length for pT res
1493 Float_t ndedx = nTrackTPC->GetTPCsignal();
1494 Float_t ntpcchi = 0;
1495 if(nTrackTPC->GetTPCNcls() > 0)
1496 ntpcchi = nTrackTPC->GetTPCchi2()/Float_t(nTrackTPC->GetTPCNcls());
1497
1498 Float_t bp[2]={0,0};
1499 Float_t bCovp[3]={0,0,0};
1500 pTrackTPC->GetImpactParameters(bp,bCovp);
1501 Float_t pdcaxy = bp[0];
1502 Float_t pdcaz = bp[1];
1503
1504 Float_t bn[2]={0,0};
1505 Float_t bCovn[3]={0,0,0};
1506 nTrackTPC->GetImpactParameters(bn,bCovn);
1507 Float_t ndcaxy = bn[0];
1508 Float_t ndcaz = bn[1];
1509
1510
1511 Float_t alpha = esdV0->AlphaV0();
1512 Float_t ptarm = esdV0->PtArmV0();
1513 // Double_t pVtxPos= v0->PrimaryVtxPosition();
1514
1515 Double_t lV0Position[3];
1516 esdV0->GetXYZ(lV0Position[0], lV0Position[1], lV0Position[2]);
1517
1518 Double_t lV0Radius = TMath::Sqrt(lV0Position[0]*lV0Position[0]+lV0Position[1]*lV0Position[1]);
1519 Double_t lV0DecayLength = TMath::Sqrt(TMath::Power(lV0Position[0] - lPrimaryVtxPosition[0],2) +
1520 TMath::Power(lV0Position[1] - lPrimaryVtxPosition[1],2) +
1521 TMath::Power(lV0Position[2] - lPrimaryVtxPosition[2],2 ));
1522 AliKFVertex primaryVtxKF( *myPrimaryVertex );
1523 AliKFParticle::SetField(ESDevent->GetMagneticField());
1524
1525 // Also implement switch here!!!!!!
1526 AliKFParticle* negEKF = 0; // e-
1527 AliKFParticle* posEKF = 0; // e+
1528 AliKFParticle* negPiKF = 0; // pi -
1529 AliKFParticle* posPiKF = 0; // pi +
1530 AliKFParticle* posPKF = 0; // p
1531 AliKFParticle* negAPKF = 0; // p-bar
1532
1533
1534 if(!isSwitched) {
1535 /*
1536 negEKF = new AliKFParticle( *(esdV0->GetParamN()) , 11);
1537 posEKF = new AliKFParticle( *(esdV0->GetParamP()) ,-11);
1538 negPiKF = new AliKFParticle( *(esdV0->GetParamN()) ,-211);
1539 posPiKF = new AliKFParticle( *(esdV0->GetParamP()) , 211);
1540 posPKF = new AliKFParticle( *(esdV0->GetParamP()) , 2212);
1541 negAPKF = new AliKFParticle( *(esdV0->GetParamN()) ,-2212);
1542 */
1543
1544 negEKF = new AliKFParticle( *(dynamic_cast<AliVTrack*>(nTrackTPC)) , 11);
1545 posEKF = new AliKFParticle( *(dynamic_cast<AliVTrack*>(pTrackTPC)) ,-11);
1546 negPiKF = new AliKFParticle( *(dynamic_cast<AliVTrack*>(nTrackTPC)) ,-211);
1547 posPiKF = new AliKFParticle( *(dynamic_cast<AliVTrack*>(pTrackTPC)) , 211);
1548 posPKF = new AliKFParticle( *(dynamic_cast<AliVTrack*>(pTrackTPC)) , 2212);
1549 negAPKF = new AliKFParticle( *(dynamic_cast<AliVTrack*>(nTrackTPC)) ,-2212);
1550
1551
1552
1553 } else { // switch + and -
1554 /*
1555 negEKF = new AliKFParticle( *(esdV0->GetParamP()) , 11);
1556 posEKF = new AliKFParticle( *(esdV0->GetParamN()) ,-11);
1557 negPiKF = new AliKFParticle( *(esdV0->GetParamP()) ,-211);
1558 posPiKF = new AliKFParticle( *(esdV0->GetParamN()) , 211);
1559 posPKF = new AliKFParticle( *(esdV0->GetParamN()) , 2212);
1560 negAPKF = new AliKFParticle( *(esdV0->GetParamP()) ,-2212);
1561 */
1562
1563 negEKF = new AliKFParticle( *(dynamic_cast<AliVTrack*>(pTrackTPC)), 11);
1564 posEKF = new AliKFParticle( *(dynamic_cast<AliVTrack*>(nTrackTPC)),-11);
1565 negPiKF = new AliKFParticle( *(dynamic_cast<AliVTrack*>(pTrackTPC)),-211);
1566 posPiKF = new AliKFParticle( *(dynamic_cast<AliVTrack*>(nTrackTPC)), 211);
1567 posPKF = new AliKFParticle( *(dynamic_cast<AliVTrack*>(nTrackTPC)), 2212);
1568 negAPKF = new AliKFParticle( *(dynamic_cast<AliVTrack*>(pTrackTPC)),-2212);
1569
1570
1571 }
1572
1573
1574
1575
1576
1577 AliKFParticle v0GKF; // Gamma e.g. from pi0
1578 v0GKF+=(*negEKF);
1579 v0GKF+=(*posEKF);
1580 v0GKF.SetProductionVertex(primaryVtxKF);
1581
1582 AliKFParticle v0K0sKF; // K0 short
1583 v0K0sKF+=(*negPiKF);
1584 v0K0sKF+=(*posPiKF);
1585 v0K0sKF.SetProductionVertex(primaryVtxKF);
1586
1587 AliKFParticle v0LambdaKF; // Lambda
1588 v0LambdaKF+=(*negPiKF);
1589 v0LambdaKF+=(*posPKF);
1590 v0LambdaKF.SetProductionVertex(primaryVtxKF);
1591
1592 AliKFParticle v0AntiLambdaKF; // Lambda-bar
1593 v0AntiLambdaKF+=(*posPiKF);
1594 v0AntiLambdaKF+=(*negAPKF);
1595 v0AntiLambdaKF.SetProductionVertex(primaryVtxKF);
1596
1597 Double_t deltaInvMassG = v0GKF.GetMass();
1598 Double_t deltaInvMassK0s = v0K0sKF.GetMass()-0.498;
1599 Double_t deltaInvMassL = v0LambdaKF.GetMass()-1.116;
1600 Double_t deltaInvMassAntiL = v0AntiLambdaKF.GetMass()-1.116;
1601
1602 if(TMath::Abs(deltaInvMassK0s) > fMassCut &&
1603 TMath::Abs(deltaInvMassL) > fMassCut &&
1604 TMath::Abs(deltaInvMassAntiL) > fMassCut)
1605 continue;
1606
1607 // TODO: Whe should these be different? Different mass hypothesis = energy loss
1608 // This is not important for us as we focus on the decay products!
1609 // Double_t ptK0s = v0K0sKF.GetPt();
1610 // Double_t ptL = v0LambdaKF.GetPt();
1611 // Double_t ptAntiL = v0AntiLambdaKF.GetPt();
1612
1613 Int_t primaryV0 = 0; // 0 means that the tracks are not both daughters of a primary particle (1 means they are)
1614 Int_t pdgV0 = 0; // 0 means that they don't have same origin for MC (1 means they have the same original mother)
1615 Float_t p_ptMC = 0;
1616 Short_t p_pidCode = 0; // 0 = real data / no mc track!
1617 Short_t p_primaryFlag = 0; // 0 = real data / not primary mc track
1618 Int_t p_pdgMother = 0;
1619 Float_t n_ptMC = 0;
1620 Short_t n_pidCode = 0; // 0 = real data / no mc track!
1621 Short_t n_primaryFlag = 0; // 0 = real data / not primary mc track
1622 Int_t n_pdgMother = 0;
1623 if(fAnalysisMC) {
1624
1625 Int_t p_mother_label = 0;
1626 Int_t p_mother_steps = 0;
1627 Int_t n_mother_label = 0;
1628 Int_t n_mother_steps = 0;
1629
1630 // positive track
1631 const Int_t p_label = TMath::Abs(pTrackTPC->GetLabel());
1632 TParticle* p_mcTrack = fMCStack->Particle(p_label);
1633 if (p_mcTrack){
1634
1635 if(fMCStack->IsPhysicalPrimary(p_label))
1636 p_primaryFlag = 1;
1637
1638 Int_t p_pdgCode = p_mcTrack->GetPdgCode();
1639 p_pidCode = GetPidCode(p_pdgCode);
1640
1641 p_ptMC = p_mcTrack->Pt();
1642
1643 p_mother_label = FindPrimaryMotherLabel(fMCStack, p_label,
1644 p_mother_steps);
1645 if(p_mother_label>0) {
1646 TParticle* p_mother = fMCStack->Particle(p_mother_label);
1647 p_pdgMother = p_mother->GetPdgCode();
1648 }
1649 }
1650
1651 // negative track
1652 const Int_t n_label = TMath::Abs(nTrackTPC->GetLabel());
1653 TParticle* n_mcTrack = fMCStack->Particle(n_label);
1654 if (n_mcTrack){
1655
1656 if(fMCStack->IsPhysicalPrimary(n_label))
1657 n_primaryFlag = 1;
1658
1659 Int_t n_pdgCode = n_mcTrack->GetPdgCode();
1660 n_pidCode = GetPidCode(n_pdgCode);
1661
1662 n_ptMC = n_mcTrack->Pt();
1663
1664 n_mother_label = FindPrimaryMotherLabel(fMCStack, n_label,
1665 n_mother_steps);
1666 if(n_mother_label>0) {
1667 TParticle* n_mother = fMCStack->Particle(n_mother_label);
1668 n_pdgMother = n_mother->GetPdgCode();
1669 }
1670 }
1671
1672 // Check if V0 is primary = first and the same mother of both partciles
1673 if(p_mother_label>0 && n_mother_label>0 && p_mother_label == n_mother_label) {
1674 pdgV0 = p_pdgMother;
1675 if(p_mother_steps == 1 && n_mother_steps == 1) {
1676 primaryV0 = 1;
1677 }
1678 }
1679 }
1680
1681
1682 if(fTreeOption) {
1683
1684 DeDxV0* v0datatpc = new((*fV0ArrayTPCPar)[nadded]) DeDxV0();
1685 nadded++;
1686
1687 // v0 data
1688 v0datatpc->p = esdV0->P();
1689 v0datatpc->pt = esdV0->Pt();
1690 v0datatpc->eta = esdV0->Eta();
1691 v0datatpc->phi = esdV0->Phi();
1692 v0datatpc->pdca = TMath::Sqrt(pdcaxy*pdcaxy + pdcaz*pdcaz);
1693 v0datatpc->ndca = TMath::Sqrt(ndcaxy*ndcaxy + ndcaz*ndcaz);
1694 v0datatpc->dmassG = deltaInvMassG;
1695 v0datatpc->dmassK0 = deltaInvMassK0s;
1696 v0datatpc->dmassL = deltaInvMassL;
1697 v0datatpc->dmassAL = deltaInvMassAntiL;
1698 v0datatpc->alpha = alpha;
1699 v0datatpc->ptarm = ptarm;
1700 v0datatpc->decayr = lV0Radius;
1701 v0datatpc->decayl = lV0DecayLength;
1702
1703 // New parameters
1704 v0datatpc->status = esdV0->GetOnFlyStatus();
1705 v0datatpc->chi2 = esdV0->GetChi2V0();
1706 v0datatpc->cospt = esdV0->GetV0CosineOfPointingAngle();
1707 // cospt: as I understand this means that the pointing to the vertex
1708 // is fine so I remove the dcaxy and dcaz for the V= class
1709 v0datatpc->dcadaughters = esdV0->GetDcaV0Daughters();
1710 v0datatpc->primary = primaryV0;
1711 v0datatpc->pdg = pdgV0;
1712
1713 // positive track
1714 v0datatpc->ptrack.p = pp;
1715 v0datatpc->ptrack.pt = ppt;
1716 // v0data->ptrack.ptcon = ppt_con;
1717 // v0data->ptrack.tpcchi = ptpcchi;
1718 v0datatpc->ptrack.eta = peta;
1719 v0datatpc->ptrack.phi = pphi;
1720 v0datatpc->ptrack.q = pcharge;
1721 v0datatpc->ptrack.ncl = pncl;
1722 v0datatpc->ptrack.neff = pneff;
1723 v0datatpc->ptrack.dedx = pdedx;
1724 v0datatpc->ptrack.dcaxy = pdcaxy;
1725 v0datatpc->ptrack.dcaz = pdcaz;
1726 v0datatpc->ptrack.pid = p_pidCode;
1727 v0datatpc->ptrack.primary = p_primaryFlag;
1728 v0datatpc->ptrack.pttrue = p_ptMC;
1729 v0datatpc->ptrack.mother = p_pdgMother;
1730 v0datatpc->ptrack.filter = filterFlag_p;
1731 v0datatpc->ptrack.filterset1 = 0;
1732 v0datatpc->ptrack.filterset2 = 0;
1733 v0datatpc->ptrack.filterset3 = 0;
1734
1735 // negative track
1736 v0datatpc->ntrack.p = np;
1737 v0datatpc->ntrack.pt = npt;
1738 // v0data->ntrack.ptcon = npt_con;
1739 // v0data->ntrack.tpcchi = ntpcchi;
1740 v0datatpc->ntrack.eta = neta;
1741 v0datatpc->ntrack.phi = nphi;
1742 v0datatpc->ntrack.q = ncharge;
1743 v0datatpc->ntrack.ncl = nncl;
1744 v0datatpc->ntrack.neff = nneff;
1745 v0datatpc->ntrack.dedx = ndedx;
1746 v0datatpc->ntrack.dcaxy = ndcaxy;
1747 v0datatpc->ntrack.dcaz = ndcaz;
1748 v0datatpc->ntrack.pid = n_pidCode;
1749 v0datatpc->ntrack.primary = n_primaryFlag;
1750 v0datatpc->ntrack.pttrue = n_ptMC;
1751 v0datatpc->ntrack.mother = n_pdgMother;
1752 v0datatpc->ntrack.filter = filterFlag_n;
1753 v0datatpc->ntrack.filterset1 = 0;
1754 v0datatpc->ntrack.filterset2 = 0;
1755 v0datatpc->ntrack.filterset3 = 0;
1756 }
1757
1758 // clean up loop over v0
1759
1760 delete negPiKF;
1761 delete posPiKF;
1762 delete posPKF;
1763 delete negAPKF;
1764
1765
1766
1767 }
1768
1769
1770
1771
1772
1773 }
1774 delete myPrimaryVertex;
1775
1776 if(fTreeOption) {
1777
1778 if( analysisMode==kGlobalTrk ){
1779
1780
1781 fEvent->trackmult = trackmult;
1782 fEvent->n = nadded;
1783
1784
1785 }
1786
1787
1788 }
1789
1790
1791
1792}
1793//_______________________________________________________________________________________________________________
1794void AliAnalysisTaskHighPtDeDxV0::ProduceArrayTrksAOD( AliAODEvent *AODevent, AnalysisMode analysisMode ){
1795 Int_t nv0s = AODevent->GetNumberOfV0s();
1796 if(nv0s<1)return;
1797 Int_t trackmult = 0; // no pt cuts
1798 Int_t nadded = 0;
1799
1800 AliAODVertex *myBestPrimaryVertex = AODevent->GetPrimaryVertex();
1801 if (!myBestPrimaryVertex) return;
1802
1803
1804 if( analysisMode == kGlobalTrk ){
1805 if(fV0ArrayGlobalPar)
1806 fV0ArrayGlobalPar->Clear();
1807 } else if( analysisMode == kTPCTrk ){
1808 if(fV0ArrayTPCPar)
1809 fV0ArrayTPCPar->Clear();
1810 }
1811
1812 if( analysisMode==kGlobalTrk ){
1813
1814
1815 // ################################
1816 // #### BEGINNING OF V0 CODE ######
1817 // ################################
1818 // This is the begining of the V0 loop
1819 for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {
1820 AliAODv0 *aodV0 = AODevent->GetV0(iV0);
1821 if (!aodV0) continue;
1822
1823 // common part
1824
1825 // AliAODTrack (V0 Daughters)
1826 AliAODVertex* vertex = aodV0->GetSecondaryVtx();
1827 if (!vertex) {
1828 Printf("ERROR: Could not retrieve vertex");
1829 continue;
1830 }
1831
1832 AliAODTrack *pTrack = (AliAODTrack*)vertex->GetDaughter(0);
1833 AliAODTrack *nTrack = (AliAODTrack*)vertex->GetDaughter(1);
1834 if (!pTrack || !nTrack) {
1835 Printf("ERROR: Could not retrieve one of the daughter track");
1836 continue;
1837 }
1838
1839 // Remove like-sign
1840 if (pTrack->Charge() == nTrack->Charge()) {
1841 //cout<< "like sign, continue"<< endl;
1842 continue;
1843 }
1844
1845 // Make sure charge ordering is ok
1846 if (pTrack->Charge() < 0) {
1847 AliAODTrack* helpTrack = pTrack;
1848 pTrack = nTrack;
1849 nTrack = helpTrack;
1850 }
1851
1852 // Eta cut on decay products
1853 if(TMath::Abs(pTrack->Eta()) > fEtaCut || TMath::Abs(nTrack->Eta()) > fEtaCut)
1854 continue;
1855
1856 // Pt cut on decay products
1857 if (aodV0->Pt() < fMinPt)
1858 // if (pTrack->Pt() < fMinPt && nTrack->Pt() < fMinPt)
1859 continue;
1860
1861 //check positive tracks
1862 UShort_t filterFlag_p = 0;
1863 Bool_t filterCut_Set1_p = kFALSE;//parameters from global tracks, with TPC cuts (filter bit =1 in AOD)
1864 Bool_t filterCut_Set2_p = kFALSE;//parameters from global tracks, cuts tpc+its 2010 W/O golden cuts
1865 Bool_t filterCut_Set3_p = kFALSE;//parameters from global tracks, cuts its+tpc 2010 WITH golden cuts
1866
1867 if (fTrackFilterGolden) {
1868 // ITSTPC2010 cuts is bit 32 according to above macro, new versions of aliroot includes the golden cuts
1869
1870 if(pTrack->TestFilterBit(32)) {
1871 filterFlag_p +=1;
1872 filterCut_Set3_p = kTRUE;
1873 }
1874 }
1875
1876
1877 if (fTrackFilterTPC) {
1878 // TPC only cuts is bit 1 according to above macro
1879 // Alex always uses 128, NOTE: FILTER 128 ARE TPC TRACKS (TPC PARAMETERS) CONTRAINED TO THE SPD VERTEX,
1880 if(pTrack->TestFilterBit(1)){
1881 filterFlag_p +=2;
1882 filterCut_Set1_p = kTRUE;
1883
1884 }
1885 }
1886
1887 if(filterFlag_p==0)
1888 continue;
1889
1890 //check negative tracks
1891 UShort_t filterFlag_n = 0;
1892 Bool_t filterCut_Set1_n = kFALSE;//parameters from global tracks, with TPC cuts (filter bit =1 in AOD)
1893 Bool_t filterCut_Set2_n = kFALSE;//parameters from global tracks, cuts tpc+its 2010 W/O golden cuts
1894 Bool_t filterCut_Set3_n = kFALSE;//parameters from global tracks, cuts its+tpc 2010 WITH golden cuts
1895
1896
1897 if (fTrackFilterGolden) {
1898
1899 // ITSTPC2010 cuts is bit 32 according to above macro, new versions of aliroot includes the golden cuts
1900 if(nTrack->TestFilterBit(32)) {
1901 filterFlag_n +=1;
1902 filterCut_Set3_n = kTRUE;
1903 }
1904 }
1905
1906
1907 if (fTrackFilterTPC) {
1908 // TPC only cuts is bit 1 according to above macro
1909 // Alex always uses 128, NOTE: FILTER 128 ARE TPC TRACKS (TPC PARAMETERS) CONTRAINED TO THE SPD VERTEX,
1910 if(nTrack->TestFilterBit(1)){
1911 filterFlag_n +=2;
1912 filterCut_Set1_n = kTRUE;
1913
1914 }
1915 }
1916
1917 if(filterFlag_n==0)
1918 continue;
1919
1920
1921 Float_t alpha = aodV0->AlphaV0();
1922 Float_t ptarm = aodV0->PtArmV0();
1923 // Double_t pVtxPos= v0->PrimaryVtxPosition();
1924
1925 Double_t lV0Radius = aodV0->RadiusV0();
1926 Double_t lV0DecayLength = aodV0->DecayLength(myBestPrimaryVertex);
1927
1928 Double_t deltaInvMassG = aodV0->InvMass2Prongs(0,1,11,11);
1929 Double_t deltaInvMassK0s = aodV0->MassK0Short()-0.498;
1930 Double_t deltaInvMassL = aodV0->MassLambda()-1.116;
1931 Double_t deltaInvMassAntiL = aodV0->MassAntiLambda()-1.116;
1932
1933 if(TMath::Abs(deltaInvMassK0s) > fMassCut &&
1934 TMath::Abs(deltaInvMassL) > fMassCut &&
1935 TMath::Abs(deltaInvMassAntiL) > fMassCut)
1936 continue;
1937
1938 // TODO: Why should these be different? Different mass hypothesis = energy loss
1939 // This is not important for us as we focus on the decay products!
1940 // Double_t ptK0s = v0K0sKF.GetPt();
1941 // Double_t ptL = v0LambdaKF.GetPt();
1942 // Double_t ptAntiL = v0AntiLambdaKF.GetPt();
1943
1944 // Extract track information
1945
1946 Double_t b[2], cov[3];
1947 if(!pTrack->PropagateToDCA(myBestPrimaryVertex, AODevent->GetMagneticField(), kVeryBig, b, cov))
1948 filterFlag_p += 32; // propagation failed!!!!!
1949
1950 Float_t pdcaxy = b[0];
1951 Float_t pdcaz = b[1];
1952 if(!nTrack->PropagateToDCA(myBestPrimaryVertex, AODevent->GetMagneticField(), kVeryBig, b, cov))
1953 filterFlag_n += 32; // propagation failed!!!!!
1954 Float_t ndcaxy = b[0];
1955 Float_t ndcaz = b[1];
1956
1957 Short_t pcharge = pTrack->Charge();
1958 Float_t ppt = pTrack->Pt();
1959 Float_t pp = pTrack->P();
1960 Float_t peta = pTrack->Eta();
1961 Float_t pphi = pTrack->Phi();
1962 // Float_t ptpcchi = pTrack->Chi2perNDF();
1963
1964 AliAODPid* pPid = pTrack->GetDetPid();
1965 Short_t pncl = -10;
1966 Short_t pneff = 0; // This is not yet there! Short_t(aodTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
1967 Float_t pdedx = -10;
1968 Float_t pbeta = -99;
1969 if(pPid) {
1970 pncl = pPid->GetTPCsignalN();
1971 pdedx = pPid->GetTPCsignal();
1972 //TOF
1973 if (pTrack->GetStatus()&AliESDtrack::kTOFpid){
1974 Double_t tof[5];
1975 pPid->GetIntegratedTimes(tof);
1976 pbeta = tof[0]/pPid->GetTOFsignal();
1977 }
1978 }
1979
1980 Short_t ncharge = nTrack->Charge();
1981 Float_t npt = nTrack->Pt();
1982 Float_t np = nTrack->P();
1983 Float_t neta = nTrack->Eta();
1984 Float_t nphi = nTrack->Phi();
1985 // Float_t ntpcchi = nTrack->Chi2perNDF();
1986
1987 AliAODPid* nPid = nTrack->GetDetPid();
1988 Short_t nncl = -10;
1989 Short_t nneff = 0; // This is not yet there! Short_t(aodTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
1990 Float_t ndedx = -10;
1991 Float_t nbeta = -99;
1992 if(pPid) {
1993 nncl = nPid->GetTPCsignalN();
1994 ndedx = nPid->GetTPCsignal();
1995 //TOF
1996 if (nTrack->GetStatus()&AliESDtrack::kTOFpid){
1997 Double_t tof[5];
1998 nPid->GetIntegratedTimes(tof);
1999 nbeta = tof[0]/nPid->GetTOFsignal();
2000 }
2001 }
2002
2003 Int_t primaryV0 = 0; // 0 means that the tracks are not both daughters of a primary particle (1 means they are)
2004 Int_t pdgV0 = 0; // 0 means that they don't have same origin for MC (1 means they have the same original mother)
2005 Float_t p_ptMC = 0;
2006 Short_t p_pidCode = 0; // 0 = real data / no mc track!
2007 Short_t p_primaryFlag = 0; // 0 = real data / not primary mc track
2008 Int_t p_pdgMother = 0;
2009 Float_t n_ptMC = 0;
2010 Short_t n_pidCode = 0; // 0 = real data / no mc track!
2011 Short_t n_primaryFlag = 0; // 0 = real data / not primary mc track
2012 Int_t n_pdgMother = 0;
2013 if(fAnalysisMC) {
2014
2015 AliAODMCParticle* p_mother = 0;
2016 Int_t p_mother_steps = 0;
2017 AliAODMCParticle* n_mother = 0;
2018 Int_t n_mother_steps = 0;
2019
2020 // positive track
2021 const Int_t p_label = TMath::Abs(pTrack->GetLabel());
2022
2023 AliAODMCParticle* p_mcTrack = dynamic_cast<AliAODMCParticle*>(fMCArray->At(p_label));
2024 if (p_mcTrack){
2025
2026 if(p_mcTrack->IsPhysicalPrimary())
2027 p_primaryFlag = 1;
2028
2029 Int_t p_pdgCode = p_mcTrack->GetPdgCode();
2030 p_pidCode = GetPidCode(p_pdgCode);
2031
2032 p_ptMC = p_mcTrack->Pt();
2033
2034 p_mother = FindPrimaryMotherAOD(p_mcTrack, p_mother_steps);
2035 if(p_mother)
2036 p_pdgMother = p_mother->GetPdgCode();
2037 }
2038
2039 // negative track
2040 const Int_t n_label = TMath::Abs(pTrack->GetLabel());
2041
2042 AliAODMCParticle* n_mcTrack = dynamic_cast<AliAODMCParticle*>(fMCArray->At(n_label));
2043 if (n_mcTrack){
2044
2045 if(n_mcTrack->IsPhysicalPrimary())
2046 n_primaryFlag = 1;
2047
2048 Int_t n_pdgCode = n_mcTrack->GetPdgCode();
2049 n_pidCode = GetPidCode(n_pdgCode);
2050
2051 n_ptMC = n_mcTrack->Pt();
2052
2053 n_mother = FindPrimaryMotherAOD(n_mcTrack, n_mother_steps);
2054 if(n_mother)
2055 n_pdgMother = n_mother->GetPdgCode();
2056 }
2057
2058 // Check if V0 is primary = first and the same mother of both partciles
2059 if(p_mother && n_mother && p_mother == n_mother) {
2060 pdgV0 = p_pdgMother;
2061 if(p_mother_steps == 1 && n_mother_steps == 1) {
2062 primaryV0 = 1;
2063 }
2064 }
2065 }
2066
2067 if(fTreeOption) {
2068
2069 DeDxV0* v0data = new((*fV0ArrayGlobalPar)[nadded]) DeDxV0();
2070 nadded++;
2071
2072 // v0 data
2073 v0data->p = aodV0->P();
2074 v0data->pt = aodV0->Pt();
2075 v0data->eta = aodV0->Eta();
2076 v0data->phi = aodV0->Phi();
2077 v0data->pdca = aodV0->DcaPosToPrimVertex();
2078 v0data->ndca = aodV0->DcaNegToPrimVertex();
2079 v0data->dmassG = deltaInvMassG;
2080 v0data->dmassK0 = deltaInvMassK0s;
2081 v0data->dmassL = deltaInvMassL;
2082 v0data->dmassAL = deltaInvMassAntiL;
2083 v0data->alpha = alpha;
2084 v0data->ptarm = ptarm;
2085 v0data->decayr = lV0Radius;
2086 v0data->decayl = lV0DecayLength;
2087 // v0data->pdca = TMath::Sqrt(pdcaxy*pdcaxy + pdcaz*pdcaz);
2088 // v0data->ndca = TMath::Sqrt(ndcaxy*ndcaxy + ndcaz*ndcaz);
2089
2090 // New parameters
2091 v0data->status = aodV0->GetOnFlyStatus();
2092 v0data->chi2 = aodV0->Chi2V0();
2093 v0data->cospt = aodV0->CosPointingAngle(myBestPrimaryVertex);
2094 // cospt: as I understand this means that the pointing to the vertex
2095 // is fine so I remove the dcaxy and dcaz for the V= class
2096 v0data->dcav0 = aodV0->DcaV0ToPrimVertex();
2097 v0data->dcadaughters = aodV0->DcaV0Daughters();
2098 v0data->primary = primaryV0;
2099 v0data->pdg = pdgV0;
2100
2101 // positive track
2102 v0data->ptrack.p = pp;
2103 v0data->ptrack.pt = ppt;
2104 // v0data->ptrack.ptcon = ppt_con;
2105 // v0data->ptrack.tpcchi = ptpcchi;
2106 v0data->ptrack.eta = peta;
2107 v0data->ptrack.phi = pphi;
2108 v0data->ptrack.q = pcharge;
2109 v0data->ptrack.ncl = pncl;
2110 v0data->ptrack.neff = pneff;
2111 v0data->ptrack.dedx = pdedx;
2112 v0data->ptrack.dcaxy = pdcaxy;
2113 v0data->ptrack.dcaz = pdcaz;
2114 v0data->ptrack.pid = p_pidCode;
2115 v0data->ptrack.primary = p_primaryFlag;
2116 v0data->ptrack.pttrue = p_ptMC;
2117 v0data->ptrack.mother = p_pdgMother;
2118 v0data->ptrack.filter = filterFlag_p;
2119 v0data->ptrack.filterset1 = filterCut_Set1_p;
2120 v0data->ptrack.filterset2 = filterCut_Set2_p;
2121 v0data->ptrack.filterset3 = filterCut_Set3_p;
2122
2123
2124 // negative track
2125 v0data->ntrack.p = np;
2126 v0data->ntrack.pt = npt;
2127 // v0data->ntrack.ptcon = npt_con;
2128 // v0data->ntrack.tpcchi = ntpcchi;
2129 v0data->ntrack.eta = neta;
2130 v0data->ntrack.phi = nphi;
2131 v0data->ntrack.q = ncharge;
2132 v0data->ntrack.ncl = nncl;
2133 v0data->ntrack.neff = nneff;
2134 v0data->ntrack.dedx = ndedx;
2135 v0data->ntrack.dcaxy = ndcaxy;
2136 v0data->ntrack.dcaz = ndcaz;
2137 v0data->ntrack.pid = n_pidCode;
2138 v0data->ntrack.primary = n_primaryFlag;
2139 v0data->ntrack.pttrue = n_ptMC;
2140 v0data->ntrack.mother = n_pdgMother;
2141 v0data->ntrack.filter = filterFlag_n;
2142 v0data->ntrack.filterset1 = filterCut_Set1_n;
2143 v0data->ntrack.filterset2 = filterCut_Set2_n;
2144 v0data->ntrack.filterset3 = filterCut_Set3_n;
2145
2146
2147 }
2148 }//end loop over v0's
2149
2150
2151 }else if( analysisMode==kTPCTrk ){
2152
2153 const AliAODVertex* vertexSPD= (AliAODVertex*)AODevent->GetPrimaryVertexSPD();//GetPrimaryVertex()
2154 cout<<"&&&&&&&&&&&&&&&&&&&&& Hello world 0"<<endl;
2155 if( vertexSPD->GetNContributors() < 1 || TMath::Abs(vertexSPD->GetZ()) > 10.0 ) return;
2156
2157 cout<<"&&&&&&&&&&&&&&&&&&&&& Hello world"<<endl;
2158 // ################################
2159 // #### BEGINNING OF V0 CODE ######
2160 // ################################
2161 // This is the begining of the V0 loop
2162
2163
2164
2165 for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {
2166 AliAODv0 *aodV0 = AODevent->GetV0(iV0);
2167 if (!aodV0) continue;
2168
2169 // common part
2170
2171 // AliAODTrack (V0 Daughters)
2172 AliAODVertex* vertex = aodV0->GetSecondaryVtx();
2173 if (!vertex) {
2174 Printf("ERROR: Could not retrieve vertex");
2175 continue;
2176 }
2177
2178 AliAODTrack *pTrack = (AliAODTrack*)vertex->GetDaughter(0);
2179 AliAODTrack *nTrack = (AliAODTrack*)vertex->GetDaughter(1);
2180 if (!pTrack || !nTrack) {
2181 Printf("ERROR: Could not retrieve one of the daughter track");
2182 continue;
2183 }
2184
2185
2186 // Remove like-sign
2187 if (pTrack->Charge() == nTrack->Charge()) {
2188 //cout<< "like sign, continue"<< endl;
2189 continue;
2190 }
2191
2192 // Make sure charge ordering is ok
2193 if (pTrack->Charge() < 0) {
2194 AliAODTrack* helpTrack = pTrack;
2195 pTrack = nTrack;
2196 nTrack = helpTrack;
2197 }
2198
2199 // Eta cut on decay products
2200 if(TMath::Abs(pTrack->Eta()) > fEtaCut || TMath::Abs(nTrack->Eta()) > fEtaCut)
2201 continue;
2202
2203 cout<<"Eta positive track:"<<pTrack->Eta()<<endl;
2204
2205
2206 // Pt cut on decay products
2207 if (aodV0->Pt() < fMinPt)
2208 // if (pTrack->Pt() < fMinPt && nTrack->Pt() < fMinPt)
2209 continue;
2210
2211 //check positive tracks
2212 UShort_t filterFlag_p = 0;
2213 Bool_t filterCut_Set1_p = kFALSE;//parameters from global tracks, with TPC cuts (filter bit =1 in AOD)
2214 Bool_t filterCut_Set2_p = kFALSE;//parameters from global tracks, cuts tpc+its 2010 W/O golden cuts
2215 Bool_t filterCut_Set3_p = kFALSE;//parameters from global tracks, cuts its+tpc 2010 WITH golden cuts
2216
2217 // TPC only cuts is bit 1 according to above macro
2218 // Alex always uses 128, NOTE: FILTER 128 ARE TPC TRACKS (TPC PARAMETERS) CONTRAINED TO THE SPD VERTEX,
2219 if(pTrack->TestFilterBit(128)) {
2220 cout<<"este track paso el corte bit 128"<<endl;
2221 filterFlag_p +=1;
2222 }
2223 cout<<"filterFlag_p="<<filterFlag_p<<endl;
2224
2225
2226
2227
2228 if(filterFlag_p==0)
2229 continue;
2230
2231
2232
2233 //check negative tracks
2234 UShort_t filterFlag_n = 0;
2235 Bool_t filterCut_Set1_n = kFALSE;//parameters from global tracks, with TPC cuts (filter bit =1 in AOD)
2236 Bool_t filterCut_Set2_n = kFALSE;//parameters from global tracks, cuts tpc+its 2010 W/O golden cuts
2237 Bool_t filterCut_Set3_n = kFALSE;//parameters from global tracks, cuts its+tpc 2010 WITH golden cuts
2238
2239 // TPC only cuts is bit 1 according to above macro
2240 // Alex always uses 128, NOTE: FILTER 128 ARE TPC TRACKS (TPC PARAMETERS) CONTRAINED TO THE SPD VERTEX,
2241 if(nTrack->TestFilterBit(128)) {
2242 filterFlag_n +=1;
2243 }
2244
2245 if(filterFlag_n==0)
2246 continue;
2247
2248
2249 Float_t alpha = aodV0->AlphaV0();
2250 Float_t ptarm = aodV0->PtArmV0();
2251 // Double_t pVtxPos= v0->PrimaryVtxPosition();
2252
2253 Double_t lV0Radius = aodV0->RadiusV0();
2254 Double_t lV0DecayLength = aodV0->DecayLength(myBestPrimaryVertex);
2255
2256 Double_t deltaInvMassG = aodV0->InvMass2Prongs(0,1,11,11);
2257 Double_t deltaInvMassK0s = aodV0->MassK0Short()-0.498;
2258 Double_t deltaInvMassL = aodV0->MassLambda()-1.116;
2259 Double_t deltaInvMassAntiL = aodV0->MassAntiLambda()-1.116;
2260
2261 if(TMath::Abs(deltaInvMassK0s) > fMassCut &&
2262 TMath::Abs(deltaInvMassL) > fMassCut &&
2263 TMath::Abs(deltaInvMassAntiL) > fMassCut)
2264 continue;
2265
2266 // TODO: Why should these be different? Different mass hypothesis = energy loss
2267 // This is not important for us as we focus on the decay products!
2268 // Double_t ptK0s = v0K0sKF.GetPt();
2269 // Double_t ptL = v0LambdaKF.GetPt();
2270 // Double_t ptAntiL = v0AntiLambdaKF.GetPt();
2271
2272 // Extract track information
2273
2274 Double_t b[2], cov[3];
2275 if(!pTrack->PropagateToDCA(vertexSPD, AODevent->GetMagneticField(), kVeryBig, b, cov))
2276 filterFlag_p += 32; // propagation failed!!!!!
2277
2278 Float_t pdcaxy = b[0];
2279 Float_t pdcaz = b[1];
2280 if(!nTrack->PropagateToDCA(vertexSPD, AODevent->GetMagneticField(), kVeryBig, b, cov))
2281 filterFlag_n += 32; // propagation failed!!!!!
2282 Float_t ndcaxy = b[0];
2283 Float_t ndcaz = b[1];
2284
2285 Short_t pcharge = pTrack->Charge();
2286 Float_t ppt = pTrack->Pt();
2287 Float_t pp = pTrack->P();
2288 Float_t peta = pTrack->Eta();
2289 Float_t pphi = pTrack->Phi();
2290 // Float_t ptpcchi = pTrack->Chi2perNDF();
2291
2292 AliAODPid* pPid = pTrack->GetDetPid();
2293 Short_t pncl = -10;
2294 Short_t pneff = 0; // This is not yet there! Short_t(aodTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
2295 Float_t pdedx = -10;
2296 Float_t pbeta = -99;
2297 if(pPid) {
2298 pncl = pPid->GetTPCsignalN();
2299 pdedx = pPid->GetTPCsignal();
2300 //TOF
2301 if (pTrack->GetStatus()&AliESDtrack::kTOFpid){
2302 Double_t tof[5];
2303 pPid->GetIntegratedTimes(tof);
2304 pbeta = tof[0]/pPid->GetTOFsignal();
2305 }
2306 }
2307
2308 Short_t ncharge = nTrack->Charge();
2309 Float_t npt = nTrack->Pt();
2310 Float_t np = nTrack->P();
2311 Float_t neta = nTrack->Eta();
2312 Float_t nphi = nTrack->Phi();
2313 // Float_t ntpcchi = nTrack->Chi2perNDF();
2314
2315 AliAODPid* nPid = nTrack->GetDetPid();
2316 Short_t nncl = -10;
2317 Short_t nneff = 0; // This is not yet there! Short_t(aodTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
2318 Float_t ndedx = -10;
2319 Float_t nbeta = -99;
2320 if(pPid) {
2321 nncl = nPid->GetTPCsignalN();
2322 ndedx = nPid->GetTPCsignal();
2323 //TOF
2324 if (nTrack->GetStatus()&AliESDtrack::kTOFpid){
2325 Double_t tof[5];
2326 nPid->GetIntegratedTimes(tof);
2327 nbeta = tof[0]/nPid->GetTOFsignal();
2328 }
2329 }
2330
2331 Int_t primaryV0 = 0; // 0 means that the tracks are not both daughters of a primary particle (1 means they are)
2332 Int_t pdgV0 = 0; // 0 means that they don't have same origin for MC (1 means they have the same original mother)
2333 Float_t p_ptMC = 0;
2334 Short_t p_pidCode = 0; // 0 = real data / no mc track!
2335 Short_t p_primaryFlag = 0; // 0 = real data / not primary mc track
2336 Int_t p_pdgMother = 0;
2337 Float_t n_ptMC = 0;
2338 Short_t n_pidCode = 0; // 0 = real data / no mc track!
2339 Short_t n_primaryFlag = 0; // 0 = real data / not primary mc track
2340 Int_t n_pdgMother = 0;
2341 if(fAnalysisMC) {
2342
2343 AliAODMCParticle* p_mother = 0;
2344 Int_t p_mother_steps = 0;
2345 AliAODMCParticle* n_mother = 0;
2346 Int_t n_mother_steps = 0;
2347
2348 // positive track
2349 const Int_t p_label = TMath::Abs(pTrack->GetLabel());
2350
2351 AliAODMCParticle* p_mcTrack = dynamic_cast<AliAODMCParticle*>(fMCArray->At(p_label));
2352 if (p_mcTrack){
2353
2354 if(p_mcTrack->IsPhysicalPrimary())
2355 p_primaryFlag = 1;
2356
2357 Int_t p_pdgCode = p_mcTrack->GetPdgCode();
2358 p_pidCode = GetPidCode(p_pdgCode);
2359
2360 p_ptMC = p_mcTrack->Pt();
2361
2362 p_mother = FindPrimaryMotherAOD(p_mcTrack, p_mother_steps);
2363 if(p_mother)
2364 p_pdgMother = p_mother->GetPdgCode();
2365 }
2366
2367 // negative track
2368 const Int_t n_label = TMath::Abs(pTrack->GetLabel());
2369
2370 AliAODMCParticle* n_mcTrack = dynamic_cast<AliAODMCParticle*>(fMCArray->At(n_label));
2371 if (n_mcTrack){
2372
2373 if(n_mcTrack->IsPhysicalPrimary())
2374 n_primaryFlag = 1;
2375
2376 Int_t n_pdgCode = n_mcTrack->GetPdgCode();
2377 n_pidCode = GetPidCode(n_pdgCode);
2378
2379 n_ptMC = n_mcTrack->Pt();
2380
2381 n_mother = FindPrimaryMotherAOD(n_mcTrack, n_mother_steps);
2382 if(n_mother)
2383 n_pdgMother = n_mother->GetPdgCode();
2384 }
2385
2386 // Check if V0 is primary = first and the same mother of both partciles
2387 if(p_mother && n_mother && p_mother == n_mother) {
2388 pdgV0 = p_pdgMother;
2389 if(p_mother_steps == 1 && n_mother_steps == 1) {
2390 primaryV0 = 1;
2391 }
2392 }
2393 }
2394
2395 if(fTreeOption) {
2396
2397 DeDxV0* v0datatpc = new((*fV0ArrayTPCPar)[nadded]) DeDxV0();
2398 nadded++;
2399
2400 // v0 data
2401 v0datatpc->p = aodV0->P();
2402 v0datatpc->pt = aodV0->Pt();
2403 v0datatpc->eta = aodV0->Eta();
2404 v0datatpc->phi = aodV0->Phi();
2405 v0datatpc->pdca = aodV0->DcaPosToPrimVertex();
2406 v0datatpc->ndca = aodV0->DcaNegToPrimVertex();
2407 v0datatpc->dmassG = deltaInvMassG;
2408 v0datatpc->dmassK0 = deltaInvMassK0s;
2409 v0datatpc->dmassL = deltaInvMassL;
2410 v0datatpc->dmassAL = deltaInvMassAntiL;
2411 v0datatpc->alpha = alpha;
2412 v0datatpc->ptarm = ptarm;
2413 v0datatpc->decayr = lV0Radius;
2414 v0datatpc->decayl = lV0DecayLength;
2415 // v0data->pdca = TMath::Sqrt(pdcaxy*pdcaxy + pdcaz*pdcaz);
2416 // v0data->ndca = TMath::Sqrt(ndcaxy*ndcaxy + ndcaz*ndcaz);
2417
2418 // New parameters
2419 v0datatpc->status = aodV0->GetOnFlyStatus();
2420 v0datatpc->chi2 = aodV0->Chi2V0();
2421 v0datatpc->cospt = aodV0->CosPointingAngle(myBestPrimaryVertex);
2422 // cospt: as I understand this means that the pointing to the vertex
2423 // is fine so I remove the dcaxy and dcaz for the V= class
2424 v0datatpc->dcav0 = aodV0->DcaV0ToPrimVertex();
2425 v0datatpc->dcadaughters = aodV0->DcaV0Daughters();
2426 v0datatpc->primary = primaryV0;
2427 v0datatpc->pdg = pdgV0;
2428
2429 // positive track
2430 v0datatpc->ptrack.p = pp;
2431 v0datatpc->ptrack.pt = ppt;
2432 // v0data->ptrack.ptcon = ppt_con;
2433 // v0data->ptrack.tpcchi = ptpcchi;
2434 v0datatpc->ptrack.eta = peta;
2435 v0datatpc->ptrack.phi = pphi;
2436 v0datatpc->ptrack.q = pcharge;
2437 v0datatpc->ptrack.ncl = pncl;
2438 v0datatpc->ptrack.neff = pneff;
2439 v0datatpc->ptrack.dedx = pdedx;
2440 v0datatpc->ptrack.dcaxy = pdcaxy;
2441 v0datatpc->ptrack.dcaz = pdcaz;
2442 v0datatpc->ptrack.pid = p_pidCode;
2443 v0datatpc->ptrack.primary = p_primaryFlag;
2444 v0datatpc->ptrack.pttrue = p_ptMC;
2445 v0datatpc->ptrack.mother = p_pdgMother;
2446 v0datatpc->ptrack.filter = filterFlag_p;
2447 v0datatpc->ptrack.filterset1 = 0;
2448 v0datatpc->ptrack.filterset2 = 0;
2449 v0datatpc->ptrack.filterset3 = 0;
2450
2451
2452 // negative track
2453 v0datatpc->ntrack.p = np;
2454 v0datatpc->ntrack.pt = npt;
2455 // v0data->ntrack.ptcon = npt_con;
2456 // v0data->ntrack.tpcchi = ntpcchi;
2457 v0datatpc->ntrack.eta = neta;
2458 v0datatpc->ntrack.phi = nphi;
2459 v0datatpc->ntrack.q = ncharge;
2460 v0datatpc->ntrack.ncl = nncl;
2461 v0datatpc->ntrack.neff = nneff;
2462 v0datatpc->ntrack.dedx = ndedx;
2463 v0datatpc->ntrack.dcaxy = ndcaxy;
2464 v0datatpc->ntrack.dcaz = ndcaz;
2465 v0datatpc->ntrack.pid = n_pidCode;
2466 v0datatpc->ntrack.primary = n_primaryFlag;
2467 v0datatpc->ntrack.pttrue = n_ptMC;
2468 v0datatpc->ntrack.mother = n_pdgMother;
2469 v0datatpc->ntrack.filter = filterFlag_n;
2470 v0datatpc->ntrack.filterset1 = 0;
2471 v0datatpc->ntrack.filterset2 = 0;
2472 v0datatpc->ntrack.filterset3 = 0;
2473
2474
2475 }
2476 }//end v0's loop
2477
2478 }
2479
2480 if(fTreeOption) {
2481
2482 if( analysisMode==kGlobalTrk ){
2483
2484
2485 fEvent->trackmult = trackmult;
2486 fEvent->n = nadded;
2487
2488
2489 }
2490
2491
2492 }
2493
2494
2495
2496}
2497
2498
2499