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