]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/AliPWG4HighPtSpectra.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
don't sort clusters after local reco, do this in AliITSUTrackerGlo
[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
f15c1f69 422 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(iTrack));
423 if(!aodtrack) AliFatal("Not a standard AOD");
9477b466 424 if(!aodtrack)
425 continue;
5ef6e458 426 if((!aodtrack->TestFilterBit(fFilterMask)) ||
427 ((!fCFManagerPos->CheckParticleCuts(kStepReconstructed,aodtrack)) && (aodtrack->Charge()>0.)) ||
428 ((!fCFManagerNeg->CheckParticleCuts(kStepReconstructed,aodtrack)) && (aodtrack->Charge()<0)) )
9477b466 429 continue;
430 else {
3abf7501 431 track = static_cast<AliVTrack*>(aodtrack);
9477b466 432 }
433 //fill the container
434 containerInputRec[0] = track->Pt();
435 containerInputRec[1] = track->Phi();
436 containerInputRec[2] = track->Eta();
3abf7501 437 containerInputRec[3] = track->GetTPCNcls();
9477b466 438
6e349fdb 439 if(track->Charge()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
440 if(track->Charge()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
9477b466 441
3abf7501 442 if(fArrayMCAOD) {
9477b466 443 Int_t label = TMath::Abs(track->GetLabel());
3abf7501 444 if(label>fArrayMCAOD->GetEntries()) {
9477b466 445 continue;
446 }
3abf7501 447 AliAODMCParticle *particle = (AliAODMCParticle*) fArrayMCAOD->At(label);
9477b466 448 if(!particle) {
9477b466 449 continue;
450 }
451 //Only select particles generated by HIJING if requested
452 if(fbSelectHIJING) {
453 if(!IsHIJINGParticle(label)) {
9477b466 454 continue;
455 }
456 }
457 containerInputRecMC[0] = particle->Pt();
458 containerInputRecMC[1] = particle->Phi();
459 containerInputRecMC[2] = particle->Eta();
460 containerInputRecMC[3] = track->GetTPCNcls();
461
462 //Container with primaries
463 if(particle->IsPhysicalPrimary()) {
5ef6e458 464 if(particle->Charge()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
465 if(particle->Charge()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
9477b466 466 //Fill pT resolution plots for primaries
3abf7501 467 //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 468 }
469
470 //Container with secondaries
471 if (!particle->IsPhysicalPrimary() ) {
6e349fdb 472 if(particle->Charge()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputRecMC,kStepSecondaries);
473 if(particle->Charge()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputRecMC,kStepSecondaries);
9477b466 474 //Fill pT resolution plots for primaries
3abf7501 475 //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 476 }
477 }
478 }//track loop
479
480 //Fill MC containers if particles are findable
3abf7501 481 if(fArrayMCAOD) {
482 int noPart = fArrayMCAOD->GetEntriesFast();
6e349fdb 483
3abf7501 484 for(int iPart = 1; iPart<noPart; iPart++) {
485 AliAODMCParticle *mcPart = (AliAODMCParticle*) fArrayMCAOD->At(iPart);
9477b466 486 if(!mcPart) continue;
487
488 //Only select particles generated by HIJING if requested
489 if(fbSelectHIJING) {
490 if(!IsHIJINGParticle(iPart))
491 continue;
492 }
493
494 Int_t pdg = mcPart->PdgCode();
495
496 // select charged pions, protons, kaons , electrons, muons
497 if(TMath::Abs(pdg) == 211 || TMath::Abs(pdg) == 2212 || TMath::Abs(pdg) ==
498 321 || TMath::Abs(pdg) == 11 || TMath::Abs(pdg) == 13){
499
500 //fill the container
501 containerInputMC[0] = mcPart->Pt();
502 containerInputMC[1] = mcPart->Phi();
503 containerInputMC[2] = mcPart->Eta();
504 // AliESDtrack *esdtrack = fESD->GetTrack(mcPart->GetLabel());
505 containerInputMC[3] = 159.;
506
507 if(mcPart->IsPhysicalPrimary()) {
6e349fdb 508 if(mcPart->Charge()>0. && fCFManagerPos->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
509 if(mcPart->Charge()<0. && fCFManagerNeg->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
9477b466 510 }
511 }
512 }
513 }
514
515 }//end of AOD analysis
516 else {//ESD analysis
517 for (Int_t iTrack = 0; iTrack<nTracks; iTrack++)
518 {
2b553e6f 519 //Get track for analysis
520 AliESDtrack *track = 0x0;
521 AliESDtrack *esdtrack = fESD->GetTrack(iTrack);
2461511f 522 if(!esdtrack)
523 continue;
524
2b553e6f 525
2461511f 526 if(fTrackType==4) {
527 if (!(fTrackCuts->AcceptTrack(esdtrack)))
528 continue;
8da84c2c 529 }
2b553e6f 530
2461511f 531 if(fTrackType==1)
532 track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
533 else if(fTrackType==2 || fTrackType==4) {
534 track = AliESDtrackCuts::GetTPCOnlyTrack(const_cast<AliESDEvent*>(fESD),esdtrack->GetID());
535 if(!track)
536 continue;
537
2b553e6f 538 AliExternalTrackParam exParam;
9477b466 539 Bool_t relate = track->RelateToVertexTPC(fESD->GetPrimaryVertexSPD(),fESD->GetMagneticField(),kVeryBig,&exParam);
2b553e6f 540 if( !relate ) {
2461511f 541 if(track) delete track;
34fc6450 542 continue;
543 }
2b553e6f 544 track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
bcd60a94 545 }
2461511f 546 else if(fTrackType==5 || fTrackType==6) {
547 if(fTrackCuts->AcceptTrack(esdtrack)) {
548 continue;
549 }
550 else {
551 if( !(fTrackCutsReject->AcceptTrack(esdtrack)) && fTrackCuts->AcceptTrack(esdtrack) ) {
552
553 if(fTrackType==5) {
554 //use TPConly constrained track
555 track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
556 if(!track)
557 continue;
558
559 AliExternalTrackParam exParam;
9477b466 560 Bool_t relate = track->RelateToVertexTPC(fESD->GetPrimaryVertexSPD(),fESD->GetMagneticField(),kVeryBig,&exParam);
2461511f 561 if( !relate ) {
562 if(track) delete track;
563 continue;
564 }
565 track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
566 }
567 else if(fTrackType==6) {
568 //use global constrained track
569 track = new AliESDtrack(*esdtrack);
570 track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
2461511f 571 }
572 }
573 }
574 }
d3a3f33d 575 else if(fTrackType==7) {
576 //use global constrained track
42881dab 577 track = new AliESDtrack(*esdtrack);
d3a3f33d 578 }
bcd60a94 579 else
2b553e6f 580 track = esdtrack;
bcd60a94 581
2461511f 582 if(!track)
583 continue;
584
585
586 if(fTrackType==2 || fTrackType==4 || fTrackType==5) {
2b553e6f 587 //Cut on chi2 of constrained fit
2461511f 588 if(track->GetConstrainedChi2TPC() > fSigmaConstrainedMax*fSigmaConstrainedMax && fSigmaConstrainedMax>0.) {
589 if(track) delete track;
2b553e6f 590 continue;
591 }
592 }
2461511f 593 if (!(fTrackCuts->AcceptTrack(track)) && fTrackType!=4 && fTrackType!=5 && fTrackType!=6) {
594 if(fTrackType==1 || fTrackType==2 || fTrackType==7) {
595 if(track) delete track;
596 }
597 continue;
598 }
fdceab34 599
2461511f 600 if(fTrackType==7) {
601 if(fTrackCutsReject ) {
602 if(fTrackCutsReject->AcceptTrack(track) ) {
603 if(track) delete track;
604 continue;
327d12da 605 }
031aac8e 606 }
2461511f 607 if(esdtrack->GetConstrainedParam())
608 track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
609 }
031aac8e 610
2461511f 611 if(!track) {
612 if(fTrackType==1 || fTrackType==2 || fTrackType==4 || fTrackType==5 || fTrackType==6 || fTrackType==7) {
613 if(track) delete track;
614 }
615 continue;
616 }
617
618 //fill the container
619 containerInputRec[0] = track->Pt();
620 containerInputRec[1] = track->Phi();
621 containerInputRec[2] = track->Eta();
622 containerInputRec[3] = track->GetTPCNclsIter1();
031aac8e 623
2461511f 624 if(track->GetSign()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
625 if(track->GetSign()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
1f329128 626
67ebd013 627
2461511f 628 //Only fill the MC containers if MC information is available
629 if(fMC) {
630 Int_t label = TMath::Abs(track->GetLabel());
631 if(label>fStack->GetNtrack()) {
632 if(fTrackType==1 || fTrackType==2 || fTrackType==7)
633 delete track;
634 continue;
635 }
636 //Only select particles generated by HIJING if requested
637 if(fbSelectHIJING) {
638 if(!IsHIJINGParticle(label)) {
5d87a047 639 if(fTrackType==1 || fTrackType==2 || fTrackType==7)
640 delete track;
641 continue;
642 }
2461511f 643 }
644 TParticle *particle = fStack->Particle(label) ;
645 if(!particle) {
646 if(fTrackType==1 || fTrackType==2 || fTrackType==7)
647 delete track;
648 continue;
649 }
650 containerInputRecMC[0] = particle->Pt();
651 containerInputRecMC[1] = particle->Phi();
652 containerInputRecMC[2] = particle->Eta();
653 containerInputRecMC[3] = track->GetTPCNclsIter1();
654
655 //Container with primaries
656 if(fStack->IsPhysicalPrimary(label)) {
6e349fdb 657 if(particle->GetPDG()->Charge()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
658 if(particle->GetPDG()->Charge()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
2461511f 659 //Fill pT resolution plots for primaries
660 fPtRelUncertainty1PtPrim->Fill(containerInputRec[0],containerInputRec[0]*TMath::Sqrt(track->GetSigma1Pt2()));
2461511f 661 }
67ebd013 662
2461511f 663 //Container with secondaries
664 if (!fStack->IsPhysicalPrimary(label) ) {
6e349fdb 665 if(particle->GetPDG()->Charge()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepSecondaries);
666 if(particle->GetPDG()->Charge()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepSecondaries);
2461511f 667 //Fill pT resolution plots for primaries
668 fPtRelUncertainty1PtSec->Fill(containerInputRec[0],containerInputRec[0]*TMath::Sqrt(track->GetSigma1Pt2()));
0cc50ca6 669 }
2461511f 670 }
60829498 671
2461511f 672 if(fTrackType==1 || fTrackType==2 || fTrackType==4 || fTrackType==5 || fTrackType==6 || fTrackType==7) {
42881dab 673 if(track) delete track;
674 }
1f329128 675 }//track loop
03372fd1 676
9477b466 677 //Fill MC containers if particles are findable
678 if(fMC) {
679 for(int iPart = 1; iPart<(fMC->GetNumberOfPrimaries()); iPart++) {
680 AliMCParticle *mcPart = (AliMCParticle*)fMC->GetTrack(iPart);
681 if(!mcPart) continue;
682
683 //Only select particles generated by HIJING if requested
684 if(fbSelectHIJING) {
685 if(!IsHIJINGParticle(iPart))
686 continue;
687 }
688
689 Int_t pdg = mcPart->PdgCode();
690
691 // select charged pions, protons, kaons , electrons, muons
692 if(TMath::Abs(pdg) == 211 || TMath::Abs(pdg) == 2212 || TMath::Abs(pdg) ==
693 321 || TMath::Abs(pdg) == 11 || TMath::Abs(pdg) == 13){
694
695 //fill the container
696 containerInputMC[0] = mcPart->Pt();
697 containerInputMC[1] = mcPart->Phi();
698 containerInputMC[2] = mcPart->Eta();
699 // AliESDtrack *esdtrack = fESD->GetTrack(mcPart->GetLabel());
700 containerInputMC[3] = 159.;
701
702 if(fStack->IsPhysicalPrimary(iPart)) {
703 if(mcPart->Charge()>0. && fCFManagerPos->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
704 if(mcPart->Charge()<0. && fCFManagerNeg->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
705 }
706 }
0cc50ca6 707 }
cce400da 708 }
9477b466 709 }//end of ESD analysis
710
67ebd013 711 PostData(0,fHistList);
712 PostData(1,fCFManagerPos->GetParticleContainer());
713 PostData(2,fCFManagerNeg->GetParticleContainer());
714
fdceab34 715}
9477b466 716
b4691ee7 717//________________________________________________________________________
9477b466 718Bool_t AliPWG4HighPtSpectra::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials)
719{
b4691ee7 720 //
721 // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
722 // This is to called in Notify and should provide the path to the AOD/ESD file
723 // Copied from AliAnalysisTaskJetSpectrum2
724 //
725
726 TString file(currFile);
727 fXsec = 0;
728 fTrials = 1;
729
730 if(file.Contains("root_archive.zip#")){
731 Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
732 Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
733 file.Replace(pos+1,20,"");
734 }
735 else {
736 // not an archive take the basename....
737 file.ReplaceAll(gSystem->BaseName(file.Data()),"");
738 }
cb76764e 739 // Printf("%s",file.Data());
b4691ee7 740
741 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
742 if(!fxsec){
743 // next trial fetch the histgram file
744 fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
745 if(!fxsec){
2461511f 746 // not a severe condition but indicates that we have no information
b4691ee7 747 return kFALSE;
748 }
749 else{
750 // find the tlist we want to be independtent of the name so use the Tkey
751 TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
752 if(!key){
753 fxsec->Close();
754 return kFALSE;
755 }
756 TList *list = dynamic_cast<TList*>(key->ReadObj());
757 if(!list){
758 fxsec->Close();
759 return kFALSE;
760 }
761 fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
762 fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
763 fxsec->Close();
764 }
765 } // no tree pyxsec.root
766 else {
767 TTree *xtree = (TTree*)fxsec->Get("Xsection");
768 if(!xtree){
769 fxsec->Close();
770 return kFALSE;
771 }
772 UInt_t ntrials = 0;
773 Double_t xsection = 0;
774 xtree->SetBranchAddress("xsection",&xsection);
775 xtree->SetBranchAddress("ntrials",&ntrials);
776 xtree->GetEntry(0);
777 fTrials = ntrials;
778 fXsec = xsection;
779 fxsec->Close();
780 }
781 return kTRUE;
782}
9477b466 783
b4691ee7 784//________________________________________________________________________
785Bool_t AliPWG4HighPtSpectra::Notify()
786{
787 //
788 // Implemented Notify() to read the cross sections
789 // and number of trials from pyxsec.root
790 // Copied from AliAnalysisTaskJetSpectrum2
791 //
792
9477b466 793 if(fNoPythiaInfo)
794 return kTRUE;
795
b4691ee7 796 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
797 Float_t xsection = 0;
798 Float_t ftrials = 1;
799
b4691ee7 800 if(tree){
801 TFile *curfile = tree->GetCurrentFile();
802 if (!curfile) {
803 Error("Notify","No current file");
804 return kFALSE;
805 }
806 if(!fh1Xsec||!fh1Trials){
cb76764e 807 // Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
b4691ee7 808 return kFALSE;
809 }
adf3920d 810 PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
b4691ee7 811 fh1Xsec->Fill("<#sigma>",xsection);
adf3920d 812 fh1Trials->Fill("#sum{ntrials}",ftrials);
813 }
b4691ee7 814 return kTRUE;
815}
816
817//________________________________________________________________________
9477b466 818AliGenPythiaEventHeader* AliPWG4HighPtSpectra::GetPythiaEventHeader(AliMCEvent *mcEvent)
819{
b4691ee7 820
821 if(!mcEvent)return 0;
822 AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
823 AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
824 if(!pythiaGenHeader){
825 // cocktail ??
826 AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
827
828 if (!genCocktailHeader) {
cb76764e 829 // AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
b4691ee7 830 // AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
831 return 0;
832 }
833 TList* headerList = genCocktailHeader->GetHeaders();
834 for (Int_t i=0; i<headerList->GetEntries(); i++) {
835 pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
836 if (pythiaGenHeader)
837 break;
838 }
839 if(!pythiaGenHeader){
840 AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
841 return 0;
842 }
843 }
844 return pythiaGenHeader;
845
846}
fdceab34 847
03372fd1 848//________________________________________________________________________
9477b466 849AliGenHijingEventHeader* AliPWG4HighPtSpectra::GetHijingEventHeader(AliMCEvent *mcEvent)
850{
03372fd1 851
852 if(!mcEvent)return 0;
853 AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
854 AliGenHijingEventHeader* hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(genHeader);
855 if(!hijingGenHeader){
856 // cocktail ??
857 AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
858
859 if (!genCocktailHeader) {
3abf7501 860 AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
03372fd1 861 // AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
862 return 0;
863 }
864 TList* headerList = genCocktailHeader->GetHeaders();
865 for (Int_t i=0; i<headerList->GetEntries(); i++) {
866 hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(headerList->At(i));
867 if (hijingGenHeader)
868 break;
869 }
870 if(!hijingGenHeader){
871 AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Hijing event header not found");
872 return 0;
873 }
874 }
875 return hijingGenHeader;
876
877}
878
3abf7501 879//________________________________________________________________________
880AliGenHijingEventHeader* AliPWG4HighPtSpectra::GetHijingEventHeaderAOD(AliAODEvent *aodEvent)
881{
882 AliGenHijingEventHeader* hijingGenHeader = 0x0;
883 if(aodEvent) {
884 AliAODMCHeader* mcHeader = (AliAODMCHeader*) aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
ba02ea81 885 TList* headerList = mcHeader->GetCocktailHeaders();
886 for (Int_t i=0; i<headerList->GetEntries(); i++) {
887 hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(headerList->At(i));
888 if (hijingGenHeader)
889 break;
890 }
3abf7501 891 }
892 return hijingGenHeader;
893}
894
fdceab34 895//___________________________________________________________________________
896void AliPWG4HighPtSpectra::Terminate(Option_t*)
897{
898 // The Terminate() function is the last function to be called during
899 // a query. It always runs on the client, it can be used to present
900 // the results graphically or save the results to file.
9b58297c 901
fdceab34 902}
903
fdceab34 904//___________________________________________________________________________
9477b466 905void AliPWG4HighPtSpectra::CreateOutputObjects()
906{
fdceab34 907 //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
908 //TO BE SET BEFORE THE EXECUTION OF THE TASK
909 //
f51451be 910 AliDebug(2,Form("CreateOutputObjects CreateOutputObjects of task %s", GetName()));
df943115 911
67ebd013 912 Bool_t oldStatus = TH1::AddDirectoryStatus();
913 TH1::AddDirectory(kFALSE);
914
fdceab34 915 //slot #1
b5cc0c6d 916 OpenFile(0);
917 fHistList = new TList();
b4691ee7 918 fHistList->SetOwner(kTRUE);
b5cc0c6d 919 fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
920 fHistList->Add(fNEventAll);
921 fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
922 fHistList->Add(fNEventSel);
fdceab34 923
cb76764e 924 fNEventReject = new TH1F("fNEventReject","Reason events are rejectected for analysis",20,0,20);
2b553e6f 925 //Set labels
9477b466 926 fNEventReject->Fill("noESD/AOD",0);
2b553e6f 927 fNEventReject->Fill("Trigger",0);
928 fNEventReject->Fill("NTracks<2",0);
929 fNEventReject->Fill("noVTX",0);
930 fNEventReject->Fill("VtxStatus",0);
931 fNEventReject->Fill("NCont<2",0);
932 fNEventReject->Fill("ZVTX>10",0);
933 fNEventReject->Fill("cent",0);
934 fNEventReject->Fill("cent>90",0);
cb76764e 935 fHistList->Add(fNEventReject);
936
2b553e6f 937 fh1Centrality = new TH1F("fh1Centrality","fh1Centrality; Centrality %",100,0,100);
938 fHistList->Add(fh1Centrality);
939
b4691ee7 940 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
941 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
942 fHistList->Add(fh1Xsec);
943
944 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
945 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
946 fHistList->Add(fh1Trials);
947
948 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
949 fHistList->Add(fh1PtHard);
950 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
951 fHistList->Add(fh1PtHardTrials);
952
a337a5a9 953 Int_t fgkNPtBins = 100;
954 Float_t kMinPt = 0.;
955 Float_t kMaxPt = 100.;
956 Double_t *binsPt = new Double_t[fgkNPtBins+1];
957 for(Int_t i=0; i<=fgkNPtBins; i++) binsPt[i]=(Double_t)kMinPt + (kMaxPt-kMinPt)/fgkNPtBins*(Double_t)i ;
958
959 Int_t fgkNRel1PtUncertaintyBins=50;
960 Float_t fgkRel1PtUncertaintyMin = 0.;
961 Float_t fgkRel1PtUncertaintyMax = 1.;
962
963 Double_t *binsRel1PtUncertainty=new Double_t[fgkNRel1PtUncertaintyBins+1];
964 for(Int_t i=0; i<=fgkNRel1PtUncertaintyBins; i++) binsRel1PtUncertainty[i]=(Double_t)fgkRel1PtUncertaintyMin + (fgkRel1PtUncertaintyMax-fgkRel1PtUncertaintyMin)/fgkNRel1PtUncertaintyBins*(Double_t)i ;
965
966 fPtRelUncertainty1PtPrim = new TH2F("fPtRelUncertainty1PtPrim","fPtRelUncertainty1PtPrim",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty);
967 fHistList->Add(fPtRelUncertainty1PtPrim);
968
969 fPtRelUncertainty1PtSec = new TH2F("fPtRelUncertainty1PtSec","fPtRelUncertainty1PtSec",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty);
970 fHistList->Add(fPtRelUncertainty1PtSec);
971
67ebd013 972 TH1::AddDirectory(oldStatus);
973
f6d670e1 974 OpenFile(1);
975 OpenFile(2);
976
11245558 977 PostData(0,fHistList);
978 PostData(1,fCFManagerPos->GetParticleContainer());
979 PostData(2,fCFManagerNeg->GetParticleContainer());
980
a337a5a9 981 if(binsPt) delete [] binsPt;
982 if(binsRel1PtUncertainty) delete [] binsRel1PtUncertainty;
983
fdceab34 984}
985
03372fd1 986//________________________________________________________________________
3abf7501 987Bool_t AliPWG4HighPtSpectra::IsHIJINGParticle(Int_t label)
9477b466 988{
03372fd1 989 //
990 // Return kTRUE in case particle is from HIJING event
991 //
992
3abf7501 993 AliGenHijingEventHeader* hijingHeader;
994 if(fReadAODData)
995 hijingHeader = GetHijingEventHeaderAOD(fAOD);
996 else
997 hijingHeader = GetHijingEventHeader(fMC);
03372fd1 998
999 Int_t nproduced = hijingHeader->NProduced();
1000
3abf7501 1001 if(fReadAODData) {
1002 AliAODMCParticle* mom = (AliAODMCParticle*) fArrayMCAOD->At(label);
1003 Int_t iMom = label;
1004 Int_t iParent = mom->GetMother();
1005 while(iParent!=-1){
1006 iMom = iParent;
1007 mom = (AliAODMCParticle*) fArrayMCAOD->At(iMom);
1008 iParent = mom->GetMother();
1009 }
03372fd1 1010
3abf7501 1011 if(iMom<nproduced) {
1012 return kTRUE;
1013 } else {
1014 return kFALSE;
1015 }
1016 }
1017 else { //ESD analysis
1018 TParticle * mom = fStack->Particle(label);
1019 Int_t iMom = label;
1020 Int_t iParent = mom->GetFirstMother();
1021 while(iParent!=-1){
1022 iMom = iParent;
1023 mom = fStack->Particle(iMom);
1024 iParent = mom->GetFirstMother();
1025 }
1026
1027 if(iMom<nproduced) {
1028 return kTRUE;
1029 } else {
1030 return kFALSE;
1031 }
03372fd1 1032 }
03372fd1 1033}
1034
fdceab34 1035#endif