[u/mrichter/AliRoot.git] / EMCAL / AliEMCALJetFinderInputSimPrep.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$ */

//_________________________________________________________________________
//  Class for JetFinder Input preparation from simulated data 
//
//*-- Author: Mark Horner (LBL/UCT)
//


#include "AliEMCALJetFinderInputSimPrep.h"

#include <TParticle.h>
#include <TParticlePDG.h>
#include <TPDGCode.h>
#include <TFile.h>
#include <TTree.h>
#include <TObjectTable.h>
#include <TMath.h>


#include "AliRun.h"
#include "AliMagF.h"
#include "AliEMCALFast.h"
#include "AliEMCAL.h"
#include "AliEMCALHit.h"
#include "AliEMCALGeometry.h"
#include "AliEMCALGetter.h"
#include "AliEMCALLoader.h"
#include "AliGenEventHeader.h"
#include "AliGenPythiaEventHeader.h"
#include "AliGenerator.h"
#include "AliHeader.h"
#include "AliMC.h"


ClassImp(AliEMCALJetFinderInputSimPrep)
	
AliEMCALJetFinderInputSimPrep::AliEMCALJetFinderInputSimPrep()
{
	// Default constructor
if (fDebug > 0) Info("AliEMCALJetFinderInputSimPrep","Beginning Constructor");	
  
  fDebug = 0;
  fSmearType = kSmearEffic;
  fEMCALType = kHits;
  fTrackType = kCharged;
  fEfficiency = 0.90;
  
}
AliEMCALJetFinderInputSimPrep::~AliEMCALJetFinderInputSimPrep()
{
if (fDebug > 0) Info("~AliEMCALJetFinderInputSimPrep","Beginning Destructor");	
	
}

void AliEMCALJetFinderInputSimPrep::Reset(AliEMCALJetFinderResetType_t resettype)
{
	// Method to reset data
if (fDebug > 1) Info("Reset","Beginning Reset");
	switch (resettype){

	case kResetData:
		fInputObject.Reset(resettype);
		break;
	case kResetTracks:
		fInputObject.Reset(resettype);
		break;
	case kResetDigits:
		fInputObject.Reset(resettype);
		break;
	case kResetParameters:
		fSmearType = kSmearEffic;
	       	fEMCALType = kHits;
		fTrackType = kCharged;
		break;
	case kResetAll:
		fSmearType = kSmearEffic;
	       	fEMCALType = kHits;
		fTrackType = kCharged;
		fInputObject.Reset(resettype);
		break;
	case kResetPartons:
	  Warning("FillFromFile", "kResetPartons not implemented") ; 
	  break;
	case kResetParticles:
	  Warning("FillFromFile", "kResetParticles not implemented") ; 
	  break;
	case kResetJets:
	  Warning("FillFromFile", "kResetJets not implemented") ; 
	  break;
	}// end switch

}

Int_t AliEMCALJetFinderInputSimPrep::FillFromFile(TString *filename, AliEMCALJetFinderFileType_t filetype,Int_t EventNumber)
{
if (fDebug > 1) Info("FillFromFile","Beginning FillFromFile");
   fFileType = filetype;	
   AliEMCALGetter *gime = AliEMCALGetter::Instance(filename->Data());
   if (fDebug > 1) Info("FillFromFile","Instantiated Getter with %s as the file",filename->Data());
   gime->Event(EventNumber,"XH") ;
   if (fDebug > 1) Info("FillFromFile","Got event %i with option \"XH\"",EventNumber);

   if (	fEMCALType == kHits ||
	fEMCALType == kTimeCut )  FillHits();
   if ( fTrackType != kNoTracks  )  FillTracks();	
   if ( fFileType  != kData){
	   FillPartons();
   }
/*   gime->EmcalLoader()->UnloadRecParticles();   
   gime->EmcalLoader()->UnloadTracks();   
   gime->Reset();*/
   return 0;	
}

void AliEMCALJetFinderInputSimPrep::FillHits()		
{
	// Fill from the hits to input object from simulation
if (fDebug > 1) Info("FillHits","Beginning FillHits");
	
// Access hit information
    AliEMCALHit *mHit;
//  AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
//  Info("AliEMCALJetFinderInputSimPrep","Title of the geometry is %s",pEMCAL->GetTitle());
    AliEMCALGeometry* geom =  AliEMCALGetter::Instance()->EMCALGeometry();
    AliEMCALGetter *gime = AliEMCALGetter::Instance();
    
    //    TTree *treeH = AliEMCALGetter::Instance()->TreeH();
//    Int_t ntracks = (Int_t) treeH->GetEntries();
//
//   Loop over tracks
//
    Int_t nbytes = 0;
    Double_t etH = 0.0;

/*    for (Int_t track=0; track<ntracks;track++) {
        gAlice->ResetHits();
        nbytes += treeH->GetEvent(track);*/
//
//  Loop over hits
//
        for(Int_t i =0; i < gime->Hits()->GetEntries() ;i++)
        {
	    mHit = gime->Hit(i);	
                if (fEMCALType == kTimeCut && 
		    (mHit->GetTime()>fTimeCut) ) continue;
            Float_t x      =    mHit->X();         // x-pos of hit
            Float_t y      =    mHit->Y();         // y-pos
            Float_t z      =    mHit->Z();         // z-pos
            Float_t eloss  =    mHit->GetEnergy(); // deposited energy

            Float_t r      =    TMath::Sqrt(x*x+y*y);
            Float_t theta  =    TMath::ATan2(r,z);
            Float_t eta    =   -TMath::Log(TMath::Tan(theta/2.));
            Float_t phi    =    TMath::ATan2(y,x);
            etH = eloss*TMath::Sin(theta);
	    if (fDebug > 10) Info("FillHits","Values of hit energy %i",Int_t(1e7*etH));
	    if (geom->TowerIndexFromEtaPhi(eta,180.0/TMath::Pi()*phi) == -1) 
	    {
		    if (fDebug >5)  
		    {
			    Error("FillHits","Hit not inside EMCAL!");
			    mHit->Dump();
		    }
		    continue;
	    }
            fInputObject.AddEnergyToDigit(geom->TowerIndexFromEtaPhi(eta,180.0/TMath::Pi()*phi)-1,Int_t(1e7*etH));
        } // Hit Loop
//    } // Track Loop


}
void AliEMCALJetFinderInputSimPrep::FillTracks()	
{
	// Fill from particles simulating a TPC to input object from simulation

    if (fDebug > 1) Info("FillTracks","Beginning FillTracks");
	
    TParticlePDG* pdgP = 0;
    TParticle *mPart;
    Int_t npart = (gAlice->GetHeader())->GetNprimary();
    if (fDebug > 1) Info("FillTracks","Header says there are %i primaries",npart);
    AliEMCALGetter *gime = AliEMCALGetter::Instance();
    Float_t bfield,rEMCAL;		 
/*	for (Int_t i =0; i<gime->NPrimaries(); i++)
	{
		cout<<gime->Primary(i)->GetFirstMother()<<endl;
		if (gime->Primary(i)->GetFirstMother() == -1) {
			counterforus++;
			gime->Primary(i)->Dump();
		}
	}*/
    
    if (fDebug > 1) Info("FillTracks","Defining the geometry");
    
    AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
    AliEMCALGeometry* geom =  AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
    fEtaMax = geom->GetArm1EtaMax();
    fEtaMin = geom->GetArm1EtaMin();
    fPhiMax = TMath::Pi()/180.0*geom->GetArm1PhiMax();
    fPhiMin = TMath::Pi()/180.0*geom->GetArm1PhiMin();

    if (fDebug > 1) Info("FillTracks","Starting particle loop");
	    
    if (fDebug > 1) Info("FillTracks","Particle loop of %i",npart);
    for (Int_t part = 0; part < npart; part++) {
	mPart = gAlice->GetMCApp()->Particle(part);
	//if (part%10) gObjectTable->Print();
	pdgP = mPart->GetPDG();

	if (fDebug > 5) Info("FillTracks","Checking if track is a primary");
	
	if (fFileType == kPythia) {
	    if (mPart->GetStatusCode() != 1) continue;
	} else if (fFileType == kHijing) {
	    if (mPart->GetFirstDaughter() >= 0 && mPart->GetFirstDaughter() < npart) continue;
	}
	
	if (fDebug > 15) Info("FillTracks","Checking if track (eta - %f, phi - %f) is in acceptance",mPart->Eta(),mPart->Phi());
	if (fDebug > 10) Info("FillTracks","Checking if EMCAL acceptance  ( %f < eta < %f, %f < phi < %f) is in acceptance",fEtaMin,fEtaMax,fPhiMin,fPhiMax);

	if ((!fPythiaComparison)&&(mPart->Eta() > fEtaMax || mPart->Eta() < fEtaMin))   continue;
	if ((!fPythiaComparison)&&(mPart->Phi() > fPhiMax || mPart->Phi() < fPhiMin))   continue;
	
/*
	{kProton, kProtonBar, kElectron, kPositron,
         kNuE, kNuEBar, kGamma, kNeutron, kNeutronBar,
         kMuonPlus, kMuonMinus, kK0Long , kPiPlus, kPiMinus,
         kKPlus, kKMinus, kLambda0, kLambda0Bar, kK0Short,
         kSigmaMinus, kSigmaPlus, kSigma0, kPi0, kK0, kK0Bar,
         0,kNuMu,kNuMuBar
*/

	if (fDebug > 15) Info("FillTracks","Checking if track is rejected");

	if ((fSmearType == kEfficiency 	||
             fSmearType == kSmearEffic)	&& 
 	     pdgP->Charge()!=0) {
	    if (AliEMCALFast::RandomReject(fEfficiency)) {
		continue;
            }
	}

   	bfield = gAlice->Field()->SolenoidField();
	rEMCAL = AliEMCALGeometry::GetInstance()->GetIPDistance();
	Float_t rB = 3335.6 * mPart->Pt() / bfield;  // [cm]  (case of |charge|=1)
	if (2.*rB < rEMCAL) continue;  // track curls away
	
	//if (part%10) gObjectTable->Print();
        switch(fTrackType)
        {

	   case kAll:  // All Stable particles to be included
		if (fDebug > 5) Info("FillTracks","Storing track");
		if (fSmearType == kSmear ||
		    fSmearType == kSmearEffic ){
			Smear(mPart);/*
			TParticle *tmp = Smear(MPart);
			fInputObject.AddTrack(Smear(MPart));
			delete tmp;*/
		}else{
			fInputObject.AddTrack(*mPart);
		}
	   break;
	   case kEM:   // All Electromagnetic particles
		if (mPart->GetPdgCode() == kElectron  ||
                    mPart->GetPdgCode() == kMuonPlus  || 
                    mPart->GetPdgCode() == kMuonMinus ||
		    mPart->GetPdgCode() == kPositron ){
		      if (fDebug > 5) Info("FillTracks","Storing electron or positron");
	 	      if (fSmearType == kSmear ||
		            fSmearType == kSmearEffic ){
			      Smear(mPart);/*
 			      TParticle *tmp = Smear(MPart); 
 			      fInputObject.AddTrack(tmp);
 			      delete tmp;*/
		      }else{
		          fInputObject.AddTrack(*mPart);
		      }
		}
		if ( mPart->GetPdgCode() == kGamma ){ 
			fInputObject.AddTrack(*mPart);
			if (fDebug > 5) Info("FillTracks","Storing gamma");
		}
			
	   break;
           case kCharged: // All Particles with non-zero charge
		if (pdgP->Charge() != 0) {
			if (fDebug > 5) Info("FillTracks","Storing charged track");
			if (fSmearType == kSmear ||
	                    fSmearType == kSmearEffic ){
				Smear(mPart);/*
				TParticle *tmp = Smear(MPart);
    				fInputObject.AddTrack(tmp);
				delete tmp;*/
	                }else{
	                 fInputObject.AddTrack(*mPart);
	                }
		}
	   break;
	   case kNeutral: // All particles with no charge
		if (pdgP->Charge() == 0){ 
			fInputObject.AddTrack(*mPart);
			if (fDebug > 5) Info("FillTracks","Storing neutral");
		}
	   break;
	   case kHadron: //All hadrons
		if (mPart->GetPdgCode() != kElectron  &&
                    mPart->GetPdgCode() != kPositron  &&
                    mPart->GetPdgCode() != kMuonPlus  &&
                    mPart->GetPdgCode() != kMuonMinus &&
                    mPart->GetPdgCode() != kGamma ) 
		{
			if (fDebug > 5) Info("FillTracks","Storing hadron");
			if (pdgP->Charge() == 0){
		    	    fInputObject.AddTrack(*mPart);
			}else{
				if (fSmearType == kSmear ||
			            fSmearType == kSmearEffic ){
					Smear(mPart);/*
					TParticle *tmp = Smear(MPart);	
    					fInputObject.AddTrack(tmp);
					delete tmp;*/
			        }else{
			            fInputObject.AddTrack(*mPart);
			        }
			}
		}
	   break;
	   case kChargedHadron:  // only charged hadrons
		if (mPart->GetPdgCode() != kElectron  &&
                    mPart->GetPdgCode() != kPositron  &&
                    mPart->GetPdgCode() != kGamma     &&
                    mPart->GetPdgCode() != kMuonPlus  &&
                    mPart->GetPdgCode() != kMuonMinus &&
		    pdgP->Charge() 	!= 0   	   ){
			if (fDebug > 5) Info("FillTracks","Storing charged hadron");
		       	if (fSmearType == kSmear ||
		            fSmearType == kSmearEffic ){
				Smear(mPart);/*
				TParticle *tmp = Smear(MPart);
 				fInputObject.AddTrack(tmp);
				delete tmp;*/
		        }else{
		               fInputObject.AddTrack(*mPart);
		        }
		}
	   break;
	   case kNoTracks:
	   break;
	   default:
	   break;
	   delete mPart;
        }	//end of switch
//	Info("FillTracks","After Particle Storage");
	//if (part%10) gObjectTable->Print();

   }	//end of particle loop
   //gObjectTable->Print();	
}

void AliEMCALJetFinderInputSimPrep::FillPartons()	
{
	// Fill partons to
	// input object from simulation
if (fDebug > 1) Info("FillPartons","Beginning FillPartons");

  AliGenEventHeader* evHeader = ((AliHeader*)(gAlice->GetHeader()))->GenEventHeader();
  Float_t triggerJetValues[4];
  AliEMCALParton tempParton;
  
  if (fFileType == kPythia)
  {
	Int_t ntriggerjets = ((AliGenPythiaEventHeader*)evHeader)->NTriggerJets();
	if (fDebug > 1) Info("FillPartons","Number of trigger jets --> %i",ntriggerjets);
       	for (Int_t jetcounter = 0 ; jetcounter < ntriggerjets; jetcounter++){
      	((AliGenPythiaEventHeader*)evHeader)->TriggerJet(jetcounter,triggerJetValues);
       	TLorentzVector tempLParton;
       	tempLParton.SetPxPyPzE(triggerJetValues[0],triggerJetValues[1],triggerJetValues[2],triggerJetValues[3]);
	tempParton.SetEnergy(tempLParton.Energy());
	tempParton.SetEta(tempLParton.Eta());
	tempParton.SetPhi(tempLParton.Phi());
	FillPartonTracks(&tempParton);
	fInputObject.AddParton(&tempParton);
  	}
 	  
  }
}

void AliEMCALJetFinderInputSimPrep::FillParticles()		
{
	// Fill particles to input object from simulation
if (fDebug > 1) Info("FillParticles","Beginning FillParticles");

    Int_t npart = (gAlice->GetHeader())->GetNprimary();
    TParticlePDG* pdgP = 0;
 
    for (Int_t part = 0; part < npart; part++) {
	TParticle *mPart = gAlice->GetMCApp()->Particle(part);
	pdgP = mPart->GetPDG();
	
	if (fDebug > 10) Info("FillParticles","Checking if particle is a primary");
	
	if (fFileType == kPythia) {
	    if (mPart->GetStatusCode() != 1) continue;
	} else if (fFileType == kHijing) {
	    if (mPart->GetFirstDaughter() >= 0 && mPart->GetFirstDaughter() < npart) continue;
	}

	if (fDebug > 10) Info("FillParticles","Checking if particle is in acceptance");
	
	if ((!fPythiaComparison)&&(mPart->Eta() > fEtaMax || mPart->Eta() < fEtaMin))   continue;
	if ((!fPythiaComparison)&&(mPart->Phi() > fPhiMax || mPart->Phi() < fPhiMin))   continue;
	

/*
	{kProton, kProtonBar, kElectron, kPositron,
         kNuE, kNuEBar, kGamma, kNeutron, kNeutronBar,
         kMuonPlus, kMuonMinus, kK0Long , kPiPlus, kPiMinus,
         kKPlus, kKMinus, kLambda0, kLambda0Bar, kK0Short,
         kSigmaMinus, kSigmaPlus, kSigma0, kPi0, kK0, kK0Bar,
         0,kNuMu,kNuMuBar
*/
        switch(fTrackType)
        {

	   case kAll:  // All Stable particles to be included
		if (fDebug > 5) Info("FillParticles","Storing particle");
		fInputObject.AddParticle(mPart);
	   break;
	   case kEM:   // All Electromagnetic particles
		if (mPart->GetPdgCode() == kElectron ||
		    mPart->GetPdgCode() == kPositron ||
		    mPart->GetPdgCode() == kGamma){
			if (fDebug > 5) Info("FillParticles","Storing electromagnetic particle");
			fInputObject.AddParticle(mPart);
		}
	   break;
           case kCharged: // All Particles with non-zero charge
		if (pdgP->Charge() != 0) {
			if (fDebug > 5) Info("FillParticles","Storing charged particle");
   			fInputObject.AddParticle(mPart);
		}
	   break;
	   case kNeutral: // All particles with no charge
		if (pdgP->Charge() == 0){
			if (fDebug > 5) Info("FillParticles","Storing neutral particle");
		       	fInputObject.AddParticle(mPart);
		}
	   break;
	   case kHadron: //All hadrons
		if (mPart->GetPdgCode() != kElectron &&
                    mPart->GetPdgCode() != kPositron &&
                    mPart->GetPdgCode() != kGamma ) 
		{

		if (fDebug > 5) Info("FillParticles","Storing hadron");
		    fInputObject.AddParticle(mPart);
		}
	   break;
	   case kChargedHadron:  // only charged hadrons
		if (mPart->GetPdgCode() != kElectron &&
                    mPart->GetPdgCode() != kPositron &&
                    mPart->GetPdgCode() != kGamma    &&
		    pdgP->Charge() 	!= 0   	   ){
		if (fDebug > 5) Info("FillParticles","Storing charged hadron");
	               fInputObject.AddParticle(mPart);
		}
	   break;
	   case kNoTracks:
	   break;
	   default:
	   break;
        }	//end of switch
    }// end of loop over particles

}
void AliEMCALJetFinderInputSimPrep::FillDigits()  
{
	// Fill digits to input object

}

void AliEMCALJetFinderInputSimPrep::Smear(TParticle *particle)
{
	// Smear particle momentum
if (fDebug > 5) Info("Smear","Beginning Smear");

 	Float_t tmp = AliEMCALFast::SmearMomentum(1,particle->P());
	Float_t tmpt = (tmp/particle->P()) * particle->Pt();
	if (fDebug > 10) Info("Smear","Smearing changed P from %f to %f",particle->Pt(),tmpt);
	TLorentzVector tmpvector;
	tmpvector.SetPtEtaPhiM(tmpt,particle->Eta(), particle->Phi(), particle->GetMass());
// create a new particle with smearing momentum - and no daughter or parent information and the same
// vertex	

	TParticle tmparticle(particle->GetPdgCode(), 1, 0, 0,0,0, 
			tmpvector.Px()	, tmpvector.Py(), tmpvector.Pz(), tmpvector.Energy(), 
			particle->Vx(), particle->Vy(), particle->Vz(), particle->T());
	fInputObject.AddTrack(tmparticle);	
	return;
}

void AliEMCALJetFinderInputSimPrep::FillPartonTracks(AliEMCALParton *parton)
{
	// Populate parton tracks for input distributions
 	if (fDebug>1) Info("FillPartonTracks","Beginning FillPartonTracks");    
	Int_t npart = (gAlice->GetHeader())->GetNprimary();
	Int_t ntracks =0;
	TParticlePDG *getpdg;
	TParticle *tempPart;
	for (Int_t part = 0; part < npart; part++)
	{
		tempPart = gAlice->GetMCApp()->Particle(part);
		if (tempPart->GetStatusCode() != 1) continue;
		if ((!fPythiaComparison)&&(tempPart->Eta() > fEtaMax || tempPart->Eta() < fEtaMin || 
		    tempPart->Phi() > fPhiMax || tempPart->Phi() < fPhiMin) ){
		 	if (fDebug>10) Info("FillPartonTracks","Excluding particle not pointing at the EMCAL");    
			continue;
		}
		getpdg = tempPart->GetPDG();
		if (getpdg->Charge() == 0.0  ) { 
			if (fDebug>10) Info("FillPartonTracks","Excluding a neutral particle");
			continue;
		}
		if ( (parton->Eta() - tempPart->Eta())*(parton->Eta() - tempPart->Eta()) +
	     	(parton->Phi() - tempPart->Phi())*(parton->Phi() - tempPart->Phi()) < 1.0 ) ntracks++;
	}
	Float_t *energy = new Float_t[ntracks];
	Float_t *eta    = new Float_t[ntracks];
	Float_t *phi    = new Float_t[ntracks];
	Int_t   *pdg    = new Int_t[ntracks];
	ntracks=0;
	for (Int_t part = 0; part < npart; part++)
	{
		tempPart = gAlice->GetMCApp()->Particle(part);
		if (tempPart->GetStatusCode() != 1) continue;
		if ((!fPythiaComparison)&&(tempPart->Eta() > fEtaMax || tempPart->Eta() < fEtaMin ||
 		    tempPart->Phi() > fPhiMax || tempPart->Phi() < fPhiMin) ){ 
			if (fDebug>10) Info("FillPartonTracks","Excluding particle not pointing at the EMCAL");       
			continue;
		}
	    	if (tempPart->GetStatusCode() != 1) continue;
  		getpdg = tempPart->GetPDG();
		if (getpdg->Charge() == 0.0  ) {
			if (fDebug>10) Info("FillPartonTracks","Excluding a neutral particle");
			continue;
		}
		if (   (parton->Eta() - tempPart->Eta())*(parton->Eta() - tempPart->Eta()) +
	       	(parton->Phi() - tempPart->Phi())*(parton->Phi() - tempPart->Phi()) < 1.0 )
	       	{
	  		energy[ntracks] = tempPart->Pt();
	      		eta[ntracks] = tempPart->Eta();
		  	phi[ntracks] = tempPart->Phi();
	       		pdg[ntracks] = tempPart->GetPdgCode();
	  		ntracks++;
		}
	}
	parton->SetTrackList(ntracks,energy,eta,phi,pdg);
	delete[] energy;
	delete[] eta;
	delete[] phi;
	delete[] pdg;
}
Commit	Line	Data
f7d5860b	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
ee6b678f	17	/* $Id$ */
f7d5860b	18
	19	//_________________________________________________________________________
	20	// Class for JetFinder Input preparation from simulated data
	21	//
	22	//*-- Author: Mark Horner (LBL/UCT)
	23	//
	24
	25
	26
	27	#include "AliEMCALJetFinderInputSimPrep.h"
	28
	29	#include <TParticle.h>
	30	#include <TParticlePDG.h>
	31	#include <TPDGCode.h>
	32	#include <TFile.h>
	33	#include <TTree.h>
	34	#include <TObjectTable.h>
	35	#include <TMath.h>
	36
	37
	38	#include "AliRun.h"
	39	#include "AliMagF.h"
	40	#include "AliEMCALFast.h"
	41	#include "AliEMCAL.h"
	42	#include "AliEMCALHit.h"
	43	#include "AliEMCALGeometry.h"
	44	#include "AliEMCALGetter.h"
b3c7d6bc	45	#include "AliEMCALLoader.h"
f7d5860b	46	#include "AliGenEventHeader.h"
	47	#include "AliGenPythiaEventHeader.h"
	48	#include "AliGenerator.h"
	49	#include "AliHeader.h"
5d12ce38	50	#include "AliMC.h"
f7d5860b	51
f7d5860b	52
b3c7d6bc	53
f7d5860b	54	ClassImp(AliEMCALJetFinderInputSimPrep)
	55
	56	AliEMCALJetFinderInputSimPrep::AliEMCALJetFinderInputSimPrep()
	57	{
44f59d68	58	// Default constructor
f7d5860b	59	if (fDebug > 0) Info("AliEMCALJetFinderInputSimPrep","Beginning Constructor");
	60
	61	fDebug = 0;
	62	fSmearType = kSmearEffic;
	63	fEMCALType = kHits;
	64	fTrackType = kCharged;
	65	fEfficiency = 0.90;
	66
	67	}
	68	AliEMCALJetFinderInputSimPrep::~AliEMCALJetFinderInputSimPrep()
	69	{
	70	if (fDebug > 0) Info("~AliEMCALJetFinderInputSimPrep","Beginning Destructor");
	71
	72	}
	73
	74	void AliEMCALJetFinderInputSimPrep::Reset(AliEMCALJetFinderResetType_t resettype)
	75	{
44f59d68	76	// Method to reset data
f7d5860b	77	if (fDebug > 1) Info("Reset","Beginning Reset");
	78	switch (resettype){
	79
	80	case kResetData:
	81	fInputObject.Reset(resettype);
	82	break;
	83	case kResetTracks:
	84	fInputObject.Reset(resettype);
	85	break;
	86	case kResetDigits:
	87	fInputObject.Reset(resettype);
	88	break;
	89	case kResetParameters:
	90	fSmearType = kSmearEffic;
	91	fEMCALType = kHits;
	92	fTrackType = kCharged;
	93	break;
	94	case kResetAll:
	95	fSmearType = kSmearEffic;
	96	fEMCALType = kHits;
	97	fTrackType = kCharged;
	98	fInputObject.Reset(resettype);
	99	break;
	100	case kResetPartons:
	101	Warning("FillFromFile", "kResetPartons not implemented") ;
	102	break;
	103	case kResetParticles:
	104	Warning("FillFromFile", "kResetParticles not implemented") ;
	105	break;
	106	case kResetJets:
	107	Warning("FillFromFile", "kResetJets not implemented") ;
	108	break;
	109	}// end switch
	110
	111	}
	112
	113	Int_t AliEMCALJetFinderInputSimPrep::FillFromFile(TString *filename, AliEMCALJetFinderFileType_t filetype,Int_t EventNumber)
	114	{
f7d5860b	115	if (fDebug > 1) Info("FillFromFile","Beginning FillFromFile");
f7d5860b	116	fFileType = filetype;
f7d5860b	117	AliEMCALGetter *gime = AliEMCALGetter::Instance(filename->Data());
b3c7d6bc	118	if (fDebug > 1) Info("FillFromFile","Instantiated Getter with %s as the file",filename->Data());
	119	gime->Event(EventNumber,"XH") ;
	120	if (fDebug > 1) Info("FillFromFile","Got event %i with option \"XH\"",EventNumber);
	121
f7d5860b	122	if ( fEMCALType == kHits \|\|
	123	fEMCALType == kTimeCut ) FillHits();
	124	if ( fTrackType != kNoTracks ) FillTracks();
	125	if ( fFileType != kData){
	126	FillPartons();
f7d5860b	127	}
b3c7d6bc	128	/* gime->EmcalLoader()->UnloadRecParticles();
	129	gime->EmcalLoader()->UnloadTracks();
	130	gime->Reset();*/
f7d5860b	131	return 0;
	132	}
	133
44f59d68	134	void AliEMCALJetFinderInputSimPrep::FillHits()
f7d5860b	135	{
44f59d68	136	// Fill from the hits to input object from simulation
f7d5860b	137	if (fDebug > 1) Info("FillHits","Beginning FillHits");
	138
	139	// Access hit information
	140	AliEMCALHit *mHit;
b3c7d6bc	141	// AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
b3c7d6bc	142	// Info("AliEMCALJetFinderInputSimPrep","Title of the geometry is %s",pEMCAL->GetTitle());
f7d5860b	143	AliEMCALGeometry* geom = AliEMCALGetter::Instance()->EMCALGeometry();
b3c7d6bc	144	AliEMCALGetter *gime = AliEMCALGetter::Instance();
	145
	146	// TTree *treeH = AliEMCALGetter::Instance()->TreeH();
	147	// Int_t ntracks = (Int_t) treeH->GetEntries();
f7d5860b	148	//
	149	// Loop over tracks
	150	//
	151	Int_t nbytes = 0;
	152	Double_t etH = 0.0;
	153
b3c7d6bc	154	/* for (Int_t track=0; track<ntracks;track++) {
f7d5860b	155	gAlice->ResetHits();
b3c7d6bc	156	nbytes += treeH->GetEvent(track);*/
f7d5860b	157	//
	158	// Loop over hits
	159	//
b3c7d6bc	160	for(Int_t i =0; i < gime->Hits()->GetEntries() ;i++)
f7d5860b	161	{
b3c7d6bc	162	mHit = gime->Hit(i);
f7d5860b	163	if (fEMCALType == kTimeCut &&
	164	(mHit->GetTime()>fTimeCut) ) continue;
	165	Float_t x = mHit->X(); // x-pos of hit
	166	Float_t y = mHit->Y(); // y-pos
	167	Float_t z = mHit->Z(); // z-pos
	168	Float_t eloss = mHit->GetEnergy(); // deposited energy
	169
	170	Float_t r = TMath::Sqrt(xx+yy);
	171	Float_t theta = TMath::ATan2(r,z);
	172	Float_t eta = -TMath::Log(TMath::Tan(theta/2.));
	173	Float_t phi = TMath::ATan2(y,x);
5255a1ce	174	etH = eloss*TMath::Sin(theta);
f7d5860b	175	if (fDebug > 10) Info("FillHits","Values of hit energy %i",Int_t(1e7*etH));
	176	if (geom->TowerIndexFromEtaPhi(eta,180.0/TMath::Pi()*phi) == -1)
	177	{
	178	if (fDebug >5)
	179	{
	180	Error("FillHits","Hit not inside EMCAL!");
	181	mHit->Dump();
	182	}
	183	continue;
	184	}
	185	fInputObject.AddEnergyToDigit(geom->TowerIndexFromEtaPhi(eta,180.0/TMath::Pi()phi)-1,Int_t(1e7etH));
	186	} // Hit Loop
b3c7d6bc	187	// } // Track Loop
f7d5860b	188
	189
	190	}
44f59d68	191	void AliEMCALJetFinderInputSimPrep::FillTracks()
f7d5860b	192	{
44f59d68	193	// Fill from particles simulating a TPC to input object from simulation
f7d5860b	194
	195	if (fDebug > 1) Info("FillTracks","Beginning FillTracks");
	196
	197	TParticlePDG* pdgP = 0;
44f59d68	198	TParticle *mPart;
f7d5860b	199	Int_t npart = (gAlice->GetHeader())->GetNprimary();
b3c7d6bc	200	if (fDebug > 1) Info("FillTracks","Header says there are %i primaries",npart);
b3c7d6bc	201	AliEMCALGetter *gime = AliEMCALGetter::Instance();
f7d5860b	202	Float_t bfield,rEMCAL;
b3c7d6bc	203	/* for (Int_t i =0; i<gime->NPrimaries(); i++)
	204	{
	205	cout<<gime->Primary(i)->GetFirstMother()<<endl;
	206	if (gime->Primary(i)->GetFirstMother() == -1) {
	207	counterforus++;
	208	gime->Primary(i)->Dump();
	209	}
	210	}*/
f7d5860b	211
	212	if (fDebug > 1) Info("FillTracks","Defining the geometry");
	213
	214	AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
	215	AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
	216	fEtaMax = geom->GetArm1EtaMax();
	217	fEtaMin = geom->GetArm1EtaMin();
	218	fPhiMax = TMath::Pi()/180.0*geom->GetArm1PhiMax();
	219	fPhiMin = TMath::Pi()/180.0*geom->GetArm1PhiMin();
	220
	221	if (fDebug > 1) Info("FillTracks","Starting particle loop");
	222
b3c7d6bc	223	if (fDebug > 1) Info("FillTracks","Particle loop of %i",npart);
f7d5860b	224	for (Int_t part = 0; part < npart; part++) {
44f59d68	225	mPart = gAlice->GetMCApp()->Particle(part);
f7d5860b	226	//if (part%10) gObjectTable->Print();
44f59d68	227	pdgP = mPart->GetPDG();
f7d5860b	228
	229	if (fDebug > 5) Info("FillTracks","Checking if track is a primary");
	230
	231	if (fFileType == kPythia) {
44f59d68	232	if (mPart->GetStatusCode() != 1) continue;
f7d5860b	233	} else if (fFileType == kHijing) {
44f59d68	234	if (mPart->GetFirstDaughter() >= 0 && mPart->GetFirstDaughter() < npart) continue;
f7d5860b	235	}
f7d5860b	236
44f59d68	237	if (fDebug > 15) Info("FillTracks","Checking if track (eta - %f, phi - %f) is in acceptance",mPart->Eta(),mPart->Phi());
f7d5860b	238	if (fDebug > 10) Info("FillTracks","Checking if EMCAL acceptance ( %f < eta < %f, %f < phi < %f) is in acceptance",fEtaMin,fEtaMax,fPhiMin,fPhiMax);
f7d5860b	239
b3c7d6bc	240	if ((!fPythiaComparison)&&(mPart->Eta() > fEtaMax \|\| mPart->Eta() < fEtaMin)) continue;
b3c7d6bc	241	if ((!fPythiaComparison)&&(mPart->Phi() > fPhiMax \|\| mPart->Phi() < fPhiMin)) continue;
f7d5860b	242
	243	/*
	244	{kProton, kProtonBar, kElectron, kPositron,
	245	kNuE, kNuEBar, kGamma, kNeutron, kNeutronBar,
	246	kMuonPlus, kMuonMinus, kK0Long , kPiPlus, kPiMinus,
	247	kKPlus, kKMinus, kLambda0, kLambda0Bar, kK0Short,
	248	kSigmaMinus, kSigmaPlus, kSigma0, kPi0, kK0, kK0Bar,
	249	0,kNuMu,kNuMuBar
	250	*/
	251
	252	if (fDebug > 15) Info("FillTracks","Checking if track is rejected");
	253
	254	if ((fSmearType == kEfficiency \|\|
	255	fSmearType == kSmearEffic) &&
	256	pdgP->Charge()!=0) {
	257	if (AliEMCALFast::RandomReject(fEfficiency)) {
	258	continue;
	259	}
	260	}
	261
	262	bfield = gAlice->Field()->SolenoidField();
	263	rEMCAL = AliEMCALGeometry::GetInstance()->GetIPDistance();
44f59d68	264	Float_t rB = 3335.6 * mPart->Pt() / bfield; // [cm] (case of \|charge\|=1)
f7d5860b	265	if (2.*rB < rEMCAL) continue; // track curls away
	266
	267	//if (part%10) gObjectTable->Print();
	268	switch(fTrackType)
	269	{
	270
	271	case kAll: // All Stable particles to be included
	272	if (fDebug > 5) Info("FillTracks","Storing track");
	273	if (fSmearType == kSmear \|\|
	274	fSmearType == kSmearEffic ){
44f59d68	275	Smear(mPart);/*
f7d5860b	276	TParticle *tmp = Smear(MPart);
	277	fInputObject.AddTrack(Smear(MPart));
	278	delete tmp;*/
	279	}else{
44f59d68	280	fInputObject.AddTrack(*mPart);
f7d5860b	281	}
	282	break;
	283	case kEM: // All Electromagnetic particles
44f59d68	284	if (mPart->GetPdgCode() == kElectron \|\|
	285	mPart->GetPdgCode() == kMuonPlus \|\|
	286	mPart->GetPdgCode() == kMuonMinus \|\|
	287	mPart->GetPdgCode() == kPositron ){
f7d5860b	288	if (fDebug > 5) Info("FillTracks","Storing electron or positron");
	289	if (fSmearType == kSmear \|\|
	290	fSmearType == kSmearEffic ){
44f59d68	291	Smear(mPart);/*
f7d5860b	292	TParticle *tmp = Smear(MPart);
	293	fInputObject.AddTrack(tmp);
	294	delete tmp;*/
	295	}else{
44f59d68	296	fInputObject.AddTrack(*mPart);
f7d5860b	297	}
f7d5860b	298	}
44f59d68	299	if ( mPart->GetPdgCode() == kGamma ){
44f59d68	300	fInputObject.AddTrack(*mPart);
f7d5860b	301	if (fDebug > 5) Info("FillTracks","Storing gamma");
	302	}
	303
	304	break;
	305	case kCharged: // All Particles with non-zero charge
	306	if (pdgP->Charge() != 0) {
	307	if (fDebug > 5) Info("FillTracks","Storing charged track");
	308	if (fSmearType == kSmear \|\|
	309	fSmearType == kSmearEffic ){
44f59d68	310	Smear(mPart);/*
f7d5860b	311	TParticle *tmp = Smear(MPart);
	312	fInputObject.AddTrack(tmp);
	313	delete tmp;*/
	314	}else{
44f59d68	315	fInputObject.AddTrack(*mPart);
f7d5860b	316	}
	317	}
	318	break;
	319	case kNeutral: // All particles with no charge
	320	if (pdgP->Charge() == 0){
44f59d68	321	fInputObject.AddTrack(*mPart);
f7d5860b	322	if (fDebug > 5) Info("FillTracks","Storing neutral");
	323	}
	324	break;
	325	case kHadron: //All hadrons
44f59d68	326	if (mPart->GetPdgCode() != kElectron &&
	327	mPart->GetPdgCode() != kPositron &&
	328	mPart->GetPdgCode() != kMuonPlus &&
	329	mPart->GetPdgCode() != kMuonMinus &&
	330	mPart->GetPdgCode() != kGamma )
f7d5860b	331	{
	332	if (fDebug > 5) Info("FillTracks","Storing hadron");
	333	if (pdgP->Charge() == 0){
44f59d68	334	fInputObject.AddTrack(*mPart);
f7d5860b	335	}else{
	336	if (fSmearType == kSmear \|\|
	337	fSmearType == kSmearEffic ){
44f59d68	338	Smear(mPart);/*
f7d5860b	339	TParticle *tmp = Smear(MPart);
	340	fInputObject.AddTrack(tmp);
	341	delete tmp;*/
	342	}else{
44f59d68	343	fInputObject.AddTrack(*mPart);
f7d5860b	344	}
	345	}
	346	}
	347	break;
	348	case kChargedHadron: // only charged hadrons
44f59d68	349	if (mPart->GetPdgCode() != kElectron &&
	350	mPart->GetPdgCode() != kPositron &&
	351	mPart->GetPdgCode() != kGamma &&
	352	mPart->GetPdgCode() != kMuonPlus &&
	353	mPart->GetPdgCode() != kMuonMinus &&
f7d5860b	354	pdgP->Charge() != 0 ){
	355	if (fDebug > 5) Info("FillTracks","Storing charged hadron");
	356	if (fSmearType == kSmear \|\|
	357	fSmearType == kSmearEffic ){
44f59d68	358	Smear(mPart);/*
f7d5860b	359	TParticle *tmp = Smear(MPart);
	360	fInputObject.AddTrack(tmp);
	361	delete tmp;*/
	362	}else{
44f59d68	363	fInputObject.AddTrack(*mPart);
f7d5860b	364	}
	365	}
	366	break;
	367	case kNoTracks:
	368	break;
	369	default:
	370	break;
44f59d68	371	delete mPart;
f7d5860b	372	} //end of switch
	373	// Info("FillTracks","After Particle Storage");
	374	//if (part%10) gObjectTable->Print();
	375
	376	} //end of particle loop
	377	//gObjectTable->Print();
	378	}
	379
44f59d68	380	void AliEMCALJetFinderInputSimPrep::FillPartons()
f7d5860b	381	{
44f59d68	382	// Fill partons to
44f59d68	383	// input object from simulation
f7d5860b	384	if (fDebug > 1) Info("FillPartons","Beginning FillPartons");
	385
	386	AliGenEventHeader* evHeader = ((AliHeader*)(gAlice->GetHeader()))->GenEventHeader();
	387	Float_t triggerJetValues[4];
	388	AliEMCALParton tempParton;
	389
	390	if (fFileType == kPythia)
	391	{
	392	Int_t ntriggerjets = ((AliGenPythiaEventHeader*)evHeader)->NTriggerJets();
	393	if (fDebug > 1) Info("FillPartons","Number of trigger jets --> %i",ntriggerjets);
	394	for (Int_t jetcounter = 0 ; jetcounter < ntriggerjets; jetcounter++){
	395	((AliGenPythiaEventHeader*)evHeader)->TriggerJet(jetcounter,triggerJetValues);
	396	TLorentzVector tempLParton;
	397	tempLParton.SetPxPyPzE(triggerJetValues[0],triggerJetValues[1],triggerJetValues[2],triggerJetValues[3]);
	398	tempParton.SetEnergy(tempLParton.Energy());
	399	tempParton.SetEta(tempLParton.Eta());
	400	tempParton.SetPhi(tempLParton.Phi());
	401	FillPartonTracks(&tempParton);
	402	fInputObject.AddParton(&tempParton);
	403	}
	404
	405	}
	406	}
	407
44f59d68	408	void AliEMCALJetFinderInputSimPrep::FillParticles()
f7d5860b	409	{
44f59d68	410	// Fill particles to input object from simulation
f7d5860b	411	if (fDebug > 1) Info("FillParticles","Beginning FillParticles");
	412
	413	Int_t npart = (gAlice->GetHeader())->GetNprimary();
	414	TParticlePDG* pdgP = 0;
	415
	416	for (Int_t part = 0; part < npart; part++) {
44f59d68	417	TParticle *mPart = gAlice->GetMCApp()->Particle(part);
44f59d68	418	pdgP = mPart->GetPDG();
f7d5860b	419
	420	if (fDebug > 10) Info("FillParticles","Checking if particle is a primary");
	421
	422	if (fFileType == kPythia) {
44f59d68	423	if (mPart->GetStatusCode() != 1) continue;
f7d5860b	424	} else if (fFileType == kHijing) {
44f59d68	425	if (mPart->GetFirstDaughter() >= 0 && mPart->GetFirstDaughter() < npart) continue;
f7d5860b	426	}
	427
	428	if (fDebug > 10) Info("FillParticles","Checking if particle is in acceptance");
	429
b3c7d6bc	430	if ((!fPythiaComparison)&&(mPart->Eta() > fEtaMax \|\| mPart->Eta() < fEtaMin)) continue;
b3c7d6bc	431	if ((!fPythiaComparison)&&(mPart->Phi() > fPhiMax \|\| mPart->Phi() < fPhiMin)) continue;
f7d5860b	432
	433
	434	/*
	435	{kProton, kProtonBar, kElectron, kPositron,
	436	kNuE, kNuEBar, kGamma, kNeutron, kNeutronBar,
	437	kMuonPlus, kMuonMinus, kK0Long , kPiPlus, kPiMinus,
	438	kKPlus, kKMinus, kLambda0, kLambda0Bar, kK0Short,
	439	kSigmaMinus, kSigmaPlus, kSigma0, kPi0, kK0, kK0Bar,
	440	0,kNuMu,kNuMuBar
	441	*/
	442	switch(fTrackType)
	443	{
	444
	445	case kAll: // All Stable particles to be included
	446	if (fDebug > 5) Info("FillParticles","Storing particle");
44f59d68	447	fInputObject.AddParticle(mPart);
f7d5860b	448	break;
f7d5860b	449	case kEM: // All Electromagnetic particles
44f59d68	450	if (mPart->GetPdgCode() == kElectron \|\|
	451	mPart->GetPdgCode() == kPositron \|\|
	452	mPart->GetPdgCode() == kGamma){
f7d5860b	453	if (fDebug > 5) Info("FillParticles","Storing electromagnetic particle");
44f59d68	454	fInputObject.AddParticle(mPart);
f7d5860b	455	}
	456	break;
	457	case kCharged: // All Particles with non-zero charge
	458	if (pdgP->Charge() != 0) {
	459	if (fDebug > 5) Info("FillParticles","Storing charged particle");
44f59d68	460	fInputObject.AddParticle(mPart);
f7d5860b	461	}
	462	break;
	463	case kNeutral: // All particles with no charge
	464	if (pdgP->Charge() == 0){
	465	if (fDebug > 5) Info("FillParticles","Storing neutral particle");
44f59d68	466	fInputObject.AddParticle(mPart);
f7d5860b	467	}
	468	break;
	469	case kHadron: //All hadrons
44f59d68	470	if (mPart->GetPdgCode() != kElectron &&
	471	mPart->GetPdgCode() != kPositron &&
	472	mPart->GetPdgCode() != kGamma )
f7d5860b	473	{
	474
	475	if (fDebug > 5) Info("FillParticles","Storing hadron");
44f59d68	476	fInputObject.AddParticle(mPart);
f7d5860b	477	}
	478	break;
	479	case kChargedHadron: // only charged hadrons
44f59d68	480	if (mPart->GetPdgCode() != kElectron &&
	481	mPart->GetPdgCode() != kPositron &&
	482	mPart->GetPdgCode() != kGamma &&
f7d5860b	483	pdgP->Charge() != 0 ){
f7d5860b	484	if (fDebug > 5) Info("FillParticles","Storing charged hadron");
44f59d68	485	fInputObject.AddParticle(mPart);
f7d5860b	486	}
	487	break;
	488	case kNoTracks:
	489	break;
	490	default:
	491	break;
	492	} //end of switch
	493	}// end of loop over particles
	494
	495	}
44f59d68	496	void AliEMCALJetFinderInputSimPrep::FillDigits()
f7d5860b	497	{
44f59d68	498	// Fill digits to input object
f7d5860b	499
	500	}
	501
	502	void AliEMCALJetFinderInputSimPrep::Smear(TParticle *particle)
	503	{
44f59d68	504	// Smear particle momentum
f7d5860b	505	if (fDebug > 5) Info("Smear","Beginning Smear");
	506
	507	Float_t tmp = AliEMCALFast::SmearMomentum(1,particle->P());
	508	Float_t tmpt = (tmp/particle->P()) * particle->Pt();
	509	if (fDebug > 10) Info("Smear","Smearing changed P from %f to %f",particle->Pt(),tmpt);
	510	TLorentzVector tmpvector;
	511	tmpvector.SetPtEtaPhiM(tmpt,particle->Eta(), particle->Phi(), particle->GetMass());
	512	// create a new particle with smearing momentum - and no daughter or parent information and the same
	513	// vertex
	514
	515	TParticle tmparticle(particle->GetPdgCode(), 1, 0, 0,0,0,
	516	tmpvector.Px() , tmpvector.Py(), tmpvector.Pz(), tmpvector.Energy(),
	517	particle->Vx(), particle->Vy(), particle->Vz(), particle->T());
	518	fInputObject.AddTrack(tmparticle);
	519	return;
	520	}
	521
	522	void AliEMCALJetFinderInputSimPrep::FillPartonTracks(AliEMCALParton *parton)
	523	{
44f59d68	524	// Populate parton tracks for input distributions
f7d5860b	525	if (fDebug>1) Info("FillPartonTracks","Beginning FillPartonTracks");
	526	Int_t npart = (gAlice->GetHeader())->GetNprimary();
	527	Int_t ntracks =0;
	528	TParticlePDG *getpdg;
	529	TParticle *tempPart;
	530	for (Int_t part = 0; part < npart; part++)
	531	{
5d12ce38	532	tempPart = gAlice->GetMCApp()->Particle(part);
f7d5860b	533	if (tempPart->GetStatusCode() != 1) continue;
b3c7d6bc	534	if ((!fPythiaComparison)&&(tempPart->Eta() > fEtaMax \|\| tempPart->Eta() < fEtaMin \|\|
b3c7d6bc	535	tempPart->Phi() > fPhiMax \|\| tempPart->Phi() < fPhiMin) ){
f7d5860b	536	if (fDebug>10) Info("FillPartonTracks","Excluding particle not pointing at the EMCAL");
	537	continue;
	538	}
	539	getpdg = tempPart->GetPDG();
	540	if (getpdg->Charge() == 0.0 ) {
	541	if (fDebug>10) Info("FillPartonTracks","Excluding a neutral particle");
	542	continue;
	543	}
	544	if ( (parton->Eta() - tempPart->Eta())*(parton->Eta() - tempPart->Eta()) +
	545	(parton->Phi() - tempPart->Phi())*(parton->Phi() - tempPart->Phi()) < 1.0 ) ntracks++;
	546	}
44f59d68	547	Float_t *energy = new Float_t[ntracks];
	548	Float_t *eta = new Float_t[ntracks];
	549	Float_t *phi = new Float_t[ntracks];
	550	Int_t *pdg = new Int_t[ntracks];
f7d5860b	551	ntracks=0;
	552	for (Int_t part = 0; part < npart; part++)
	553	{
5d12ce38	554	tempPart = gAlice->GetMCApp()->Particle(part);
f7d5860b	555	if (tempPart->GetStatusCode() != 1) continue;
b3c7d6bc	556	if ((!fPythiaComparison)&&(tempPart->Eta() > fEtaMax \|\| tempPart->Eta() < fEtaMin \|\|
b3c7d6bc	557	tempPart->Phi() > fPhiMax \|\| tempPart->Phi() < fPhiMin) ){
f7d5860b	558	if (fDebug>10) Info("FillPartonTracks","Excluding particle not pointing at the EMCAL");
	559	continue;
	560	}
	561	if (tempPart->GetStatusCode() != 1) continue;
	562	getpdg = tempPart->GetPDG();
	563	if (getpdg->Charge() == 0.0 ) {
	564	if (fDebug>10) Info("FillPartonTracks","Excluding a neutral particle");
	565	continue;
	566	}
	567	if ( (parton->Eta() - tempPart->Eta())*(parton->Eta() - tempPart->Eta()) +
	568	(parton->Phi() - tempPart->Phi())*(parton->Phi() - tempPart->Phi()) < 1.0 )
	569	{
44f59d68	570	energy[ntracks] = tempPart->Pt();
	571	eta[ntracks] = tempPart->Eta();
	572	phi[ntracks] = tempPart->Phi();
	573	pdg[ntracks] = tempPart->GetPdgCode();
f7d5860b	574	ntracks++;
	575	}
	576	}
44f59d68	577	parton->SetTrackList(ntracks,energy,eta,phi,pdg);
b3c7d6bc	578	delete[] energy;
	579	delete[] eta;
	580	delete[] phi;
	581	delete[] pdg;
f7d5860b	582	}
f7d5860b	583