]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/SPECTRA/IdentifiedHighPt/grid/AliAnalysisTaskHighPtDeDx.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / IdentifiedHighPt / grid / AliAnalysisTaskHighPtDeDx.cxx
CommitLineData
4ebdd20e 1#include "AliAnalysisTaskHighPtDeDx.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#include <AliESDtrackCuts.h>
20#include <AliESDVZERO.h>
21#include <AliAODVZERO.h>
22
23#include <AliMCEventHandler.h>
24#include <AliMCEvent.h>
25#include <AliStack.h>
26
27#include <TTreeStream.h>
28
29#include <AliHeader.h>
30#include <AliGenPythiaEventHeader.h>
31#include <AliGenDPMjetEventHeader.h>
32
33#include "AliCentrality.h"
34
35#include <AliAODTrack.h>
36#include <AliAODPid.h>
37#include <AliAODMCHeader.h>
38
39
40// STL includes
41#include <iostream>
42using namespace std;
43
44
45//
46// Responsible:
47// Alexandru Dobrin (Wayne State)
48// Peter Christiansen (Lund)
49//
50
51/*
52To do:
53
54Separate the code into two
55
56*/
57
58
59
60
61ClassImp(AliAnalysisTaskHighPtDeDx)
62
63const Double_t AliAnalysisTaskHighPtDeDx::fgkClight = 2.99792458e-2;
64
65//_____________________________________________________________________________
66AliAnalysisTaskHighPtDeDx::AliAnalysisTaskHighPtDeDx():
67 AliAnalysisTaskSE(),
68 fESD(0x0),
69 fAOD(0x0),
70 fMC(0x0),
71 fMCStack(0x0),
72 fMCArray(0x0),
73 fTrackFilter(0x0),
74 fTrackFilterGolden(0x0),
75 fTrackFilterTPC(0x0),
76 fAnalysisType("ESD"),
77 fAnalysisMC(kFALSE),
78 fAnalysisPbPb(kFALSE),
79 fTPCBranch(kFALSE),
80 ftrigBit1(0x0),
81 ftrigBit2(0x0),
82 fRandom(0x0),
83 fEvent(0x0),
84 fTrackArrayGlobalPar(0x0),
85 fTrackArrayTPCPar(0x0),
86 fTrackArrayMC(0x0),
87 fVZEROArray(0x0),
88 fVtxCut(10.0),
89 fEtaCut(0.9),
90 fMinPt(0.1),
91 fMinCent(0.0),
92 fMaxCent(100.0),
93 fLowPtFraction(0.01),
94 fTreeOption(0),
95 fMcProcessType(-999),
96 fTriggeredEventMB(-999),
97 fVtxStatus(-999),
98 fZvtx(-999),
99 fZvtxMC(-999),
100 fRun(-999),
101 fEventId(-999),
102 fListOfObjects(0),
103 fEvents(0x0), fVtx(0x0), fVtxMC(0x0), fVtxBeforeCuts(0x0), fVtxAfterCuts(0x0),
104 fTree(0x0),
105 fn1(0),
106 fn2(0),
107 fcent(0)
108
109
110{
111 // Default constructor (should not be used)
112
113 // fRandom = new TRandom(0); // 0 means random seed
114}
115
116//______________________________________________________________________________
117AliAnalysisTaskHighPtDeDx::AliAnalysisTaskHighPtDeDx(const char *name):
118 AliAnalysisTaskSE(name),
119 fESD(0x0),
120 fAOD(0x0),
121 fMC(0x0),
122 fMCStack(0x0),
123 fMCArray(0x0),
124 fTrackFilter(0x0),
125 fTrackFilterGolden(0x0),
126 fTrackFilterTPC(0x0),
127 fAnalysisType("ESD"),
128 fAnalysisMC(kFALSE),
129 fAnalysisPbPb(kFALSE),
130 fTPCBranch(kFALSE),
131 ftrigBit1(0x0),
132 ftrigBit2(0x0),
133 fRandom(0x0),
134 fEvent(0x0),
135 fTrackArrayGlobalPar(0x0),
136 fTrackArrayMC(0x0),
137 fVtxCut(10.0),
138 fEtaCut(0.9),
139 fMinPt(0.1),
140 fLowPtFraction(0.01),
141 fTreeOption(0),
142 fMcProcessType(-999),
143 fTriggeredEventMB(-999),
144 fVtxStatus(-999),
145 fZvtx(-999),
146 fZvtxMC(-999),
147 fRun(-999),
148 fEventId(-999),
149 fListOfObjects(0),
150 fEvents(0x0), fVtx(0x0), fVtxMC(0x0), fVtxBeforeCuts(0x0), fVtxAfterCuts(0x0),
151 fTree(0x0),
152 fn1(0),
153 fn2(0),
154 fcent(0)
155
156{
157
158 if(fTree)fTree=0;
159 // Output slot #1 writes into a TList
160 DefineOutput(1, TList::Class());
161}
162
163//_____________________________________________________________________________
164AliAnalysisTaskHighPtDeDx::~AliAnalysisTaskHighPtDeDx()
165{
166 // Destructor
167 // histograms are in the output list and deleted when the output
168 // list is deleted by the TSelector dtor
169 if (fListOfObjects && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
170 delete fListOfObjects;
171 fListOfObjects = 0;
172 }
173 if (fRandom) delete fRandom;
174 fRandom=0;
175
176 // //for proof running; I cannot create tree do to memory limitations -> work with THnSparse
177 // if (fListOfObjects && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fOutputList;
178
179
180
181}
182
183//______________________________________________________________________________
184void AliAnalysisTaskHighPtDeDx::UserCreateOutputObjects()
185{
186 // This method is called once per worker node
187 // Here we define the output: histograms and debug tree if requested
188 // We also create the random generator here so it might get different seeds...
189 fRandom = new TRandom(0); // 0 means random seed
190
191 OpenFile(1);
192 fListOfObjects = new TList();
193 fListOfObjects->SetOwner();
194
195 //
196 // Histograms
197 //
198 fEvents = new TH1I("fEvents","Number of analyzed events; Events; Counts", 3, 0, 3);
199 fListOfObjects->Add(fEvents);
200
201 fn1=new TH1F("fn1","fn1",5001,-1,5000);
202 fListOfObjects->Add(fn1);
203
204 fn2=new TH1F("fn2","fn2",5001,-1,5000);
205 fListOfObjects->Add(fn2);
206
207 fcent=new TH1F("fcent","fcent",104,-2,102);
208 fListOfObjects->Add(fcent);
209
210 fVtx = new TH1I("fVtx","Vtx info (0=no, 1=yes); Vtx; Counts", 2, -0.5, 1.5);
211 fListOfObjects->Add(fVtx);
212
213 fVtxBeforeCuts = new TH1F("fVtxBeforeCuts", "Vtx distribution (before cuts); Vtx z [cm]; Counts", 120, -30, 30);
214 fListOfObjects->Add(fVtxBeforeCuts);
215
216 fVtxAfterCuts = new TH1F("fVtxAfterCuts", "Vtx distribution (before cuts); Vtx z [cm]; Counts", 120, -30, 30);
217 fListOfObjects->Add(fVtxAfterCuts);
218
219 if (fAnalysisMC) {
220 fVtxMC = new TH1I("fVtxMC","Vtx info - no trigger cut (0=no, 1=yes); Vtx; Counts", 2, -0.5, 1.5);
221 fListOfObjects->Add(fVtxMC);
222 }
223
224 if (fTreeOption) {
225
226 fTree = new TTree("tree","Event data");
227 fEvent = new DeDxEvent();
228 fTree->Branch("event", &fEvent);
229
230 fTrackArrayGlobalPar = new TClonesArray("DeDxTrack", 1000);
231 fTree->Bronch("trackGlobalPar", "TClonesArray", &fTrackArrayGlobalPar);
232 if(fTPCBranch){
233 fTrackArrayTPCPar = new TClonesArray("DeDxTrack", 1000);
234 fTree->Bronch("trackTPCPar", "TClonesArray", &fTrackArrayTPCPar);
235 }
236
237 fVZEROArray = new TClonesArray("VZEROCell", 1000);
238 fTree->Bronch("cellVZERO", "TClonesArray", &fVZEROArray);
239
240 if (fAnalysisMC) {
241 fTrackArrayMC = new TClonesArray("DeDxTrackMC", 1000);
242 fTree->Bronch("trackMC", "TClonesArray", &fTrackArrayMC);
243 }
244
245 fTree->SetDirectory(0);
246
247 fListOfObjects->Add(fTree);
248
249 }
250
251 // Post output data.
252 PostData(1, fListOfObjects);
253}
254
255//______________________________________________________________________________
256void AliAnalysisTaskHighPtDeDx::UserExec(Option_t *)
257{
258 // Main loop
259
260 //
261 // First we make sure that we have valid input(s)!
262 //
263 AliVEvent *event = InputEvent();
264 if (!event) {
265 Error("UserExec", "Could not retrieve event");
266 return;
267 }
268
269
270 if (fAnalysisType == "ESD"){
271 fESD = dynamic_cast<AliESDEvent*>(event);
272 if(!fESD){
273 Printf("%s:%d ESDEvent not found in Input Manager",(char*)__FILE__,__LINE__);
274 this->Dump();
275 return;
276 }
277 } else {
278 fAOD = dynamic_cast<AliAODEvent*>(event);
279 if(!fAOD){
280 Printf("%s:%d AODEvent not found in Input Manager",(char*)__FILE__,__LINE__);
281 this->Dump();
282 return;
283 }
284 }
285
286
287
288 if (fAnalysisMC) {
289
290 if (fAnalysisType == "ESD"){
291 fMC = dynamic_cast<AliMCEvent*>(MCEvent());
292 if(!fMC){
293 Printf("%s:%d MCEvent not found in Input Manager",(char*)__FILE__,__LINE__);
294 this->Dump();
295 return;
296 }
297
298 fMCStack = fMC->Stack();
299
300 if(!fMCStack){
301 Printf("%s:%d MCStack not found in Input Manager",(char*)__FILE__,__LINE__);
302 this->Dump();
303 return;
304 }
305 } else { // AOD
306
307 fMC = dynamic_cast<AliMCEvent*>(MCEvent());
308 if(fMC)
309 fMC->Dump();
310
311 fMCArray = (TClonesArray*)fAOD->FindListObject("mcparticles");
312 if(!fMCArray){
313 Printf("%s:%d AOD MC array not found in Input Manager",(char*)__FILE__,__LINE__);
314 this->Dump();
315 return;
316 }
317 }
318 }
319
320
321 // Get trigger decision
322 fTriggeredEventMB = 0; //init
323
324
325 if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
326 ->IsEventSelected() & ftrigBit1 ){
327 fn1->Fill(1);
328 fTriggeredEventMB = 1; //event triggered as minimum bias
329 }
330 if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
331 ->IsEventSelected() & ftrigBit2 ){
332 // From AliVEvent:
333 // kINT7 = BIT(1), // V0AND trigger, offline V0 selection
334 fTriggeredEventMB += 2;
335 fn2->Fill(1);
336 }
337
338
339 // Get process type for MC
340 fMcProcessType = 0; // -1=invalid, 0=data, 1=ND, 2=SD, 3=DD
341
342 // real data that are not triggered we skip
343 if(!fAnalysisMC && !fTriggeredEventMB)
344 return;
345
346 if (fAnalysisMC) {
347
348 if (fAnalysisType == "ESD"){
349
350 AliHeader* headerMC = fMC->Header();
351 if (headerMC) {
352
353 AliGenEventHeader* genHeader = headerMC->GenEventHeader();
354 TArrayF vtxMC(3); // primary vertex MC
355 vtxMC[0]=9999; vtxMC[1]=9999; vtxMC[2]=9999; //initialize with dummy
356 if (genHeader) {
357 genHeader->PrimaryVertex(vtxMC);
358 }
359 fZvtxMC = vtxMC[2];
360
361 // PYTHIA:
362 AliGenPythiaEventHeader* pythiaGenHeader =
363 dynamic_cast<AliGenPythiaEventHeader*>(headerMC->GenEventHeader());
364 if (pythiaGenHeader) { //works only for pythia
365 fMcProcessType = GetPythiaEventProcessType(pythiaGenHeader->ProcessType());
366 }
367 // PHOJET:
368 AliGenDPMjetEventHeader* dpmJetGenHeader =
369 dynamic_cast<AliGenDPMjetEventHeader*>(headerMC->GenEventHeader());
370 if (dpmJetGenHeader) {
371 fMcProcessType = GetDPMjetEventProcessType(dpmJetGenHeader->ProcessType());
372 }
373 }
374 } else { // AOD
375
376 AliAODMCHeader* mcHeader = dynamic_cast<AliAODMCHeader*>(fAOD->FindListObject("mcHeader"));
377 if(mcHeader) {
378 fZvtxMC = mcHeader->GetVtxZ();
379
380 if(strstr(mcHeader->GetGeneratorName(), "Pythia")) {
381 fMcProcessType = GetPythiaEventProcessType(mcHeader->GetEventType());
382 } else {
383 fMcProcessType = GetDPMjetEventProcessType(mcHeader->GetEventType());
384 }
385 }
386 }
387 }
388
389 // There are 3 cases
390 // Vertex: NO - status -1
391 // Vertex: YES : outside cut - status 0
392 // : inside cut - status 1
393 // We have to be careful how we normalize because we probably want to
394 // normalize to:
395 // Nevents=(No vertex + outside + inside)/(out + in)*in
396
397
398 if (fAnalysisType == "ESD")
399 fZvtx = GetVertex(fESD);
400 else // AOD
401 fZvtx = GetVertex(fAOD);
402
403 fVtxStatus = -999;
404
405 if(fZvtx<-990) {
406
407 fVtxStatus = -1;
408 if(fTriggeredEventMB)
409 fVtx->Fill(0);
410 if(fAnalysisMC)
411 fVtxMC->Fill(0);
412 } else {
413
414 if(fTriggeredEventMB)
415 fVtx->Fill(1);
416 if(fAnalysisMC)
417 fVtxMC->Fill(1);
418 fVtxBeforeCuts->Fill(fZvtx);
419 fVtxStatus = 0;
420 if (TMath::Abs(fZvtx) < fVtxCut) {
421 fVtxAfterCuts->Fill(fZvtx);
422 fVtxStatus = 1;
423 }
424 }
425
426 // store MC event data no matter what
427 if(fAnalysisMC) {
428
429 if (fAnalysisType == "ESD") {
430 ProcessMCTruthESD();
431 } else { // AOD
432 ProcessMCTruthAOD();
433 }
434 }
435
436
437
438 Float_t centrality = -10;
439 // only analyze triggered events
440 if(fTriggeredEventMB) {
441
442 if (fAnalysisType == "ESD"){
443 if(fAnalysisPbPb){
444 AliCentrality *centObject = fESD->GetCentrality();
445 centrality = centObject->GetCentralityPercentile("V0M");
446 if((centrality>fMaxCent)||(centrality<fMinCent))return;
447 }
448 fcent->Fill(centrality);
449 AnalyzeESD(fESD);
450 } else { // AOD
451 if(fAnalysisPbPb){
452 AliCentrality *centObject = fAOD->GetCentrality();
453 if(centObject)
454 centrality = centObject->GetCentralityPercentile("V0M");
455 if((centrality>fMaxCent)||(centrality<fMinCent))return;
456 }
457 fcent->Fill(centrality);
458 AnalyzeAOD(fAOD);
459 }
460 }
461
462 if( fTreeOption) {
463
464 fEvent->process = fMcProcessType;
465 fEvent->trig = fTriggeredEventMB;
466 fEvent->zvtxMC = fZvtxMC;
467 fEvent->cent = centrality;
468
469 fTree->Fill();
470 fTrackArrayGlobalPar->Clear();
471 if(fTPCBranch)
472 fTrackArrayTPCPar->Clear();
473 fVZEROArray->Clear();
474
475 if (fAnalysisMC)
476 fTrackArrayMC->Clear();
477 }
478
479 // Post output data.
480 PostData(1, fListOfObjects);
481}
482
483//________________________________________________________________________
484void AliAnalysisTaskHighPtDeDx::AnalyzeESD(AliESDEvent* esdEvent)
485{
486 fRun = esdEvent->GetRunNumber();
487 fEventId = 0;
488 if(esdEvent->GetHeader())
489 fEventId = GetEventIdAsLong(esdEvent->GetHeader());
490
491 Short_t isPileup = esdEvent->IsPileupFromSPD();
492
493 // Int_t event = esdEvent->GetEventNumberInFile();
494 UInt_t time = esdEvent->GetTimeStamp();
495 // ULong64_t trigger = esdEvent->GetTriggerMask();
496 Float_t magf = esdEvent->GetMagneticField();
497
498
499
500
501
502 if(fTriggeredEventMB) {// Only MC case can we have not triggered events
503
504 // accepted event
505 fEvents->Fill(0);
506
507
508 if(fVtxStatus!=1) return; // accepted vertex
509 Int_t nESDTracks = esdEvent->GetNumberOfTracks();
510
511 if(nESDTracks<1)return;
512 cout<<"Multiplicity="<<nESDTracks<<endl;
513
514 ProduceArrayTrksESD( esdEvent, kGlobalTrk );//produce array with global track parameters
515 if(fTPCBranch)
516 ProduceArrayTrksESD( esdEvent, kTPCTrk );//array with tpc track parametes
517
518 fEvents->Fill(1);
519 AliESDVZERO *esdV0 = esdEvent->GetVZEROData();// loop sobre canales del V0 para obtener las multiplicidad
520
521 for (Int_t iCh=0; iCh<64; ++iCh) {
522 Float_t multv=esdV0->GetMultiplicity(iCh);
523 Int_t intexv=iCh;
524 VZEROCell* cellv0 = new((*fVZEROArray)[iCh]) VZEROCell();
525 cellv0->cellindex=intexv;
526 cellv0->cellmult= multv;
527 }
528
529
530
531 } // end if triggered
532
533 if(fTreeOption) {
534
535 fEvent->run = fRun;
536 fEvent->eventid = fEventId;
537 fEvent->time = time;
538 //fEvent->cent = centrality;
539 fEvent->mag = magf;
540 fEvent->zvtx = fZvtx;
541 fEvent->vtxstatus = fVtxStatus;
542 fEvent->pileup = isPileup;
543
544 }
545
546
547
548
549}
550
551//________________________________________________________________________
552void AliAnalysisTaskHighPtDeDx::AnalyzeAOD(AliAODEvent* aodEvent)
553{
554 fRun = aodEvent->GetRunNumber();
555 fEventId = 0;
556 if(aodEvent->GetHeader())
557 fEventId = GetEventIdAsLong(aodEvent->GetHeader());
558
559 UInt_t time = 0; // Missing AOD info? aodEvent->GetTimeStamp();
560 Float_t magf = aodEvent->GetMagneticField();
561
562 //Int_t trackmult = 0; // no pt cuts
563 //Int_t nadded = 0;
564
565 Short_t isPileup = aodEvent->IsPileupFromSPD();
566
567
568
569
570 if(fTriggeredEventMB) {// Only MC case can we have not triggered events
571
572 // accepted event
573 fEvents->Fill(0);
574
575 if(fVtxStatus!=1) return; // accepted vertex
576 Int_t nAODTracks = aodEvent->GetNumberOfTracks();
577 if(nAODTracks<1)return;
578
579 ProduceArrayTrksAOD( aodEvent, kGlobalTrk );
580 if(fTPCBranch)
581 ProduceArrayTrksAOD( aodEvent, kTPCTrk );
582 fEvents->Fill(1);
583
584 AliAODVZERO *esdV0 = aodEvent->GetVZEROData();// loop sobre canales del V0 para obtener las multiplicidad
585
586 for (Int_t iCh=0; iCh<64; ++iCh) {
587 Float_t multv=esdV0->GetMultiplicity(iCh);
588 Int_t intexv=iCh;
589 VZEROCell* cellv0 = new((*fVZEROArray)[iCh]) VZEROCell();
590 cellv0->cellindex=intexv;
591 cellv0->cellmult= multv;
592 }
593
594
595
596
597 } // end if triggered
598
599 if(fTreeOption) {
600
601 //Sort(fTrackArrayGlobalPar, kFALSE);
602
603 fEvent->run = fRun;
604 fEvent->eventid = fEventId;
605 fEvent->time = time;
606 //fEvent->cent = centrality;
607 fEvent->mag = magf;
608 fEvent->zvtx = fZvtx;
609 fEvent->vtxstatus = fVtxStatus;
610 //fEvent->trackmult = trackmult;
611 //fEvent->n = nadded;
612 fEvent->pileup = isPileup;
613 }
614}
615
616//_____________________________________________________________________________
617Float_t AliAnalysisTaskHighPtDeDx::GetVertex(const AliVEvent* event) const
618{
619 Float_t zvtx = -999;
620
621 const AliVVertex* primaryVertex = event->GetPrimaryVertex();
622
623 if(primaryVertex->GetNContributors()>0)
624 zvtx = primaryVertex->GetZ();
625
626 return zvtx;
627}
628
629//_____________________________________________________________________________
630Short_t AliAnalysisTaskHighPtDeDx::GetPidCode(Int_t pdgCode) const
631{
632 // return our internal code for pions, kaons, and protons
633
634 Short_t pidCode = 6;
635
636 switch (TMath::Abs(pdgCode)) {
637 case 211:
638 pidCode = 1; // pion
639 break;
640 case 321:
641 pidCode = 2; // kaon
642 break;
643 case 2212:
644 pidCode = 3; // proton
645 break;
646 case 11:
647 pidCode = 4; // electron
648 break;
649 case 13:
650 pidCode = 5; // muon
651 break;
652 default:
653 pidCode = 6; // something else?
654 };
655
656 return pidCode;
657}
658
659
660//_____________________________________________________________________________
661void AliAnalysisTaskHighPtDeDx::ProcessMCTruthESD()
662{
663 // Fill the special MC histogram with the MC truth info
664
665 Short_t trackmult = 0;
666 Short_t nadded = 0;
667 const Int_t nTracksMC = fMCStack->GetNtrack();
668
669 for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
670
671 //Cuts
672 if(!(fMCStack->IsPhysicalPrimary(iTracks)))
673 continue;
674
675 TParticle* trackMC = fMCStack->Particle(iTracks);
676
677 Double_t chargeMC = trackMC->GetPDG()->Charge();
678 if (chargeMC == 0)
679 continue;
680
681 if (TMath::Abs(trackMC->Eta()) > fEtaCut )
682 continue;
683
684 trackmult++;
685
686 // "charge:pt:p:eta:phi:pidCode"
687 Float_t ptMC = trackMC->Pt();
688 Float_t pMC = trackMC->P();
689 Float_t etaMC = trackMC->Eta();
690 Float_t phiMC = trackMC->Phi();
691
692 Int_t pdgCode = trackMC->GetPdgCode();
693 Short_t pidCodeMC = 0;
694 pidCodeMC = GetPidCode(pdgCode);
695
696 // Here we want to add some of the MC histograms!
697
698 // And therefore we first cut here!
699 if (trackMC->Pt() < fMinPt) {
700
701 // Keep small fraction of low pT tracks
702 if(fRandom->Rndm() > fLowPtFraction)
703 continue;
704 } // else {
705 // Here we want to add the high pt part of the MC histograms!
706 // }
707
708 if(fTreeOption) {
709
710 DeDxTrackMC* track = new((*fTrackArrayMC)[nadded]) DeDxTrackMC();
711 nadded++;
712
713 track->pMC = pMC;
714 track->ptMC = ptMC;
715 track->etaMC = etaMC;
716 track->phiMC = phiMC;
717 track->qMC = Short_t(chargeMC);
718 track->pidMC = pidCodeMC;
719 track->pdgMC = pdgCode;
720 }
721
722 }//MC track loop
723
724 if(fTreeOption) {
725
726 Sort(fTrackArrayMC, kTRUE);
727
728 fEvent->trackmultMC = trackmult;
729 fEvent->nMC = nadded;
730 }
731
732}
733
734//_____________________________________________________________________________
735void AliAnalysisTaskHighPtDeDx::ProcessMCTruthAOD()
736{
737 // Fill the special MC histogram with the MC truth info
738
739 Short_t trackmult = 0;
740 Short_t nadded = 0;
741 const Int_t nTracksMC = fMCArray->GetEntriesFast();
742
743 for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
744
745 AliAODMCParticle* trackMC = dynamic_cast<AliAODMCParticle*>(fMCArray->At(iTracks));
746
747 //Cuts
748 if(!(trackMC->IsPhysicalPrimary()))
749 continue;
750
751
752 Double_t chargeMC = trackMC->Charge();
753 if (chargeMC == 0)
754 continue;
755
756 if (TMath::Abs(trackMC->Eta()) > fEtaCut )
757 continue;
758
759 trackmult++;
760
761 // "charge:pt:p:eta:phi:pidCode"
762 Float_t ptMC = trackMC->Pt();
763 Float_t pMC = trackMC->P();
764 Float_t etaMC = trackMC->Eta();
765 Float_t phiMC = trackMC->Phi();
766
767 Int_t pdgCode = trackMC->PdgCode();
768 Short_t pidCodeMC = 0;
769 pidCodeMC = GetPidCode(pdgCode);
770
771 // Here we want to add some of the MC histograms!
772
773 // And therefore we first cut here!
774 if (trackMC->Pt() < fMinPt) {
775
776 // Keep small fraction of low pT tracks
777 if(fRandom->Rndm() > fLowPtFraction)
778 continue;
779 } // else {
780 // Here we want to add the high pt part of the MC histograms!
781 // }
782
783 if(fTreeOption) {
784
785 DeDxTrackMC* track = new((*fTrackArrayMC)[nadded]) DeDxTrackMC();
786 nadded++;
787
788 track->pMC = pMC;
789 track->ptMC = ptMC;
790 track->etaMC = etaMC;
791 track->phiMC = phiMC;
792 track->qMC = Short_t(chargeMC);
793 track->pidMC = pidCodeMC;
794 track->pdgMC = pdgCode;
795 }
796
797 }//MC track loop
798
799 if(fTreeOption) {
800
801 Sort(fTrackArrayMC, kTRUE);
802
803 fEvent->trackmultMC = trackmult;
804 fEvent->nMC = nadded;
805 }
806
807}
808
809//_____________________________________________________________________________
810Short_t AliAnalysisTaskHighPtDeDx::GetPythiaEventProcessType(Int_t pythiaType) {
811 //
812 // Get the process type of the event. PYTHIA
813 //
814 // source PWG0 dNdpt
815
816 Short_t globalType = -1; //init
817
818 if(pythiaType==92||pythiaType==93){
819 globalType = 2; //single diffractive
820 }
821 else if(pythiaType==94){
822 globalType = 3; //double diffractive
823 }
824 //else if(pythiaType != 91){ // also exclude elastic to be sure... CKB??
825 else {
826 globalType = 1; //non diffractive
827 }
828 return globalType;
829}
830
831//_____________________________________________________________________________
832Short_t AliAnalysisTaskHighPtDeDx::GetDPMjetEventProcessType(Int_t dpmJetType) {
833 //
834 // get the process type of the event. PHOJET
835 //
836 //source PWG0 dNdpt
837 // can only read pythia headers, either directly or from cocktalil header
838 Short_t globalType = -1;
839
840 if (dpmJetType == 1 || dpmJetType == 4) { // explicitly inelastic plus central diffraction
841 globalType = 1;
842 }
843 else if (dpmJetType==5 || dpmJetType==6) {
844 globalType = 2;
845 }
846 else if (dpmJetType==7) {
847 globalType = 3;
848 }
849 return globalType;
850}
851
852//_____________________________________________________________________________
853ULong64_t AliAnalysisTaskHighPtDeDx::GetEventIdAsLong(AliVHeader* header) const
854{
855 // To have a unique id for each event in a run!
856 // Modified from AliRawReader.h
857 return ((ULong64_t)header->GetBunchCrossNumber()+
858 (ULong64_t)header->GetOrbitNumber()*3564+
859 (ULong64_t)header->GetPeriodNumber()*16777215*3564);
860}
861
862
863//____________________________________________________________________
864TParticle* AliAnalysisTaskHighPtDeDx::FindPrimaryMother(AliStack* stack, Int_t label)
865{
866 //
867 // Finds the first mother among the primary particles of the particle identified by <label>,
868 // i.e. the primary that "caused" this particle
869 //
870 // Taken from AliPWG0Helper class
871 //
872
873 Int_t motherLabel = FindPrimaryMotherLabel(stack, label);
874 if (motherLabel < 0)
875 return 0;
876
877 return stack->Particle(motherLabel);
878}
879
880//____________________________________________________________________
881Int_t AliAnalysisTaskHighPtDeDx::FindPrimaryMotherLabel(AliStack* stack, Int_t label)
882{
883 //
884 // Finds the first mother among the primary particles of the particle identified by <label>,
885 // i.e. the primary that "caused" this particle
886 //
887 // returns its label
888 //
889 // Taken from AliPWG0Helper class
890 //
891 const Int_t nPrim = stack->GetNprimary();
892
893 while (label >= nPrim) {
894
895 //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0));
896
897 TParticle* particle = stack->Particle(label);
898 if (!particle) {
899
900 AliDebugGeneral("FindPrimaryMotherLabel", AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack.", label));
901 return -1;
902 }
903
904 // find mother
905 if (particle->GetMother(0) < 0) {
906
907 AliDebugGeneral("FindPrimaryMotherLabel", AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label));
908 return -1;
909 }
910
911 label = particle->GetMother(0);
912 }
913
914 return label;
915}
916
917//____________________________________________________________________
918AliAODMCParticle* AliAnalysisTaskHighPtDeDx::FindPrimaryMotherAOD(AliAODMCParticle* startParticle)
919{
920 //
921 // Finds the first mother among the primary particles of the particle identified by <label>,
922 // i.e. the primary that "caused" this particle
923 //
924 // Taken from AliPWG0Helper class
925 //
926
927 AliAODMCParticle* mcPart = startParticle;
928
929 while (mcPart)
930 {
931
932 if(mcPart->IsPrimary())
933 return mcPart;
934
935 Int_t mother = mcPart->GetMother();
936
937 mcPart = dynamic_cast<AliAODMCParticle*>(fMCArray->At(mother));
938 }
939
940 return 0;
941}
942
943//_____________________________________________________________________________
944void AliAnalysisTaskHighPtDeDx::Sort(TClonesArray* array, Bool_t isMC)
945{
946 const Int_t n = array->GetEntriesFast();
947 if(n==0) {
948
949 if(isMC)
950 fEvent->ptmaxMC = 0;
951 else
952 fEvent->ptmax = 0;
953
954 return;
955 }
956
957 Float_t ptArray[n];
958 Int_t indexArray[n];
959
960 for(Int_t i = 0; i < n; i++) {
961
962 Float_t pt = 0;
963 if(isMC) {
964 DeDxTrackMC* track = (DeDxTrackMC*)array->At(i);
965 pt = track->ptMC;
966 } else {
967 DeDxTrack* track = (DeDxTrack*)array->At(i);
968 pt = track->pt;
969 }
970 ptArray[i] = pt;
971 indexArray[i] = i;
972 }
973
974 TMath::Sort(n, ptArray, indexArray);
975
976 // set max pt
977 if(isMC) {
978 DeDxTrackMC* track = (DeDxTrackMC*)array->At(indexArray[0]);
979 fEvent->ptmaxMC = track->ptMC;
980 } else {
981 DeDxTrack* track = (DeDxTrack*)array->At(indexArray[0]);
982 fEvent->ptmax = track->pt;
983 }
984
985 // set order of each track
986 for(Int_t i = 0; i < n; i++) {
987
988 if(isMC) {
989 DeDxTrackMC* track = (DeDxTrackMC*)array->At(indexArray[i]);
990 track->orderMC = i;
991 } else {
992 DeDxTrack* track = (DeDxTrack*)array->At(indexArray[i]);
993 track->order = i;
994 }
995 }
996}
997//__________________________________________________________________
998void AliAnalysisTaskHighPtDeDx::ProduceArrayTrksESD( AliESDEvent *ESDevent, AnalysisMode analysisMode ){
999
1000 const Int_t nESDTracks = ESDevent->GetNumberOfTracks();
1001 Int_t trackmult=0;
1002 Int_t nadded=0;
1003 if( analysisMode == kGlobalTrk ){
1004 if(fTrackArrayGlobalPar)
1005 fTrackArrayGlobalPar->Clear();
1006 } else if( analysisMode == kTPCTrk ){
1007 if(fTrackArrayTPCPar)
1008 fTrackArrayTPCPar->Clear();
1009 }
1010
1011 if( analysisMode==kGlobalTrk ){
1012
1013 for(Int_t iT = 0; iT < nESDTracks; iT++) {
1014
1015 AliESDtrack* esdTrack = ESDevent->GetTrack(iT);
1016
1017
1018 if(TMath::Abs(esdTrack->Eta()) > fEtaCut)
1019 continue;
1020
1021 UShort_t filterFlag = 0;
1022 Bool_t filterCut_Set1 = kFALSE;//parameters from global tracks, with TPC cuts (filter bit =1 in AOD)
1023 Bool_t filterCut_Set2 = kFALSE;//parameters from global tracks, cuts tpc+its 2010 W/O golden cuts
1024 Bool_t filterCut_Set3 = kFALSE;//parameters from global tracks, cuts its+tpc 2010 WITH golden cuts
1025
1026 UInt_t selectDebug = 0;
1027 if (fTrackFilterGolden) {
1028 selectDebug = fTrackFilterGolden->IsSelected(esdTrack);
1029 if (selectDebug) {
1030 filterFlag +=1;
1031 filterCut_Set3=kTRUE;
1032 }
1033 }
1034
1035 if (fTrackFilterTPC) {
1036
1037 selectDebug = fTrackFilterTPC->IsSelected(esdTrack);
1038 if (selectDebug){//only tracks which pass the TPC-only track cuts
1039 trackmult++;
1040 filterFlag +=2;
1041 filterCut_Set1=kTRUE;
1042
1043 }
1044
1045 }
1046
1047 if (fTrackFilter) {
1048 selectDebug = fTrackFilter->IsSelected(esdTrack);
1049 if (selectDebug) {
1050 filterCut_Set2=kTRUE;
1051 }
1052 }
1053
1054
1055 if(filterFlag==0)
1056 continue;
1057
1058
1059 //
1060 // Track was accepted
1061 //
1062
1063 // Here we want to add histograms!
1064
1065 if (esdTrack->Pt() < fMinPt) {
1066
1067 // Keep small fraction of low pT tracks
1068 if(fRandom->Rndm() > fLowPtFraction)
1069 continue;
1070 } // else {
1071 // Here we want to add the high pt part of the histograms!
1072 // }
1073
1074 Short_t charge = esdTrack->Charge();
1075 Float_t pt = esdTrack->Pt();
1076 Float_t p = esdTrack->P();
1077 Float_t eta = esdTrack->Eta();
1078 Float_t phi = esdTrack->Phi();
1079 Short_t ncl = esdTrack->GetTPCsignalN();
1080 Short_t neff = Short_t(esdTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
1081 // Short_t nclf = esdTrack->GetTPCNclsF();
1082 Float_t dedx = esdTrack->GetTPCsignal();
1083 Float_t tpcchi = 0;
1084 if(esdTrack->GetTPCNcls() > 0)
1085 tpcchi = esdTrack->GetTPCchi2()/Float_t(esdTrack->GetTPCNcls());
1086 Float_t b[2];
1087 Float_t bCov[3];
1088 esdTrack->GetImpactParameters(b,bCov);
1089 Float_t dcaxy = b[0];
1090 Float_t dcaz = b[1];
1091 Double_t p_con[3] = {0, 0, 0};
1092 esdTrack->GetConstrainedPxPyPz(p_con);
1093
1094
1095 // Float_t pt_con = (Float_t)TMath::Sqrt(p_con[0]*p_con[0] + p_con[1]*p_con[1]);
1096 // const AliExternalTrackParam* tpcParam = esdTrack->GetTPCInnerParam();
1097 // Float_t pttpc = tpcParam->Pt();
1098 // Float_t ptpc = tpcParam->P();
1099
1100 Float_t ptMC = 0;
1101 Short_t pidCode = 0; // 0 = real data / no mc track!
1102 Short_t primaryFlag = 0; // 0 = real data / not primary mc track
1103 Int_t pdgMother = 0;
1104
1105
1106 //with Globaltrack parameters????
1107 if(fAnalysisMC) {
1108
1109 const Int_t label = TMath::Abs(esdTrack->GetLabel());
1110 TParticle* mcTrack = fMCStack->Particle(label);
1111 if (mcTrack){
1112
1113 if(fMCStack->IsPhysicalPrimary(label))
1114 primaryFlag = 1;
1115
1116 Int_t pdgCode = mcTrack->GetPdgCode();
1117 pidCode = GetPidCode(pdgCode);
1118
1119 ptMC = mcTrack->Pt();
1120
1121 TParticle* mother = FindPrimaryMother(fMCStack, label);
1122 pdgMother = mother->GetPdgCode();
1123 }
1124 }
1125
1126
1127 //TOF
1128 Float_t beta = -99;
1129 if (esdTrack->GetStatus()&AliESDtrack::kTOFpid){
1130 if ((esdTrack->GetIntegratedLength() != 0) && (esdTrack->GetTOFsignal() != 0))
1131 beta = esdTrack->GetIntegratedLength()/esdTrack->GetTOFsignal()/fgkClight;
1132 }
1133
1134 if(fTreeOption) {
1135
1136 DeDxTrack* track = new((*fTrackArrayGlobalPar)[nadded]) DeDxTrack();
1137 nadded++;
1138
1139 track->p = p;
1140 track->pt = pt;
1141 // track->ptcon = pt_con;
1142 // track->tpcchi = tpcchi;
1143 track->eta = eta;
1144 track->phi = phi;
1145 track->q = charge;
1146 track->filter = filterFlag;
1147 track->ncl = ncl;
1148 track->neff = neff;
1149 track->dedx = dedx;
1150 track->beta = beta;
1151 track->dcaxy = dcaxy;
1152 track->dcaz = dcaz;
1153 track->pid = pidCode;
1154 track->primary = primaryFlag;
1155 track->pttrue = ptMC;
1156 track->mother = pdgMother;
1157 track->filterset1 = filterCut_Set1;
1158 track->filterset2 = filterCut_Set2;
1159 track->filterset3 = filterCut_Set3;
1160
1161
1162 }
1163 }//end of track loop
1164
1165 }//end global
1166
1167 else if( analysisMode==kTPCTrk ){
1168 const AliESDVertex *vtxSPD = ESDevent->GetPrimaryVertexSPD();
1169 if( vtxSPD->GetNContributors() < 1 || TMath::Abs(vtxSPD->GetZ()) > 10.0 ) return;
1170
1171
1172 for(Int_t iT = 0; iT < nESDTracks; iT++) {
1173
1174 AliESDtrack* esdTrack = ESDevent->GetTrack(iT);
1175
1176 AliESDtrack *trackTPC = AliESDtrackCuts::GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(ESDevent),esdTrack->GetID());
1177
1178 if(!trackTPC) continue;
1179
1180 UInt_t selectDebug = 0;
1181 selectDebug = fTrackFilterTPC->IsSelected(trackTPC);
1182 if(selectDebug==0) continue;
1183
1184 if(trackTPC->Pt()>0.){
1185 // only constrain tracks above threshold
1186 AliExternalTrackParam exParam;
1187 // take the B-field from the ESD, no 3D fieldMap available at this point
1188 Bool_t relate = false;
1189 relate = trackTPC->RelateToVertexTPC(vtxSPD,ESDevent->GetMagneticField(),
1190 kVeryBig,&exParam);
1191 if(!relate){
1192 delete trackTPC;
1193 continue;
1194 }
1195 trackTPC->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),
1196 exParam.GetCovariance());
1197 }
1198 else continue;
1199
1200
1201
1202
1203 if(TMath::Abs(trackTPC->Eta()) > fEtaCut)
1204 continue;
1205
1206 //
1207 // Track was accepted
1208 //
1209
1210 // Here we want to add histograms!
1211
1212 if (trackTPC->Pt() < fMinPt) {
1213
1214 // Keep small fraction of low pT tracks
1215 if(fRandom->Rndm() > fLowPtFraction)
1216 continue;
1217 } // else {
1218 // Here we want to add the high pt part of the histograms!
1219 // }
1220
1221 Short_t charge = trackTPC->Charge();
1222 Float_t pt = trackTPC->Pt();
1223 Float_t p = trackTPC->P();
1224 Float_t eta = trackTPC->Eta();
1225 Float_t phi = trackTPC->Phi();
1226 Short_t ncl = trackTPC->GetTPCsignalN();
1227 Short_t neff = Short_t(trackTPC->GetTPCClusterInfo(2, 1)); // effective track length for pT res
1228 // Short_t nclf = esdTrack->GetTPCNclsF();
1229 Float_t dedx = trackTPC->GetTPCsignal();
1230 Float_t tpcchi = 0;
1231 if(trackTPC->GetTPCNcls() > 0)
1232 tpcchi = trackTPC->GetTPCchi2()/Float_t(trackTPC->GetTPCNcls());
1233 //Float_t b[2];
1234 //Float_t bCov[3];
1235 //trackTPC->GetImpactParameters(b,bCov);
1236 //Float_t dcaxy = b[0];
1237 //Float_t dcaz = b[1];
1238
1239 Float_t dcaxy = 0.;
1240 Float_t dcaz = 0.;
1241 trackTPC->GetImpactParameters(dcaxy,dcaz);
1242
1243
1244 Double_t p_con[3] = {0, 0, 0};
1245 trackTPC->GetConstrainedPxPyPz(p_con);
1246
1247
1248 // Float_t pt_con = (Float_t)TMath::Sqrt(p_con[0]*p_con[0] + p_con[1]*p_con[1]);
1249 // const AliExternalTrackParam* tpcParam = esdTrack->GetTPCInnerParam();
1250 // Float_t pttpc = tpcParam->Pt();
1251 // Float_t ptpc = tpcParam->P();
1252
1253 Float_t ptMC = 0;
1254 Short_t pidCode = 0; // 0 = real data / no mc track!
1255 Short_t primaryFlag = 0; // 0 = real data / not primary mc track
1256 Int_t pdgMother = 0;
1257
1258
1259 //with Globaltrack parameters????
1260 if(fAnalysisMC) {
1261
1262 const Int_t label = TMath::Abs(esdTrack->GetLabel());
1263 TParticle* mcTrack = fMCStack->Particle(label);
1264 if (mcTrack){
1265
1266 if(fMCStack->IsPhysicalPrimary(label))
1267 primaryFlag = 1;
1268
1269 Int_t pdgCode = mcTrack->GetPdgCode();
1270 pidCode = GetPidCode(pdgCode);
1271
1272 ptMC = mcTrack->Pt();
1273
1274 TParticle* mother = FindPrimaryMother(fMCStack, label);
1275 pdgMother = mother->GetPdgCode();
1276 }
1277 }
1278
1279
1280 //TOF
1281 Float_t beta = -99;
1282 if (esdTrack->GetStatus()&AliESDtrack::kTOFpid){
1283 if ((esdTrack->GetIntegratedLength() != 0) && (esdTrack->GetTOFsignal() != 0))
1284 beta = esdTrack->GetIntegratedLength()/esdTrack->GetTOFsignal()/fgkClight;
1285 }
1286
1287 if(fTreeOption) {
1288
1289 DeDxTrack* tracktpc = new((*fTrackArrayTPCPar)[nadded]) DeDxTrack();
1290 nadded++;
1291
1292 tracktpc->p = p;
1293 tracktpc->pt = pt;
1294 // track->ptcon = pt_con;
1295 // track->tpcchi = tpcchi;
1296 tracktpc->eta = eta;
1297 tracktpc->phi = phi;
1298 tracktpc->q = charge;
1299 tracktpc->filter = 1;
1300 tracktpc->ncl = ncl;
1301 tracktpc->neff = neff;
1302 tracktpc->dedx = dedx;
1303 tracktpc->beta = beta;//computed with Global tracks
1304 tracktpc->dcaxy = dcaxy;
1305 tracktpc->dcaz = dcaz;
1306 tracktpc->pid = pidCode;
1307 tracktpc->primary = primaryFlag;
1308 tracktpc->pttrue = ptMC;
1309 tracktpc->mother = pdgMother;
1310 tracktpc->filterset1 = 0;
1311 tracktpc->filterset2 = 0;
1312 tracktpc->filterset3 = 0;
1313
1314 }
1315
1316
1317 if(trackTPC)
1318 delete trackTPC;
1319
1320
1321 }//end of track loop
1322
1323 }//end of: if isglobal==kFALSE
1324
1325
1326 if(fTreeOption) {
1327
1328 if( analysisMode==kGlobalTrk ){
1329 Sort(fTrackArrayGlobalPar, kFALSE);
1330
1331 fEvent->trackmult = trackmult;
1332 fEvent->n = nadded;
1333 }
1334 else if( analysisMode==kTPCTrk ){
1335 Sort(fTrackArrayTPCPar, kFALSE);
1336 }
1337
1338 }
1339
1340
1341}
1342//__________________________________________________________________
1343void AliAnalysisTaskHighPtDeDx::ProduceArrayTrksAOD( AliAODEvent *AODevent, AnalysisMode analysisMode ){
1344
1345
1346 Int_t trackmult = 0; // no pt cuts
1347 Int_t nadded = 0;
1348 Int_t nAODTracks = AODevent->GetNumberOfTracks();
1349 if( analysisMode == kGlobalTrk ){
1350
1351 if(fTrackArrayGlobalPar)
1352 fTrackArrayGlobalPar->Clear();
1353
1354 }
1355 if( analysisMode == kTPCTrk ){
1356 if(fTrackArrayTPCPar)
1357 fTrackArrayTPCPar->Clear();
1358 }
1359
1360
1361
1362
1363 //const AliAODVertex* vertexSPD= (AliAODVertex*)AODevent->GetPrimaryVertexSPD();//GetPrimaryVertex()
1364 //if( vertexSPD->GetNContributors() < 1 || TMath::Abs( vertexSPD->GetZ()) > 10.0 ) return;
1365
1366
1367
1368 if( analysisMode==kGlobalTrk ){
1369
1370
1371 const AliAODVertex* vertex = AODevent->GetPrimaryVertex();
1372 for(Int_t iT = 0; iT < nAODTracks; iT++) {
1373
1374 AliAODTrack* aodTrack = AODevent->GetTrack(iT);
1375
1376
1377
1378 UShort_t filterFlag = 0;
1379 Bool_t filterCut_Set1 = kFALSE;//parameters from global tracks, with TPC cuts (filter bit =1 in AOD)
1380 Bool_t filterCut_Set2 = kFALSE;//parameters from global tracks, cuts tpc+its 2010 W/O golden cuts
1381 Bool_t filterCut_Set3 = kFALSE;//parameters from global tracks, cuts its+tpc 2010 WITH golden cuts
1382
1383
1384
1385 if (fTrackFilterGolden) {
1386
1387 // ITSTPC2010 cuts is bit 32 according to above macro, new versions of aliroot includes the golden cuts
1388 if(aodTrack->TestFilterBit(32)) {
1389 filterFlag +=1;
1390 filterCut_Set3 = kTRUE;
1391 }
1392 }
1393
1394
1395 if (fTrackFilterTPC) {
1396 // TPC only cuts is bit 1 according to above macro
1397 // Alex always uses 128, NOTE: FILTER 128 ARE TPC TRACKS (TPC PARAMETERS) CONTRAINED TO THE SPD VERTEX,
1398 if(aodTrack->TestFilterBit(1)){
1399 filterFlag +=2;
1400 filterCut_Set1 = kTRUE;
1401 trackmult++;
1402
1403 }
1404 }
1405
1406
1407
1408 if(filterFlag==0)
1409 continue;
1410
1411
1412 Double_t b[2], cov[3];
1413 if (!(aodTrack->PropagateToDCA(vertex, AODevent->GetMagneticField(), kVeryBig, b, cov)))
1414 filterFlag = 32; // propagation failed!!!!!;
1415 Float_t dcaxy = b[0];
1416 Float_t dcaz = b[1];
1417
1418
1419
1420
1421
1422 // As I understand this routine recalculates the momentum so it should be called first!
1423 //Double_t b[2], cov[3];
1424
1425
1426 //if(!aodTrack->PropagateToDCA(vertex, AODevent->GetMagneticField(), kVeryBig, b, cov))
1427 // filterFlag = 32; // propagation failed!!!!!
1428
1429 if(TMath::Abs(aodTrack->Eta()) > fEtaCut)
1430 continue;
1431 if (aodTrack->Pt() < fMinPt) {
1432
1433 // Keep small fraction of low pT tracks
1434 if(fRandom->Rndm() > fLowPtFraction)
1435 continue;
1436 } // else {
1437 // Here we want to add the high pt part of the histograms!
1438 // }
1439
1440
1441
1442
1443 Short_t charge = aodTrack->Charge();
1444 Float_t pt = aodTrack->Pt();
1445 Float_t p = aodTrack->P();
1446 Float_t eta = aodTrack->Eta();
1447 Float_t phi = aodTrack->Phi();
1448 AliAODPid* aodPid = aodTrack->GetDetPid();
1449 Short_t ncl = -10;
1450 Short_t neff = 0; // This is not yet there! Short_t(aodTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
1451 // Short_t nclf = aodTrack->GetTPCNclsF();
1452 Float_t dedx = -10;
1453 Float_t beta = -99;
1454 if(aodPid) {
1455 ncl = aodPid->GetTPCsignalN();
1456 dedx = aodPid->GetTPCsignal();
1457 //TOF
1458 if (aodTrack->GetStatus()&AliESDtrack::kTOFpid){
1459 Double_t tof[5];
1460 aodPid->GetIntegratedTimes(tof);
1461 beta = tof[0]/aodPid->GetTOFsignal();
1462 }
1463 }
1464 // Float_t tpcchi = aodTrack->Chi2perNDF();
1465
1466 // Double_t p_con[3] = {0, 0, 0};
1467 // aodTrack->GetConstrainedPxPyPz(p_con);
1468 // Float_t pt_con = 0; // This is not there! (Float_t)TMath::Sqrt(p_con[0]*p_con[0] + p_con[1]*p_con[1]);
1469 // const AliExternalTrackParam* tpcParam = aodTrack->GetTPCInnerParam();
1470 // Float_t pttpc = tpcParam->Pt();
1471 // Float_t ptpc = tpcParam->P();
1472
1473 Float_t ptMC = 0;
1474 Short_t pidCode = 0; // 0 = real data / no mc track!
1475 Short_t primaryFlag = 0; // 0 = real data / not primary mc track
1476 Int_t pdgMother = 0;
1477
1478 if(fAnalysisMC) {
1479
1480 const Int_t label = TMath::Abs(aodTrack->GetLabel());
1481
1482 AliAODMCParticle* mcTrack = dynamic_cast<AliAODMCParticle*>(fMCArray->At(label));
1483 if (mcTrack){
1484
1485 if(mcTrack->IsPhysicalPrimary())
1486 primaryFlag = 1;
1487
1488 Int_t pdgCode = mcTrack->GetPdgCode();
1489 pidCode = GetPidCode(pdgCode);
1490
1491 ptMC = mcTrack->Pt();
1492
1493 AliAODMCParticle* mother = FindPrimaryMotherAOD(mcTrack);
1494 pdgMother = mother->GetPdgCode();
1495 }
1496 }
1497
1498
1499 if(fTreeOption) {
1500
1501 DeDxTrack* track = new((*fTrackArrayGlobalPar)[nadded]) DeDxTrack();
1502 nadded++;
1503
1504 track->p = p;
1505 track->pt = pt;
1506 // track->ptcon = pt_con;
1507 // track->tpcchi = tpcchi;
1508 track->eta = eta;
1509 track->phi = phi;
1510 track->q = charge;
1511 track->filter = filterFlag;
1512 track->ncl = ncl;
1513 track->neff = neff;
1514 track->dedx = dedx;
1515 track->beta = beta;
1516 track->dcaxy = dcaxy;
1517 track->dcaz = dcaz;
1518 track->pid = pidCode;
1519 track->primary = primaryFlag;
1520 track->pttrue = ptMC;
1521 track->mother = pdgMother;
1522 track->filterset1 = filterCut_Set1;
1523 track->filterset2 = filterCut_Set2;
1524 track->filterset3 = filterCut_Set3;
1525 }
1526 }//end of track loop
1527
1528
1529
1530 }//end of global track analysis
1531
1532
1533
1534 else if( analysisMode==kTPCTrk ){
1535
1536 const AliAODVertex* vertexSPD= (AliAODVertex*)AODevent->GetPrimaryVertexSPD();//GetPrimaryVertex()
1537 if( vertexSPD->GetNContributors() < 1 || TMath::Abs(vertexSPD->GetZ()) > 10.0 ) return;
1538
1539
1540 //const AliAODVertex* vertex = AODevent->GetPrimaryVertex();
1541 for(Int_t iT = 0; iT < nAODTracks; iT++) {
1542
1543 AliAODTrack* aodTrack = AODevent->GetTrack(iT);
1544
1545
1546
1547 UShort_t filterFlag = 0;
1548 Bool_t filterCut_Set1 = kFALSE;//parameters from global tracks, with TPC cuts (filter bit =1 in AOD)
1549 Bool_t filterCut_Set2 = kFALSE;//parameters from global tracks, cuts tpc+its 2010 W/O golden cuts
1550 Bool_t filterCut_Set3 = kFALSE;//parameters from global tracks, cuts its+tpc 2010 WITH golden cuts
1551
1552
1553
1554 // TPC only cuts is bit 1 according to above macro
1555 // Alex always uses 128, NOTE: FILTER 128 ARE TPC TRACKS (TPC PARAMETERS) CONTRAINED TO THE SPD VERTEX,
1556 if(aodTrack->TestFilterBit(128)){
1557 filterFlag +=2;
1558 trackmult++;
1559 }
1560
1561
1562 if(filterFlag==0)
1563 continue;
1564
1565
1566 Double_t b[2], cov[3];
1567 //AliAODVertex* vertex = AODevent->GetPrimaryVertex();
1568 //Double_t b[2] = {-99., -99.};
1569 //Double_t bCov[3] = {-99., -99., -99.};
1570 //aodTrack->PropagateToDCA(aod->GetPrimaryVertex(), aod->GetMagneticField(), 100., b, bCov);
1571 if (!(aodTrack->PropagateToDCA(vertexSPD, AODevent->GetMagneticField(), kVeryBig, b, cov)))
1572 filterFlag = 32; // propagation failed!!!!!;
1573 Float_t dcaxy = b[0];
1574 Float_t dcaz = b[1];
1575
1576
1577
1578
1579 // As I understand this routine recalculates the momentum so it should be called first!
1580 //Double_t b[2], cov[3];
1581
1582
1583 //if(!aodTrack->PropagateToDCA(vertex, AODevent->GetMagneticField(), kVeryBig, b, cov))
1584 // filterFlag = 32; // propagation failed!!!!!
1585
1586 if(TMath::Abs(aodTrack->Eta()) > fEtaCut)
1587 continue;
1588 if (aodTrack->Pt() < fMinPt) {
1589
1590 // Keep small fraction of low pT tracks
1591 if(fRandom->Rndm() > fLowPtFraction)
1592 continue;
1593 } // else {
1594 // Here we want to add the high pt part of the histograms!
1595 // }
1596
1597
1598
1599
1600 Short_t charge = aodTrack->Charge();
1601 Float_t pt = aodTrack->Pt();
1602 Float_t p = aodTrack->P();
1603 Float_t eta = aodTrack->Eta();
1604 Float_t phi = aodTrack->Phi();
1605 AliAODPid* aodPid = aodTrack->GetDetPid();
1606 Short_t ncl = -10;
1607 Short_t neff = 0; // This is not yet there! Short_t(aodTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
1608 // Short_t nclf = aodTrack->GetTPCNclsF();
1609 Float_t dedx = -10;
1610 Float_t beta = -99;
1611 if(aodPid) {
1612 ncl = aodPid->GetTPCsignalN();
1613 dedx = aodPid->GetTPCsignal();
1614 //TOF
1615 if (aodTrack->GetStatus()&AliESDtrack::kTOFpid){
1616 Double_t tof[5];
1617 aodPid->GetIntegratedTimes(tof);
1618 beta = tof[0]/aodPid->GetTOFsignal();
1619 }
1620 }
1621 // Float_t tpcchi = aodTrack->Chi2perNDF();
1622
1623 // Double_t p_con[3] = {0, 0, 0};
1624 // aodTrack->GetConstrainedPxPyPz(p_con);
1625 // Float_t pt_con = 0; // This is not there! (Float_t)TMath::Sqrt(p_con[0]*p_con[0] + p_con[1]*p_con[1]);
1626 // const AliExternalTrackParam* tpcParam = aodTrack->GetTPCInnerParam();
1627 // Float_t pttpc = tpcParam->Pt();
1628 // Float_t ptpc = tpcParam->P();
1629
1630 Float_t ptMC = 0;
1631 Short_t pidCode = 0; // 0 = real data / no mc track!
1632 Short_t primaryFlag = 0; // 0 = real data / not primary mc track
1633 Int_t pdgMother = 0;
1634
1635 if(fAnalysisMC) {
1636
1637 const Int_t label = TMath::Abs(aodTrack->GetLabel());
1638
1639 AliAODMCParticle* mcTrack = dynamic_cast<AliAODMCParticle*>(fMCArray->At(label));
1640 if (mcTrack){
1641
1642 if(mcTrack->IsPhysicalPrimary())
1643 primaryFlag = 1;
1644
1645 Int_t pdgCode = mcTrack->GetPdgCode();
1646 pidCode = GetPidCode(pdgCode);
1647
1648 ptMC = mcTrack->Pt();
1649
1650 AliAODMCParticle* mother = FindPrimaryMotherAOD(mcTrack);
1651 pdgMother = mother->GetPdgCode();
1652 }
1653 }
1654
1655
1656 if(fTreeOption) {
1657
1658 DeDxTrack* tracktpc = new((*fTrackArrayTPCPar)[nadded]) DeDxTrack();
1659 nadded++;
1660
1661 tracktpc->p = p;
1662 tracktpc->pt = pt;
1663 // track->ptcon = pt_con;
1664 // track->tpcchi = tpcchi;
1665 tracktpc->eta = eta;
1666 tracktpc->phi = phi;
1667 tracktpc->q = charge;
1668 tracktpc->filter = filterFlag;
1669 tracktpc->ncl = ncl;
1670 tracktpc->neff = neff;
1671 tracktpc->dedx = dedx;
1672 tracktpc->beta = beta;
1673 tracktpc->dcaxy = dcaxy;
1674 tracktpc->dcaz = dcaz;
1675 tracktpc->pid = pidCode;
1676 tracktpc->primary = primaryFlag;
1677 tracktpc->pttrue = ptMC;
1678 tracktpc->mother = pdgMother;
1679 tracktpc->filterset1 = filterCut_Set1;
1680 tracktpc->filterset2 = filterCut_Set2;
1681 tracktpc->filterset3 = filterCut_Set3;
1682 }
1683 }//end of track loop
1684
1685
1686 }//end of global track analysis
1687
1688
1689 if(fTreeOption) {
1690
1691 if( analysisMode==kGlobalTrk ){
1692 Sort(fTrackArrayGlobalPar, kFALSE);
1693
1694 fEvent->trackmult = trackmult;
1695 fEvent->n = nadded;
1696
1697
1698 }
1699 if( analysisMode==kTPCTrk ){
1700 Sort(fTrackArrayTPCPar, kFALSE);
1701 }
1702
1703 }
1704
1705}