move EMCALJetTasks from PWGGA to PWGJE
[u/mrichter/AliRoot.git] / PWGJE / AliPWG4HighPtSpectra.cxx
CommitLineData
fdceab34 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//-----------------------------------------------------------------------
2b553e6f 17// Efficiency between different steps of the procedure.
fdceab34 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
67ebd013 25#ifndef ALIPWG4HIGHPTSPECTRA_CXX
26#define ALIPWG4HIGHPTSPECTRA_CXX
fdceab34 27
28#include "AliPWG4HighPtSpectra.h"
29
67ebd013 30#include "TVector3.h"
31#include <iostream>
32#include "TH1.h"
33#include "TH2.h"
34#include "TH3.h"
b4691ee7 35#include "TProfile.h"
67ebd013 36#include "TList.h"
37#include "TChain.h"
b4691ee7 38#include "TKey.h"
39#include "TSystem.h"
40#include "TFile.h"
67ebd013 41
42#include "AliAnalysisManager.h"
43#include "AliESDInputHandler.h"
44#include "AliESDtrack.h"
45#include "AliESDtrackCuts.h"
46#include "AliExternalTrackParam.h"
2b553e6f 47#include "AliCentrality.h"
65e8ecdd 48
67ebd013 49#include "AliLog.h"
50
fdceab34 51#include "AliStack.h"
52#include "TParticle.h"
fdceab34 53#include "AliMCEvent.h"
54#include "AliMCEventHandler.h"
fdceab34 55#include "AliCFContainer.h"
b4691ee7 56#include "AliGenPythiaEventHeader.h"
57#include "AliGenCocktailEventHeader.h"
67ebd013 58//#include "$ALICE_ROOT/PWG4/JetTasks/AliAnalysisHelperJetTasks.h"
59
60//#include <iostream>
fdceab34 61using namespace std; //required for resolving the 'cout' symbol
62
63ClassImp(AliPWG4HighPtSpectra)
64
65//__________________________________________________________________________
66AliPWG4HighPtSpectra::AliPWG4HighPtSpectra() : AliAnalysisTask("AliPWG4HighPtSpectra", ""),
67 fReadAODData(0),
67ebd013 68 fCFManagerPos(0x0),
69 fCFManagerNeg(0x0),
cce400da 70 fESD(0x0),
71 fMC(0x0),
72 fStack(0x0),
34fc6450 73 fVtx(0x0),
2b553e6f 74 fIsPbPb(0),
75 fCentClass(10),
34fc6450 76 fTrackType(0),
cce400da 77 fTrackCuts(0x0),
327d12da 78 fTrackCutsReject(0x0),
dae7dd67 79 fSigmaConstrainedMax(100.),
b4691ee7 80 fAvgTrials(1),
fdceab34 81 fHistList(0),
b5cc0c6d 82 fNEventAll(0),
b4691ee7 83 fNEventSel(0),
cb76764e 84 fNEventReject(0),
2b553e6f 85 fh1Centrality(0x0),
b4691ee7 86 fh1Xsec(0),
87 fh1Trials(0),
88 fh1PtHard(0),
a337a5a9 89 fh1PtHardTrials(0),
90 fPtRelUncertainty1PtPrim(0x0),
91 fPtRelUncertainty1PtSec(0x0)
fdceab34 92{
93 //
94 //Default ctor
95 //
96}
97//___________________________________________________________________________
98AliPWG4HighPtSpectra::AliPWG4HighPtSpectra(const Char_t* name) :
99 AliAnalysisTask(name,""),
100 fReadAODData(0),
67ebd013 101 fCFManagerPos(0x0),
102 fCFManagerNeg(0x0),
cce400da 103 fESD(0x0),
104 fMC(0x0),
105 fStack(0x0),
34fc6450 106 fVtx(0x0),
2b553e6f 107 fIsPbPb(0),
108 fCentClass(10),
34fc6450 109 fTrackType(0),
cce400da 110 fTrackCuts(0x0),
327d12da 111 fTrackCutsReject(0x0),
dae7dd67 112 fSigmaConstrainedMax(100.),
b4691ee7 113 fAvgTrials(1),
fdceab34 114 fHistList(0),
b5cc0c6d 115 fNEventAll(0),
b4691ee7 116 fNEventSel(0),
cb76764e 117 fNEventReject(0),
2b553e6f 118 fh1Centrality(0x0),
b4691ee7 119 fh1Xsec(0),
120 fh1Trials(0),
121 fh1PtHard(0),
a337a5a9 122 fh1PtHardTrials(0),
123 fPtRelUncertainty1PtPrim(0x0),
124 fPtRelUncertainty1PtSec(0x0)
fdceab34 125{
126 //
127 // Constructor. Initialization of Inputs and Outputs
128 //
f51451be 129 AliDebug(2,Form("AliPWG4HighPtSpectra Calling Constructor"));
fdceab34 130 // Input slot #0 works with a TChain ESD
131 DefineInput(0, TChain::Class());
e5abcde9 132 // Output slot #0 writes into a TList
fdceab34 133 DefineOutput(0,TList::Class());
e5abcde9 134 // Output slot #1, #2 writes into a AliCFContainer
fdceab34 135 DefineOutput(1,AliCFContainer::Class());
67ebd013 136 DefineOutput(2,AliCFContainer::Class());
e5abcde9 137 // Output slot #3 writes into a AliESDtrackCuts
138 DefineOutput(3, AliESDtrackCuts::Class());
1f329128 139 DefineOutput(4, AliESDtrackCuts::Class());
e5abcde9 140}
141
142//________________________________________________________________________
143void AliPWG4HighPtSpectra::LocalInit()
144{
145 //
146 // Only called once at beginning
147 //
148 PostData(3,fTrackCuts);
fdceab34 149}
150
fdceab34 151//________________________________________________________________________
152void AliPWG4HighPtSpectra::ConnectInputData(Option_t *)
153{
154 // Connect ESD here
155 // Called once
df943115 156 AliDebug(2,Form(">> AliPWG4HighPtSpectra::ConnectInputData \n"));
fdceab34 157
158 TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
159 if (!tree) {
cce400da 160 AliDebug(2,Form( "ERROR: Could not read chain from input slot 0 \n"));
161 return;
fdceab34 162 }
cce400da 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();
67ebd013 171
cce400da 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
fdceab34 179}
cce400da 180
181//________________________________________________________________________
182Bool_t AliPWG4HighPtSpectra::SelectEvent() {
fdceab34 183 //
cce400da 184 // Decide if event should be selected for analysis
fdceab34 185 //
fdceab34 186
cce400da 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;
b5cc0c6d 195
cce400da 196 //fESD object available?
fdceab34 197 if (!fESD) {
cce400da 198 AliDebug(2,Form("ERROR: fInputEvent not available\n"));
cb76764e 199 fNEventReject->Fill("noESD",1);
cce400da 200 selectEvent = kFALSE;
201 return selectEvent;
fdceab34 202 }
203
cce400da 204 //Trigger
bfd58011 205 UInt_t isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
206 if(!(isSelected&AliVEvent::kMB)) { //Select collison candidates
b5cc0c6d 207 AliDebug(2,Form(" Trigger Selection: event REJECTED ... "));
cb76764e 208 fNEventReject->Fill("Trigger",1);
cce400da 209 selectEvent = kFALSE;
210 return selectEvent;
67ebd013 211 }
b5cc0c6d 212
cce400da 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;
fdceab34 218 }
b4691ee7 219
cce400da 220 //Check if vertex is reconstructed
34fc6450 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
fdceab34 235 // Need vertex cut
34fc6450 236 // TString vtxName(fVtx->GetName());
237 if(fVtx->GetNContributors()<2) {
238 fNEventReject->Fill("NCont<2",1);
239 selectEvent = kFALSE;
240 return selectEvent;
cd9a6fa2 241 }
cce400da 242
243 //Check if z-vertex < 10 cm
b5cc0c6d 244 double primVtx[3];
34fc6450 245 fVtx->GetXYZ(primVtx);
67ebd013 246 if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
cb76764e 247 fNEventReject->Fill("ZVTX>10",1);
cce400da 248 selectEvent = kFALSE;
249 return selectEvent;
67ebd013 250 }
cce400da 251
2b553e6f 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
cce400da 273 return selectEvent;
274
275}
276
2b553e6f 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
b1041e3b 289 if(cent<0) return 5;
2b553e6f 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
cce400da 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()) {
cb76764e 310 fNEventReject->Fill("NTracks<2",1);
67ebd013 311 // Post output data
312 PostData(0,fHistList);
313 PostData(1,fCFManagerPos->GetParticleContainer());
314 PostData(2,fCFManagerNeg->GetParticleContainer());
cd9a6fa2 315 return;
316 }
cce400da 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
fdceab34 343 Int_t nTracks = fESD->GetNumberOfTracks();
344 AliDebug(2,Form("nTracks %d", nTracks));
fdceab34 345
67ebd013 346 if(!fTrackCuts) {
cb76764e 347 fNEventReject->Fill("noTrackCuts",1);
67ebd013 348 // Post output data
349 PostData(0,fHistList);
350 PostData(1,fCFManagerPos->GetParticleContainer());
351 PostData(2,fCFManagerNeg->GetParticleContainer());
352 return;
353 }
354
b5cc0c6d 355 // Selected events for analysis
356 fNEventSel->Fill(0.);
67ebd013 357
dae7dd67 358 const Int_t nvar = 4;
67ebd013 359
dae7dd67 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
ec555a3c 363
fdceab34 364 //Now go to rec level
365 for (Int_t iTrack = 0; iTrack<nTracks; iTrack++)
366 {
2b553e6f 367 //Get track for analysis
368 AliESDtrack *track = 0x0;
369 AliESDtrack *esdtrack = fESD->GetTrack(iTrack);
370 if(!esdtrack) continue;
371
8da84c2c 372 if(fTrackType==1) {
2b553e6f 373 track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
8da84c2c 374 if(!track) continue;
375 }
2b553e6f 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;
34fc6450 384 continue;
385 }
2b553e6f 386 track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
bcd60a94 387 }
d3a3f33d 388 else if(fTrackType==7) {
389 //use global constrained track
42881dab 390 track = new AliESDtrack(*esdtrack);
391 // track = esdtrack;
031aac8e 392 // track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
d3a3f33d 393
394 }
bcd60a94 395 else
2b553e6f 396 track = esdtrack;
bcd60a94 397
83ff93cd 398 if(!track) continue;
2b553e6f 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 }
fdceab34 407
1f329128 408 if (fTrackCuts->AcceptTrack(track)) {
031aac8e 409
410 if(fTrackType==7) {
327d12da 411 if(fTrackCutsReject ) {
42881dab 412 if(fTrackCutsReject->AcceptTrack(track) ) {
413 if(track) delete track;
327d12da 414 continue;
42881dab 415 }
327d12da 416 }
417
031aac8e 418 if(esdtrack->GetConstrainedParam())
419 track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
031aac8e 420 }
421
422 //fill the container
423 containerInputRec[0] = track->Pt();
424 containerInputRec[1] = track->Phi();
425 containerInputRec[2] = track->Eta();
dae7dd67 426 containerInputRec[3] = track->GetTPCNclsIter1();
031aac8e 427
1f329128 428 if(track->GetSign()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
429 if(track->GetSign()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
430
67ebd013 431
1f329128 432 //Only fill the MC containers if MC information is available
cce400da 433 if(fMC) {
0cc50ca6 434 Int_t label = TMath::Abs(track->GetLabel());
cce400da 435 TParticle *particle = fStack->Particle(label) ;
bcd60a94 436 if(!particle) {
42881dab 437 if(fTrackType==1 || fTrackType==2 || fTrackType==7)
2b553e6f 438 delete track;
bcd60a94 439 continue;
440 }
ec555a3c 441 containerInputRecMC[0] = particle->Pt();
442 containerInputRecMC[1] = particle->Phi();
443 containerInputRecMC[2] = particle->Eta();
dae7dd67 444 containerInputRecMC[3] = track->GetTPCNclsIter1();
67ebd013 445
1f329128 446 //Container with primaries
cce400da 447 if(fStack->IsPhysicalPrimary(label)) {
1f329128 448 if(particle->GetPDG()->Charge()>0.) {
ec555a3c 449 fCFManagerPos->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
1f329128 450 }
451 if(particle->GetPDG()->Charge()<0.) {
ec555a3c 452 fCFManagerNeg->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
1f329128 453 }
a337a5a9 454
455 //Fill pT resolution plots for primaries
b8f6bad5 456 fPtRelUncertainty1PtPrim->Fill(containerInputRec[0],containerInputRec[0]*TMath::Sqrt(track->GetSigma1Pt2()));
a337a5a9 457
67ebd013 458 }
459
1f329128 460 //Container with secondaries
cce400da 461 if (!fStack->IsPhysicalPrimary(label) ) {
67ebd013 462 if(particle->GetPDG()->Charge()>0.) {
aa3ba8d2 463 fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepSecondaries);
67ebd013 464 }
465 if(particle->GetPDG()->Charge()<0.) {
aa3ba8d2 466 fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepSecondaries);
67ebd013 467 }
a337a5a9 468 //Fill pT resolution plots for primaries
b8f6bad5 469 fPtRelUncertainty1PtSec->Fill(containerInputRec[0],containerInputRec[0]*TMath::Sqrt(track->GetSigma1Pt2()));
0cc50ca6 470 }
471 }
67ebd013 472
2b553e6f 473 }//trackCuts global tracks
60829498 474
42881dab 475 if(fTrackType==1 || fTrackType==2 || fTrackType==7) {
476 if(track) delete track;
477 }
1f329128 478 }//track loop
67ebd013 479
fdceab34 480
a337a5a9 481 //Fill MC containers if particles are findable
cce400da 482 if(fMC) {
483 for(int iPart = 1; iPart<(fMC->GetNumberOfPrimaries()); iPart++) {
484 AliMCParticle *mcPart = (AliMCParticle*)fMC->GetTrack(iPart);
485 if(!mcPart) continue;
4f6d2d78 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 }
0cc50ca6 504 }
cce400da 505 }
0cc50ca6 506 }
67ebd013 507
67ebd013 508 PostData(0,fHistList);
509 PostData(1,fCFManagerPos->GetParticleContainer());
510 PostData(2,fCFManagerNeg->GetParticleContainer());
511
fdceab34 512}
b4691ee7 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 }
cb76764e 534 // Printf("%s",file.Data());
b4691ee7 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){
cb76764e 599 // Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
b4691ee7 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) {
cb76764e 622 // AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
b4691ee7 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}
fdceab34 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.
9b58297c 648
fdceab34 649}
650
fdceab34 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 //
f51451be 656 AliDebug(2,Form("CreateOutputObjects CreateOutputObjects of task %s", GetName()));
df943115 657
67ebd013 658 Bool_t oldStatus = TH1::AddDirectoryStatus();
659 TH1::AddDirectory(kFALSE);
660
fdceab34 661 //slot #1
b5cc0c6d 662 OpenFile(0);
663 fHistList = new TList();
b4691ee7 664 fHistList->SetOwner(kTRUE);
b5cc0c6d 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);
fdceab34 669
cb76764e 670 fNEventReject = new TH1F("fNEventReject","Reason events are rejectected for analysis",20,0,20);
2b553e6f 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);
cb76764e 681 fHistList->Add(fNEventReject);
682
2b553e6f 683 fh1Centrality = new TH1F("fh1Centrality","fh1Centrality; Centrality %",100,0,100);
684 fHistList->Add(fh1Centrality);
685
b4691ee7 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
a337a5a9 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
67ebd013 718 TH1::AddDirectory(oldStatus);
719
11245558 720 PostData(0,fHistList);
721 PostData(1,fCFManagerPos->GetParticleContainer());
722 PostData(2,fCFManagerNeg->GetParticleContainer());
723
a337a5a9 724 if(binsPt) delete [] binsPt;
725 if(binsRel1PtUncertainty) delete [] binsRel1PtUncertainty;
726
fdceab34 727}
728
729#endif