]> 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
Merge branch 'AD'
[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 "AliGenHijingEventHeader.h"
58#include "AliGenCocktailEventHeader.h"
59//#include "$ALICE_ROOT/PWG4/JetTasks/AliAnalysisHelperJetTasks.h"
60
61//#include <iostream>
62using namespace std; //required for resolving the 'cout' symbol
63
64ClassImp(AliPWG4HighPtSpectra)
65
66//__________________________________________________________________________
67AliPWG4HighPtSpectra::AliPWG4HighPtSpectra() : AliAnalysisTask("AliPWG4HighPtSpectra", ""),
68 fReadAODData(0),
69 fCFManagerPos(0x0),
70 fCFManagerNeg(0x0),
71 fESD(0x0),
72 fMC(0x0),
73 fStack(0x0),
74 fVtx(0x0),
75 fTriggerMask(AliVEvent::kMB),
76 fIsPbPb(0),
77 fCentClass(10),
78 fTrackType(0),
79 fTrackCuts(0x0),
80 fTrackCutsReject(0x0),
81 fbSelectHIJING(kFALSE),
82 fSigmaConstrainedMax(100.),
83 fAvgTrials(1),
84 fHistList(0),
85 fNEventAll(0),
86 fNEventSel(0),
87 fNEventReject(0),
88 fh1Centrality(0x0),
89 fh1Xsec(0),
90 fh1Trials(0),
91 fh1PtHard(0),
92 fh1PtHardTrials(0),
93 fPtRelUncertainty1PtPrim(0x0),
94 fPtRelUncertainty1PtSec(0x0)
95{
96 //
97 //Default ctor
98 //
99}
100//___________________________________________________________________________
101AliPWG4HighPtSpectra::AliPWG4HighPtSpectra(const Char_t* name) :
102 AliAnalysisTask(name,""),
103 fReadAODData(0),
104 fCFManagerPos(0x0),
105 fCFManagerNeg(0x0),
106 fESD(0x0),
107 fMC(0x0),
108 fStack(0x0),
109 fVtx(0x0),
110 fTriggerMask(AliVEvent::kMB),
111 fIsPbPb(0),
112 fCentClass(10),
113 fTrackType(0),
114 fTrackCuts(0x0),
115 fTrackCutsReject(0x0),
116 fbSelectHIJING(kFALSE),
117 fSigmaConstrainedMax(100.),
118 fAvgTrials(1),
119 fHistList(0),
120 fNEventAll(0),
121 fNEventSel(0),
122 fNEventReject(0),
123 fh1Centrality(0x0),
124 fh1Xsec(0),
125 fh1Trials(0),
126 fh1PtHard(0),
127 fh1PtHardTrials(0),
128 fPtRelUncertainty1PtPrim(0x0),
129 fPtRelUncertainty1PtSec(0x0)
130{
131 //
132 // Constructor. Initialization of Inputs and Outputs
133 //
134 AliDebug(2,Form("AliPWG4HighPtSpectra Calling Constructor"));
135 // Input slot #0 works with a TChain ESD
136 DefineInput(0, TChain::Class());
137 // Output slot #0 writes into a TList
138 DefineOutput(0,TList::Class());
139 // Output slot #1, #2 writes into a AliCFContainer
140 DefineOutput(1,AliCFContainer::Class());
141 DefineOutput(2,AliCFContainer::Class());
142 // Output slot #3 writes into a AliESDtrackCuts
143 DefineOutput(3, AliESDtrackCuts::Class());
144 DefineOutput(4, AliESDtrackCuts::Class());
145}
146
147//________________________________________________________________________
148void AliPWG4HighPtSpectra::LocalInit()
149{
150 //
151 // Only called once at beginning
152 //
153 PostData(3,fTrackCuts);
154}
155
156//________________________________________________________________________
157void AliPWG4HighPtSpectra::ConnectInputData(Option_t *)
158{
159 // Connect ESD here
160 // Called once
161 AliDebug(2,Form(">> AliPWG4HighPtSpectra::ConnectInputData \n"));
162
163 TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
164 if (!tree) {
165 AliDebug(2,Form( "ERROR: Could not read chain from input slot 0 \n"));
166 return;
167 }
168
169 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
170
171 if (!esdH) {
172 AliDebug(2,Form("ERROR: Could not get ESDInputHandler"));
173 return;
174 } else
175 fESD = esdH->GetEvent();
176
177 AliMCEventHandler *eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
178 if (!eventHandler) {
179 AliDebug(2,Form( "ERROR: Could not retrieve MC event handler \n"));
180 }
181 else
182 fMC = eventHandler->MCEvent();
183
184}
185
186//________________________________________________________________________
187Bool_t AliPWG4HighPtSpectra::SelectEvent() {
188 //
189 // Decide if event should be selected for analysis
190 //
191
192 // Checks following requirements:
193 // - fESD available
194 // - trigger info from AliPhysicsSelection
195 // - number of reconstructed tracks > 1
196 // - primary vertex reconstructed
197 // - z-vertex < 10 cm
198
199 Bool_t selectEvent = kTRUE;
200
201 //fESD object available?
202 if (!fESD) {
203 AliDebug(2,Form("ERROR: fInputEvent not available\n"));
204 fNEventReject->Fill("noESD",1);
205 selectEvent = kFALSE;
206 return selectEvent;
207 }
208
209 //Trigger
210 UInt_t isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
211 if(fTriggerMask != AliVEvent::kAny && !(isSelected&fTriggerMask)) { //Select collison candidates
212 AliDebug(2,Form(" Trigger Selection: event REJECTED ... "));
213 fNEventReject->Fill("Trigger",1);
214 selectEvent = kFALSE;
215 return selectEvent;
216 }
217
218 //Check if number of reconstructed tracks is larger than 1
219 if(!fESD->GetNumberOfTracks() || fESD->GetNumberOfTracks()<2) {
220 fNEventReject->Fill("NTracks<2",1);
221 selectEvent = kFALSE;
222 return selectEvent;
223 }
224
225 //Check if vertex is reconstructed
226 fVtx = fESD->GetPrimaryVertexSPD();
227
228 if(!fVtx) {
229 fNEventReject->Fill("noVTX",1);
230 selectEvent = kFALSE;
231 return selectEvent;
232 }
233
234 if(!fVtx->GetStatus()) {
235 fNEventReject->Fill("VtxStatus",1);
236 selectEvent = kFALSE;
237 return selectEvent;
238 }
239
240 // Need vertex cut
241 // TString vtxName(fVtx->GetName());
242 if(fVtx->GetNContributors()<2) {
243 fNEventReject->Fill("NCont<2",1);
244 selectEvent = kFALSE;
245 return selectEvent;
246 }
247
248 //Check if z-vertex < 10 cm
249 double primVtx[3];
250 fVtx->GetXYZ(primVtx);
251 if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
252 fNEventReject->Fill("ZVTX>10",1);
253 selectEvent = kFALSE;
254 return selectEvent;
255 }
256
257 //Centrality selection should only be done in case of PbPb
258 if(IsPbPb()) {
259 Float_t cent = 0.;
260 if(fCentClass!=CalculateCentrality(fESD) && fCentClass!=10) {
261 fNEventReject->Fill("cent",1);
262 selectEvent = kFALSE;
263 return selectEvent;
264 }
265 else {
266 if(dynamic_cast<AliESDEvent*>(fESD)->GetCentrality()) {
267 cent = dynamic_cast<AliESDEvent*>(fESD)->GetCentrality()->GetCentralityPercentile("V0M");
268 }
269 if(cent>90.) {
270 fNEventReject->Fill("cent>90",1);
271 selectEvent = kFALSE;
272 return selectEvent;
273 }
274 fh1Centrality->Fill(cent);
275 }
276 }
277
278 return selectEvent;
279
280}
281
282//________________________________________________________________________
283Int_t AliPWG4HighPtSpectra::CalculateCentrality(AliESDEvent *esd){
284
285
286 Float_t cent = 999;
287
288 if(esd){
289 if(esd->GetCentrality()){
290 cent = esd->GetCentrality()->GetCentralityPercentile("V0M");
291 }
292 }
293
294 if(cent<0) return 5;
295 if(cent>80)return 4;
296 if(cent>50)return 3;
297 if(cent>30)return 2;
298 if(cent>10)return 1;
299 return 0;
300
301}
302
303//_________________________________________________
304void AliPWG4HighPtSpectra::Exec(Option_t *)
305{
306 //
307 // Main loop function
308 //
309 AliDebug(2,Form(">> AliPWG4HighPtSpectra::Exec \n"));
310
311 if(!fMC) {
312 AliMCEventHandler *eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
313 if (!eventHandler) {
314 AliDebug(2,Form( "ERROR: Could not retrieve MC event handler \n"));
315 }
316 else
317 fMC = eventHandler->MCEvent();
318 }
319
320 // All events without selection
321 fNEventAll->Fill(0.);
322
323 if(!SelectEvent()) {
324 // Post output data
325 PostData(0,fHistList);
326 PostData(1,fCFManagerPos->GetParticleContainer());
327 PostData(2,fCFManagerNeg->GetParticleContainer());
328 return;
329 }
330
331 //MCEvent available?
332 //if yes: get stack
333 if(fMC) {
334 AliDebug(2,Form("MC particles: %d", fMC->GetNumberOfTracks()));
335 fStack = fMC->Stack(); //Particles Stack
336 AliDebug(2,Form("MC particles stack: %d", fStack->GetNtrack()));
337 }
338
339 Int_t nTracks = fESD->GetNumberOfTracks();
340 AliDebug(2,Form("nTracks %d", nTracks));
341
342 if(!fTrackCuts) {
343 fNEventReject->Fill("noTrackCuts",1);
344 // Post output data
345 PostData(0,fHistList);
346 PostData(1,fCFManagerPos->GetParticleContainer());
347 PostData(2,fCFManagerNeg->GetParticleContainer());
348 return;
349 }
350
351 // Selected events for analysis
352 fNEventSel->Fill(0.);
353
354 const Int_t nvar = 4;
355
356 Double_t containerInputRec[nvar] = {0.,0.,0.,0.};
357 Double_t containerInputMC[nvar] = {0.,0.,0.,0.};
358 Double_t containerInputRecMC[nvar] = {0.,0.,0.,0.}; //reconstructed yield as function of MC variable
359
360 //Now go to rec level
361 for (Int_t iTrack = 0; iTrack<nTracks; iTrack++)
362 {
363 //Get track for analysis
364 AliESDtrack *track = 0x0;
365 AliESDtrack *esdtrack = fESD->GetTrack(iTrack);
366 if(!esdtrack)
367 continue;
368
369
370 if(fTrackType==4) {
371 if (!(fTrackCuts->AcceptTrack(esdtrack)))
372 continue;
373 }
374
375 if(fTrackType==1)
376 track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
377 else if(fTrackType==2 || fTrackType==4) {
378 track = AliESDtrackCuts::GetTPCOnlyTrack(const_cast<AliESDEvent*>(fESD),esdtrack->GetID());
379 if(!track)
380 continue;
381
382 AliExternalTrackParam exParam;
383 Bool_t relate = track->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
384 if( !relate ) {
385 if(track) delete track;
386 continue;
387 }
388 track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
389 }
390 else if(fTrackType==5 || fTrackType==6) {
391 if(fTrackCuts->AcceptTrack(esdtrack)) {
392 continue;
393 }
394 else {
395 if( !(fTrackCutsReject->AcceptTrack(esdtrack)) && fTrackCuts->AcceptTrack(esdtrack) ) {
396
397 if(fTrackType==5) {
398 //use TPConly constrained track
399 track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
400 if(!track)
401 continue;
402
403 AliExternalTrackParam exParam;
404 Bool_t relate = track->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
405 if( !relate ) {
406 if(track) delete track;
407 continue;
408 }
409 track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
410 }
411 else if(fTrackType==6) {
412 //use global constrained track
413 track = new AliESDtrack(*esdtrack);
414 track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
415
416 }
417 }
418 }
419 }
420 else if(fTrackType==7) {
421 //use global constrained track
422 track = new AliESDtrack(*esdtrack);
423 }
424 else
425 track = esdtrack;
426
427 if(!track)
428 continue;
429
430
431 if(fTrackType==2 || fTrackType==4 || fTrackType==5) {
432 //Cut on chi2 of constrained fit
433 if(track->GetConstrainedChi2TPC() > fSigmaConstrainedMax*fSigmaConstrainedMax && fSigmaConstrainedMax>0.) {
434 if(track) delete track;
435 continue;
436 }
437 }
438 if (!(fTrackCuts->AcceptTrack(track)) && fTrackType!=4 && fTrackType!=5 && fTrackType!=6) {
439 if(fTrackType==1 || fTrackType==2 || fTrackType==7) {
440 if(track) delete track;
441 }
442 continue;
443 }
444
445 if(fTrackType==7) {
446 if(fTrackCutsReject ) {
447 if(fTrackCutsReject->AcceptTrack(track) ) {
448 if(track) delete track;
449 continue;
450 }
451 }
452
453 if(esdtrack->GetConstrainedParam())
454 track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
455 }
456
457 if(!track) {
458 if(fTrackType==1 || fTrackType==2 || fTrackType==4 || fTrackType==5 || fTrackType==6 || fTrackType==7) {
459 if(track) delete track;
460 }
461 continue;
462 }
463
464 //fill the container
465 containerInputRec[0] = track->Pt();
466 containerInputRec[1] = track->Phi();
467 containerInputRec[2] = track->Eta();
468 containerInputRec[3] = track->GetTPCNclsIter1();
469
470 if(track->GetSign()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
471 if(track->GetSign()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
472
473
474 //Only fill the MC containers if MC information is available
475 if(fMC) {
476 Int_t label = TMath::Abs(track->GetLabel());
477 if(label>fStack->GetNtrack()) {
478 if(fTrackType==1 || fTrackType==2 || fTrackType==7)
479 delete track;
480 continue;
481 }
482 //Only select particles generated by HIJING if requested
483 if(fbSelectHIJING) {
484 if(!IsHIJINGParticle(label)) {
485 if(fTrackType==1 || fTrackType==2 || fTrackType==7)
486 delete track;
487 continue;
488 }
489 }
490 TParticle *particle = fStack->Particle(label) ;
491 if(!particle) {
492 if(fTrackType==1 || fTrackType==2 || fTrackType==7)
493 delete track;
494 continue;
495 }
496 containerInputRecMC[0] = particle->Pt();
497 containerInputRecMC[1] = particle->Phi();
498 containerInputRecMC[2] = particle->Eta();
499 containerInputRecMC[3] = track->GetTPCNclsIter1();
500
501 //Container with primaries
502 if(fStack->IsPhysicalPrimary(label)) {
503 if(particle->GetPDG()->Charge()>0.) {
504 fCFManagerPos->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
505 }
506 if(particle->GetPDG()->Charge()<0.) {
507 fCFManagerNeg->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
508 }
509
510 //Fill pT resolution plots for primaries
511 fPtRelUncertainty1PtPrim->Fill(containerInputRec[0],containerInputRec[0]*TMath::Sqrt(track->GetSigma1Pt2()));
512
513 }
514
515 //Container with secondaries
516 if (!fStack->IsPhysicalPrimary(label) ) {
517 if(particle->GetPDG()->Charge()>0.) {
518 fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepSecondaries);
519 }
520 if(particle->GetPDG()->Charge()<0.) {
521 fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepSecondaries);
522 }
523 //Fill pT resolution plots for primaries
524 fPtRelUncertainty1PtSec->Fill(containerInputRec[0],containerInputRec[0]*TMath::Sqrt(track->GetSigma1Pt2()));
525 }
526 }
527
528 if(fTrackType==1 || fTrackType==2 || fTrackType==4 || fTrackType==5 || fTrackType==6 || fTrackType==7) {
529 if(track) delete track;
530 }
531
532
533 }//track loop
534
535
536 //Fill MC containers if particles are findable
537 if(fMC) {
538 for(int iPart = 1; iPart<(fMC->GetNumberOfPrimaries()); iPart++) {
539 AliMCParticle *mcPart = (AliMCParticle*)fMC->GetTrack(iPart);
540 if(!mcPart) continue;
541
542 //Only select particles generated by HIJING if requested
543 if(fbSelectHIJING) {
544 if(!IsHIJINGParticle(iPart))
545 continue;
546 }
547
548 Int_t pdg = mcPart->PdgCode();
549
550 // select charged pions, protons, kaons , electrons, muons
551 if(TMath::Abs(pdg) == 211 || TMath::Abs(pdg) == 2212 || TMath::Abs(pdg) ==
552 321 || TMath::Abs(pdg) == 11 || TMath::Abs(pdg) == 13){
553
554 //fill the container
555 containerInputMC[0] = mcPart->Pt();
556 containerInputMC[1] = mcPart->Phi();
557 containerInputMC[2] = mcPart->Eta();
558 // AliESDtrack *esdtrack = fESD->GetTrack(mcPart->GetLabel());
559 containerInputMC[3] = 159.;
560
561 if(fStack->IsPhysicalPrimary(iPart)) {
562 if(mcPart->Charge()>0. && fCFManagerPos->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
563 if(mcPart->Charge()<0. && fCFManagerNeg->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
564 }
565 }
566 }
567 }
568
569 PostData(0,fHistList);
570 PostData(1,fCFManagerPos->GetParticleContainer());
571 PostData(2,fCFManagerNeg->GetParticleContainer());
572
573}
574//________________________________________________________________________
575Bool_t AliPWG4HighPtSpectra::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
576 //
577 // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
578 // This is to called in Notify and should provide the path to the AOD/ESD file
579 // Copied from AliAnalysisTaskJetSpectrum2
580 //
581
582 TString file(currFile);
583 fXsec = 0;
584 fTrials = 1;
585
586 if(file.Contains("root_archive.zip#")){
587 Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
588 Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
589 file.Replace(pos+1,20,"");
590 }
591 else {
592 // not an archive take the basename....
593 file.ReplaceAll(gSystem->BaseName(file.Data()),"");
594 }
595 // Printf("%s",file.Data());
596
597 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
598 if(!fxsec){
599 // next trial fetch the histgram file
600 fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
601 if(!fxsec){
602 // not a severe condition but indicates that we have no information
603 return kFALSE;
604 }
605 else{
606 // find the tlist we want to be independtent of the name so use the Tkey
607 TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
608 if(!key){
609 fxsec->Close();
610 return kFALSE;
611 }
612 TList *list = dynamic_cast<TList*>(key->ReadObj());
613 if(!list){
614 fxsec->Close();
615 return kFALSE;
616 }
617 fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
618 fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
619 fxsec->Close();
620 }
621 } // no tree pyxsec.root
622 else {
623 TTree *xtree = (TTree*)fxsec->Get("Xsection");
624 if(!xtree){
625 fxsec->Close();
626 return kFALSE;
627 }
628 UInt_t ntrials = 0;
629 Double_t xsection = 0;
630 xtree->SetBranchAddress("xsection",&xsection);
631 xtree->SetBranchAddress("ntrials",&ntrials);
632 xtree->GetEntry(0);
633 fTrials = ntrials;
634 fXsec = xsection;
635 fxsec->Close();
636 }
637 return kTRUE;
638}
639//________________________________________________________________________
640Bool_t AliPWG4HighPtSpectra::Notify()
641{
642 //
643 // Implemented Notify() to read the cross sections
644 // and number of trials from pyxsec.root
645 // Copied from AliAnalysisTaskJetSpectrum2
646 //
647
648 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
649 Float_t xsection = 0;
650 Float_t ftrials = 1;
651
652 if(tree){
653 TFile *curfile = tree->GetCurrentFile();
654 if (!curfile) {
655 Error("Notify","No current file");
656 return kFALSE;
657 }
658 if(!fh1Xsec||!fh1Trials){
659 // Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
660 return kFALSE;
661 }
662 PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
663 fh1Xsec->Fill("<#sigma>",xsection);
664 fh1Trials->Fill("#sum{ntrials}",ftrials);
665 }
666 return kTRUE;
667}
668
669//________________________________________________________________________
670AliGenPythiaEventHeader* AliPWG4HighPtSpectra::GetPythiaEventHeader(AliMCEvent *mcEvent){
671
672 if(!mcEvent)return 0;
673 AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
674 AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
675 if(!pythiaGenHeader){
676 // cocktail ??
677 AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
678
679 if (!genCocktailHeader) {
680 // AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
681 // AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
682 return 0;
683 }
684 TList* headerList = genCocktailHeader->GetHeaders();
685 for (Int_t i=0; i<headerList->GetEntries(); i++) {
686 pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
687 if (pythiaGenHeader)
688 break;
689 }
690 if(!pythiaGenHeader){
691 AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
692 return 0;
693 }
694 }
695 return pythiaGenHeader;
696
697}
698
699//________________________________________________________________________
700AliGenHijingEventHeader* AliPWG4HighPtSpectra::GetHijingEventHeader(AliMCEvent *mcEvent){
701
702 if(!mcEvent)return 0;
703 AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
704 AliGenHijingEventHeader* hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(genHeader);
705 if(!hijingGenHeader){
706 // cocktail ??
707 AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
708
709 if (!genCocktailHeader) {
710 // AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
711 // AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
712 return 0;
713 }
714 TList* headerList = genCocktailHeader->GetHeaders();
715 for (Int_t i=0; i<headerList->GetEntries(); i++) {
716 hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(headerList->At(i));
717 if (hijingGenHeader)
718 break;
719 }
720 if(!hijingGenHeader){
721 AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Hijing event header not found");
722 return 0;
723 }
724 }
725 return hijingGenHeader;
726
727}
728
729
730//___________________________________________________________________________
731void AliPWG4HighPtSpectra::Terminate(Option_t*)
732{
733 // The Terminate() function is the last function to be called during
734 // a query. It always runs on the client, it can be used to present
735 // the results graphically or save the results to file.
736
737}
738
739//___________________________________________________________________________
740void AliPWG4HighPtSpectra::CreateOutputObjects() {
741 //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
742 //TO BE SET BEFORE THE EXECUTION OF THE TASK
743 //
744 AliDebug(2,Form("CreateOutputObjects CreateOutputObjects of task %s", GetName()));
745
746 Bool_t oldStatus = TH1::AddDirectoryStatus();
747 TH1::AddDirectory(kFALSE);
748
749 //slot #1
750 OpenFile(0);
751 fHistList = new TList();
752 fHistList->SetOwner(kTRUE);
753 fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
754 fHistList->Add(fNEventAll);
755 fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
756 fHistList->Add(fNEventSel);
757
758 fNEventReject = new TH1F("fNEventReject","Reason events are rejectected for analysis",20,0,20);
759 //Set labels
760 fNEventReject->Fill("noESD",0);
761 fNEventReject->Fill("Trigger",0);
762 fNEventReject->Fill("NTracks<2",0);
763 fNEventReject->Fill("noVTX",0);
764 fNEventReject->Fill("VtxStatus",0);
765 fNEventReject->Fill("NCont<2",0);
766 fNEventReject->Fill("ZVTX>10",0);
767 fNEventReject->Fill("cent",0);
768 fNEventReject->Fill("cent>90",0);
769 fHistList->Add(fNEventReject);
770
771 fh1Centrality = new TH1F("fh1Centrality","fh1Centrality; Centrality %",100,0,100);
772 fHistList->Add(fh1Centrality);
773
774 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
775 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
776 fHistList->Add(fh1Xsec);
777
778 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
779 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
780 fHistList->Add(fh1Trials);
781
782 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
783 fHistList->Add(fh1PtHard);
784 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
785 fHistList->Add(fh1PtHardTrials);
786
787 Int_t fgkNPtBins = 100;
788 Float_t kMinPt = 0.;
789 Float_t kMaxPt = 100.;
790 Double_t *binsPt = new Double_t[fgkNPtBins+1];
791 for(Int_t i=0; i<=fgkNPtBins; i++) binsPt[i]=(Double_t)kMinPt + (kMaxPt-kMinPt)/fgkNPtBins*(Double_t)i ;
792
793 Int_t fgkNRel1PtUncertaintyBins=50;
794 Float_t fgkRel1PtUncertaintyMin = 0.;
795 Float_t fgkRel1PtUncertaintyMax = 1.;
796
797 Double_t *binsRel1PtUncertainty=new Double_t[fgkNRel1PtUncertaintyBins+1];
798 for(Int_t i=0; i<=fgkNRel1PtUncertaintyBins; i++) binsRel1PtUncertainty[i]=(Double_t)fgkRel1PtUncertaintyMin + (fgkRel1PtUncertaintyMax-fgkRel1PtUncertaintyMin)/fgkNRel1PtUncertaintyBins*(Double_t)i ;
799
800 fPtRelUncertainty1PtPrim = new TH2F("fPtRelUncertainty1PtPrim","fPtRelUncertainty1PtPrim",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty);
801 fHistList->Add(fPtRelUncertainty1PtPrim);
802
803 fPtRelUncertainty1PtSec = new TH2F("fPtRelUncertainty1PtSec","fPtRelUncertainty1PtSec",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty);
804 fHistList->Add(fPtRelUncertainty1PtSec);
805
806 TH1::AddDirectory(oldStatus);
807
808 PostData(0,fHistList);
809 PostData(1,fCFManagerPos->GetParticleContainer());
810 PostData(2,fCFManagerNeg->GetParticleContainer());
811
812 if(binsPt) delete [] binsPt;
813 if(binsRel1PtUncertainty) delete [] binsRel1PtUncertainty;
814
815}
816
817//________________________________________________________________________
818Bool_t AliPWG4HighPtSpectra::IsHIJINGParticle(Int_t label) {
819 //
820 // Return kTRUE in case particle is from HIJING event
821 //
822
823 AliGenHijingEventHeader* hijingHeader = GetHijingEventHeader(fMC);
824
825 Int_t nproduced = hijingHeader->NProduced();
826
827 TParticle * mom = fStack->Particle(label);
828 Int_t iMom = label;
829 Int_t iParent = mom->GetFirstMother();
830 while(iParent!=-1){
831 iMom = iParent;
832 mom = fStack->Particle(iMom);
833 iParent = mom->GetFirstMother();
834 }
835
836 if(iMom<nproduced) {
837 return kTRUE;
838 } else {
839 return kFALSE;
840 }
841
842}
843
844#endif