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"
41 //__________________________________________________________________________
42 AliCFRsnTask::AliCFRsnTask() :
46 fHistEventsProcessed(0x0)
52 //___________________________________________________________________________
53 AliCFRsnTask::AliCFRsnTask(const Char_t* name) :
54 AliAnalysisTaskSE(name),
57 fHistEventsProcessed(0x0)
60 // Constructor. Initialization of Inputs and Outputs
62 Info("AliCFRsnTask","Calling Constructor");
64 DefineInput(0) and DefineOutput(0)
65 are taken care of by AliAnalysisTaskSE constructor
67 DefineOutput(1,TH1I::Class());
68 DefineOutput(2,AliCFContainer::Class());
69 // DefineOutput(3,TList::Class());
72 //___________________________________________________________________________
73 AliCFRsnTask& AliCFRsnTask::operator=(const AliCFRsnTask& c)
76 // Assignment operator
79 AliAnalysisTaskSE::operator=(c) ;
81 fCFManager = c.fCFManager;
82 fHistEventsProcessed = c.fHistEventsProcessed;
87 //___________________________________________________________________________
88 AliCFRsnTask::AliCFRsnTask(const AliCFRsnTask& c) :
91 fCFManager(c.fCFManager),
92 fHistEventsProcessed(c.fHistEventsProcessed)
99 //___________________________________________________________________________
100 AliCFRsnTask::~AliCFRsnTask() {
104 Info("~AliCFRsnTask","Calling Destructor");
105 if (fCFManager) delete fCFManager ;
106 if (fHistEventsProcessed) delete fHistEventsProcessed ;
109 //_________________________________________________
110 void AliCFRsnTask::UserExec(Option_t *)
113 // Main loop function
115 Info("UserExec","") ;
118 Error("UserExec","NO EVENT FOUND!");
122 if (!fMCEvent) Error("UserExec","NO MC INFO FOUND!");
123 fCFManager->SetEventInfo(fMCEvent);
125 AliStack* stack = fMCEvent->Stack();
127 Bool_t isESDEvent = strcmp(fInputEvent->ClassName(),"AliESDEvent") == 0 ? kTRUE : kFALSE ;
128 Bool_t isAODEvent = strcmp(fInputEvent->ClassName(),"AliAODEvent") == 0 ? kTRUE : kFALSE ;
130 // MC-event selection
131 Double_t containerInput[2] ;
133 //loop on the MC event
134 Info("UserExec","Looping on MC event");
135 for (Int_t ipart=0; ipart<stack->GetNprimary(); ipart++) {
136 AliMCParticle *mcPart = fMCEvent->GetTrack(ipart);
138 //check the MC-level cuts
139 if (!fCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,mcPart)) continue;
140 containerInput[0] = mcPart->Pt();
141 containerInput[1] = mcPart->Y() ;
142 //fill the container for Gen-level selection
143 fCFManager->GetParticleContainer()->Fill(containerInput,kStepGenerated);
145 //check the Acceptance-level cuts
146 if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mcPart)) continue;
147 //fill the container for Acceptance-level selection
148 fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructible);
152 //Now go to rec level
153 Info("UserExec","Looping on %s",fInputEvent->ClassName());
155 // Loop on negative tracks
156 for (Int_t iTrack1 = 0; iTrack1<fInputEvent->GetNumberOfTracks(); iTrack1++) {
157 AliVParticle* track1 = fInputEvent->GetTrack(iTrack1);
159 if (track1->Charge()>=0) continue;
160 Int_t label1 = track1->GetLabel();
161 if (label1<0) continue;
163 //Loop on positive tracks
164 for (Int_t iTrack2 = 0; iTrack2<fInputEvent->GetNumberOfTracks(); iTrack2++) {
165 AliVParticle* track2 = fInputEvent->GetTrack(iTrack2);
167 if (track2->Charge()<=0) continue;
168 Int_t label2 = track2->GetLabel();
169 if (label2<0) continue;
171 //Create Resonance daughter objects
172 AliRsnDaughter* daughter1tmp = 0x0 ;
173 AliRsnDaughter* daughter2tmp = 0x0 ;
175 daughter1tmp = new AliRsnDaughter((AliESDtrack*)track1) ;
176 daughter2tmp = new AliRsnDaughter((AliESDtrack*)track2) ;
178 else if (isAODEvent) {
179 daughter1tmp = new AliRsnDaughter((AliAODTrack*)track1) ;
180 daughter2tmp = new AliRsnDaughter((AliAODTrack*)track2) ;
183 Error("UserExec","Error: input data file is not an ESD nor an AOD");
187 AliRsnDaughter daughter1(*daughter1tmp);
188 AliRsnDaughter daughter2(*daughter2tmp);
192 AliCFPair pair(track1,track2); // This object is used for cuts
193 // (to be replaced when AliRsnPairParticle
194 // inherits from AliVParticle)
196 //Set MC PDG information to resonance daughters
197 TParticle *part1 = stack->Particle(label1);
198 daughter1.InitMCInfo(part1);
199 daughter1.GetMCInfo()->SetPDG(part1->GetPdgCode());
200 daughter1.SetM(part1->GetCalcMass()); // assign true mass
202 Int_t mother1 = part1->GetFirstMother();
203 daughter1.GetMCInfo()->SetMother(mother1);
205 TParticle *mum = stack->Particle(mother1);
206 daughter1.GetMCInfo()->SetMotherPDG(mum->GetPdgCode());
209 TParticle *part2 = stack->Particle(label2);
210 daughter2.InitMCInfo(part2);
211 daughter2.GetMCInfo()->SetPDG(part2->GetPdgCode());
212 daughter2.SetM(part2->GetCalcMass()); // assign true mass
214 Int_t mother2 = part2->GetFirstMother();
215 daughter2.GetMCInfo()->SetMother(mother2);
217 TParticle *mum = stack->Particle(mother2);
218 daughter2.GetMCInfo()->SetMotherPDG(mum->GetPdgCode());
221 //make a mother resonance from the 2 candidate daughters
222 AliRsnPairParticle rsn;
223 rsn.SetPair(&daughter1,&daughter2);
225 //check if true resonance
226 if (!rsn.IsTruePair(fRsnPDG)) continue;
227 if (!fCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,&pair)) continue;
229 //check if associated MC resonance passes the cuts
230 Int_t motherLabel = rsn.GetDaughter(0)->GetMCInfo()->Mother() ;
231 if (motherLabel<0) continue ;
233 AliMCParticle* mcRsn = fMCEvent->GetTrack(motherLabel);
234 if (!mcRsn) continue;
235 if (!fCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,mcRsn)) continue;
237 // fill the container
238 Double_t rsnEnergy = rsn.GetDaughter(0)->E() + rsn.GetDaughter(1)->E() ;
239 Double_t rsnPz = rsn.GetDaughter(0)->Pz() + rsn.GetDaughter(1)->Pz() ;
241 containerInput[0] = rsn.GetPt() ;
242 containerInput[1] = GetRapidity(rsnEnergy,rsnPz);
243 fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed) ;
245 if (!fCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,&pair)) continue ;
246 fCFManager->GetParticleContainer()->Fill(containerInput,kStepSelected);
250 fHistEventsProcessed->Fill(0);
251 /* PostData(0) is taken care of by AliAnalysisTaskSE */
252 PostData(1,fHistEventsProcessed) ;
253 PostData(2,fCFManager->GetParticleContainer()) ;
255 // TList * list = new TList();
256 // fCFManager->AddQAHistosToList(list);
257 // PostData(2,list) ;
261 //___________________________________________________________________________
262 void AliCFRsnTask::Terminate(Option_t*)
264 // The Terminate() function is the last function to be called during
265 // a query. It always runs on the client, it can be used to present
266 // the results graphically or save the results to file.
268 Info("Terminate","");
269 AliAnalysisTaskSE::Terminate();
272 //draw some example plots....
274 AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
276 TH1D* h00 = cont->ShowProjection(0,0) ;
277 TH1D* h01 = cont->ShowProjection(0,1) ;
278 TH1D* h02 = cont->ShowProjection(0,2) ;
279 TH1D* h03 = cont->ShowProjection(0,3) ;
281 TH1D* h10 = cont->ShowProjection(1,0) ;
282 TH1D* h11 = cont->ShowProjection(1,1) ;
283 TH1D* h12 = cont->ShowProjection(1,2) ;
284 TH1D* h13 = cont->ShowProjection(1,3) ;
286 Double_t max1 = h00->GetMaximum();
287 Double_t max2 = h10->GetMaximum();
289 h00->GetYaxis()->SetRangeUser(0,max1*1.2);
290 h01->GetYaxis()->SetRangeUser(0,max1*1.2);
291 h02->GetYaxis()->SetRangeUser(0,max1*1.2);
292 h03->GetYaxis()->SetRangeUser(0,max1*1.2);
294 h10->GetYaxis()->SetRangeUser(0,max2*1.2);
295 h11->GetYaxis()->SetRangeUser(0,max2*1.2);
296 h12->GetYaxis()->SetRangeUser(0,max2*1.2);
297 h13->GetYaxis()->SetRangeUser(0,max2*1.2);
299 h00->SetMarkerStyle(23) ;
300 h01->SetMarkerStyle(24) ;
301 h02->SetMarkerStyle(25) ;
302 h03->SetMarkerStyle(26) ;
304 h10->SetMarkerStyle(23) ;
305 h11->SetMarkerStyle(24) ;
306 h12->SetMarkerStyle(25) ;
307 h13->SetMarkerStyle(26) ;
309 TCanvas * c =new TCanvas("c","",1400,800);
329 c->SaveAs("plots.eps");
332 //___________________________________________________________________________
333 void AliCFRsnTask::UserCreateOutputObjects() {
334 //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
335 //TO BE SET BEFORE THE EXECUTION OF THE TASK
337 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
341 fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ;
344 //___________________________________________________________________________
345 Double_t AliCFRsnTask::GetRapidity(Double_t energy, Double_t pz) {
346 if (energy == pz || energy == -pz) {
347 printf("GetRapidity : ERROR : rapidity for 4-vector with E = Pz -- infinite result");
351 printf("GetRapidity : ERROR : rapidity for 4-vector with E = Pz -- infinite result");
354 Double_t y = 0.5 * log((energy + pz) / (energy - pz));