new classes for resonance and V0 analysis (from R. Vernet)
[u/mrichter/AliRoot.git] / CORRFW / test / AliCFRsnTask.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 // Author : R. Vernet, Consorzio Cometa - Catania (it)
18 //-----------------------------------------------------------------------
19
20
21 #include <TROOT.h>
22 #include <TInterpreter.h>
23
24 #include "AliCFRsnTask.h"
25 #include "TCanvas.h"
26 #include "AliStack.h"
27 #include "TParticle.h"
28 #include "TH1I.h"
29 #include "TChain.h"
30 #include "AliMCEventHandler.h"
31 #include "AliMCEvent.h"
32 #include "AliAnalysisManager.h"
33 #include "AliESDEvent.h"
34 #include "AliCFManager.h"
35 #include "AliCFCutBase.h"
36 #include "AliCFContainer.h"
37 #include "TChain.h"
38 #include "AliESDtrack.h"
39 #include "AliLog.h"
40 #include "AliRsnDaughter.h"
41 #include "AliCFPair.h"
42 #include <memory>
43
44 //__________________________________________________________________________
45 AliCFRsnTask::AliCFRsnTask() :
46   fRsnPDG(313),
47   fChain(0x0),
48   fESD(0x0),
49   fCFManager(0x0),
50   fHistEventsProcessed(0x0)
51 {
52   //Defual ctor
53 }
54 //___________________________________________________________________________
55 AliCFRsnTask::AliCFRsnTask(const Char_t* name) :
56   AliAnalysisTask(name,"AliCFRsnTask"),
57   fRsnPDG(313),
58   fChain(0x0),
59   fESD(0x0),
60   fCFManager(0x0),
61   fHistEventsProcessed(0x0)
62 {
63   //
64   // Constructor. Initialization of Inputs and Outputs
65   //
66   Info("AliCFRsnTask","Calling Constructor");
67   DefineInput (0,TChain::Class());
68   DefineOutput(0,TH1I::Class());
69   DefineOutput(1,AliCFContainer::Class());
70   //   DefineOutput(2,TList::Class());
71 }
72
73 //___________________________________________________________________________
74 AliCFRsnTask& AliCFRsnTask::operator=(const AliCFRsnTask& c) 
75 {
76   //
77   // Assignment operator
78   //
79   if (this!=&c) {
80     AliAnalysisTask::operator=(c) ;
81     fRsnPDG     = c.fRsnPDG;
82     fChain      = c.fChain;
83     fESD        = c.fESD;
84     fCFManager  = c.fCFManager;
85     fHistEventsProcessed = c.fHistEventsProcessed;
86   }
87   return *this;
88 }
89
90 //___________________________________________________________________________
91 AliCFRsnTask::AliCFRsnTask(const AliCFRsnTask& c) :
92   AliAnalysisTask(c),
93   fRsnPDG(c.fRsnPDG),
94   fChain(c.fChain),
95   fESD(c.fESD),
96   fCFManager(c.fCFManager),
97   fHistEventsProcessed(c.fHistEventsProcessed)
98 {
99   //
100   // Copy Constructor
101   //
102 }
103
104 //___________________________________________________________________________
105 AliCFRsnTask::~AliCFRsnTask() {
106   //
107   //destructor
108   //
109   Info("~AliCFRsnTask","Calling Destructor");
110   if (fChain)               delete fChain ;
111   if (fESD)                 delete fESD ;
112   if (fCFManager)           delete fCFManager ;
113   if (fHistEventsProcessed) delete fHistEventsProcessed ;
114 }
115
116 //___________________________________________________________________________
117
118 void AliCFRsnTask::Init()
119 {
120
121 }
122 //_________________________________________________
123 void AliCFRsnTask::Exec(Option_t *)
124 {
125   //
126   // Main loop function
127   //
128   Info("Exec","") ;
129   // Get the mc truth
130   AliMCEventHandler* mcTruth = (AliMCEventHandler*)((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
131
132   if (!mcTruth) Error("Exec","NO MC INFO FOUND... EXITING");
133
134   // transform possible old AliESD into AliESDEvent
135
136   if (fESD->GetAliESDOld()) fESD->CopyFromOldESD(); //transition to new ESD format
137
138   //pass the MC evt handler to the cuts that need it 
139
140   fCFManager->SetEventInfo(mcTruth);
141
142   // Get the MC event 
143   AliMCEvent* mcEvent = mcTruth->MCEvent();
144   AliStack*   stack   = mcEvent->Stack();
145
146   // MC-event selection
147   Double_t containerInput[2] ;
148         
149   //loop on the MC event
150   Info("Exec","Looping on MC event");
151   for (Int_t ipart=0; ipart<stack->GetNprimary(); ipart++) { 
152     AliMCParticle *mcPart  = mcEvent->GetTrack(ipart);
153
154     //check the MC-level cuts
155     if (!fCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,mcPart)) continue;
156     containerInput[0] = mcPart->Pt();
157     containerInput[1] = mcPart->Eta() ;
158     //fill the container for Gen-level selection
159     fCFManager->GetParticleContainer()->Fill(containerInput,kStepGenerated);
160     
161     //check the Acceptance-level cuts
162     if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mcPart)) continue;
163     //fill the container for Acceptance-level selection
164     fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructible);
165   }    
166
167
168   //Now go to rec level
169   Info("Exec","Looping on ESD event");
170
171   //SET THE ESD AS EVENT INFO IN RECONSTRUCTION CUTS
172   TObjArray* fCutsReco = fCFManager->GetParticleCutsList(AliCFManager::kPartRecCuts);
173   TObjArray* fCutsPID = fCFManager->GetParticleCutsList(AliCFManager::kPartSelCuts);
174   TObjArrayIter iter1(fCutsReco);
175   TObjArrayIter iter2(fCutsPID);
176   AliCFCutBase *cut = 0;
177   while ( (cut = (AliCFCutBase*)iter1.Next()) ) {
178     cut->SetEvtInfo(fESD);
179   }
180   while ( (cut = (AliCFCutBase*)iter2.Next()) ) {
181     cut->SetEvtInfo(fESD);
182   }
183   
184   for (Int_t iTrack1 = 0; iTrack1<fESD->GetNumberOfTracks(); iTrack1++) {
185     AliESDtrack* esdTrack1 = fESD->GetTrack(iTrack1);
186     //track1 is negative
187     if (esdTrack1->Charge()>=0) continue;
188     Int_t esdLabel1 = esdTrack1->GetLabel();
189     if (esdLabel1<0) continue;
190
191     for (Int_t iTrack2 = 0; iTrack2<fESD->GetNumberOfTracks(); iTrack2++) {
192       AliESDtrack* esdTrack2 = fESD->GetTrack(iTrack2);
193       //track2 is positive
194       if (esdTrack2->Charge()<=0) continue;
195       Int_t esdLabel2 = esdTrack2->GetLabel();
196       if (esdLabel2<0) continue;
197         
198       // copy ESDtrack data into RsnDaughter (and make Bayesian PID)
199       //if problem with copy constructor, use the following special command to destroy the pointer when exiting scope
200       //std::auto_ptr<AliRsnDaughter> track1(AliRsnDaughter::Adopt(esdTrack1,iTrack1));
201
202       AliRsnDaughter* tmp1=AliRsnDaughter::Adopt(esdTrack1,iTrack1);
203       AliRsnDaughter* tmp2=AliRsnDaughter::Adopt(esdTrack2,iTrack2);
204       AliRsnDaughter track1(*tmp1);
205       AliRsnDaughter track2(*tmp2);
206       delete tmp1;
207       delete tmp2;
208
209       TParticle *part1 = stack->Particle(esdLabel1);
210       track1.SetTruePDG(part1->GetPdgCode());
211       Int_t mother1 = part1->GetFirstMother();
212       track1.SetMother(mother1);
213       if (mother1 >= 0) {
214         TParticle *mum = stack->Particle(mother1);
215         track1.SetMotherPDG(mum->GetPdgCode());
216       }
217       
218       TParticle *part2 = stack->Particle(esdLabel2);
219       track2.SetTruePDG(part2->GetPdgCode());
220       Int_t mother2 = part2->GetFirstMother();
221       track2.SetMother(mother2);
222       if (mother2 >= 0) {
223         TParticle *mum = stack->Particle(mother2);
224         track2.SetMotherPDG(mum->GetPdgCode());
225       }
226         
227       //make a mother resonance from the 2 candidate daughters
228       AliRsnDaughter rsn = AliRsnDaughter::Sum(track1,track2) ;
229       AliCFPair pair(esdTrack1,esdTrack2);
230
231       //check if true resonance
232       if (rsn.GetMotherPDG() != fRsnPDG) continue;
233       if (!fCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,&pair)) continue;
234
235       //check if associated MC resonance passes the cuts
236       Int_t motherLabel=rsn.GetLabel();
237       if (motherLabel<0) continue ;
238       AliMCParticle* mcRsn = mcEvent->GetTrack(motherLabel);
239       if (!mcRsn) continue;
240       if (!fCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,mcRsn)) continue; 
241     
242       //fill the container
243       containerInput[0] = rsn.GetPt() ;
244       containerInput[1] = GetRapidity(rsn.GetEnergy(),rsn.GetPz());
245       fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed) ;   
246
247       if (!fCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,&pair)) continue ;
248       fCFManager->GetParticleContainer()->Fill(containerInput,kStepSelected);
249     }
250   }
251     
252   fHistEventsProcessed->Fill(0);
253   PostData(0,fHistEventsProcessed) ;
254   PostData(1,fCFManager->GetParticleContainer()) ;
255   
256   //   TList * list = new TList();
257   //   fCFManager->AddQAHistosToList(list);
258   //   PostData(2,list) ;
259 }
260
261
262 //___________________________________________________________________________
263 void AliCFRsnTask::Terminate(Option_t*)
264 {
265   // The Terminate() function is the last function to be called during
266   // a query. It always runs on the client, it can be used to present
267   // the results graphically or save the results to file.
268
269   Info("Terminate","");
270   AliAnalysisTask::Terminate();
271
272
273   Double_t max1 = fCFManager->GetParticleContainer()->ShowProjection(0,0)->GetMaximum();
274   Double_t max2 = fCFManager->GetParticleContainer()->ShowProjection(1,0)->GetMaximum();
275
276   fCFManager->GetParticleContainer()->ShowProjection(0,0)->GetYaxis()->SetRangeUser(0,max1*1.2);
277   fCFManager->GetParticleContainer()->ShowProjection(0,1)->GetYaxis()->SetRangeUser(0,max1*1.2);
278   fCFManager->GetParticleContainer()->ShowProjection(0,2)->GetYaxis()->SetRangeUser(0,max1*1.2);
279   fCFManager->GetParticleContainer()->ShowProjection(0,3)->GetYaxis()->SetRangeUser(0,max1*1.2);
280
281   fCFManager->GetParticleContainer()->ShowProjection(1,0)->GetYaxis()->SetRangeUser(0,max2*1.2);
282   fCFManager->GetParticleContainer()->ShowProjection(1,1)->GetYaxis()->SetRangeUser(0,max2*1.2);
283   fCFManager->GetParticleContainer()->ShowProjection(1,2)->GetYaxis()->SetRangeUser(0,max2*1.2);
284   fCFManager->GetParticleContainer()->ShowProjection(1,3)->GetYaxis()->SetRangeUser(0,max2*1.2);
285
286   fCFManager->GetParticleContainer()->ShowProjection(0,0)->SetMarkerStyle(23) ;
287   fCFManager->GetParticleContainer()->ShowProjection(0,1)->SetMarkerStyle(24) ;
288   fCFManager->GetParticleContainer()->ShowProjection(0,2)->SetMarkerStyle(25) ;
289   fCFManager->GetParticleContainer()->ShowProjection(0,3)->SetMarkerStyle(26) ;
290
291   fCFManager->GetParticleContainer()->ShowProjection(1,0)->SetMarkerStyle(23) ;
292   fCFManager->GetParticleContainer()->ShowProjection(1,1)->SetMarkerStyle(24) ;
293   fCFManager->GetParticleContainer()->ShowProjection(1,2)->SetMarkerStyle(25) ;
294   fCFManager->GetParticleContainer()->ShowProjection(1,3)->SetMarkerStyle(26) ;
295
296   TCanvas * c =new TCanvas("c","",1400,800);
297   c->Divide(4,2);
298
299   c->cd(1);
300   fCFManager->GetParticleContainer()->ShowProjection(0,0)->Draw("p");
301   c->cd(2);
302   fCFManager->GetParticleContainer()->ShowProjection(0,1)->Draw("p");
303   c->cd(3);
304   fCFManager->GetParticleContainer()->ShowProjection(0,2)->Draw("p");
305   c->cd(4);
306   fCFManager->GetParticleContainer()->ShowProjection(0,3)->Draw("p");
307   c->cd(5);
308   fCFManager->GetParticleContainer()->ShowProjection(1,0)->Draw("p");
309   c->cd(6);
310   fCFManager->GetParticleContainer()->ShowProjection(1,1)->Draw("p");
311   c->cd(7);
312   fCFManager->GetParticleContainer()->ShowProjection(1,2)->Draw("p");
313   c->cd(8);
314   fCFManager->GetParticleContainer()->ShowProjection(1,3)->Draw("p");
315
316   c->SaveAs("plots.eps");
317
318   delete fHistEventsProcessed ;
319 }
320
321 //___________________________________________________________________________
322 void AliCFRsnTask::ConnectInputData(Option_t *) {
323   //
324   // Initialize branches.
325   //
326   Info("ConnectInputData","ConnectInputData of task %s\n",GetName());
327
328   fChain = (TChain*)GetInputData(0);
329   fChain->SetBranchStatus("*FMD*",0);
330   fChain->SetBranchStatus("*CaloClusters*",0);
331   fESD = new AliESDEvent();
332   fESD->ReadFromTree(fChain);
333 }
334
335 //___________________________________________________________________________
336 void AliCFRsnTask::CreateOutputObjects() {
337   //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
338   //TO BE SET BEFORE THE EXECUTION OF THE TASK
339   //
340   Info("CreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
341
342   //slot #0
343   fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ;
344 }
345
346 //___________________________________________________________________________
347 Double_t AliCFRsnTask::GetRapidity(Double_t energy, Double_t pz) {
348   if (energy == pz || energy == -pz) {
349     printf("GetRapidity : ERROR : rapidity for 4-vector with E = Pz -- infinite result");
350     return 999;
351   }
352   if (energy < pz) {
353     printf("GetRapidity : ERROR : rapidity for 4-vector with E = Pz -- infinite result");
354     return 999;
355   }
356   Double_t y = 0.5 * log((energy + pz) / (energy - pz));
357   return y;
358
359
360