]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PWGJE/AliPWG4HighPtSpectra.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Now also the QA task only considers charged primary pions, kaons, etc.(M. Verweij)
[u/mrichter/AliRoot.git] / PWGJE / AliPWG4HighPtSpectra.cxx
... / ...
CommitLineData
1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16//-----------------------------------------------------------------------
17// Efficiency between different steps of the procedure.
18// The ouptut of the task is a AliCFContainer from which the efficiencies
19// can be calculated
20//-----------------------------------------------------------------------
21// Author : Marta Verweij - UU
22//-----------------------------------------------------------------------
23
24
25#ifndef ALIPWG4HIGHPTSPECTRA_CXX
26#define ALIPWG4HIGHPTSPECTRA_CXX
27
28#include "AliPWG4HighPtSpectra.h"
29
30#include "TVector3.h"
31#include <iostream>
32#include "TH1.h"
33#include "TH2.h"
34#include "TH3.h"
35#include "TProfile.h"
36#include "TList.h"
37#include "TChain.h"
38#include "TKey.h"
39#include "TSystem.h"
40#include "TFile.h"
41
42#include "AliAnalysisManager.h"
43#include "AliESDInputHandler.h"
44#include "AliESDtrack.h"
45#include "AliESDtrackCuts.h"
46#include "AliExternalTrackParam.h"
47#include "AliCentrality.h"
48
49#include "AliLog.h"
50
51#include "AliStack.h"
52#include "TParticle.h"
53#include "AliMCEvent.h"
54#include "AliMCEventHandler.h"
55#include "AliCFContainer.h"
56#include "AliGenPythiaEventHeader.h"
57#include "AliGenCocktailEventHeader.h"
58//#include "$ALICE_ROOT/PWG4/JetTasks/AliAnalysisHelperJetTasks.h"
59
60//#include <iostream>
61using namespace std; //required for resolving the 'cout' symbol
62
63ClassImp(AliPWG4HighPtSpectra)
64
65//__________________________________________________________________________
66AliPWG4HighPtSpectra::AliPWG4HighPtSpectra() : AliAnalysisTask("AliPWG4HighPtSpectra", ""),
67 fReadAODData(0),
68 fCFManagerPos(0x0),
69 fCFManagerNeg(0x0),
70 fESD(0x0),
71 fMC(0x0),
72 fStack(0x0),
73 fVtx(0x0),
74 fIsPbPb(0),
75 fCentClass(10),
76 fTrackType(0),
77 fTrackCuts(0x0),
78 fTrackCutsReject(0x0),
79 fSigmaConstrainedMax(100.),
80 fAvgTrials(1),
81 fHistList(0),
82 fNEventAll(0),
83 fNEventSel(0),
84 fNEventReject(0),
85 fh1Centrality(0x0),
86 fh1Xsec(0),
87 fh1Trials(0),
88 fh1PtHard(0),
89 fh1PtHardTrials(0),
90 fPtRelUncertainty1PtPrim(0x0),
91 fPtRelUncertainty1PtSec(0x0)
92{
93 //
94 //Default ctor
95 //
96}
97//___________________________________________________________________________
98AliPWG4HighPtSpectra::AliPWG4HighPtSpectra(const Char_t* name) :
99 AliAnalysisTask(name,""),
100 fReadAODData(0),
101 fCFManagerPos(0x0),
102 fCFManagerNeg(0x0),
103 fESD(0x0),
104 fMC(0x0),
105 fStack(0x0),
106 fVtx(0x0),
107 fIsPbPb(0),
108 fCentClass(10),
109 fTrackType(0),
110 fTrackCuts(0x0),
111 fTrackCutsReject(0x0),
112 fSigmaConstrainedMax(100.),
113 fAvgTrials(1),
114 fHistList(0),
115 fNEventAll(0),
116 fNEventSel(0),
117 fNEventReject(0),
118 fh1Centrality(0x0),
119 fh1Xsec(0),
120 fh1Trials(0),
121 fh1PtHard(0),
122 fh1PtHardTrials(0),
123 fPtRelUncertainty1PtPrim(0x0),
124 fPtRelUncertainty1PtSec(0x0)
125{
126 //
127 // Constructor. Initialization of Inputs and Outputs
128 //
129 AliDebug(2,Form("AliPWG4HighPtSpectra Calling Constructor"));
130 // Input slot #0 works with a TChain ESD
131 DefineInput(0, TChain::Class());
132 // Output slot #0 writes into a TList
133 DefineOutput(0,TList::Class());
134 // Output slot #1, #2 writes into a AliCFContainer
135 DefineOutput(1,AliCFContainer::Class());
136 DefineOutput(2,AliCFContainer::Class());
137 // Output slot #3 writes into a AliESDtrackCuts
138 DefineOutput(3, AliESDtrackCuts::Class());
139 DefineOutput(4, AliESDtrackCuts::Class());
140}
141
142//________________________________________________________________________
143void AliPWG4HighPtSpectra::LocalInit()
144{
145 //
146 // Only called once at beginning
147 //
148 PostData(3,fTrackCuts);
149}
150
151//________________________________________________________________________
152void AliPWG4HighPtSpectra::ConnectInputData(Option_t *)
153{
154 // Connect ESD here
155 // Called once
156 AliDebug(2,Form(">> AliPWG4HighPtSpectra::ConnectInputData \n"));
157
158 TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
159 if (!tree) {
160 AliDebug(2,Form( "ERROR: Could not read chain from input slot 0 \n"));
161 return;
162 }
163
164 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
165
166 if (!esdH) {
167 AliDebug(2,Form("ERROR: Could not get ESDInputHandler"));
168 return;
169 } else
170 fESD = esdH->GetEvent();
171
172 AliMCEventHandler *eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
173 if (!eventHandler) {
174 AliDebug(2,Form( "ERROR: Could not retrieve MC event handler \n"));
175 }
176 else
177 fMC = eventHandler->MCEvent();
178
179}
180
181//________________________________________________________________________
182Bool_t AliPWG4HighPtSpectra::SelectEvent() {
183 //
184 // Decide if event should be selected for analysis
185 //
186
187 // Checks following requirements:
188 // - fESD available
189 // - trigger info from AliPhysicsSelection
190 // - number of reconstructed tracks > 1
191 // - primary vertex reconstructed
192 // - z-vertex < 10 cm
193
194 Bool_t selectEvent = kTRUE;
195
196 //fESD object available?
197 if (!fESD) {
198 AliDebug(2,Form("ERROR: fInputEvent not available\n"));
199 fNEventReject->Fill("noESD",1);
200 selectEvent = kFALSE;
201 return selectEvent;
202 }
203
204 //Trigger
205 UInt_t isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
206 if(!(isSelected&AliVEvent::kMB)) { //Select collison candidates
207 AliDebug(2,Form(" Trigger Selection: event REJECTED ... "));
208 fNEventReject->Fill("Trigger",1);
209 selectEvent = kFALSE;
210 return selectEvent;
211 }
212
213 //Check if number of reconstructed tracks is larger than 1
214 if(!fESD->GetNumberOfTracks() || fESD->GetNumberOfTracks()<2) {
215 fNEventReject->Fill("NTracks<2",1);
216 selectEvent = kFALSE;
217 return selectEvent;
218 }
219
220 //Check if vertex is reconstructed
221 fVtx = fESD->GetPrimaryVertexSPD();
222
223 if(!fVtx) {
224 fNEventReject->Fill("noVTX",1);
225 selectEvent = kFALSE;
226 return selectEvent;
227 }
228
229 if(!fVtx->GetStatus()) {
230 fNEventReject->Fill("VtxStatus",1);
231 selectEvent = kFALSE;
232 return selectEvent;
233 }
234
235 // Need vertex cut
236 // TString vtxName(fVtx->GetName());
237 if(fVtx->GetNContributors()<2) {
238 fNEventReject->Fill("NCont<2",1);
239 selectEvent = kFALSE;
240 return selectEvent;
241 }
242
243 //Check if z-vertex < 10 cm
244 double primVtx[3];
245 fVtx->GetXYZ(primVtx);
246 if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
247 fNEventReject->Fill("ZVTX>10",1);
248 selectEvent = kFALSE;
249 return selectEvent;
250 }
251
252 //Centrality selection should only be done in case of PbPb
253 if(IsPbPb()) {
254 Float_t cent = 0.;
255 if(fCentClass!=CalculateCentrality(fESD) && fCentClass!=10) {
256 fNEventReject->Fill("cent",1);
257 selectEvent = kFALSE;
258 return selectEvent;
259 }
260 else {
261 if(dynamic_cast<AliESDEvent*>(fESD)->GetCentrality()) {
262 cent = dynamic_cast<AliESDEvent*>(fESD)->GetCentrality()->GetCentralityPercentile("V0M");
263 }
264 if(cent>90.) {
265 fNEventReject->Fill("cent>90",1);
266 selectEvent = kFALSE;
267 return selectEvent;
268 }
269 fh1Centrality->Fill(cent);
270 }
271 }
272
273 return selectEvent;
274
275}
276
277//________________________________________________________________________
278Int_t AliPWG4HighPtSpectra::CalculateCentrality(AliESDEvent *esd){
279
280
281 Float_t cent = 999;
282
283 if(esd){
284 if(esd->GetCentrality()){
285 cent = esd->GetCentrality()->GetCentralityPercentile("V0M");
286 }
287 }
288
289 if(cent<0) return 5;
290 if(cent>80)return 4;
291 if(cent>50)return 3;
292 if(cent>30)return 2;
293 if(cent>10)return 1;
294 return 0;
295
296}
297
298//_________________________________________________
299void AliPWG4HighPtSpectra::Exec(Option_t *)
300{
301 //
302 // Main loop function
303 //
304 AliDebug(2,Form(">> AliPWG4HighPtSpectra::Exec \n"));
305
306 // All events without selection
307 fNEventAll->Fill(0.);
308
309 if(!SelectEvent()) {
310 fNEventReject->Fill("NTracks<2",1);
311 // Post output data
312 PostData(0,fHistList);
313 PostData(1,fCFManagerPos->GetParticleContainer());
314 PostData(2,fCFManagerNeg->GetParticleContainer());
315 return;
316 }
317
318 //MCEvent available?
319 //if yes: get stack
320 if(fMC) {
321 AliDebug(2,Form("MC particles: %d", fMC->GetNumberOfTracks()));
322 fStack = fMC->Stack(); //Particles Stack
323 AliDebug(2,Form("MC particles stack: %d", fStack->GetNtrack()));
324 }
325
326 // ---- Get MC Header information (for MC productions in pThard bins) ----
327 Double_t ptHard = 0.;
328 Double_t nTrials = 1; // trials for MC trigger weight for real data
329
330 if(fMC){
331 AliGenPythiaEventHeader* pythiaGenHeader = GetPythiaEventHeader(fMC);
332 if(pythiaGenHeader){
333 nTrials = pythiaGenHeader->Trials();
334 ptHard = pythiaGenHeader->GetPtHard();
335
336 fh1PtHard->Fill(ptHard);
337 fh1PtHardTrials->Fill(ptHard,nTrials);
338
339 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
340 }
341 }
342
343 Int_t nTracks = fESD->GetNumberOfTracks();
344 AliDebug(2,Form("nTracks %d", nTracks));
345
346 if(!fTrackCuts) {
347 fNEventReject->Fill("noTrackCuts",1);
348 // Post output data
349 PostData(0,fHistList);
350 PostData(1,fCFManagerPos->GetParticleContainer());
351 PostData(2,fCFManagerNeg->GetParticleContainer());
352 return;
353 }
354
355 // Selected events for analysis
356 fNEventSel->Fill(0.);
357
358 const Int_t nvar = 4;
359
360 Double_t containerInputRec[nvar] = {0.,0.,0.,0.};
361 Double_t containerInputMC[nvar] = {0.,0.,0.,0.};
362 Double_t containerInputRecMC[nvar] = {0.,0.,0.,0.}; //reconstructed yield as function of MC variable
363
364 //Now go to rec level
365 for (Int_t iTrack = 0; iTrack<nTracks; iTrack++)
366 {
367 //Get track for analysis
368 AliESDtrack *track = 0x0;
369 AliESDtrack *esdtrack = fESD->GetTrack(iTrack);
370 if(!esdtrack) continue;
371
372 if(fTrackType==1) {
373 track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
374 if(!track) continue;
375 }
376 else if(fTrackType==2) {
377 track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
378 if(!track) continue;
379
380 AliExternalTrackParam exParam;
381 Bool_t relate = track->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
382 if( !relate ) {
383 delete track;
384 continue;
385 }
386 track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
387 }
388 else if(fTrackType==7) {
389 //use global constrained track
390 track = new AliESDtrack(*esdtrack);
391 // track = esdtrack;
392 // track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
393
394 }
395 else
396 track = esdtrack;
397
398 if(!track) continue;
399
400 if(fTrackType==2) {
401 //Cut on chi2 of constrained fit
402 if(track->GetConstrainedChi2TPC() > fSigmaConstrainedMax*fSigmaConstrainedMax) {
403 delete track;
404 continue;
405 }
406 }
407
408 if (fTrackCuts->AcceptTrack(track)) {
409
410 if(fTrackType==7) {
411 if(fTrackCutsReject ) {
412 if(fTrackCutsReject->AcceptTrack(track) ) {
413 if(track) delete track;
414 continue;
415 }
416 }
417
418 if(esdtrack->GetConstrainedParam())
419 track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
420 }
421
422 //fill the container
423 containerInputRec[0] = track->Pt();
424 containerInputRec[1] = track->Phi();
425 containerInputRec[2] = track->Eta();
426 containerInputRec[3] = track->GetTPCNclsIter1();
427
428 if(track->GetSign()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
429 if(track->GetSign()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
430
431
432 //Only fill the MC containers if MC information is available
433 if(fMC) {
434 Int_t label = TMath::Abs(track->GetLabel());
435 TParticle *particle = fStack->Particle(label) ;
436 if(!particle) {
437 if(fTrackType==1 || fTrackType==2 || fTrackType==7)
438 delete track;
439 continue;
440 }
441 containerInputRecMC[0] = particle->Pt();
442 containerInputRecMC[1] = particle->Phi();
443 containerInputRecMC[2] = particle->Eta();
444 containerInputRecMC[3] = track->GetTPCNclsIter1();
445
446 //Container with primaries
447 if(fStack->IsPhysicalPrimary(label)) {
448 if(particle->GetPDG()->Charge()>0.) {
449 fCFManagerPos->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
450 }
451 if(particle->GetPDG()->Charge()<0.) {
452 fCFManagerNeg->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
453 }
454
455 //Fill pT resolution plots for primaries
456 fPtRelUncertainty1PtPrim->Fill(containerInputRec[0],containerInputRec[0]*TMath::Sqrt(track->GetSigma1Pt2()));
457
458 }
459
460 //Container with secondaries
461 if (!fStack->IsPhysicalPrimary(label) ) {
462 if(particle->GetPDG()->Charge()>0.) {
463 fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepSecondaries);
464 }
465 if(particle->GetPDG()->Charge()<0.) {
466 fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepSecondaries);
467 }
468 //Fill pT resolution plots for primaries
469 fPtRelUncertainty1PtSec->Fill(containerInputRec[0],containerInputRec[0]*TMath::Sqrt(track->GetSigma1Pt2()));
470 }
471 }
472
473 }//trackCuts global tracks
474
475 if(fTrackType==1 || fTrackType==2 || fTrackType==7) {
476 if(track) delete track;
477 }
478 }//track loop
479
480
481 //Fill MC containers if particles are findable
482 if(fMC) {
483 for(int iPart = 1; iPart<(fMC->GetNumberOfPrimaries()); iPart++) {
484 AliMCParticle *mcPart = (AliMCParticle*)fMC->GetTrack(iPart);
485 if(!mcPart) continue;
486
487 Int_t pdg = mcPart->PdgCode();
488
489 // select charged pions, protons, kaons , electrons, muons
490 if(TMath::Abs(pdg) == 211 || TMath::Abs(pdg) == 2212 || TMath::Abs(pdg) ==
491 321 || TMath::Abs(pdg) == 11 || TMath::Abs(pdg) == 13){
492
493 //fill the container
494 containerInputMC[0] = mcPart->Pt();
495 containerInputMC[1] = mcPart->Phi();
496 containerInputMC[2] = mcPart->Eta();
497 // AliESDtrack *esdtrack = fESD->GetTrack(mcPart->GetLabel());
498 containerInputMC[3] = 159.;
499
500 if(fStack->IsPhysicalPrimary(iPart)) {
501 if(mcPart->Charge()>0. && fCFManagerPos->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
502 if(mcPart->Charge()<0. && fCFManagerNeg->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
503 }
504 }
505 }
506 }
507
508 PostData(0,fHistList);
509 PostData(1,fCFManagerPos->GetParticleContainer());
510 PostData(2,fCFManagerNeg->GetParticleContainer());
511
512}
513//________________________________________________________________________
514Bool_t AliPWG4HighPtSpectra::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
515 //
516 // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
517 // This is to called in Notify and should provide the path to the AOD/ESD file
518 // Copied from AliAnalysisTaskJetSpectrum2
519 //
520
521 TString file(currFile);
522 fXsec = 0;
523 fTrials = 1;
524
525 if(file.Contains("root_archive.zip#")){
526 Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
527 Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
528 file.Replace(pos+1,20,"");
529 }
530 else {
531 // not an archive take the basename....
532 file.ReplaceAll(gSystem->BaseName(file.Data()),"");
533 }
534 // Printf("%s",file.Data());
535
536 TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")); // problem that we cannot really test the existance of a file in a archive so we have to lvie with open error message from root
537 if(!fxsec){
538 // next trial fetch the histgram file
539 fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
540 if(!fxsec){
541 // not a severe condition but inciate that we have no information
542 return kFALSE;
543 }
544 else{
545 // find the tlist we want to be independtent of the name so use the Tkey
546 TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
547 if(!key){
548 fxsec->Close();
549 return kFALSE;
550 }
551 TList *list = dynamic_cast<TList*>(key->ReadObj());
552 if(!list){
553 fxsec->Close();
554 return kFALSE;
555 }
556 fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
557 fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
558 fxsec->Close();
559 }
560 } // no tree pyxsec.root
561 else {
562 TTree *xtree = (TTree*)fxsec->Get("Xsection");
563 if(!xtree){
564 fxsec->Close();
565 return kFALSE;
566 }
567 UInt_t ntrials = 0;
568 Double_t xsection = 0;
569 xtree->SetBranchAddress("xsection",&xsection);
570 xtree->SetBranchAddress("ntrials",&ntrials);
571 xtree->GetEntry(0);
572 fTrials = ntrials;
573 fXsec = xsection;
574 fxsec->Close();
575 }
576 return kTRUE;
577}
578//________________________________________________________________________
579Bool_t AliPWG4HighPtSpectra::Notify()
580{
581 //
582 // Implemented Notify() to read the cross sections
583 // and number of trials from pyxsec.root
584 // Copied from AliAnalysisTaskJetSpectrum2
585 //
586
587 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
588 Float_t xsection = 0;
589 Float_t ftrials = 1;
590
591 fAvgTrials = 1;
592 if(tree){
593 TFile *curfile = tree->GetCurrentFile();
594 if (!curfile) {
595 Error("Notify","No current file");
596 return kFALSE;
597 }
598 if(!fh1Xsec||!fh1Trials){
599 // Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
600 return kFALSE;
601 }
602 PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
603 fh1Xsec->Fill("<#sigma>",xsection);
604 // construct a poor man average trials
605 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
606 if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
607 }
608 return kTRUE;
609}
610
611//________________________________________________________________________
612AliGenPythiaEventHeader* AliPWG4HighPtSpectra::GetPythiaEventHeader(AliMCEvent *mcEvent){
613
614 if(!mcEvent)return 0;
615 AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
616 AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
617 if(!pythiaGenHeader){
618 // cocktail ??
619 AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
620
621 if (!genCocktailHeader) {
622 // AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
623 // AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
624 return 0;
625 }
626 TList* headerList = genCocktailHeader->GetHeaders();
627 for (Int_t i=0; i<headerList->GetEntries(); i++) {
628 pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
629 if (pythiaGenHeader)
630 break;
631 }
632 if(!pythiaGenHeader){
633 AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
634 return 0;
635 }
636 }
637 return pythiaGenHeader;
638
639}
640
641
642//___________________________________________________________________________
643void AliPWG4HighPtSpectra::Terminate(Option_t*)
644{
645 // The Terminate() function is the last function to be called during
646 // a query. It always runs on the client, it can be used to present
647 // the results graphically or save the results to file.
648
649}
650
651//___________________________________________________________________________
652void AliPWG4HighPtSpectra::CreateOutputObjects() {
653 //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
654 //TO BE SET BEFORE THE EXECUTION OF THE TASK
655 //
656 AliDebug(2,Form("CreateOutputObjects CreateOutputObjects of task %s", GetName()));
657
658 Bool_t oldStatus = TH1::AddDirectoryStatus();
659 TH1::AddDirectory(kFALSE);
660
661 //slot #1
662 OpenFile(0);
663 fHistList = new TList();
664 fHistList->SetOwner(kTRUE);
665 fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
666 fHistList->Add(fNEventAll);
667 fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
668 fHistList->Add(fNEventSel);
669
670 fNEventReject = new TH1F("fNEventReject","Reason events are rejectected for analysis",20,0,20);
671 //Set labels
672 fNEventReject->Fill("noESD",0);
673 fNEventReject->Fill("Trigger",0);
674 fNEventReject->Fill("NTracks<2",0);
675 fNEventReject->Fill("noVTX",0);
676 fNEventReject->Fill("VtxStatus",0);
677 fNEventReject->Fill("NCont<2",0);
678 fNEventReject->Fill("ZVTX>10",0);
679 fNEventReject->Fill("cent",0);
680 fNEventReject->Fill("cent>90",0);
681 fHistList->Add(fNEventReject);
682
683 fh1Centrality = new TH1F("fh1Centrality","fh1Centrality; Centrality %",100,0,100);
684 fHistList->Add(fh1Centrality);
685
686 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
687 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
688 fHistList->Add(fh1Xsec);
689
690 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
691 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
692 fHistList->Add(fh1Trials);
693
694 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
695 fHistList->Add(fh1PtHard);
696 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
697 fHistList->Add(fh1PtHardTrials);
698
699 Int_t fgkNPtBins = 100;
700 Float_t kMinPt = 0.;
701 Float_t kMaxPt = 100.;
702 Double_t *binsPt = new Double_t[fgkNPtBins+1];
703 for(Int_t i=0; i<=fgkNPtBins; i++) binsPt[i]=(Double_t)kMinPt + (kMaxPt-kMinPt)/fgkNPtBins*(Double_t)i ;
704
705 Int_t fgkNRel1PtUncertaintyBins=50;
706 Float_t fgkRel1PtUncertaintyMin = 0.;
707 Float_t fgkRel1PtUncertaintyMax = 1.;
708
709 Double_t *binsRel1PtUncertainty=new Double_t[fgkNRel1PtUncertaintyBins+1];
710 for(Int_t i=0; i<=fgkNRel1PtUncertaintyBins; i++) binsRel1PtUncertainty[i]=(Double_t)fgkRel1PtUncertaintyMin + (fgkRel1PtUncertaintyMax-fgkRel1PtUncertaintyMin)/fgkNRel1PtUncertaintyBins*(Double_t)i ;
711
712 fPtRelUncertainty1PtPrim = new TH2F("fPtRelUncertainty1PtPrim","fPtRelUncertainty1PtPrim",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty);
713 fHistList->Add(fPtRelUncertainty1PtPrim);
714
715 fPtRelUncertainty1PtSec = new TH2F("fPtRelUncertainty1PtSec","fPtRelUncertainty1PtSec",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty);
716 fHistList->Add(fPtRelUncertainty1PtSec);
717
718 TH1::AddDirectory(oldStatus);
719
720 PostData(0,fHistList);
721 PostData(1,fCFManagerPos->GetParticleContainer());
722 PostData(2,fCFManagerNeg->GetParticleContainer());
723
724 if(binsPt) delete [] binsPt;
725 if(binsRel1PtUncertainty) delete [] binsRel1PtUncertainty;
726
727}
728
729#endif