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