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