Updating track cuts and removing functionality that is now in separate cosmics tasks
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliPWG4HighPtSpectra.cxx
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
27 #ifndef ALIPWG4HIGHPTSPECTRA_CXX
28 #define ALIPWG4HIGHPTSPECTRA_CXX
29
30 #include "AliPWG4HighPtSpectra.h"
31
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"
39
40 #include "AliAnalysisManager.h"
41 #include "AliESDInputHandler.h"
42 #include "AliESDtrack.h"
43 #include "AliESDtrackCuts.h"
44 #include "AliExternalTrackParam.h"
45
46 #include "AliLog.h"
47
48 #include "AliStack.h"
49 #include "TParticle.h"
50 #include "AliMCEvent.h"
51 #include "AliMCEventHandler.h"
52 #include "AliCFContainer.h"
53
54 //#include "$ALICE_ROOT/PWG4/JetTasks/AliAnalysisHelperJetTasks.h"
55
56 //#include <iostream>
57 using namespace std; //required for resolving the 'cout' symbol
58
59 ClassImp(AliPWG4HighPtSpectra)
60
61 //__________________________________________________________________________
62 AliPWG4HighPtSpectra::AliPWG4HighPtSpectra() : AliAnalysisTask("AliPWG4HighPtSpectra", ""), 
63   fReadAODData(0),
64   fCFManagerPos(0x0),
65   fCFManagerNeg(0x0),
66   fESD(0),
67   fTrackCuts(0),
68   fTrackCutsTPConly(0),
69   fHistList(0),
70   fNEventAll(0),
71   fNEventSel(0)
72 {
73   //
74   //Default ctor
75   //
76 }
77 //___________________________________________________________________________
78 AliPWG4HighPtSpectra::AliPWG4HighPtSpectra(const Char_t* name) :
79   AliAnalysisTask(name,""),
80   fReadAODData(0),
81   fCFManagerPos(0x0),
82   fCFManagerNeg(0x0),
83   fESD(0),
84   fTrackCuts(),
85   fTrackCutsTPConly(0),
86   fHistList(0),
87   fNEventAll(0),
88   fNEventSel(0)
89 {
90   //
91   // Constructor. Initialization of Inputs and Outputs
92   //
93   AliDebug(2,Form("AliPWG4HighPtSpectra","Calling Constructor"));
94   // Input slot #0 works with a TChain ESD
95   DefineInput(0, TChain::Class());
96   // Output slot #0 writes into a TList
97   DefineOutput(0,TList::Class());
98   // Output slot #1, #2 writes into a AliCFContainer
99   DefineOutput(1,AliCFContainer::Class());
100   DefineOutput(2,AliCFContainer::Class());
101   // Output slot #3 writes into a AliESDtrackCuts
102   DefineOutput(3, AliESDtrackCuts::Class());
103   DefineOutput(4, AliESDtrackCuts::Class());
104 }
105
106 //________________________________________________________________________
107 void AliPWG4HighPtSpectra::LocalInit()
108 {
109   //
110   // Only called once at beginning
111   //
112   PostData(3,fTrackCuts);
113   PostData(4,fTrackCutsTPConly);
114 }
115
116 //________________________________________________________________________
117 void AliPWG4HighPtSpectra::ConnectInputData(Option_t *) 
118 {
119   // Connect ESD here
120   // Called once
121   AliDebug(2,Form(">> AliPWG4HighPtSpectra::ConnectInputData \n"));
122   //  cout << "cout >> AliPWG4HighPtSpectra::ConnectInputData" << endl;
123   printf(">> AliPWG4HighPtSpectra::ConnectInputData \n");
124
125   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
126   if (!tree) {
127     AliDebug(2,Form("ERROR: Could not read chain from input slot 0"));
128   } else {
129     
130     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
131     
132     if (!esdH) {
133       AliDebug(2,Form("ERROR: Could not get ESDInputHandler"));
134     } else {
135       fESD = esdH->GetEvent();
136     }
137   }
138   
139 }
140 //_________________________________________________
141 void AliPWG4HighPtSpectra::Exec(Option_t *)
142 {
143   //
144   // Main loop function
145   //
146   AliDebug(2,Form(">> AliPWG4HighPtSpectra::Exec \n"));  
147
148   // All events without selection
149   fNEventAll->Fill(0.);
150
151   if (!fESD) {
152     AliDebug(2,Form("ERROR: fESD not available"));
153     PostData(0,fHistList);
154     PostData(1,fCFManagerPos->GetParticleContainer());
155     PostData(2,fCFManagerNeg->GetParticleContainer());
156     return;
157   }
158
159   Bool_t isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
160   if(!isSelected) { //Select collison candidates
161     AliDebug(2,Form(" Trigger Selection: event REJECTED ... "));
162     PostData(0,fHistList);
163     PostData(1,fCFManagerPos->GetParticleContainer());
164     PostData(2,fCFManagerNeg->GetParticleContainer());
165     return;
166   }
167
168   // Process MC truth, therefore we receive the AliAnalysisManager and ask it for the AliMCEventHandler
169   // This handler can return the current MC event
170   
171   AliMCEventHandler *eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
172   //  AliMCEventHandler* eventHandler = (AliMCEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
173   
174   AliStack* stack = 0x0;
175   AliMCEvent* mcEvent = 0x0;
176   
177   if(eventHandler) {
178     mcEvent = eventHandler->MCEvent();
179     if (!mcEvent) {
180       AliDebug(2,Form("ERROR: Could not retrieve MC event"));
181       PostData(0,fHistList);
182       PostData(1,fCFManagerPos->GetParticleContainer());
183     PostData(2,fCFManagerNeg->GetParticleContainer());
184     return;
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()));
192   }
193   
194   const AliESDVertex *vtx = fESD->GetPrimaryVertex();
195   AliDebug(2,Form("Vertex title %s, status %d, nCont %d\n",vtx->GetTitle(), vtx->GetStatus(), vtx->GetNContributors()));
196   // Need vertex cut
197   if (vtx->GetNContributors() < 2) {
198     PostData(0,fHistList);
199     PostData(1,fCFManagerPos->GetParticleContainer());
200     PostData(2,fCFManagerNeg->GetParticleContainer());
201     return;
202   }
203   
204   double primVtx[3];
205   vtx->GetXYZ(primVtx);
206   if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
207     PostData(0,fHistList);
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());
218     return;
219   }
220   Int_t nTracks = fESD->GetNumberOfTracks();
221   AliDebug(2,Form("nTracks %d", nTracks));
222
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
231   // Selected events for analysis
232   fNEventSel->Fill(0.);
233   
234   
235   Double_t containerInputRec[5] ;
236   Double_t containerInputTPConly[5];
237   Double_t containerInputMC[5];
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
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);
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.);
257
258       //fill the container
259       containerInputRec[0] = track->Pt();
260       containerInputRec[1] = track->Phi();
261       containerInputRec[2] = track->Eta();
262       containerInputRec[3] = dca2D;
263       containerInputRec[4] = chi2PerClusterTPC;
264
265       //Store TPC Inner Params for TPConly tracks
266       containerInputTPConly[0] = trackTPC->Pt();
267       containerInputTPConly[1] = trackTPC->Phi();
268       containerInputTPConly[2] = trackTPC->Eta();
269       containerInputTPConly[3] = dca2DTPC/10.; //Divide by 10 in order to store in same container. Should be corrected back when looking at output.
270       containerInputTPConly[4] = chi2PerClusterTPCIter1;//TPC;
271
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);
277         }
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
284         
285         //Only fill the MC containers if MC information is available
286         if(eventHandler) {
287           Int_t label = TMath::Abs(track->GetLabel());
288           TParticle *particle = stack->Particle(label) ;
289           if(!particle) continue;
290
291           containerInputMC[0] = particle->Pt();      
292           containerInputMC[1] = particle->Phi();      
293           containerInputMC[2] = particle->Eta();  
294           containerInputMC[3] = 0.0;      
295           containerInputMC[4] = 0.0;  
296
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             }
305           }
306
307           //Container with secondaries
308           if (!stack->IsPhysicalPrimary(label) ) {
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             }
315           }
316         }
317         
318       }//trackCuts
319
320       delete trackTPCESD;
321     }//track loop
322   
323
324   //Fill MC containters if particles are findable
325   if(eventHandler) {
326     for(int iPart = 1; iPart<(mcEvent->GetNumberOfPrimaries()); iPart++)//stack->GetNprimary();
327       {
328         AliMCParticle *mcPart  = (AliMCParticle*)mcEvent->GetTrack(iPart);
329         if(!mcPart) continue;
330         
331         //fill the container
332         containerInputMC[0] = mcPart->Pt();
333         containerInputMC[1] = mcPart->Phi();      
334         containerInputMC[2] = mcPart->Eta();  
335         containerInputMC[3] = 0.0;
336         containerInputMC[4] = 0.0;
337
338         if(stack->IsPhysicalPrimary(iPart)) {
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);
341         }
342       }
343   }
344   
345   PostData(0,fHistList);
346   PostData(1,fCFManagerPos->GetParticleContainer());
347   PostData(2,fCFManagerNeg->GetParticleContainer());
348   
349 }
350
351
352 //___________________________________________________________________________
353 void 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.
358
359 }
360
361 //___________________________________________________________________________
362 void 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   //
366   AliDebug(2,Form("CreateOutputObjects","CreateOutputObjects of task %s", GetName()));
367
368   Bool_t oldStatus = TH1::AddDirectoryStatus();
369   TH1::AddDirectory(kFALSE); 
370
371   //slot #1
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);
378
379   TH1::AddDirectory(oldStatus);   
380
381 }
382
383 #endif