#include <TDatabasePDG.h>
#include <TRandom.h>
// From AliRoot ...
-#include "AliJet.h"
+#include "AliAODJet.h"
+#include "AliPDG.h"
#include "AliJetKineReaderHeader.h"
#include "AliJetKineReader.h"
-#include "AliRunLoader.h"
+#include "AliMCEventHandler.h"
+#include "AliMCEvent.h"
#include "AliStack.h"
#include "AliHeader.h"
#include "AliGenPythiaEventHeader.h"
+#include "AliLog.h"
ClassImp(AliJetKineReader)
AliJetKineReader::AliJetKineReader():
AliJetReader(),
- fRunLoader(0x0),
- fAliHeader(0x0)
+ fAliHeader(0),
+ fMCEvent(0),
+ fGenJets(0)
{
// Default constructor
}
AliJetKineReader::~AliJetKineReader()
{
// Destructor
- delete fAliHeader;
+ if (fAliHeader) {
+ delete fAliHeader;
+ fAliHeader = 0;
+ }
}
//____________________________________________________________________________
void AliJetKineReader::OpenInputFiles()
{
// Opens the input file using the run loader
- const char* dirName = fReaderHeader->GetDirectory();
- char path[256];
- sprintf(path, "%s/galice.root",dirName);
- fRunLoader = AliRunLoader::Open(path);
- fRunLoader->LoadKinematics();
- fRunLoader->LoadHeader();
-
- // get number of events
- Int_t nMax = fRunLoader->GetNumberOfEvents();
- printf("\nTotal number of events = %d", nMax);
-
- // set number of events in header
- if (fReaderHeader->GetLastEvent() == -1)
- fReaderHeader->SetLastEvent(nMax);
- else {
- Int_t nUsr = fReaderHeader->GetLastEvent();
- fReaderHeader->SetLastEvent(TMath::Min(nMax, nUsr));
- }
+ // Not used anymore
}
//____________________________________________________________________________
-Bool_t AliJetKineReader::FillMomentumArray(Int_t event)
+Bool_t AliJetKineReader::FillMomentumArray()
{
// Fill momentum array for event
- char path[256];
- const char* dirName = fReaderHeader->GetDirectory();
- const char* bgdirName = fReaderHeader->GetBgDirectory();
- Int_t goodTrack = 0;
- // Clear array
- // temporary storage of signal and cut flags
- Int_t* sflag = new Int_t[50000];
- Int_t* cflag = new Int_t[50000];
- Int_t spb = 0;
+ Int_t goodTrack = 0;
+ // Clear array
+ // temporary storage of signal and cut flags
+ Int_t* sflag = new Int_t[50000];
+ Int_t* cflag = new Int_t[50000];
- ClearArray();
- for (Int_t i = 0; i < 2; i++) {
- if (i == 1) {
-// Pythia Events
- if (fRunLoader) delete fRunLoader;
- sprintf(path, "%s/galice.root", dirName);
- fRunLoader = AliRunLoader::Open(path);
- fRunLoader->LoadKinematics();
- fRunLoader->LoadHeader();
- // Get event from runloader
- fRunLoader->GetEvent(event);
- } else {
-// Background events
- if ((spb = fReaderHeader->GetSignalPerBg()) > 0) {
- if (fRunLoader) delete fRunLoader;
- sprintf(path, "%s/galice.root", bgdirName);
- fRunLoader = AliRunLoader::Open(path);
- fRunLoader->LoadKinematics();
- fRunLoader->LoadHeader();
- // Get event from runloader
- fRunLoader->GetEvent(event / spb);
- } else {
- continue;
- }
- }
-
- // Get the stack
- AliStack* stack = fRunLoader->Stack();
- // Number of primaries
- Int_t nt = stack->GetNprimary();
-
- // Get cuts set by user and header
- Double_t ptMin = ((AliJetKineReaderHeader*) fReaderHeader)->GetPtCut();
- Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
- Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
- fAliHeader = fRunLoader->GetHeader();
+ ClearArray();
+ // Get the stack
+ AliStack* stack = fMCEvent->Stack();
+ // Number of primaries
+ Int_t nt = stack->GetNprimary();
+ // Get cuts set by user and header
+ Double_t ptMin = ((AliJetKineReaderHeader*) fReaderHeader)->GetPtCut();
+ Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
+ Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
+ fAliHeader = fMCEvent->Header();
- TLorentzVector p4;
- // Loop over particles
- for (Int_t it = 0; it < nt; it++) {
- TParticle *part = stack->Particle(it);
- Int_t status = part->GetStatusCode();
- Int_t pdg = TMath::Abs(part->GetPdgCode());
- Float_t pt = part->Pt();
+
+ TLorentzVector p4;
+ // Loop over particles
+ for (Int_t it = 0; it < nt; it++) {
+ TParticle *part = stack->Particle(it);
+ Int_t status = part->GetStatusCode();
+ Int_t pdg = TMath::Abs(part->GetPdgCode());
+ Float_t pt = part->Pt();
- // Skip non-final state particles, neutrinos
- // and particles with pt < pt_min
- if ((status != 1)
- || (pdg == 12 || pdg == 14 || pdg == 16)) continue;
-
- Float_t p = part->P();
- Float_t p0 = p;
- Float_t eta = part->Eta();
- Float_t phi = part->Phi();
-
- // Fast simulation of TPC if requested
- if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimTPC()) {
- // Charged particles only
- Float_t charge =
- TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
- if (charge == 0) continue;
- // Simulate efficiency
- if (!Efficiency(p0, eta, phi)) continue;
- // Simulate resolution
- p = SmearMomentum(4, p0);
- } // End fast TPC
+ // Skip non-final state particles, neutrinos
+ // and particles with pt < pt_min
+ if ((status != 1)
+ || (pdg == 12 || pdg == 14 || pdg == 16)) continue;
+
+ Float_t p = part->P();
+ Float_t p0 = p;
+ Float_t eta = part->Eta();
+ Float_t phi = part->Phi();
+
+
+ if (((AliJetKineReaderHeader*)fReaderHeader)->ChargedOnly()) {
+ // Charged particles only
+ Float_t charge =
+ TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
+ if (charge == 0) continue;
+ } // End charged only
+
+
+ // Fast simulation of TPC if requested
+ if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimTPC()) {
+ // Charged particles only
+ Float_t charge =
+ TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
+ if (charge == 0) continue;
+ // Simulate efficiency
+ if (!Efficiency(p0, eta, phi)) continue;
+ // Simulate resolution
+ p = SmearMomentum(4, p0);
+ } // End fast TPC
// Fast simulation of EMCAL if requested
- if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimEMCAL()) {
- Float_t charge =
- TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
- // Charged particles only
- if (charge != 0){
- // Simulate efficiency
- if (!Efficiency(p0, eta, phi)) continue;
- // Simulate resolution
- p = SmearMomentum(4, p0);
- } // end "if" charged particles
+ if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimEMCAL()) {
+ Float_t charge =
+ TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
+ // Charged particles only
+ if (charge != 0){
+ // Simulate efficiency
+ if (!Efficiency(p0, eta, phi)) continue;
+ // Simulate resolution
+ p = SmearMomentum(4, p0);
+ } // end "if" charged particles
// Neutral particles (exclude K0L, n, nbar)
- if (pdg == kNeutron || pdg == kK0Long) continue;
- } // End fast EMCAL
-
+ if (pdg == kNeutron || pdg == kK0Long) continue;
+ } // End fast EMCAL
+
// Fill momentum array
- Float_t r = p/p0;
- Float_t px = r * part->Px();
- Float_t py = r * part->Py();
- Float_t pz = r * part->Pz();
- Float_t m = part->GetMass();
- Float_t e = TMath::Sqrt(px * px + py * py + pz * pz + m * m);
- p4 = TLorentzVector(px, py, pz, e);
- if ( (p4.Eta()>etaMax) || (p4.Eta()<etaMin)) continue;
- new ((*fMomentumArray)[goodTrack]) TLorentzVector(p4);
-
- // all tracks are signal ... ?
- sflag[goodTrack]=1;
- cflag[goodTrack]=0;
- if (pt > ptMin) cflag[goodTrack]=1; // track surviving pt cut
- goodTrack++;
- }
- }
+ Float_t r = p/p0;
+ Float_t px = r * part->Px();
+ Float_t py = r * part->Py();
+ Float_t pz = r * part->Pz();
+ TParticlePDG *pPDG = part->GetPDG();
+ Float_t m = 0;
+ if(!pPDG){
+ // this is very rare...
+ // Is avoided by AliPDG::AddParticlesToPdgDataBase();
+ // but this should be called only once... (how in proof?)
+ // Calucluate mass with unsmeared momentum values
+ m = part->Energy()*part->Energy() -
+ (px * px + py * py + pz * pz)/r/r;
+ if(m>0)m = TMath::Sqrt(m);
+ else m = 0;
+ AliInfo(Form("Unknown Particle using %d calculated mass m = %3.3f",part->GetPdgCode(),m));
+
+ }
+ else{
+ m = pPDG->Mass();
+ }
+ Float_t e = TMath::Sqrt(px * px + py * py + pz * pz + m * m);
+ p4 = TLorentzVector(px, py, pz, e);
+ if ( (p4.Eta()>etaMax) || (p4.Eta()<etaMin)) continue;
+ new ((*fMomentumArray)[goodTrack]) TLorentzVector(p4);
+
+ // all tracks are signal ... ?
+ sflag[goodTrack]=1;
+ cflag[goodTrack]=0;
+ if (pt > ptMin) cflag[goodTrack] = 1; // track surviving pt cut
+ goodTrack++;
+ } // track loop
- // set the signal flags
- fSignalFlag.Set(goodTrack,sflag);
- fCutFlag.Set(goodTrack,cflag);
- printf("\nIn event %d, number of good tracks %d \n", event, goodTrack);
- return kTRUE;
+// set the signal flags
+ fSignalFlag.Set(goodTrack, sflag);
+ fCutFlag.Set(goodTrack, cflag);
+ AliInfo(Form(" Number of good tracks %d \n", goodTrack));
+
+ delete[] sflag;
+ delete[] cflag;
+
+ return kTRUE;
}
}
-Bool_t AliJetKineReader::GetGenJets(AliJet* genJets)
+TClonesArray* AliJetKineReader::GetGeneratedJets()
{
+ //
// Get generated jets from mc header
+ //
+ if (!fGenJets) {
+ fGenJets = new TClonesArray("AliAODJets", 100);
+ } else {
+ fGenJets->Clear();
+ }
+
+
AliHeader* alih = GetAliHeader();
- if (alih == 0) return kFALSE;
+ if (alih == 0) return 0;
AliGenEventHeader * genh = alih->GenEventHeader();
- if (genh == 0) return kFALSE;
+ if (genh == 0) return 0;
+
Int_t nj =((AliGenPythiaEventHeader*)genh)->NTriggerJets();
- Int_t* m = new Int_t[nj];
- Int_t* k = new Int_t[nj];
- for (Int_t i=0; i< nj; i++) {
+ for (Int_t i = 0; i < nj; i++) {
Float_t p[4];
((AliGenPythiaEventHeader*)genh)->TriggerJet(i,p);
- genJets->AddJet(p[0],p[1],p[2],p[3]);
- m[i]=1;
- k[i]=i;
+ new ((*fGenJets)[i]) AliAODJet(p[0], p[1], p[2], p[3]);
}
- genJets->SetNinput(nj);
- genJets->SetMultiplicities(m);
- genJets->SetInJet(k);
- return kTRUE;
+
+ return fGenJets;
+}
+
+void AliJetKineReader::SetInputEvent(const TObject* /*esd*/, const TObject* /*aod*/, const TObject* mc)
+{
+ // Connect the input event
+ fMCEvent = (AliMCEvent*) mc;
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff