d73d3bfbe6d67665ec853cd1a946e8865551c661
[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   fTrigger(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   fTrigger(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 }
104
105 //________________________________________________________________________
106 void AliPWG4HighPtSpectra::LocalInit()
107 {
108   //
109   // Only called once at beginning
110   //
111   PostData(3,fTrackCuts);
112 }
113
114 // //___________________________________________________________________________
115 // AliPWG4HighPtSpectra& AliPWG4HighPtSpectra::operator=(const AliPWG4HighPtSpectra& c) 
116 // {
117 //   //
118 //   // Assignment operator
119 //   //
120 //   if (this!=&c) {
121 //     AliAnalysisTask::operator=(c) ;
122 //     fReadAODData = c.fReadAODData ;
123 //     fCFManagerPos  = c.fCFManagerPos;
124 //     fCFManagerNeg  = c.fCFManagerNeg;
125 //     fHistList = c.fHistList;
126 //     fNEventAll = c.fNEventAll;
127 //     fNEventSel = c.fNEventSel;
128 //   }
129 //   return *this;
130 // }
131
132 // //___________________________________________________________________________
133 // AliPWG4HighPtSpectra::AliPWG4HighPtSpectra(const AliPWG4HighPtSpectra& c) :
134 //   AliAnalysisTask(c),
135 //   fReadAODData(c.fReadAODData),
136 //   fCFManagerPos(c.fCFManagerPos),
137 //   fCFManagerNeg(c.fCFManagerNeg),
138 //   fESD(c.fESD),
139 //   fTrackCuts(c.fTrackCuts),
140 //   fTrigger(c.fTrigger),
141 //   fHistList(c.fHistList),
142 //   fNEventAll(c.fNEventAll),
143 //   fNEventSel(c.fNEventSel)
144 // {
145 //   //
146 //   // Copy Constructor
147 //   //
148 // }
149
150 // //___________________________________________________________________________
151 // AliPWG4HighPtSpectra::~AliPWG4HighPtSpectra() {
152 //   //
153 //   //destructor
154 //   //
155 //   Info("~AliPWG4HighPtSpectra","Calling Destructor");
156 //   if (fCFManagerPos)           delete fCFManagerPos ;
157 //   if (fCFManagerNeg)           delete fCFManagerNeg ;
158 //   if (fNEventAll)              delete fNEventAll ;
159 //   if (fNEventSel)              delete fNEventSel ;
160 // }
161 //________________________________________________________________________
162 void AliPWG4HighPtSpectra::ConnectInputData(Option_t *) 
163 {
164   // Connect ESD here
165   // Called once
166   AliDebug(2,Form(">> AliPWG4HighPtSpectra::ConnectInputData \n"));
167   //  cout << "cout >> AliPWG4HighPtSpectra::ConnectInputData" << endl;
168   printf(">> AliPWG4HighPtSpectra::ConnectInputData \n");
169
170   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
171   if (!tree) {
172     AliDebug(2,Form("ERROR: Could not read chain from input slot 0"));
173   } else {
174     
175     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
176     
177     if (!esdH) {
178       AliDebug(2,Form("ERROR: Could not get ESDInputHandler"));
179     } else {
180       fESD = esdH->GetEvent();
181     }
182   }
183   
184 }
185 //_________________________________________________
186 void AliPWG4HighPtSpectra::Exec(Option_t *)
187 {
188   //
189   // Main loop function
190   //
191   AliDebug(2,Form(">> AliPWG4HighPtSpectra::Exec \n"));  
192
193   // All events without selection
194   fNEventAll->Fill(0.);
195
196   if (!fESD) {
197     AliDebug(2,Form("ERROR: fESD not available"));
198     PostData(0,fHistList);
199     PostData(1,fCFManagerPos->GetParticleContainer());
200     PostData(2,fCFManagerNeg->GetParticleContainer());
201     return;
202   }
203
204   Bool_t isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
205   if(!isSelected) { //Select collison candidates
206     AliDebug(2,Form(" Trigger Selection: event REJECTED ... "));
207     PostData(0,fHistList);
208     PostData(1,fCFManagerPos->GetParticleContainer());
209     PostData(2,fCFManagerNeg->GetParticleContainer());
210     return;
211   }
212
213   // Process MC truth, therefore we receive the AliAnalysisManager and ask it for the AliMCEventHandler
214   // This handler can return the current MC event
215   
216   AliMCEventHandler *eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
217   //  AliMCEventHandler* eventHandler = (AliMCEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
218   
219   AliStack* stack = 0x0;
220   AliMCEvent* mcEvent = 0x0;
221   
222   if(eventHandler) {
223     mcEvent = eventHandler->MCEvent();
224     if (!mcEvent) {
225       AliDebug(2,Form("ERROR: Could not retrieve MC event"));
226       PostData(0,fHistList);
227       PostData(1,fCFManagerPos->GetParticleContainer());
228     PostData(2,fCFManagerNeg->GetParticleContainer());
229     return;
230     }
231     
232     AliDebug(2,Form("MC particles: %d", mcEvent->GetNumberOfTracks()));
233     
234     stack = mcEvent->Stack();                //Particles Stack
235     
236     AliDebug(2,Form("MC particles stack: %d", stack->GetNtrack()));
237   }
238   
239   const AliESDVertex *vtx = fESD->GetPrimaryVertex();
240   AliDebug(2,Form("Vertex title %s, status %d, nCont %d\n",vtx->GetTitle(), vtx->GetStatus(), vtx->GetNContributors()));
241   // Need vertex cut
242   if (vtx->GetNContributors() < 2) {
243     PostData(0,fHistList);
244     PostData(1,fCFManagerPos->GetParticleContainer());
245     PostData(2,fCFManagerNeg->GetParticleContainer());
246     return;
247   }
248   
249   double primVtx[3];
250   vtx->GetXYZ(primVtx);
251   if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
252     PostData(0,fHistList);
253     PostData(1,fCFManagerPos->GetParticleContainer());
254     PostData(2,fCFManagerNeg->GetParticleContainer());
255     return;
256   }
257   
258   if(!fESD->GetNumberOfTracks() || fESD->GetNumberOfTracks()<2){ 
259     // Post output data
260     PostData(0,fHistList);
261     PostData(1,fCFManagerPos->GetParticleContainer());
262     PostData(2,fCFManagerNeg->GetParticleContainer());
263     return;
264   }
265   Int_t nTracks = fESD->GetNumberOfTracks();
266   AliDebug(2,Form("nTracks %d", nTracks));
267
268   if(!fTrackCuts) { 
269     // Post output data
270     PostData(0,fHistList);
271     PostData(1,fCFManagerPos->GetParticleContainer());
272     PostData(2,fCFManagerNeg->GetParticleContainer());
273     return;
274   }
275
276   // Selected events for analysis
277   fNEventSel->Fill(0.);
278   
279   
280   Double_t containerInputRec[5] ;
281   Double_t containerInputTPConly[5];
282   Double_t containerInputMC[5];
283   //Now go to rec level
284   for (Int_t iTrack = 0; iTrack<nTracks; iTrack++) 
285     {   
286       if(!fESD->GetTrack(iTrack) ) continue;
287       AliESDtrack* track = fESD->GetTrack(iTrack);
288       if(!(AliExternalTrackParam *)track->GetTPCInnerParam()) continue;
289       AliExternalTrackParam *trackTPC = (AliExternalTrackParam *)track->GetTPCInnerParam();
290       if(!track || !trackTPC) continue;
291
292       Float_t dca2D, dcaZ;
293       track->GetImpactParameters(dca2D,dcaZ);
294       Float_t dca2DTPC, dcaZTPC;
295       track->GetImpactParametersTPC(dca2DTPC,dcaZTPC); 
296       Float_t chi2PerClusterTPC = -1.;
297       Float_t nClustersTPC = track->GetTPCNcls();//track->GetTPCclusters(0);
298       if(nClustersTPC>0.) chi2PerClusterTPC = track->GetTPCchi2()/(2.*nClustersTPC-5.);
299       Float_t chi2PerClusterTPCIter1 = -1.;
300       Float_t nClustersTPCIter1 = track->GetTPCNclsIter1();   
301       if(nClustersTPCIter1>0.) chi2PerClusterTPCIter1 = track->GetTPCchi2Iter1()/(2.*nClustersTPCIter1-5.);
302
303       //fill the container
304       containerInputRec[0] = track->Pt();
305       containerInputRec[1] = track->Phi();
306       containerInputRec[2] = track->Eta();
307       containerInputRec[3] = dca2D;
308       containerInputRec[4] = chi2PerClusterTPC;
309
310       containerInputTPConly[0] = trackTPC->Pt();
311       containerInputTPConly[1] = trackTPC->Phi();
312       containerInputTPConly[2] = trackTPC->Eta();
313       containerInputTPConly[3] = dca2DTPC/10.; //Divide by 10 in order to store in same container. Should be corrected back when looking at output.
314       containerInputTPConly[4] = chi2PerClusterTPCIter1;//TPC;
315
316       if (fTrackCuts->AcceptTrack(track)) {
317         if(track->GetSign()>0.) {
318           fCFManagerPos->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
319           fCFManagerPos->GetParticleContainer()->Fill(containerInputTPConly,kStepReconstructedTPCOnly);
320         }
321         if(track->GetSign()<0.) {
322           fCFManagerNeg->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
323           fCFManagerNeg->GetParticleContainer()->Fill(containerInputTPConly,kStepReconstructedTPCOnly);
324         }
325         
326         
327         //Only fill the secondary particle container if MC information is available
328         if(eventHandler) {
329           Int_t label = TMath::Abs(track->GetLabel());
330           TParticle *particle = stack->Particle(label) ;
331           if(!particle) continue;
332
333           containerInputMC[0] = particle->Pt();      
334           containerInputMC[1] = particle->Phi();      
335           containerInputMC[2] = particle->Eta();  
336           containerInputMC[3] = 0.0;      
337           containerInputMC[4] = 0.0;  
338
339           if(particle->GetPDG()->Charge()>0.) {
340             fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepReconstructedMC);
341           }
342           if(particle->GetPDG()->Charge()<0.) {
343             fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepReconstructedMC);
344           }
345
346           if (!stack->IsPhysicalPrimary(label) ) {
347             if(particle->GetPDG()->Charge()>0.) {
348               fCFManagerPos->GetParticleContainer()->Fill(containerInputRec,kStepSecondaries);
349             }
350             if(particle->GetPDG()->Charge()<0.) {
351               fCFManagerNeg->GetParticleContainer()->Fill(containerInputRec,kStepSecondaries);
352             }
353           }
354         }
355         
356       }
357       
358     }
359   
360
361   //Fill MC containters if particles are findable
362   if(eventHandler) {
363     for(int iPart = 1; iPart<(mcEvent->GetNumberOfTracks()); iPart++)//stack->GetNprimary();
364       {
365         AliMCParticle *mcPart  = (AliMCParticle*)mcEvent->GetTrack(iPart);
366         if(!mcPart) continue;
367         
368         //fill the container
369         containerInputMC[0] = mcPart->Pt();
370         containerInputMC[1] = mcPart->Phi();      
371         containerInputMC[2] = mcPart->Eta();  
372         containerInputMC[3] = 0.0;
373         containerInputMC[4] = 0.0;
374
375         int counter;
376         Float_t trackLengthTPC = 0.;
377         if(mcPart->Charge()>0. && fCFManagerPos->CheckParticleCuts(3,mcPart)) {
378           fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance) ;
379           trackLengthTPC = mcPart->GetTPCTrackLength(fESD->GetMagneticField(),0.1,counter,3.0);
380           if(trackLengthTPC > (float)fTrackCuts->GetMinNClusterTPC()) fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepMCtrackable) ;
381         }
382         if(mcPart->Charge()<0. && fCFManagerNeg->CheckParticleCuts(3,mcPart)) {
383           fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance) ;
384           trackLengthTPC = mcPart->GetTPCTrackLength(fESD->GetMagneticField(),0.1,counter,3.0);
385           if(trackLengthTPC > (float)fTrackCuts->GetMinNClusterTPC()) fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepMCtrackable) ;
386         }
387       }
388   }
389   
390   PostData(0,fHistList);
391   PostData(1,fCFManagerPos->GetParticleContainer());
392   PostData(2,fCFManagerNeg->GetParticleContainer());
393   
394 }
395
396
397 //___________________________________________________________________________
398 void AliPWG4HighPtSpectra::Terminate(Option_t*)
399 {
400   // The Terminate() function is the last function to be called during
401   // a query. It always runs on the client, it can be used to present
402   // the results graphically or save the results to file.
403 }
404
405 //___________________________________________________________________________
406 void AliPWG4HighPtSpectra::CreateOutputObjects() {
407   //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
408   //TO BE SET BEFORE THE EXECUTION OF THE TASK
409   //
410   AliDebug(2,Form("CreateOutputObjects","CreateOutputObjects of task %s", GetName()));
411
412   Bool_t oldStatus = TH1::AddDirectoryStatus();
413   TH1::AddDirectory(kFALSE); 
414
415   //slot #1
416   OpenFile(0);
417   fHistList = new TList();
418   fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
419   fHistList->Add(fNEventAll);
420   fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
421   fHistList->Add(fNEventSel);
422
423   TH1::AddDirectory(oldStatus);   
424
425 }
426
427 #endif