2 /**************************************************************************
3 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
5 * Author: The ALICE Off-line Project. *
6 * Contributors are mentioned in the code where appropriate. *
8 * Permission to use, copy, modify and distribute this software and its *
9 * documentation strictly for non-commercial purposes is hereby granted *
10 * without fee, provided that the above copyright notice appears in all *
11 * copies and that both the copyright notice and this permission notice *
12 * appear in the supporting documentation. The authors make no claims *
13 * about the suitability of this software for any purpose. It is *
14 * provided "as is" without express or implied warranty. *
15 **************************************************************************/
17 //-----------------------------------------------------------------------
18 // Author : R. Vernet, Consorzio Cometa - Catania (it)
19 //-----------------------------------------------------------------------
22 #include "AliCFRsnTask.h"
25 #include "TParticle.h"
28 #include "AliMCEvent.h"
29 #include "AliESDEvent.h"
30 #include "AliCFManager.h"
31 #include "AliCFCutBase.h"
32 #include "AliCFContainer.h"
33 #include "AliESDtrack.h"
35 #include "AliRsnDaughter.h"
36 #include "AliCFPair.h"
37 #include "AliRsnMCInfo.h"
38 #include "AliRsnPairParticle.h"
39 #include "AliAODEvent.h"
40 #include "AliRsnReader.h"
42 //__________________________________________________________________________
43 AliCFRsnTask::AliCFRsnTask() :
47 fHistEventsProcessed(0x0)
53 //___________________________________________________________________________
54 AliCFRsnTask::AliCFRsnTask(const Char_t* name) :
55 AliAnalysisTaskSE(name),
58 fHistEventsProcessed(0x0)
61 // Constructor. Initialization of Inputs and Outputs
63 Info("AliCFRsnTask","Calling Constructor");
65 DefineInput(0) and DefineOutput(0)
66 are taken care of by AliAnalysisTaskSE constructor
68 DefineOutput(1,TH1I::Class());
69 DefineOutput(2,AliCFContainer::Class());
70 // DefineOutput(3,TList::Class());
73 //___________________________________________________________________________
74 AliCFRsnTask& AliCFRsnTask::operator=(const AliCFRsnTask& c)
77 // Assignment operator
80 AliAnalysisTaskSE::operator=(c) ;
82 fCFManager = c.fCFManager;
83 fHistEventsProcessed = c.fHistEventsProcessed;
88 //___________________________________________________________________________
89 AliCFRsnTask::AliCFRsnTask(const AliCFRsnTask& c) :
92 fCFManager(c.fCFManager),
93 fHistEventsProcessed(c.fHistEventsProcessed)
100 //___________________________________________________________________________
101 AliCFRsnTask::~AliCFRsnTask() {
105 Info("~AliCFRsnTask","Calling Destructor");
106 if (fCFManager) delete fCFManager ;
107 if (fHistEventsProcessed) delete fHistEventsProcessed ;
110 //_________________________________________________
111 void AliCFRsnTask::UserExec(Option_t *)
114 // Main loop function
116 Info("UserExec","") ;
119 Error("UserExec","NO EVENT FOUND!");
123 if (!fMCEvent) Error("UserExec","NO MC INFO FOUND!");
124 fCFManager->SetEventInfo(fMCEvent);
126 AliStack* stack = fMCEvent->Stack();
128 Bool_t isESDEvent = strcmp(fInputEvent->ClassName(),"AliESDEvent") == 0 ? kTRUE : kFALSE ;
129 Bool_t isAODEvent = strcmp(fInputEvent->ClassName(),"AliAODEvent") == 0 ? kTRUE : kFALSE ;
131 // MC-event selection
132 Double_t containerInput[2] ;
134 //loop on the MC event
135 Info("UserExec","Looping on MC event");
136 for (Int_t ipart=0; ipart<stack->GetNprimary(); ipart++) {
137 AliMCParticle *mcPart = fMCEvent->GetTrack(ipart);
139 //check the MC-level cuts
140 if (!fCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,mcPart)) continue;
141 containerInput[0] = mcPart->Pt();
142 containerInput[1] = mcPart->Y() ;
143 //fill the container for Gen-level selection
144 fCFManager->GetParticleContainer()->Fill(containerInput,kStepGenerated);
146 //check the Acceptance-level cuts
147 if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mcPart)) continue;
148 //fill the container for Acceptance-level selection
149 fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructible);
153 //Now go to rec level
154 Info("UserExec","Looping on %s",fInputEvent->ClassName());
156 // Loop on negative tracks
157 for (Int_t iTrack1 = 0; iTrack1<fInputEvent->GetNumberOfTracks(); iTrack1++) {
158 AliVParticle* track1 = fInputEvent->GetTrack(iTrack1);
160 if (track1->Charge()>=0) continue;
161 Int_t label1 = track1->GetLabel();
162 if (label1<0) continue;
164 //Loop on positive tracks
165 for (Int_t iTrack2 = 0; iTrack2<fInputEvent->GetNumberOfTracks(); iTrack2++) {
166 AliVParticle* track2 = fInputEvent->GetTrack(iTrack2);
168 if (track2->Charge()<=0) continue;
169 Int_t label2 = track2->GetLabel();
170 if (label2<0) continue;
172 AliRsnDaughter daughter1;
173 AliRsnDaughter daughter2;
175 AliRsnReader reader ; //needed for conversion ESD/AOD->AliRsnDaughter
178 if (!reader.ConvertTrack(&daughter1,(AliESDtrack*)track1)) Error("UserExec","Could not convert ESD track") ;
179 if (!reader.ConvertTrack(&daughter2,(AliESDtrack*)track2)) Error("UserExec","Could not convert ESD track") ;
181 else if (isAODEvent) {
182 if (!reader.ConvertTrack(&daughter1,(AliAODTrack*)track1)) Error("UserExec","Could not convert AOD track") ;
183 if (!reader.ConvertTrack(&daughter2,(AliAODTrack*)track2)) Error("UserExec","Could not convert AOD track") ;
186 Error("UserExec","Error: input data file is not an ESD nor an AOD");
190 AliCFPair pair(track1,track2); // This object is used for cuts
191 // (to be replaced when AliRsnPairParticle
192 // inherits from AliVParticle)
194 //Set MC PDG information to resonance daughters
195 TParticle *part1 = stack->Particle(label1);
196 daughter1.InitMCInfo(part1);
197 daughter1.GetMCInfo()->SetPDG(part1->GetPdgCode());
198 daughter1.SetM(part1->GetCalcMass()); // assign true mass
200 Int_t mother1 = part1->GetFirstMother();
201 daughter1.GetMCInfo()->SetMother(mother1);
203 TParticle *mum = stack->Particle(mother1);
204 daughter1.GetMCInfo()->SetMotherPDG(mum->GetPdgCode());
207 TParticle *part2 = stack->Particle(label2);
208 daughter2.InitMCInfo(part2);
209 daughter2.GetMCInfo()->SetPDG(part2->GetPdgCode());
210 daughter2.SetM(part2->GetCalcMass()); // assign true mass
212 Int_t mother2 = part2->GetFirstMother();
213 daughter2.GetMCInfo()->SetMother(mother2);
215 TParticle *mum = stack->Particle(mother2);
216 daughter2.GetMCInfo()->SetMotherPDG(mum->GetPdgCode());
219 //make a mother resonance from the 2 candidate daughters
220 AliRsnPairParticle rsn;
221 rsn.SetPair(&daughter1,&daughter2);
223 //check if true resonance
224 if (!rsn.IsTruePair(fRsnPDG)) continue;
225 if (!fCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,&pair)) continue;
227 //check if associated MC resonance passes the cuts
228 Int_t motherLabel = rsn.GetDaughter(0)->GetMCInfo()->Mother() ;
229 if (motherLabel<0) continue ;
231 AliMCParticle* mcRsn = fMCEvent->GetTrack(motherLabel);
232 if (!mcRsn) continue;
233 if (!fCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,mcRsn)) continue;
235 // fill the container
236 Double_t rsnEnergy = rsn.GetDaughter(0)->E() + rsn.GetDaughter(1)->E() ;
237 Double_t rsnPz = rsn.GetDaughter(0)->Pz() + rsn.GetDaughter(1)->Pz() ;
239 containerInput[0] = rsn.GetPt() ;
240 containerInput[1] = GetRapidity(rsnEnergy,rsnPz);
241 fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed) ;
243 if (!fCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,&pair)) continue ;
244 fCFManager->GetParticleContainer()->Fill(containerInput,kStepSelected);
248 fHistEventsProcessed->Fill(0);
249 /* PostData(0) is taken care of by AliAnalysisTaskSE */
250 PostData(1,fHistEventsProcessed) ;
251 PostData(2,fCFManager->GetParticleContainer()) ;
253 // TList * list = new TList();
254 // fCFManager->AddQAHistosToList(list);
255 // PostData(2,list) ;
259 //___________________________________________________________________________
260 void AliCFRsnTask::Terminate(Option_t*)
262 // The Terminate() function is the last function to be called during
263 // a query. It always runs on the client, it can be used to present
264 // the results graphically or save the results to file.
266 Info("Terminate","");
267 AliAnalysisTaskSE::Terminate();
270 //draw some example plots....
272 AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
274 TH1D* h00 = cont->ShowProjection(0,0) ;
275 TH1D* h01 = cont->ShowProjection(0,1) ;
276 TH1D* h02 = cont->ShowProjection(0,2) ;
277 TH1D* h03 = cont->ShowProjection(0,3) ;
279 TH1D* h10 = cont->ShowProjection(1,0) ;
280 TH1D* h11 = cont->ShowProjection(1,1) ;
281 TH1D* h12 = cont->ShowProjection(1,2) ;
282 TH1D* h13 = cont->ShowProjection(1,3) ;
284 Double_t max1 = h00->GetMaximum();
285 Double_t max2 = h10->GetMaximum();
287 h00->GetYaxis()->SetRangeUser(0,max1*1.2);
288 h01->GetYaxis()->SetRangeUser(0,max1*1.2);
289 h02->GetYaxis()->SetRangeUser(0,max1*1.2);
290 h03->GetYaxis()->SetRangeUser(0,max1*1.2);
292 h10->GetYaxis()->SetRangeUser(0,max2*1.2);
293 h11->GetYaxis()->SetRangeUser(0,max2*1.2);
294 h12->GetYaxis()->SetRangeUser(0,max2*1.2);
295 h13->GetYaxis()->SetRangeUser(0,max2*1.2);
297 h00->SetMarkerStyle(23) ;
298 h01->SetMarkerStyle(24) ;
299 h02->SetMarkerStyle(25) ;
300 h03->SetMarkerStyle(26) ;
302 h10->SetMarkerStyle(23) ;
303 h11->SetMarkerStyle(24) ;
304 h12->SetMarkerStyle(25) ;
305 h13->SetMarkerStyle(26) ;
307 TCanvas * c =new TCanvas("c","",1400,800);
327 c->SaveAs("plots.eps");
330 //___________________________________________________________________________
331 void AliCFRsnTask::UserCreateOutputObjects() {
332 //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
333 //TO BE SET BEFORE THE EXECUTION OF THE TASK
335 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
339 fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ;
342 //___________________________________________________________________________
343 Double_t AliCFRsnTask::GetRapidity(Double_t energy, Double_t pz) {
344 if (energy == pz || energy == -pz) {
345 printf("GetRapidity : ERROR : rapidity for 4-vector with E = Pz -- infinite result");
349 printf("GetRapidity : ERROR : rapidity for 4-vector with E = Pz -- infinite result");
352 Double_t y = 0.5 * log((energy + pz) / (energy - pz));