[u/mrichter/AliRoot.git] / PWG4 / PartCorrBase / AliMCAnalysisUtils.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/
/* $Id: AliMCAnalysisUtils.cxx 21839 2007-10-29 13:49:42Z gustavo $ */

//_________________________________________________________________________
// Class for analysis utils for MC data
// stored in stack or event header.
// Contains:
//  - method to check the origin of a given track/cluster
//  - method to obtain the generated jets
//                
//*-- Author: Gustavo Conesa (LNF-INFN) 
//////////////////////////////////////////////////////////////////////////////


// --- ROOT system ---
#include <TMath.h>
#include <TString.h>
#include <TList.h>

//---- ANALYSIS system ----
#include "AliLog.h"
#include "AliMCAnalysisUtils.h"
#include "AliStack.h"
#include "TParticle.h"
#include "AliGenPythiaEventHeader.h"

ClassImp(AliMCAnalysisUtils)


//________________________________________________
AliMCAnalysisUtils::AliMCAnalysisUtils() : 
TObject(), fCurrentEvent(-1), fDebug(-1), 
fJetsList(new TList), fMCGenerator("PYTHIA")
{
	//Ctor
}

//____________________________________________________________________________
AliMCAnalysisUtils::AliMCAnalysisUtils(const AliMCAnalysisUtils & mcutils) :   
TObject(mcutils), fCurrentEvent(mcutils.fCurrentEvent), fDebug(mcutils.fDebug),
fJetsList(mcutils.fJetsList), fMCGenerator(mcutils.fMCGenerator)
{
	// cpy ctor
	
}

//_________________________________________________________________________
AliMCAnalysisUtils & AliMCAnalysisUtils::operator = (const AliMCAnalysisUtils & mcutils)
{
	// assignment operator
	
	if(&mcutils == this) return *this;
	fCurrentEvent = mcutils.fCurrentEvent ;
	fDebug        = mcutils.fDebug;
	fJetsList     = mcutils.fJetsList;
	fMCGenerator  = mcutils.fMCGenerator;

	return *this;
	
}

//____________________________________________________________________________
AliMCAnalysisUtils::~AliMCAnalysisUtils() 
{
	// Remove all pointers.
	
	if (fJetsList) {
		fJetsList->Clear();
		delete fJetsList ;
	}   

}

//_________________________________________________________________________
Int_t AliMCAnalysisUtils::CheckOrigin(const Int_t label, AliStack * stack) const {
  //Play with the MC stack if available
  //Check origin of the candidates, good for PYTHIA
  
  if(!stack) AliFatal("Stack is not available, check analysis settings in configuration file, STOP!!");
  //  printf("label %d, ntrack %d, nprim %d\n",label, stack->GetNtrack(), stack->GetNprimary());
//   for(Int_t i = 0; i< stack->GetNprimary(); i++){
//      TParticle *particle =   stack->Particle(i);
// 			//particle->Print();
//   }
  if(label >= 0 && label <  stack->GetNtrack()){
    //Mother
    TParticle * mom = stack->Particle(label);
    Int_t mPdg = TMath::Abs(mom->GetPdgCode());
    Int_t mStatus =  mom->GetStatusCode() ;
    Int_t iParent =  mom->GetFirstMother() ;
    if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOrigin: Mother is parton %d\n",iParent);
    
    //GrandParent
    TParticle * parent = new TParticle ;
    Int_t pPdg = -1;
    Int_t pStatus =-1;
    if(iParent > 0){
      parent = stack->Particle(iParent);
      pPdg = TMath::Abs(parent->GetPdgCode());
      pStatus = parent->GetStatusCode();  
    }
    else if(fDebug > 0 ) printf("AliMCAnalysisUtils::CheckOrigin: Parent with label %d\n",iParent);
    
    //return tag
    if(mPdg == 22){
      if(mStatus == 1){
	if(fMCGenerator == "PYTHIA"){
	  if(iParent < 8 && iParent > 5) {//outgoing partons
	    if(pPdg == 22) return kMCPrompt;
	    else  return kMCFragmentation;
	  }//Outgoing partons
	  else if(pStatus == 11){//Decay
	    if(pPdg == 111) return kMCPi0Decay ;
	    else if (pPdg == 321)  return kMCEtaDecay ;
	    else  return kMCOtherDecay ;
	  }//Decay
	  else return kMCISR; //Initial state radiation
	}//PYTHIA

	else if(fMCGenerator == "HERWIG"){	  
	  if(pStatus < 197){//Not decay
 	    while(1){
	      if(parent->GetFirstMother()<=5) break;
	      iParent = parent->GetFirstMother();
	      parent=stack->Particle(iParent);
	      pStatus= parent->GetStatusCode();
	      pPdg = parent->GetPdgCode();
	    }//Look for the parton
	    
	    if(iParent < 8 && iParent > 5) {
	      if(pPdg == 22) return kMCPrompt;
	      else  return kMCFragmentation;
	    }
	    return kMCISR;//Initial state radiation
	  }//Not decay
	  else{//Decay
	    if(pPdg == 111) return kMCPi0Decay ;
	    else if (pPdg == 321)  return kMCEtaDecay ;
	    else  return kMCOtherDecay ;
	  }//Decay
	}//HERWIG
	else return  kMCUnknown;
      }//Status 1 : Pythia generated
      else if(mStatus == 0){
	if(pPdg ==22 || pPdg ==11|| pPdg == 2112 ||  pPdg == 211 ||  
	   pPdg == 321 ||  pPdg == 2212  ||  pPdg == 130  ||  pPdg == 13 ) 
	  return kMCConversion ;
	if(pPdg == 111) return kMCPi0Decay ;
	else if (pPdg == 221)  return kMCEtaDecay ;
	else  return kMCOtherDecay ;
      }//status 0 : geant generated
    }//Mother Photon
    else if(mPdg == 111)  return kMCPi0 ;
    else if(mPdg == 221)  return kMCEta ;
    else if(mPdg ==11){
//       if(pPdg !=22 &&mStatus == 0) {
// 	printf("Origin electron, pT %f, status %d, parent %s, pstatus %d, vx %f, vy %f, vz %f\n",
// 	       mom->Pt(),mStatus, parent->GetName(),pStatus,mom->Vx(),mom->Vy(), mom->Vz());

//       }

      if(mStatus == 0) {
	if(pPdg ==22) return kMCConversion ;
	else if (pStatus == 0) return kMCConversion;
	else return kMCElectron ;
      }
      else return kMCElectron ;
    }
    else return kMCUnknown;
  }//Good label value
  else{
    if(label < 0 ) printf("AliMCAnalysisUtils::CheckOrigin: *** bad label or no stack ***:  label %d \n", label);
    if(label >=  stack->GetNtrack()) printf("AliMCAnalysisUtils::CheckOrigin: *** large label ***:  label %d, n tracks %d \n", label, stack->GetNtrack());
    return kMCUnknown;
  }//Bad label
	
  return kMCUnknown;
  
}

//_________________________________________________________________________
TList * AliMCAnalysisUtils::GetJets(Int_t iEvent, AliStack * stack, AliGenEventHeader * geh) {
 //Return list of jets (TParticles) and index of most likely parton that originated it.
	
	if(fCurrentEvent!=iEvent){
		fCurrentEvent = iEvent;
		fJetsList = new TList;
		Int_t nTriggerJets = 0;
		Float_t tmpjet[]={0,0,0,0};
		
		//printf("Event %d %d\n",fCurrentEvent,iEvent);
		//Get outgoing partons
		if(stack->GetNtrack() < 8) return fJetsList;
		TParticle * parton1 =  stack->Particle(6);
		TParticle * parton2 =  stack->Particle(7);
		if(fDebug > 2){
			printf("parton 6 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
				parton1->GetName(),parton1->Pt(),parton1->Energy(),parton1->Phi()*TMath::RadToDeg(),parton1->Eta());
			printf("parton 7 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
				parton2->GetName(),parton2->Pt(),parton2->Energy(),parton2->Phi()*TMath::RadToDeg(),parton2->Eta());
		}
// 		//Trace the jet from the mother parton
// 		Float_t pt  = 0;
// 		Float_t pt1 = 0;
// 		Float_t pt2 = 0;
// 		Float_t e   = 0;
// 		Float_t e1  = 0;
// 		Float_t e2  = 0;
// 		TParticle * tmptmp = new TParticle;
// 		for(Int_t i = 0; i< stack->GetNprimary(); i++){
// 			tmptmp = stack->Particle(i);
		
// 			if(tmptmp->GetStatusCode() == 1){
// 				pt = tmptmp->Pt();
// 				e =  tmptmp->Energy();			
// 				Int_t imom = tmptmp->GetFirstMother();
// 				Int_t imom1 = 0;
// 				//printf("1st imom %d\n",imom);
// 				while(imom > 5){
// 					imom1=imom;
// 					tmptmp = stack->Particle(imom);
// 					imom = tmptmp->GetFirstMother();
// 					//printf("imom %d	\n",imom);
// 				}
// 				//printf("Last imom %d %d\n",imom1, imom);
// 				if(imom1 == 6) {
// 					pt1+=pt;
// 					e1+=e;				
// 				}
// 				else if (imom1 == 7){
// 					pt2+=pt;
// 					e2+=e;					}
// 			}// status 1
				
// 		}// for
		
// 		printf("JET 1, pt %2.2f, e %2.2f; JET 2, pt %2.2f, e %2.2f \n",pt1,e1,pt2,e2);
		
		//Get the jet, different way for different generator
		//PYTHIA
		if(fMCGenerator == "PYTHIA"){
			TParticle * jet =  new TParticle;
			AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) geh;
			nTriggerJets =  pygeh->NTriggerJets();
			//if(fDebug > 1)
				printf("PythiaEventHeader: Njets: %d\n",nTriggerJets);
			Int_t iparton = -1;
			for(Int_t i = 0; i< nTriggerJets; i++){
				iparton=-1;
				pygeh->TriggerJet(i, tmpjet);
				jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
				//Assign an outgoing parton as mother
				Float_t phidiff1 = TMath::Abs(jet->Phi()-parton1->Phi());		
				Float_t phidiff2 = TMath::Abs(jet->Phi()-parton2->Phi());
				if(phidiff1 > phidiff2) jet->SetFirstMother(7);
				else  jet->SetFirstMother(6);
				//jet->Print();
				if(fDebug > 1)
				  printf("PYTHIA Jet %d: mother %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
					 i, jet->GetFirstMother(),jet->Pt(),jet->Energy(),jet->Phi()*TMath::RadToDeg(),jet->Eta());
				fJetsList->Add(jet);			
			}
		}//Pythia triggered jets
		//HERWIG
		else if (fMCGenerator=="HERWIG"){
			Int_t pdg = -1;		
			//Check parton 1
			TParticle * tmp = parton1;
			if(parton1->GetPdgCode()!=22){
				while(pdg != 94){
					if(tmp->GetFirstDaughter()==-1) return fJetsList;
					tmp = stack->Particle(tmp->GetFirstDaughter());
					pdg = tmp->GetPdgCode();
				}//while
				
				//Add found jet to list
				TParticle *jet1 = new TParticle(*tmp);
				jet1->SetFirstMother(6);
				fJetsList->Add(jet1);
				//printf("jet 1:  first daughter %d, last daughter %d\n", tmp->GetFirstDaughter(), tmp->GetLastDaughter());
				//tmp = stack->Particle(tmp->GetFirstDaughter());
				//tmp->Print();
				//jet1->Print();
				if(fDebug > 1)			
					printf("HERWIG Jet 1: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
						jet1->GetFirstMother(),jet1->GetStatusCode(),jet1->Pt(),jet1->Energy(),jet1->Phi()*TMath::RadToDeg(),jet1->Eta());
			}//not photon
				
			//Check parton 2
			pdg = -1;
			tmp = parton2;
			Int_t i = -1;
			if(parton2->GetPdgCode()!=22){
				while(pdg != 94){
					if(tmp->GetFirstDaughter()==-1) return fJetsList;
					i = tmp->GetFirstDaughter();
					tmp = stack->Particle(tmp->GetFirstDaughter());
					pdg = tmp->GetPdgCode();
				}//while
				//Add found jet to list
				TParticle *jet2 = new TParticle(*tmp);
				jet2->SetFirstMother(7);
				fJetsList->Add(jet2);
				//jet2->Print();
				if(fDebug > 1)
					printf("HERWIG Jet 2: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
					jet2->GetFirstMother(),jet2->GetStatusCode(),jet2->Pt(),jet2->Energy(),jet2->Phi()*TMath::RadToDeg(),jet2->Eta());
					Int_t first =  tmp->GetFirstDaughter();
					Int_t last  =  tmp->GetLastDaughter();
					printf("jet 2:  first daughter %d, last daughter %d, pdg %d\n",first, last, tmp->GetPdgCode());
				//	for(Int_t d = first ; d < last+1; d++){
//						tmp = stack->Particle(d);
//						if(i == tmp->GetFirstMother())
//							printf("Daughter n %d, Mother %d, name %s, status %d, pT %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
//							d,tmp->GetFirstMother(), tmp->GetName(), tmp->GetStatusCode(),tmp->Pt(),tmp->Energy(),tmp->Phi()*TMath::RadToDeg(),tmp->Eta());			   
//			   }
			  			   //tmp->Print();
			   }//not photon
		}//Herwig generated jets
	}
	
	return fJetsList;
}

//________________________________________________________________
void AliMCAnalysisUtils::Print(const Option_t * opt) const
{
	
	//Print some relevant parameters set for the analysis
	if(! opt)
		return;
	
	printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
	
	printf("Debug level    = %d\n",fDebug);
	printf("MC Generator   = %s\n",fMCGenerator.Data());

	printf(" \n");
	
}
Commit	Line	Data
6639984f	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	/* $Id: AliMCAnalysisUtils.cxx 21839 2007-10-29 13:49:42Z gustavo $ */
	16
	17	//_________________________________________________________________________
	18	// Class for analysis utils for MC data
	19	// stored in stack or event header.
	20	// Contains:
	21	// - method to check the origin of a given track/cluster
	22	// - method to obtain the generated jets
	23	//
	24	//*-- Author: Gustavo Conesa (LNF-INFN)
	25	//////////////////////////////////////////////////////////////////////////////
	26
	27
	28	// --- ROOT system ---
	29	#include <TMath.h>
	30	#include <TString.h>
	31	#include <TList.h>
	32
	33	//---- ANALYSIS system ----
	34	#include "AliLog.h"
	35	#include "AliMCAnalysisUtils.h"
	36	#include "AliStack.h"
	37	#include "TParticle.h"
	38	#include "AliGenPythiaEventHeader.h"
	39
	40	ClassImp(AliMCAnalysisUtils)
	41
	42
	43	//________________________________________________
	44	AliMCAnalysisUtils::AliMCAnalysisUtils() :
	45	TObject(), fCurrentEvent(-1), fDebug(-1),
	46	fJetsList(new TList), fMCGenerator("PYTHIA")
	47	{
	48	//Ctor
	49	}
	50
	51	//____________________________________________________________________________
	52	AliMCAnalysisUtils::AliMCAnalysisUtils(const AliMCAnalysisUtils & mcutils) :
	53	TObject(mcutils), fCurrentEvent(mcutils.fCurrentEvent), fDebug(mcutils.fDebug),
	54	fJetsList(mcutils.fJetsList), fMCGenerator(mcutils.fMCGenerator)
	55	{
	56	// cpy ctor
	57
	58	}
	59
	60	//_________________________________________________________________________
	61	AliMCAnalysisUtils & AliMCAnalysisUtils::operator = (const AliMCAnalysisUtils & mcutils)
	62	{
	63	// assignment operator
	64
65	if(&mcutils == this) return *this;
66	fCurrentEvent = mcutils.fCurrentEvent ;
67	fDebug = mcutils.fDebug;
68	fJetsList = mcutils.fJetsList;
69	fMCGenerator = mcutils.fMCGenerator;
70
71	return *this;
72
73	}
74
75	//____________________________________________________________________________
76	AliMCAnalysisUtils::~AliMCAnalysisUtils()
77	{
78	// Remove all pointers.
79
80	if (fJetsList) {
81	fJetsList->Clear();
82	delete fJetsList ;
83	}
84
85	}
86
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
91
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();
97	// }
98	if(label >= 0 && label < stack->GetNtrack()){
99	//Mother
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);
105
106	//GrandParent
107	TParticle * parent = new TParticle ;
108	Int_t pPdg = -1;
109	Int_t pStatus =-1;
110	if(iParent > 0){
111	parent = stack->Particle(iParent);
112	pPdg = TMath::Abs(parent->GetPdgCode());
113	pStatus = parent->GetStatusCode();
114	}
115	else if(fDebug > 0 ) printf("AliMCAnalysisUtils::CheckOrigin: Parent with label %d\n",iParent);
116
117	//return tag
118	if(mPdg == 22){
119	if(mStatus == 1){
120	if(fMCGenerator == "PYTHIA"){
121	if(iParent < 8 && iParent > 5) {//outgoing partons
122	if(pPdg == 22) return kMCPrompt;
123	else return kMCFragmentation;
124	}//Outgoing partons
125	else if(pStatus == 11){//Decay
126	if(pPdg == 111) return kMCPi0Decay ;
127	else if (pPdg == 321) return kMCEtaDecay ;
128	else return kMCOtherDecay ;
129	}//Decay
130	else return kMCISR; //Initial state radiation
131	}//PYTHIA
132
133	else if(fMCGenerator == "HERWIG"){
134	if(pStatus < 197){//Not decay
135	while(1){
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
142
143	if(iParent < 8 && iParent > 5) {
144	if(pPdg == 22) return kMCPrompt;
145	else return kMCFragmentation;
146	}
147	return kMCISR;//Initial state radiation
148	}//Not decay
149	else{//Decay
150	if(pPdg == 111) return kMCPi0Decay ;
151	else if (pPdg == 321) return kMCEtaDecay ;
152	else return kMCOtherDecay ;
153	}//Decay
154	}//HERWIG
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
165	}//Mother Photon
166	else if(mPdg == 111) return kMCPi0 ;
167	else if(mPdg == 221) return kMCEta ;
168	else if(mPdg ==11){
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());
172
173	// }
174
175	if(mStatus == 0) {
176	if(pPdg ==22) return kMCConversion ;
177	else if (pStatus == 0) return kMCConversion;
178	else return kMCElectron ;
179	}
180	else return kMCElectron ;
181	}
182	else return kMCUnknown;
183	}//Good label value
184	else{
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());
187	return kMCUnknown;
188	}//Bad label
189
190	return kMCUnknown;
191
192	}
193
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.
197
198	if(fCurrentEvent!=iEvent){
199	fCurrentEvent = iEvent;
200	fJetsList = new TList;
201	Int_t nTriggerJets = 0;
202	Float_t tmpjet[]={0,0,0,0};
203
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);
209	if(fDebug > 2){
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());
214	}
215	// //Trace the jet from the mother parton
216	// Float_t pt = 0;
217	// Float_t pt1 = 0;
218	// Float_t pt2 = 0;
219	// Float_t e = 0;
220	// Float_t e1 = 0;
221	// Float_t e2 = 0;
222	// TParticle * tmptmp = new TParticle;
223	// for(Int_t i = 0; i< stack->GetNprimary(); i++){
224	// tmptmp = stack->Particle(i);
225
226	// if(tmptmp->GetStatusCode() == 1){
227	// pt = tmptmp->Pt();
228	// e = tmptmp->Energy();
229	// Int_t imom = tmptmp->GetFirstMother();
230	// Int_t imom1 = 0;
231	// //printf("1st imom %d\n",imom);
232	// while(imom > 5){
233	// imom1=imom;
234	// tmptmp = stack->Particle(imom);
235	// imom = tmptmp->GetFirstMother();
236	// //printf("imom %d \n",imom);
237	// }
238	// //printf("Last imom %d %d\n",imom1, imom);
239	// if(imom1 == 6) {
240	// pt1+=pt;
241	// e1+=e;
242	// }
243	// else if (imom1 == 7){
244	// pt2+=pt;
245	// e2+=e; }
246	// }// status 1
247
248	// }// for
249
250	// printf("JET 1, pt %2.2f, e %2.2f; JET 2, pt %2.2f, e %2.2f \n",pt1,e1,pt2,e2);
251
252	//Get the jet, different way for different generator
253	//PYTHIA
254	if(fMCGenerator == "PYTHIA"){
255	TParticle * jet = new TParticle;
256	AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) geh;
257	nTriggerJets = pygeh->NTriggerJets();
258	//if(fDebug > 1)
259	printf("PythiaEventHeader: Njets: %d\n",nTriggerJets);
260	Int_t iparton = -1;
261	for(Int_t i = 0; i< nTriggerJets; i++){
262	iparton=-1;
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);
270	//jet->Print();
271	if(fDebug > 1)
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());
274	fJetsList->Add(jet);
275	}
276	}//Pythia triggered jets
277	//HERWIG
278	else if (fMCGenerator=="HERWIG"){
279	Int_t pdg = -1;
280	//Check parton 1
281	TParticle * tmp = parton1;
282	if(parton1->GetPdgCode()!=22){
283	while(pdg != 94){
284	if(tmp->GetFirstDaughter()==-1) return fJetsList;
285	tmp = stack->Particle(tmp->GetFirstDaughter());
286	pdg = tmp->GetPdgCode();
287	}//while
288
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());
295	//tmp->Print();
296	//jet1->Print();
297	if(fDebug > 1)
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());
300	}//not photon
301
302	//Check parton 2
303	pdg = -1;
304	tmp = parton2;
305	Int_t i = -1;
306	if(parton2->GetPdgCode()!=22){
307	while(pdg != 94){
308	if(tmp->GetFirstDaughter()==-1) return fJetsList;
309	i = tmp->GetFirstDaughter();
310	tmp = stack->Particle(tmp->GetFirstDaughter());
311	pdg = tmp->GetPdgCode();
312	}//while
313	//Add found jet to list
314	TParticle jet2 = new TParticle(tmp);
315	jet2->SetFirstMother(7);
316	fJetsList->Add(jet2);
317	//jet2->Print();
318	if(fDebug > 1)
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());
329	// }
330	//tmp->Print();
331	}//not photon
332	}//Herwig generated jets
333	}
334
335	return fJetsList;
336	}
337
338	//________________________________________________________________
339	void AliMCAnalysisUtils::Print(const Option_t * opt) const
340	{
341
342	//Print some relevant parameters set for the analysis
343	if(! opt)
344	return;
345
346	printf("*** Print: %s %s ****\n", GetName(), GetTitle() ) ;
347
348	printf("Debug level = %d\n",fDebug);
349	printf("MC Generator = %s\n",fMCGenerator.Data());
350
351	printf(" \n");
352
353	}
354
355