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