**************************************************************************/
-/*
-
-$Log$
-Revision 1.1.1.1 2003/05/29 18:56:43 horner
-Initial import - Mark
-
-
-*/
+/* $Id$ */
//_________________________________________________________________________
// Class for JetFinder Input preparation from simulated data
#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;
void AliEMCALJetFinderInputSimPrep::Reset(AliEMCALJetFinderResetType_t resettype)
{
+ // Method to reset data
if (fDebug > 1) Info("Reset","Beginning Reset");
switch (resettype){
Int_t AliEMCALJetFinderInputSimPrep::FillFromFile(TString *filename, AliEMCALJetFinderFileType_t filetype,Int_t EventNumber)
{
- // gObjectTable->Print();
- // Test that file exists - Getter doesn't like bogus filenames
if (fDebug > 1) Info("FillFromFile","Beginning FillFromFile");
fFileType = filetype;
- TFile file(filename->Data());
- if (!file.IsOpen()){
- Error("FillFromFile","Could not open file in FillFromFile");
- return -1;
- }
- file.Close();
- /*
- gAlice = static_cast<AliRun *>(file.Get("gAlice"));
- AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
- if (gAlice->GetEvent(EventNumber) < 0){
- Error("FillFromFile","Could not open event in FillFromFile");
- file.Close();
- return -1;
- }*/
AliEMCALGetter *gime = AliEMCALGetter::Instance(filename->Data());
- gime->Event(EventNumber) ;
-
+ 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();
-// FillParticles();
}
- //gime->CloseFile();
-// gObjectTable->Print();
- delete gime;
+/* gime->EmcalLoader()->UnloadRecParticles();
+ gime->EmcalLoader()->UnloadTracks();
+ gime->Reset();*/
return 0;
}
-void AliEMCALJetFinderInputSimPrep::FillHits() // Fill from the hits to input object from simulation
+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());
+// AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
+// Info("AliEMCALJetFinderInputSimPrep","Title of the geometry is %s",pEMCAL->GetTitle());
AliEMCALGeometry* geom = AliEMCALGetter::Instance()->EMCALGeometry();
-// Float_t samplingF = geom->GetSampling();
- Float_t samplingF = 11.6;
- Info("AliEMCALJetFinderInputSimPrep","Sampling fraction is %f",samplingF);
- TTree *treeH = AliEMCALGetter::Instance()->TreeH();
- Int_t ntracks = (Int_t) treeH->GetEntries();
+ 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++) {
+/* for (Int_t track=0; track<ntracks;track++) {
gAlice->ResetHits();
- nbytes += treeH->GetEvent(track);
+ nbytes += treeH->GetEvent(track);*/
//
// Loop over hits
//
- for(mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1);
- mHit;
- mHit=(AliEMCALHit*) pEMCAL->NextHit())
+ 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 theta = TMath::ATan2(r,z);
Float_t eta = -TMath::Log(TMath::Tan(theta/2.));
Float_t phi = TMath::ATan2(y,x);
- etH = samplingF*eloss*TMath::Sin(theta);
+ 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)
{
}
fInputObject.AddEnergyToDigit(geom->TowerIndexFromEtaPhi(eta,180.0/TMath::Pi()*phi)-1,Int_t(1e7*etH));
} // Hit Loop
- } // Track Loop
+// } // Track Loop
}
-void AliEMCALJetFinderInputSimPrep::FillTracks() // Fill from particles simulating a TPC to input object from simulation
+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;
+ 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");
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->Particle(part);
+ mPart = gAlice->GetMCApp()->Particle(part);
//if (part%10) gObjectTable->Print();
- pdgP = MPart->GetPDG();
+ pdgP = mPart->GetPDG();
if (fDebug > 5) Info("FillTracks","Checking if track is a primary");
if (fFileType == kPythia) {
- if (MPart->GetStatusCode() != 1) continue;
+ if (mPart->GetStatusCode() != 1) continue;
} else if (fFileType == kHijing) {
- if (MPart->GetFirstDaughter() >= 0 && MPart->GetFirstDaughter() < npart) continue;
+ 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 > 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 (MPart->Eta() > fEtaMax || MPart->Eta() < fEtaMin) continue;
- if (MPart->Phi() > fPhiMax || MPart->Phi() < fPhiMin ) continue;
+ if ((!fPythiaComparison)&&(mPart->Eta() > fEtaMax || mPart->Eta() < fEtaMin)) continue;
+ if ((!fPythiaComparison)&&(mPart->Phi() > fPhiMax || mPart->Phi() < fPhiMin)) continue;
/*
{kProton, kProtonBar, kElectron, kPositron,
bfield = gAlice->Field()->SolenoidField();
rEMCAL = AliEMCALGeometry::GetInstance()->GetIPDistance();
- Float_t rB = 3335.6 * MPart->Pt() / bfield; // [cm] (case of |charge|=1)
+ 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();
if (fDebug > 5) Info("FillTracks","Storing track");
if (fSmearType == kSmear ||
fSmearType == kSmearEffic ){
- Smear(MPart);/*
+ Smear(mPart);/*
TParticle *tmp = Smear(MPart);
fInputObject.AddTrack(Smear(MPart));
delete tmp;*/
}else{
- fInputObject.AddTrack(*MPart);
+ fInputObject.AddTrack(*mPart);
}
break;
case kEM: // All Electromagnetic particles
- if (MPart->GetPdgCode() == kElectron ||
- MPart->GetPdgCode() == kPositron ){
+ 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);/*
+ Smear(mPart);/*
TParticle *tmp = Smear(MPart);
fInputObject.AddTrack(tmp);
delete tmp;*/
}else{
- fInputObject.AddTrack(*MPart);
+ fInputObject.AddTrack(*mPart);
}
}
- if ( MPart->GetPdgCode() == kGamma ){
- fInputObject.AddTrack(*MPart);
+ if ( mPart->GetPdgCode() == kGamma ){
+ fInputObject.AddTrack(*mPart);
if (fDebug > 5) Info("FillTracks","Storing gamma");
}
if (fDebug > 5) Info("FillTracks","Storing charged track");
if (fSmearType == kSmear ||
fSmearType == kSmearEffic ){
- Smear(MPart);/*
+ Smear(mPart);/*
TParticle *tmp = Smear(MPart);
fInputObject.AddTrack(tmp);
delete tmp;*/
}else{
- fInputObject.AddTrack(*MPart);
+ fInputObject.AddTrack(*mPart);
}
}
break;
case kNeutral: // All particles with no charge
if (pdgP->Charge() == 0){
- fInputObject.AddTrack(*MPart);
+ fInputObject.AddTrack(*mPart);
if (fDebug > 5) Info("FillTracks","Storing neutral");
}
break;
case kHadron: //All hadrons
- if (MPart->GetPdgCode() != kElectron &&
- MPart->GetPdgCode() != kPositron &&
- MPart->GetPdgCode() != kGamma )
+ 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);
+ fInputObject.AddTrack(*mPart);
}else{
if (fSmearType == kSmear ||
fSmearType == kSmearEffic ){
- Smear(MPart);/*
+ Smear(mPart);/*
TParticle *tmp = Smear(MPart);
fInputObject.AddTrack(tmp);
delete tmp;*/
}else{
- fInputObject.AddTrack(*MPart);
+ fInputObject.AddTrack(*mPart);
}
}
}
break;
case kChargedHadron: // only charged hadrons
- if (MPart->GetPdgCode() != kElectron &&
- MPart->GetPdgCode() != kPositron &&
- MPart->GetPdgCode() != kGamma &&
+ 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);/*
+ Smear(mPart);/*
TParticle *tmp = Smear(MPart);
fInputObject.AddTrack(tmp);
delete tmp;*/
}else{
- fInputObject.AddTrack(*MPart);
+ fInputObject.AddTrack(*mPart);
}
}
break;
break;
default:
break;
- delete MPart;
+ delete mPart;
} //end of switch
// Info("FillTracks","After Particle Storage");
//if (part%10) gObjectTable->Print();
//gObjectTable->Print();
}
-void AliEMCALJetFinderInputSimPrep::FillPartons() // Fill partons to input object from simulation
+void AliEMCALJetFinderInputSimPrep::FillPartons()
{
+ // Fill partons to
+ // input object from simulation
if (fDebug > 1) Info("FillPartons","Beginning FillPartons");
AliGenEventHeader* evHeader = ((AliHeader*)(gAlice->GetHeader()))->GenEventHeader();
}
}
-void AliEMCALJetFinderInputSimPrep::FillParticles() // Fill particles to input object from simulation
+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->Particle(part);
- pdgP = MPart->GetPDG();
+ 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;
+ if (mPart->GetStatusCode() != 1) continue;
} else if (fFileType == kHijing) {
- if (MPart->GetFirstDaughter() >= 0 && MPart->GetFirstDaughter() < npart) continue;
+ if (mPart->GetFirstDaughter() >= 0 && mPart->GetFirstDaughter() < npart) continue;
}
if (fDebug > 10) Info("FillParticles","Checking if particle is in acceptance");
- if (MPart->Eta() > fEtaMax || MPart->Eta() < fEtaMin) continue;
- if (MPart->Phi() > fPhiMax || MPart->Phi() < fPhiMin ) continue;
+ if ((!fPythiaComparison)&&(mPart->Eta() > fEtaMax || mPart->Eta() < fEtaMin)) continue;
+ if ((!fPythiaComparison)&&(mPart->Phi() > fPhiMax || mPart->Phi() < fPhiMin)) continue;
/*
case kAll: // All Stable particles to be included
if (fDebug > 5) Info("FillParticles","Storing particle");
- fInputObject.AddParticle(MPart);
+ fInputObject.AddParticle(mPart);
break;
case kEM: // All Electromagnetic particles
- if (MPart->GetPdgCode() == kElectron ||
- MPart->GetPdgCode() == kPositron ||
- MPart->GetPdgCode() == kGamma){
+ if (mPart->GetPdgCode() == kElectron ||
+ mPart->GetPdgCode() == kPositron ||
+ mPart->GetPdgCode() == kGamma){
if (fDebug > 5) Info("FillParticles","Storing electromagnetic particle");
- fInputObject.AddParticle(MPart);
+ 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);
+ 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);
+ fInputObject.AddParticle(mPart);
}
break;
case kHadron: //All hadrons
- if (MPart->GetPdgCode() != kElectron &&
- MPart->GetPdgCode() != kPositron &&
- MPart->GetPdgCode() != kGamma )
+ if (mPart->GetPdgCode() != kElectron &&
+ mPart->GetPdgCode() != kPositron &&
+ mPart->GetPdgCode() != kGamma )
{
if (fDebug > 5) Info("FillParticles","Storing hadron");
- fInputObject.AddParticle(MPart);
+ fInputObject.AddParticle(mPart);
}
break;
case kChargedHadron: // only charged hadrons
- if (MPart->GetPdgCode() != kElectron &&
- MPart->GetPdgCode() != kPositron &&
- MPart->GetPdgCode() != kGamma &&
+ if (mPart->GetPdgCode() != kElectron &&
+ mPart->GetPdgCode() != kPositron &&
+ mPart->GetPdgCode() != kGamma &&
pdgP->Charge() != 0 ){
if (fDebug > 5) Info("FillParticles","Storing charged hadron");
- fInputObject.AddParticle(MPart);
+ fInputObject.AddParticle(mPart);
}
break;
case kNoTracks:
}// end of loop over particles
}
-void AliEMCALJetFinderInputSimPrep::FillDigits() // Fill digits to input object
+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());
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;
TParticle *tempPart;
for (Int_t part = 0; part < npart; part++)
{
- tempPart = gAlice->Particle(part);
+ tempPart = gAlice->GetMCApp()->Particle(part);
if (tempPart->GetStatusCode() != 1) continue;
- if (tempPart->Eta() > fEtaMax || tempPart->Eta() < fEtaMin ||
- tempPart->Phi() > fPhiMax || tempPart->Phi() < fPhiMin ){
+ 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 ( (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];
+ 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->Particle(part);
+ tempPart = gAlice->GetMCApp()->Particle(part);
if (tempPart->GetStatusCode() != 1) continue;
- if (tempPart->Eta() > fEtaMax || tempPart->Eta() < fEtaMin ||
- tempPart->Phi() > fPhiMax || tempPart->Phi() < fPhiMin ){
+ 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 ( (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();
+ energy[ntracks] = tempPart->Pt();
+ eta[ntracks] = tempPart->Eta();
+ phi[ntracks] = tempPart->Phi();
+ pdg[ntracks] = tempPart->GetPdgCode();
ntracks++;
}
}
- parton->SetTrackList(ntracks,Energy,Eta,Phi,Pdg);
+ parton->SetTrackList(ntracks,energy,eta,phi,pdg);
+ delete[] energy;
+ delete[] eta;
+ delete[] phi;
+ delete[] pdg;
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff