new classes for resonance and V0 analysis (from R. Vernet)
[u/mrichter/AliRoot.git] / CORRFW / test / AliCFV0Task.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 on AliEn (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 : R. Vernet, Consorzio Cometa - Catania (it)
24 //-----------------------------------------------------------------------
25
26
27 #ifndef ALICFV0TASK_CXX
28 #define ALICFV0TASK_CXX
29 #include <TROOT.h>
30 #include <TInterpreter.h>
31
32 #include "AliCFV0Task.h"
33 #include "TCanvas.h"
34 #include "AliStack.h"
35 #include "TParticle.h"
36 #include "TH1I.h"
37 #include "TChain.h"
38 #include "AliMCEventHandler.h"
39 #include "AliMCEvent.h"
40 #include "AliAnalysisManager.h"
41 #include "AliESDEvent.h"
42 #include "AliCFManager.h"
43 #include "AliCFCutBase.h"
44 #include "AliCFContainer.h"
45 #include "TChain.h"
46 #include "AliESDtrack.h"
47 #include "AliLog.h"
48 #include "AliESDv0.h"
49 #include "AliV0vertexer.h"
50 #include "AliCFPair.h"
51
52 //__________________________________________________________________________
53 AliCFV0Task::AliCFV0Task() :
54   fRebuildV0s(0),
55   fV0PDG(310),
56   fChain(0x0),
57   fESD(0x0),
58   fCFManager(0x0),
59   fHistEventsProcessed(0x0)
60 {
61 //Defual ctor
62 }
63 //___________________________________________________________________________
64 AliCFV0Task::AliCFV0Task(const Char_t* name) :
65   AliAnalysisTask(name,"AliCFV0Task"),
66   fRebuildV0s(0),
67   fV0PDG(310),
68   fChain(0x0),
69   fESD(0x0),
70   fCFManager(0x0),
71   fHistEventsProcessed(0x0)
72 {
73   //
74   // Constructor. Initialization of Inputs and Outputs
75   //
76   Info("AliCFV0Task","Calling Constructor");
77   DefineInput (0,TChain::Class());
78   DefineOutput(0,TH1I::Class());
79   DefineOutput(1,AliCFContainer::Class());
80 //   DefineOutput(2,TList::Class());
81 }
82
83 //___________________________________________________________________________
84 AliCFV0Task& AliCFV0Task::operator=(const AliCFV0Task& c) 
85 {
86   //
87   // Assignment operator
88   //
89   if (this!=&c) {
90     AliAnalysisTask::operator=(c) ;
91     fRebuildV0s = c.fRebuildV0s;
92     fV0PDG      = c.fV0PDG;
93     fChain      = c.fChain;
94     fESD        = c.fESD;
95     fCFManager  = c.fCFManager;
96     fHistEventsProcessed = c.fHistEventsProcessed;
97   }
98   return *this;
99 }
100
101 //___________________________________________________________________________
102 AliCFV0Task::AliCFV0Task(const AliCFV0Task& c) :
103   AliAnalysisTask(c),
104   fRebuildV0s(c.fRebuildV0s),
105   fV0PDG(c.fV0PDG),
106   fChain(c.fChain),
107   fESD(c.fESD),
108   fCFManager(c.fCFManager),
109   fHistEventsProcessed(c.fHistEventsProcessed)
110 {
111   //
112   // Copy Constructor
113   //
114 }
115
116 //___________________________________________________________________________
117 AliCFV0Task::~AliCFV0Task() {
118   //
119   //destructor
120   //
121   Info("~AliCFV0Task","Calling Destructor");
122   if (fChain)               delete fChain ;
123   if (fESD)                 delete fESD ;
124   if (fCFManager)           delete fCFManager ;
125   if (fHistEventsProcessed) delete fHistEventsProcessed ;
126 }
127
128 //___________________________________________________________________________
129
130 void AliCFV0Task::Init()
131 {
132
133 }
134 //_________________________________________________
135 void AliCFV0Task::Exec(Option_t *)
136 {
137   //
138   // Main loop function
139   //
140   Info("Exec","") ;
141   // Get the mc truth
142   AliMCEventHandler* mcTruth = (AliMCEventHandler*)((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
143
144   if (!mcTruth) Error("Exec","NO MC INFO FOUND... EXITING");
145
146   // transform possible old AliESD into AliESDEvent
147
148   if (fESD->GetAliESDOld()) fESD->CopyFromOldESD(); //transition to new ESD format
149
150   //pass the MC evt handler to the cuts that need it 
151
152   fCFManager->SetEventInfo(mcTruth);
153
154   // Get the MC event 
155   AliMCEvent* mcEvent = mcTruth->MCEvent();
156
157   // MC-event selection
158   Double_t containerInput[2] ;
159         
160   //loop on the MC event
161   Info("Exec","Looping on MC event");
162   for (Int_t ipart=0; ipart<mcEvent->GetNumberOfTracks(); ipart++) { 
163     AliMCParticle *mcPart  = mcEvent->GetTrack(ipart);
164
165     //check the MC-level cuts
166     if (!fCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,mcPart)) continue;
167     containerInput[0] = mcPart->Pt();
168     containerInput[1] = mcPart->Eta() ;
169     //fill the container for Gen-level selection
170     fCFManager->GetParticleContainer()->Fill(containerInput,kStepGenerated);
171     
172     //check the Acceptance-level cuts
173     if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mcPart)) continue;
174     //fill the container for Acceptance-level selection
175     fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructible);
176   }    
177
178
179   //Now go to rec level
180   Info("Exec","Looping on ESD event");
181
182   if (fRebuildV0s) RebuildV0s() ;
183   printf("There are %d V0s in event\n",fESD->GetNumberOfV0s());
184   for (Int_t iV0 = 0; iV0<fESD->GetNumberOfV0s(); iV0++) {
185
186     AliESDv0* esdV0 = fESD->GetV0(iV0);
187
188     //check if mother reconstructed V0 can be associated to a MC V0
189     Int_t labMCV0 = IsMcV0(esdV0,fESD,mcEvent->Stack()) ;
190     if (labMCV0 <0) continue;
191
192     esdV0->ChangeMassHypothesis(fV0PDG); //important to do that before entering the cut check
193
194     AliCFPair pair(esdV0,fESD);
195     if (!fCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,&pair)) continue;
196
197     //check if associated MC v0 passes the cuts
198     AliMCParticle* mcV0 = mcEvent->GetTrack(labMCV0);
199     if (!mcV0) continue;
200
201     if (!fCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,mcV0)) continue; 
202     
203     //fill the container
204     Double_t mom[3];
205     esdV0->GetPxPyPz(mom[0],mom[1],mom[2]);
206     Double_t pt2 = mom[0]*mom[0]+mom[1]*mom[1] ;
207     Double_t pt  = TMath::Sqrt(pt2);
208     Double_t energy  = TMath::Sqrt(pt2 + mom[2]*mom[2] + TMath::Power(TDatabasePDG::Instance()->GetParticle(fV0PDG)->Mass(),2));
209
210     containerInput[0] = pt ;
211     containerInput[1] = GetRapidity(energy,mom[2]);
212     fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed) ;   
213     if (!fCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,&pair)) continue ;
214     fCFManager->GetParticleContainer()->Fill(containerInput,kStepSelected);
215   }
216   
217   fHistEventsProcessed->Fill(0);
218   PostData(0,fHistEventsProcessed) ;
219   PostData(1,fCFManager->GetParticleContainer()) ;
220   
221 //   TList * list = new TList();
222 //   fCFManager->AddQAHistosToList(list);
223 //   PostData(2,list) ;
224 }
225
226
227 //___________________________________________________________________________
228 void AliCFV0Task::Terminate(Option_t*)
229 {
230   // The Terminate() function is the last function to be called during
231   // a query. It always runs on the client, it can be used to present
232   // the results graphically or save the results to file.
233
234   Info("Terminate","");
235   AliAnalysisTask::Terminate();
236
237
238   Double_t max1 = fCFManager->GetParticleContainer()->ShowProjection(0,0)->GetMaximum();
239   Double_t max2 = fCFManager->GetParticleContainer()->ShowProjection(1,0)->GetMaximum();
240
241   fCFManager->GetParticleContainer()->ShowProjection(0,0)->GetYaxis()->SetRangeUser(0,max1*1.2);
242   fCFManager->GetParticleContainer()->ShowProjection(0,1)->GetYaxis()->SetRangeUser(0,max1*1.2);
243   fCFManager->GetParticleContainer()->ShowProjection(0,2)->GetYaxis()->SetRangeUser(0,max1*1.2);
244   fCFManager->GetParticleContainer()->ShowProjection(0,3)->GetYaxis()->SetRangeUser(0,max1*1.2);
245
246   fCFManager->GetParticleContainer()->ShowProjection(1,0)->GetYaxis()->SetRangeUser(0,max2*1.2);
247   fCFManager->GetParticleContainer()->ShowProjection(1,1)->GetYaxis()->SetRangeUser(0,max2*1.2);
248   fCFManager->GetParticleContainer()->ShowProjection(1,2)->GetYaxis()->SetRangeUser(0,max2*1.2);
249   fCFManager->GetParticleContainer()->ShowProjection(1,3)->GetYaxis()->SetRangeUser(0,max2*1.2);
250
251   fCFManager->GetParticleContainer()->ShowProjection(0,0)->SetMarkerStyle(23) ;
252   fCFManager->GetParticleContainer()->ShowProjection(0,1)->SetMarkerStyle(24) ;
253   fCFManager->GetParticleContainer()->ShowProjection(0,2)->SetMarkerStyle(25) ;
254   fCFManager->GetParticleContainer()->ShowProjection(0,3)->SetMarkerStyle(26) ;
255
256   fCFManager->GetParticleContainer()->ShowProjection(1,0)->SetMarkerStyle(23) ;
257   fCFManager->GetParticleContainer()->ShowProjection(1,1)->SetMarkerStyle(24) ;
258   fCFManager->GetParticleContainer()->ShowProjection(1,2)->SetMarkerStyle(25) ;
259   fCFManager->GetParticleContainer()->ShowProjection(1,3)->SetMarkerStyle(26) ;
260
261   TCanvas * c =new TCanvas("c","",1400,800);
262   c->Divide(4,2);
263
264   c->cd(1);
265   fCFManager->GetParticleContainer()->ShowProjection(0,0)->Draw("p");
266   c->cd(2);
267   fCFManager->GetParticleContainer()->ShowProjection(0,1)->Draw("p");
268   c->cd(3);
269   fCFManager->GetParticleContainer()->ShowProjection(0,2)->Draw("p");
270   c->cd(4);
271   fCFManager->GetParticleContainer()->ShowProjection(0,3)->Draw("p");
272   c->cd(5);
273   fCFManager->GetParticleContainer()->ShowProjection(1,0)->Draw("p");
274   c->cd(6);
275   fCFManager->GetParticleContainer()->ShowProjection(1,1)->Draw("p");
276   c->cd(7);
277   fCFManager->GetParticleContainer()->ShowProjection(1,2)->Draw("p");
278   c->cd(8);
279   fCFManager->GetParticleContainer()->ShowProjection(1,3)->Draw("p");
280
281   c->SaveAs("plots.eps");
282
283   delete fHistEventsProcessed ;
284 }
285
286 //___________________________________________________________________________
287 void AliCFV0Task::ConnectInputData(Option_t *) {
288   //
289   // Initialize branches.
290   //
291   Info("ConnectInputData","ConnectInputData of task %s\n",GetName());
292
293   fChain = (TChain*)GetInputData(0);
294   fChain->SetBranchStatus("*FMD*",0);
295   fChain->SetBranchStatus("*CaloClusters*",0);
296   fESD = new AliESDEvent();
297   fESD->ReadFromTree(fChain);
298 }
299
300 //___________________________________________________________________________
301 void AliCFV0Task::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   Info("CreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
306
307   //slot #0
308   fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ;
309 }
310
311 //___________________________________________________________________________
312 Int_t AliCFV0Task::GetV0Label(UInt_t lab1, UInt_t lab2, AliStack* stack) const {
313   //
314   // returns the label of the V0, given the labels of the 2 daughter tracks
315   // returns -1 if the V0 is fake
316   //
317   
318   TParticle* part1 = stack->Particle(lab1) ;
319   TParticle* part2 = stack->Particle(lab2) ;
320   
321   Int_t part1MotherLab=part1->GetFirstMother();
322   Int_t part2MotherLab=part2->GetFirstMother();
323   
324   if (part1MotherLab==-1 || part2MotherLab==-1) return -1 ;
325   if (part1MotherLab != part2MotherLab )        return -1 ;
326   if (stack->Particle(part1MotherLab)->GetPdgCode() != fV0PDG ) return -1 ;
327
328   switch (fV0PDG) {
329   case kK0Short : 
330     if ( (part1->GetPdgCode()==  -211 && part2->GetPdgCode()==  211)  || 
331          (part1->GetPdgCode()==   211 && part2->GetPdgCode()== -211)   )   return part1MotherLab ;
332     break ;
333   case kLambda0 :
334     if ( (part1->GetPdgCode()==  -211 && part2->GetPdgCode()== 2212)  || 
335          (part1->GetPdgCode()==  2212 && part2->GetPdgCode()== -211)   )   return part1MotherLab ;
336     break ;
337   case kLambda0Bar :
338     if ( (part1->GetPdgCode()==   211 && part2->GetPdgCode()== -2212) || 
339          (part1->GetPdgCode()== -2212 && part2->GetPdgCode()==   211)  )   return part1MotherLab ;
340     break ;
341   default :
342     return -1;
343     break ;
344   }
345   
346   return -1 ;
347 }
348
349 //___________________________________________________________________________
350 Double_t AliCFV0Task::GetRapidity(Double_t energy, Double_t pz) {
351   if (energy == pz || energy == -pz) {
352     printf("GetRapidity : ERROR : rapidity for 4-vector with E = Pz -- infinite result");
353     return 999;
354   }
355   if (energy < pz) {
356     printf("GetRapidity : ERROR : rapidity for 4-vector with E = Pz -- infinite result");
357     return 999;
358   }
359   Double_t y = 0.5 * log((energy + pz) / (energy - pz));
360   return y;
361
362
363 //___________________________________________________________________________
364 Int_t AliCFV0Task::IsMcV0(AliESDv0* v0, AliESDEvent* esd, AliStack* stack) const {
365   Int_t nindex = v0->GetNindex();
366   Int_t pindex = v0->GetPindex();
367   AliESDtrack *nTrack = esd->GetTrack(nindex) ;
368   AliESDtrack *pTrack = esd->GetTrack(pindex) ;
369   
370   if (!nTrack || !pTrack) return -1 ;
371
372   Int_t nlab  = nTrack->GetLabel() ;
373   Int_t plab  = pTrack->GetLabel() ;
374   
375   if (nlab <0 || plab <0) return -1 ;
376
377   Int_t v0Label = GetV0Label((UInt_t)nlab,(UInt_t)plab,stack) ;
378   return v0Label ;
379 }
380
381
382 //___________________________________________________________________________
383 void AliCFV0Task::RebuildV0s() {
384   
385   fESD->ResetV0s();
386
387   //These are pp cuts : to change if Pb+Pb !
388   Double_t cuts[]={33,  // max. allowed chi2
389                    0.1,// min. allowed negative daughter's impact parameter 
390                    0.1,// min. allowed positive daughter's impact parameter 
391                    0.1,// max. allowed DCA between the daughter tracks
392                    0.999,// max. allowed cosine of V0's pointing angle
393                    0.9,  // min. radius of the fiducial volume
394                    100.   // max. radius of the fiducial volume
395   };
396   //   AliV0vertexer* v0Vertexer = new AliV0vertexer(cuts);
397   //   v0Vertexer->SetVertex(primVertexPos);
398   AliV0vertexer* v0Vertexer = new AliV0vertexer();
399   v0Vertexer->SetCuts(cuts) ;
400   v0Vertexer->Tracks2V0vertices(fESD);
401   delete v0Vertexer ;
402 }
403
404 #endif