Updating track cuts and removing functionality that is now in separate cosmics tasks
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / 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//-----------------------------------------------------------------------
17// Example of task (running locally, on AliEn and CAF),
18// which provides standard way of calculating acceptance and efficiency
19// between different steps of the procedure.
20// The ouptut of the task is a AliCFContainer from which the efficiencies
21// can be calculated
22//-----------------------------------------------------------------------
23// Author : Marta Verweij - UU
24//-----------------------------------------------------------------------
25
26
67ebd013 27#ifndef ALIPWG4HIGHPTSPECTRA_CXX
28#define ALIPWG4HIGHPTSPECTRA_CXX
fdceab34 29
30#include "AliPWG4HighPtSpectra.h"
31
67ebd013 32#include "TVector3.h"
33#include <iostream>
34#include "TH1.h"
35#include "TH2.h"
36#include "TH3.h"
37#include "TList.h"
38#include "TChain.h"
67ebd013 39
40#include "AliAnalysisManager.h"
41#include "AliESDInputHandler.h"
42#include "AliESDtrack.h"
43#include "AliESDtrackCuts.h"
44#include "AliExternalTrackParam.h"
65e8ecdd 45
67ebd013 46#include "AliLog.h"
47
fdceab34 48#include "AliStack.h"
49#include "TParticle.h"
fdceab34 50#include "AliMCEvent.h"
51#include "AliMCEventHandler.h"
fdceab34 52#include "AliCFContainer.h"
fdceab34 53
67ebd013 54//#include "$ALICE_ROOT/PWG4/JetTasks/AliAnalysisHelperJetTasks.h"
55
56//#include <iostream>
fdceab34 57using namespace std; //required for resolving the 'cout' symbol
58
59ClassImp(AliPWG4HighPtSpectra)
60
61//__________________________________________________________________________
62AliPWG4HighPtSpectra::AliPWG4HighPtSpectra() : AliAnalysisTask("AliPWG4HighPtSpectra", ""),
63 fReadAODData(0),
67ebd013 64 fCFManagerPos(0x0),
65 fCFManagerNeg(0x0),
fdceab34 66 fESD(0),
67 fTrackCuts(0),
1f329128 68 fTrackCutsTPConly(0),
fdceab34 69 fHistList(0),
b5cc0c6d 70 fNEventAll(0),
71 fNEventSel(0)
fdceab34 72{
73 //
74 //Default ctor
75 //
76}
77//___________________________________________________________________________
78AliPWG4HighPtSpectra::AliPWG4HighPtSpectra(const Char_t* name) :
79 AliAnalysisTask(name,""),
80 fReadAODData(0),
67ebd013 81 fCFManagerPos(0x0),
82 fCFManagerNeg(0x0),
fdceab34 83 fESD(0),
67ebd013 84 fTrackCuts(),
1f329128 85 fTrackCutsTPConly(0),
fdceab34 86 fHistList(0),
b5cc0c6d 87 fNEventAll(0),
88 fNEventSel(0)
fdceab34 89{
90 //
91 // Constructor. Initialization of Inputs and Outputs
92 //
67ebd013 93 AliDebug(2,Form("AliPWG4HighPtSpectra","Calling Constructor"));
fdceab34 94 // Input slot #0 works with a TChain ESD
95 DefineInput(0, TChain::Class());
e5abcde9 96 // Output slot #0 writes into a TList
fdceab34 97 DefineOutput(0,TList::Class());
e5abcde9 98 // Output slot #1, #2 writes into a AliCFContainer
fdceab34 99 DefineOutput(1,AliCFContainer::Class());
67ebd013 100 DefineOutput(2,AliCFContainer::Class());
e5abcde9 101 // Output slot #3 writes into a AliESDtrackCuts
102 DefineOutput(3, AliESDtrackCuts::Class());
1f329128 103 DefineOutput(4, AliESDtrackCuts::Class());
e5abcde9 104}
105
106//________________________________________________________________________
107void AliPWG4HighPtSpectra::LocalInit()
108{
109 //
110 // Only called once at beginning
111 //
112 PostData(3,fTrackCuts);
1f329128 113 PostData(4,fTrackCutsTPConly);
fdceab34 114}
115
fdceab34 116//________________________________________________________________________
117void AliPWG4HighPtSpectra::ConnectInputData(Option_t *)
118{
119 // Connect ESD here
120 // Called once
df943115 121 AliDebug(2,Form(">> AliPWG4HighPtSpectra::ConnectInputData \n"));
67ebd013 122 // cout << "cout >> AliPWG4HighPtSpectra::ConnectInputData" << endl;
123 printf(">> AliPWG4HighPtSpectra::ConnectInputData \n");
fdceab34 124
125 TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
126 if (!tree) {
df943115 127 AliDebug(2,Form("ERROR: Could not read chain from input slot 0"));
fdceab34 128 } else {
129
130 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
131
132 if (!esdH) {
df943115 133 AliDebug(2,Form("ERROR: Could not get ESDInputHandler"));
67ebd013 134 } else {
fdceab34 135 fESD = esdH->GetEvent();
67ebd013 136 }
fdceab34 137 }
67ebd013 138
fdceab34 139}
140//_________________________________________________
67ebd013 141void AliPWG4HighPtSpectra::Exec(Option_t *)
fdceab34 142{
143 //
144 // Main loop function
145 //
df943115 146 AliDebug(2,Form(">> AliPWG4HighPtSpectra::Exec \n"));
fdceab34 147
b5cc0c6d 148 // All events without selection
149 fNEventAll->Fill(0.);
150
fdceab34 151 if (!fESD) {
df943115 152 AliDebug(2,Form("ERROR: fESD not available"));
67ebd013 153 PostData(0,fHistList);
154 PostData(1,fCFManagerPos->GetParticleContainer());
155 PostData(2,fCFManagerNeg->GetParticleContainer());
fdceab34 156 return;
157 }
158
67ebd013 159 Bool_t isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
160 if(!isSelected) { //Select collison candidates
b5cc0c6d 161 AliDebug(2,Form(" Trigger Selection: event REJECTED ... "));
cd9a6fa2 162 PostData(0,fHistList);
67ebd013 163 PostData(1,fCFManagerPos->GetParticleContainer());
164 PostData(2,fCFManagerNeg->GetParticleContainer());
b5cc0c6d 165 return;
67ebd013 166 }
b5cc0c6d 167
fdceab34 168 // Process MC truth, therefore we receive the AliAnalysisManager and ask it for the AliMCEventHandler
169 // This handler can return the current MC event
67ebd013 170
171 AliMCEventHandler *eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
172 // AliMCEventHandler* eventHandler = (AliMCEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
173
0cc50ca6 174 AliStack* stack = 0x0;
175 AliMCEvent* mcEvent = 0x0;
67ebd013 176
0cc50ca6 177 if(eventHandler) {
178 mcEvent = eventHandler->MCEvent();
179 if (!mcEvent) {
180 AliDebug(2,Form("ERROR: Could not retrieve MC event"));
67ebd013 181 PostData(0,fHistList);
182 PostData(1,fCFManagerPos->GetParticleContainer());
183 PostData(2,fCFManagerNeg->GetParticleContainer());
184 return;
0cc50ca6 185 }
186
187 AliDebug(2,Form("MC particles: %d", mcEvent->GetNumberOfTracks()));
188
189 stack = mcEvent->Stack(); //Particles Stack
190
191 AliDebug(2,Form("MC particles stack: %d", stack->GetNtrack()));
fdceab34 192 }
0cc50ca6 193
fdceab34 194 const AliESDVertex *vtx = fESD->GetPrimaryVertex();
67ebd013 195 AliDebug(2,Form("Vertex title %s, status %d, nCont %d\n",vtx->GetTitle(), vtx->GetStatus(), vtx->GetNContributors()));
fdceab34 196 // Need vertex cut
67ebd013 197 if (vtx->GetNContributors() < 2) {
cd9a6fa2 198 PostData(0,fHistList);
67ebd013 199 PostData(1,fCFManagerPos->GetParticleContainer());
200 PostData(2,fCFManagerNeg->GetParticleContainer());
cd9a6fa2 201 return;
202 }
67ebd013 203
b5cc0c6d 204 double primVtx[3];
205 vtx->GetXYZ(primVtx);
67ebd013 206 if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
cd9a6fa2 207 PostData(0,fHistList);
67ebd013 208 PostData(1,fCFManagerPos->GetParticleContainer());
209 PostData(2,fCFManagerNeg->GetParticleContainer());
210 return;
211 }
212
213 if(!fESD->GetNumberOfTracks() || fESD->GetNumberOfTracks()<2){
214 // Post output data
215 PostData(0,fHistList);
216 PostData(1,fCFManagerPos->GetParticleContainer());
217 PostData(2,fCFManagerNeg->GetParticleContainer());
cd9a6fa2 218 return;
219 }
fdceab34 220 Int_t nTracks = fESD->GetNumberOfTracks();
221 AliDebug(2,Form("nTracks %d", nTracks));
fdceab34 222
67ebd013 223 if(!fTrackCuts) {
224 // Post output data
225 PostData(0,fHistList);
226 PostData(1,fCFManagerPos->GetParticleContainer());
227 PostData(2,fCFManagerNeg->GetParticleContainer());
228 return;
229 }
230
b5cc0c6d 231 // Selected events for analysis
232 fNEventSel->Fill(0.);
67ebd013 233
234
235 Double_t containerInputRec[5] ;
236 Double_t containerInputTPConly[5];
237 Double_t containerInputMC[5];
fdceab34 238 //Now go to rec level
239 for (Int_t iTrack = 0; iTrack<nTracks; iTrack++)
240 {
241 if(!fESD->GetTrack(iTrack) ) continue;
242 AliESDtrack* track = fESD->GetTrack(iTrack);
243 if(!(AliExternalTrackParam *)track->GetTPCInnerParam()) continue;
244 AliExternalTrackParam *trackTPC = (AliExternalTrackParam *)track->GetTPCInnerParam();
245 if(!track || !trackTPC) continue;
246
67ebd013 247 Float_t dca2D, dcaZ;
248 track->GetImpactParameters(dca2D,dcaZ);
249 Float_t dca2DTPC, dcaZTPC;
250 track->GetImpactParametersTPC(dca2DTPC,dcaZTPC);
251 Float_t chi2PerClusterTPC = -1.;
252 Float_t nClustersTPC = track->GetTPCNcls();//track->GetTPCclusters(0);
38ecb6a5 253 if(nClustersTPC>0.) chi2PerClusterTPC = track->GetTPCchi2()/(2.*nClustersTPC-5.);
254 Float_t chi2PerClusterTPCIter1 = -1.;
255 Float_t nClustersTPCIter1 = track->GetTPCNclsIter1();
256 if(nClustersTPCIter1>0.) chi2PerClusterTPCIter1 = track->GetTPCchi2Iter1()/(2.*nClustersTPCIter1-5.);
67ebd013 257
fdceab34 258 //fill the container
259 containerInputRec[0] = track->Pt();
67ebd013 260 containerInputRec[1] = track->Phi();
261 containerInputRec[2] = track->Eta();
262 containerInputRec[3] = dca2D;
263 containerInputRec[4] = chi2PerClusterTPC;
264
1f329128 265 //Store TPC Inner Params for TPConly tracks
fdceab34 266 containerInputTPConly[0] = trackTPC->Pt();
67ebd013 267 containerInputTPConly[1] = trackTPC->Phi();
268 containerInputTPConly[2] = trackTPC->Eta();
add06ceb 269 containerInputTPConly[3] = dca2DTPC/10.; //Divide by 10 in order to store in same container. Should be corrected back when looking at output.
38ecb6a5 270 containerInputTPConly[4] = chi2PerClusterTPCIter1;//TPC;
0cc50ca6 271
1f329128 272 AliESDtrack* trackTPCESD = fTrackCutsTPConly->GetTPCOnlyTrack(fESD, iTrack);
273 if(trackTPCESD) {
274 if (fTrackCutsTPConly->AcceptTrack(trackTPCESD)) {
275 if(trackTPC->GetSign()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputTPConly,kStepReconstructedTPCOnly);
276 if(trackTPC->GetSign()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputTPConly,kStepReconstructedTPCOnly);
67ebd013 277 }
1f329128 278 }
279
280 if (fTrackCuts->AcceptTrack(track)) {
281 if(track->GetSign()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
282 if(track->GetSign()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
283
67ebd013 284
1f329128 285 //Only fill the MC containers if MC information is available
0cc50ca6 286 if(eventHandler) {
287 Int_t label = TMath::Abs(track->GetLabel());
288 TParticle *particle = stack->Particle(label) ;
289 if(!particle) continue;
65e8ecdd 290
0cc50ca6 291 containerInputMC[0] = particle->Pt();
67ebd013 292 containerInputMC[1] = particle->Phi();
293 containerInputMC[2] = particle->Eta();
294 containerInputMC[3] = 0.0;
295 containerInputMC[4] = 0.0;
296
1f329128 297 //Container with primaries
298 if(stack->IsPhysicalPrimary(label)) {
299 if(particle->GetPDG()->Charge()>0.) {
300 fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepReconstructedMC);
301 }
302 if(particle->GetPDG()->Charge()<0.) {
303 fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepReconstructedMC);
304 }
67ebd013 305 }
306
1f329128 307 //Container with secondaries
0cc50ca6 308 if (!stack->IsPhysicalPrimary(label) ) {
67ebd013 309 if(particle->GetPDG()->Charge()>0.) {
310 fCFManagerPos->GetParticleContainer()->Fill(containerInputRec,kStepSecondaries);
311 }
312 if(particle->GetPDG()->Charge()<0.) {
313 fCFManagerNeg->GetParticleContainer()->Fill(containerInputRec,kStepSecondaries);
314 }
0cc50ca6 315 }
316 }
67ebd013 317
1f329128 318 }//trackCuts
60829498 319
320 delete trackTPCESD;
1f329128 321 }//track loop
67ebd013 322
fdceab34 323
65e8ecdd 324 //Fill MC containters if particles are findable
0cc50ca6 325 if(eventHandler) {
1f329128 326 for(int iPart = 1; iPart<(mcEvent->GetNumberOfPrimaries()); iPart++)//stack->GetNprimary();
0cc50ca6 327 {
328 AliMCParticle *mcPart = (AliMCParticle*)mcEvent->GetTrack(iPart);
329 if(!mcPart) continue;
330
331 //fill the container
332 containerInputMC[0] = mcPart->Pt();
67ebd013 333 containerInputMC[1] = mcPart->Phi();
334 containerInputMC[2] = mcPart->Eta();
335 containerInputMC[3] = 0.0;
336 containerInputMC[4] = 0.0;
337
1f329128 338 if(stack->IsPhysicalPrimary(iPart)) {
9b58297c 339 if(mcPart->Charge()>0. && fCFManagerPos->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
340 if(mcPart->Charge()<0. && fCFManagerNeg->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
67ebd013 341 }
0cc50ca6 342 }
343 }
67ebd013 344
67ebd013 345 PostData(0,fHistList);
346 PostData(1,fCFManagerPos->GetParticleContainer());
347 PostData(2,fCFManagerNeg->GetParticleContainer());
348
fdceab34 349}
350
351
352//___________________________________________________________________________
353void AliPWG4HighPtSpectra::Terminate(Option_t*)
354{
355 // The Terminate() function is the last function to be called during
356 // a query. It always runs on the client, it can be used to present
357 // the results graphically or save the results to file.
9b58297c 358
fdceab34 359}
360
fdceab34 361//___________________________________________________________________________
362void AliPWG4HighPtSpectra::CreateOutputObjects() {
363 //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
364 //TO BE SET BEFORE THE EXECUTION OF THE TASK
365 //
df943115 366 AliDebug(2,Form("CreateOutputObjects","CreateOutputObjects of task %s", GetName()));
367
67ebd013 368 Bool_t oldStatus = TH1::AddDirectoryStatus();
369 TH1::AddDirectory(kFALSE);
370
fdceab34 371 //slot #1
b5cc0c6d 372 OpenFile(0);
373 fHistList = new TList();
374 fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
375 fHistList->Add(fNEventAll);
376 fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
377 fHistList->Add(fNEventSel);
fdceab34 378
67ebd013 379 TH1::AddDirectory(oldStatus);
380
fdceab34 381}
382
383#endif