New book-keeping class
[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 "AliStack.h"
33 #include "TParticle.h"
34 #include "TH1I.h"
35 #include "AliMCEvent.h"
36 #include "AliMCEventHandler.h"
37 #include "AliAnalysisManager.h"
38 #include "AliCFContainer.h"
39 #include "TChain.h"
40 #include "AliESDtrack.h"
41 #include "AliESDtrackCuts.h"
42 #include "AliExternalTrackParam.h"
43 #include "AliESDInputHandler.h"
44 #include "AliAnalysisHelperJetTasks.h"
45
46 using namespace std; //required for resolving the 'cout' symbol
47
48 ClassImp(AliPWG4HighPtSpectra)
49
50 //__________________________________________________________________________
51 AliPWG4HighPtSpectra::AliPWG4HighPtSpectra() : AliAnalysisTask("AliPWG4HighPtSpectra", ""), 
52   fReadAODData(0),
53   fCFManager(0x0),
54   fESD(0),
55   fTrackCuts(0),
56   fTrigger(0),
57   fHistList(0),
58   fNEventAll(0),
59   fNEventSel(0)
60   //  fHistEventsProcessed(0x0)
61 {
62   //
63   //Default ctor
64   //
65 }
66 //___________________________________________________________________________
67 AliPWG4HighPtSpectra::AliPWG4HighPtSpectra(const Char_t* name) :
68   AliAnalysisTask(name,""),
69   fReadAODData(0),
70   fCFManager(0x0),
71   fESD(0),
72   fTrackCuts(),//new AliESDtrackCuts),
73   fTrigger(0),
74   fHistList(0),
75   fNEventAll(0),
76   fNEventSel(0)
77   //  fHistEventsProcessed(0x0)
78 {
79   //
80   // Constructor. Initialization of Inputs and Outputs
81   //
82   AliDebug(2,Form("AliPWG4HighPtQAMC","Calling Constructor"));
83   // Input slot #0 works with a TChain ESD
84   DefineInput(0, TChain::Class());
85   DefineOutput(0,TList::Class());
86   DefineOutput(1,AliCFContainer::Class());
87 }
88
89 //___________________________________________________________________________
90 AliPWG4HighPtSpectra& AliPWG4HighPtSpectra::operator=(const AliPWG4HighPtSpectra& c) 
91 {
92   //
93   // Assignment operator
94   //
95   if (this!=&c) {
96     AliAnalysisTask::operator=(c) ;
97     fReadAODData = c.fReadAODData ;
98     fCFManager  = c.fCFManager;
99     fHistList = c.fHistList;
100     fNEventAll = c.fNEventAll;
101     fNEventSel = c.fNEventSel;
102     //    fHistEventsProcessed = c.fHistEventsProcessed;
103   }
104   return *this;
105 }
106
107 //___________________________________________________________________________
108 AliPWG4HighPtSpectra::AliPWG4HighPtSpectra(const AliPWG4HighPtSpectra& c) :
109   AliAnalysisTask(c),
110   fReadAODData(c.fReadAODData),
111   fCFManager(c.fCFManager),
112   fESD(c.fESD),
113   fTrackCuts(c.fTrackCuts),
114   fTrigger(c.fTrigger),
115   fHistList(c.fHistList),
116   fNEventAll(c.fNEventAll),
117   fNEventSel(c.fNEventSel)
118   //  fHistEventsProcessed(c.fHistEventsProcessed)
119 {
120   //
121   // Copy Constructor
122   //
123 }
124
125 //___________________________________________________________________________
126 AliPWG4HighPtSpectra::~AliPWG4HighPtSpectra() {
127   //
128   //destructor
129   //
130   Info("~AliPWG4HighPtSpectra","Calling Destructor");
131   if (fCFManager)           delete fCFManager ;
132   if (fNEventAll) delete fNEventAll ;
133   if (fNEventSel) delete fNEventSel ;
134 }
135 //________________________________________________________________________
136 void AliPWG4HighPtSpectra::ConnectInputData(Option_t *) 
137 {
138   // Connect ESD here
139   // Called once
140   AliDebug(2,Form(">> AliPWG4HighPtSpectra::ConnectInputData \n"));
141   //  printf(">> AliPWG4HighPtSpectra::ConnectInputData \n");
142
143   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
144   if (!tree) {
145     AliDebug(2,Form("ERROR: Could not read chain from input slot 0"));
146   } else {
147     
148     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
149     
150     if (!esdH) {
151       AliDebug(2,Form("ERROR: Could not get ESDInputHandler"));
152     } else
153       fESD = esdH->GetEvent();
154   }
155  
156 }
157 //_________________________________________________
158 void AliPWG4HighPtSpectra::Exec(Option_t *)//UserExec(Option_t *)
159 {
160   //
161   // Main loop function
162   //
163   AliDebug(2,Form(">> AliPWG4HighPtSpectra::Exec \n"));  
164
165   // All events without selection
166   fNEventAll->Fill(0.);
167
168   if (!fESD) {
169     AliDebug(2,Form("ERROR: fESD not available"));
170     return;
171   }
172
173   //Trigger selection
174   AliAnalysisHelperJetTasks::Trigger trig;
175   trig = (AliAnalysisHelperJetTasks::Trigger)fTrigger;
176   if (AliAnalysisHelperJetTasks::IsTriggerFired(fESD,trig)){
177     AliDebug(2,Form(" Trigger Selection: event ACCEPTED ... "));
178   }else{
179     AliDebug(2,Form(" Trigger Selection: event REJECTED ... "));
180     PostData(0,fHistList);
181     PostData(1,fCFManager->GetParticleContainer());
182     return;
183   } 
184   //  if(!fESD->IsTriggerClassFired("CINT1B-ABCE-NOPF-ALL") || !fESD->IsTriggerClassFired("CSMBB-ABCE-NOPF-ALL")) return;
185
186   // Process MC truth, therefore we receive the AliAnalysisManager and ask it for the AliMCEventHandler
187   // This handler can return the current MC event
188   AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
189   AliStack* stack = 0x0;
190   AliMCEvent* mcEvent = 0x0;
191
192   if(eventHandler) {
193     mcEvent = eventHandler->MCEvent();
194     if (!mcEvent) {
195       AliDebug(2,Form("ERROR: Could not retrieve MC event"));
196       return;
197     }
198     
199     AliDebug(2,Form("MC particles: %d", mcEvent->GetNumberOfTracks()));
200     
201     stack = mcEvent->Stack();                //Particles Stack
202     
203     AliDebug(2,Form("MC particles stack: %d", stack->GetNtrack()));
204   }
205   
206   const AliESDVertex *vtx = fESD->GetPrimaryVertex();
207   // Need vertex cut
208   if (vtx->GetNContributors() < 2){
209     PostData(0,fHistList);
210     PostData(1,fCFManager->GetParticleContainer());
211     return;
212   }
213   double primVtx[3];
214   vtx->GetXYZ(primVtx);
215   if(TMath::Abs(primVtx[0]>1. )|| TMath::Abs(primVtx[1]>1.) || TMath::Abs(primVtx[2]>10.)){
216     PostData(0,fHistList);
217     PostData(1,fCFManager->GetParticleContainer());
218     return;
219   }
220   AliDebug(2,Form("Vertex title %s, status %d, nCont %d\n",vtx->GetTitle(), vtx->GetStatus(), vtx->GetNContributors()));
221
222   Int_t nTracks = fESD->GetNumberOfTracks();
223   AliDebug(2,Form("nTracks %d", nTracks));
224
225   // Selected events for analysis
226   fNEventSel->Fill(0.);
227
228   Double_t containerInputRec[1] ;
229   Double_t containerInputTPConly[1] ;
230   Double_t containerInputMC[1] ;
231   //Now go to rec level
232   for (Int_t iTrack = 0; iTrack<nTracks; iTrack++) 
233     {   
234       if(!fESD->GetTrack(iTrack) ) continue;
235       AliESDtrack* track = fESD->GetTrack(iTrack);
236       if(!(AliExternalTrackParam *)track->GetTPCInnerParam()) continue;
237       AliExternalTrackParam *trackTPC = (AliExternalTrackParam *)track->GetTPCInnerParam();
238       if(!track || !trackTPC) continue;
239
240     
241       //fill the container
242       containerInputRec[0] = track->Pt();
243       containerInputTPConly[0] = trackTPC->Pt();
244
245       if (fTrackCuts->AcceptTrack(track)) {
246         fCFManager->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
247         fCFManager->GetParticleContainer()->Fill(containerInputTPConly,kStepReconstructedTPCOnly);
248         
249         //Only fill the secondary particle container if MC information is available
250         if(eventHandler) {
251           Int_t label = TMath::Abs(track->GetLabel());
252           TParticle *particle = stack->Particle(label) ;
253           if(!particle) continue;
254           containerInputMC[0] = particle->Pt();      
255           fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepReconstructedMC);
256           if (!stack->IsPhysicalPrimary(label) ) {
257             fCFManager->GetParticleContainer()->Fill(containerInputRec,kStepSecondaries);
258           }
259         }
260       }
261       
262     }
263
264   if(eventHandler) {
265     for(int iPart = 1; iPart<(mcEvent->GetNumberOfTracks()); iPart++)//stack->GetNprimary();
266       {
267         AliMCParticle *mcPart  = (AliMCParticle*)mcEvent->GetTrack(iPart);
268         if(!mcPart) continue;
269         
270         //fill the container
271         containerInputMC[0] = mcPart->Pt();
272         
273         if (!fCFManager->CheckParticleCuts(3,mcPart)) continue ;
274         
275         int counter;
276         
277         Float_t trackLengthTPC = mcPart->GetTPCTrackLength(fESD->GetMagneticField(),0.1,counter,3.0);
278         
279         if(trackLengthTPC>80.) fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepMCtrackable) ;
280         
281       }
282   }
283
284    PostData(0,fHistList);
285    PostData(1,fCFManager->GetParticleContainer());
286
287 }
288
289
290 //___________________________________________________________________________
291 void AliPWG4HighPtSpectra::Terminate(Option_t*)
292 {
293   // The Terminate() function is the last function to be called during
294   // a query. It always runs on the client, it can be used to present
295   // the results graphically or save the results to file.
296
297
298 }
299
300 //___________________________________________________________________________
301 void AliPWG4HighPtSpectra::CreateOutputObjects() {
302   //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
303   //TO BE SET BEFORE THE EXECUTION OF THE TASK
304   //
305   AliDebug(2,Form("CreateOutputObjects","CreateOutputObjects of task %s", GetName()));
306
307   //slot #1
308   OpenFile(0);
309   fHistList = new TList();
310   fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
311   fHistList->Add(fNEventAll);
312   fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
313   fHistList->Add(fNEventSel);
314
315 }
316
317 #endif