1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 /* $Id: AliMCAnalysisUtils.cxx 21839 2007-10-29 13:49:42Z gustavo $ */
17 //_________________________________________________________________________
18 // Class for analysis utils for MC data
19 // stored in stack or event header.
21 // - method to check the origin of a given track/cluster
22 // - method to obtain the generated jets
24 //*-- Author: Gustavo Conesa (LNF-INFN)
25 //////////////////////////////////////////////////////////////////////////////
28 // --- ROOT system ---
33 //---- ANALYSIS system ----
35 #include "AliMCAnalysisUtils.h"
37 #include "TParticle.h"
38 #include "AliGenPythiaEventHeader.h"
40 ClassImp(AliMCAnalysisUtils)
43 //________________________________________________
44 AliMCAnalysisUtils::AliMCAnalysisUtils() :
45 TObject(), fCurrentEvent(-1), fDebug(-1),
46 fJetsList(new TList), fMCGenerator("PYTHIA")
51 //____________________________________________________________________________
52 AliMCAnalysisUtils::AliMCAnalysisUtils(const AliMCAnalysisUtils & mcutils) :
53 TObject(mcutils), fCurrentEvent(mcutils.fCurrentEvent), fDebug(mcutils.fDebug),
54 fJetsList(mcutils.fJetsList), fMCGenerator(mcutils.fMCGenerator)
60 //_________________________________________________________________________
61 AliMCAnalysisUtils & AliMCAnalysisUtils::operator = (const AliMCAnalysisUtils & mcutils)
63 // assignment operator
65 if(&mcutils == this) return *this;
66 fCurrentEvent = mcutils.fCurrentEvent ;
67 fDebug = mcutils.fDebug;
68 fJetsList = mcutils.fJetsList;
69 fMCGenerator = mcutils.fMCGenerator;
75 //____________________________________________________________________________
76 AliMCAnalysisUtils::~AliMCAnalysisUtils()
78 // Remove all pointers.
87 //_________________________________________________________________________
88 Int_t AliMCAnalysisUtils::CheckOrigin(const Int_t label, AliStack * stack) const {
89 //Play with the MC stack if available
90 //Check origin of the candidates, good for PYTHIA
92 if(!stack) AliFatal("Stack is not available, check analysis settings in configuration file, STOP!!");
93 // printf("label %d, ntrack %d, nprim %d\n",label, stack->GetNtrack(), stack->GetNprimary());
94 // for(Int_t i = 0; i< stack->GetNprimary(); i++){
95 // TParticle *particle = stack->Particle(i);
96 // //particle->Print();
98 if(label >= 0 && label < stack->GetNtrack()){
100 TParticle * mom = stack->Particle(label);
101 Int_t mPdg = TMath::Abs(mom->GetPdgCode());
102 Int_t mStatus = mom->GetStatusCode() ;
103 Int_t iParent = mom->GetFirstMother() ;
104 if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOrigin: Mother is parton %d\n",iParent);
107 TParticle * parent = new TParticle ;
111 parent = stack->Particle(iParent);
112 pPdg = TMath::Abs(parent->GetPdgCode());
113 pStatus = parent->GetStatusCode();
115 else if(fDebug > 0 ) printf("AliMCAnalysisUtils::CheckOrigin: Parent with label %d\n",iParent);
120 if(fMCGenerator == "PYTHIA"){
121 if(iParent < 8 && iParent > 5) {//outgoing partons
122 if(pPdg == 22) return kMCPrompt;
123 else return kMCFragmentation;
125 else if(pStatus == 11){//Decay
126 if(pPdg == 111) return kMCPi0Decay ;
127 else if (pPdg == 321) return kMCEtaDecay ;
128 else return kMCOtherDecay ;
130 else return kMCISR; //Initial state radiation
133 else if(fMCGenerator == "HERWIG"){
134 if(pStatus < 197){//Not decay
136 if(parent->GetFirstMother()<=5) break;
137 iParent = parent->GetFirstMother();
138 parent=stack->Particle(iParent);
139 pStatus= parent->GetStatusCode();
140 pPdg = parent->GetPdgCode();
141 }//Look for the parton
143 if(iParent < 8 && iParent > 5) {
144 if(pPdg == 22) return kMCPrompt;
145 else return kMCFragmentation;
147 return kMCISR;//Initial state radiation
150 if(pPdg == 111) return kMCPi0Decay ;
151 else if (pPdg == 321) return kMCEtaDecay ;
152 else return kMCOtherDecay ;
155 else return kMCUnknown;
156 }//Status 1 : Pythia generated
157 else if(mStatus == 0){
158 if(pPdg ==22 || pPdg ==11|| pPdg == 2112 || pPdg == 211 ||
159 pPdg == 321 || pPdg == 2212 || pPdg == 130 || pPdg == 13 )
160 return kMCConversion ;
161 if(pPdg == 111) return kMCPi0Decay ;
162 else if (pPdg == 221) return kMCEtaDecay ;
163 else return kMCOtherDecay ;
164 }//status 0 : geant generated
166 else if(mPdg == 111) return kMCPi0 ;
167 else if(mPdg == 221) return kMCEta ;
169 // if(pPdg !=22 &&mStatus == 0) {
170 // printf("Origin electron, pT %f, status %d, parent %s, pstatus %d, vx %f, vy %f, vz %f\n",
171 // mom->Pt(),mStatus, parent->GetName(),pStatus,mom->Vx(),mom->Vy(), mom->Vz());
176 if(pPdg ==22) return kMCConversion ;
177 else if (pStatus == 0) return kMCConversion;
178 else return kMCElectron ;
180 else return kMCElectron ;
182 else return kMCUnknown;
185 if(label < 0 ) printf("AliMCAnalysisUtils::CheckOrigin: *** bad label or no stack ***: label %d \n", label);
186 if(label >= stack->GetNtrack()) printf("AliMCAnalysisUtils::CheckOrigin: *** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack());
194 //_________________________________________________________________________
195 TList * AliMCAnalysisUtils::GetJets(Int_t iEvent, AliStack * stack, AliGenEventHeader * geh) {
196 //Return list of jets (TParticles) and index of most likely parton that originated it.
198 if(fCurrentEvent!=iEvent){
199 fCurrentEvent = iEvent;
200 fJetsList = new TList;
201 Int_t nTriggerJets = 0;
202 Float_t tmpjet[]={0,0,0,0};
204 //printf("Event %d %d\n",fCurrentEvent,iEvent);
205 //Get outgoing partons
206 if(stack->GetNtrack() < 8) return fJetsList;
207 TParticle * parton1 = stack->Particle(6);
208 TParticle * parton2 = stack->Particle(7);
210 printf("parton 6 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
211 parton1->GetName(),parton1->Pt(),parton1->Energy(),parton1->Phi()*TMath::RadToDeg(),parton1->Eta());
212 printf("parton 7 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
213 parton2->GetName(),parton2->Pt(),parton2->Energy(),parton2->Phi()*TMath::RadToDeg(),parton2->Eta());
215 // //Trace the jet from the mother parton
222 // TParticle * tmptmp = new TParticle;
223 // for(Int_t i = 0; i< stack->GetNprimary(); i++){
224 // tmptmp = stack->Particle(i);
226 // if(tmptmp->GetStatusCode() == 1){
227 // pt = tmptmp->Pt();
228 // e = tmptmp->Energy();
229 // Int_t imom = tmptmp->GetFirstMother();
231 // //printf("1st imom %d\n",imom);
234 // tmptmp = stack->Particle(imom);
235 // imom = tmptmp->GetFirstMother();
236 // //printf("imom %d \n",imom);
238 // //printf("Last imom %d %d\n",imom1, imom);
243 // else if (imom1 == 7){
250 // printf("JET 1, pt %2.2f, e %2.2f; JET 2, pt %2.2f, e %2.2f \n",pt1,e1,pt2,e2);
252 //Get the jet, different way for different generator
254 if(fMCGenerator == "PYTHIA"){
255 TParticle * jet = new TParticle;
256 AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) geh;
257 nTriggerJets = pygeh->NTriggerJets();
259 printf("PythiaEventHeader: Njets: %d\n",nTriggerJets);
261 for(Int_t i = 0; i< nTriggerJets; i++){
263 pygeh->TriggerJet(i, tmpjet);
264 jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
265 //Assign an outgoing parton as mother
266 Float_t phidiff1 = TMath::Abs(jet->Phi()-parton1->Phi());
267 Float_t phidiff2 = TMath::Abs(jet->Phi()-parton2->Phi());
268 if(phidiff1 > phidiff2) jet->SetFirstMother(7);
269 else jet->SetFirstMother(6);
272 printf("PYTHIA Jet %d: mother %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
273 i, jet->GetFirstMother(),jet->Pt(),jet->Energy(),jet->Phi()*TMath::RadToDeg(),jet->Eta());
276 }//Pythia triggered jets
278 else if (fMCGenerator=="HERWIG"){
281 TParticle * tmp = parton1;
282 if(parton1->GetPdgCode()!=22){
284 if(tmp->GetFirstDaughter()==-1) return fJetsList;
285 tmp = stack->Particle(tmp->GetFirstDaughter());
286 pdg = tmp->GetPdgCode();
289 //Add found jet to list
290 TParticle *jet1 = new TParticle(*tmp);
291 jet1->SetFirstMother(6);
292 fJetsList->Add(jet1);
293 //printf("jet 1: first daughter %d, last daughter %d\n", tmp->GetFirstDaughter(), tmp->GetLastDaughter());
294 //tmp = stack->Particle(tmp->GetFirstDaughter());
298 printf("HERWIG Jet 1: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
299 jet1->GetFirstMother(),jet1->GetStatusCode(),jet1->Pt(),jet1->Energy(),jet1->Phi()*TMath::RadToDeg(),jet1->Eta());
306 if(parton2->GetPdgCode()!=22){
308 if(tmp->GetFirstDaughter()==-1) return fJetsList;
309 i = tmp->GetFirstDaughter();
310 tmp = stack->Particle(tmp->GetFirstDaughter());
311 pdg = tmp->GetPdgCode();
313 //Add found jet to list
314 TParticle *jet2 = new TParticle(*tmp);
315 jet2->SetFirstMother(7);
316 fJetsList->Add(jet2);
319 printf("HERWIG Jet 2: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
320 jet2->GetFirstMother(),jet2->GetStatusCode(),jet2->Pt(),jet2->Energy(),jet2->Phi()*TMath::RadToDeg(),jet2->Eta());
321 Int_t first = tmp->GetFirstDaughter();
322 Int_t last = tmp->GetLastDaughter();
323 printf("jet 2: first daughter %d, last daughter %d, pdg %d\n",first, last, tmp->GetPdgCode());
324 // for(Int_t d = first ; d < last+1; d++){
325 // tmp = stack->Particle(d);
326 // if(i == tmp->GetFirstMother())
327 // printf("Daughter n %d, Mother %d, name %s, status %d, pT %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
328 // d,tmp->GetFirstMother(), tmp->GetName(), tmp->GetStatusCode(),tmp->Pt(),tmp->Energy(),tmp->Phi()*TMath::RadToDeg(),tmp->Eta());
332 }//Herwig generated jets
338 //________________________________________________________________
339 void AliMCAnalysisUtils::Print(const Option_t * opt) const
342 //Print some relevant parameters set for the analysis
346 printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
348 printf("Debug level = %d\n",fDebug);
349 printf("MC Generator = %s\n",fMCGenerator.Data());