* provided "as is" without express or implied warranty. *
**************************************************************************/
-//
+//-------------------------------------------------------------------------
// Jet kinematics reader
// MC reader for jet analysis
-// Author: Andreas Morsch
-// andreas.morsch@cern.ch
-//
+// Author: Andreas Morsch (andreas.morsch@cern.ch)
+//-------------------------------------------------------------------------
// From root ...
#include <TClonesArray.h>
#include <TPDGCode.h>
#include <TParticle.h>
#include <TParticlePDG.h>
-#include <TVector3.h>
#include <TLorentzVector.h>
#include <TSystem.h>
#include <TDatabasePDG.h>
#include <TRandom.h>
// From AliRoot ...
+#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()
+AliJetKineReader::AliJetKineReader():
+ AliJetReader(),
+ fAliHeader(0),
+ fMCEvent(0),
+ fGenJets(0)
{
- // Constructor
- fReaderHeader = 0x0;
- fRunLoader = 0x0;
- fMass = 0;
- fPdgC = 0;
+ // Default constructor
}
//____________________________________________________________________________
AliJetKineReader::~AliJetKineReader()
{
// Destructor
- delete fReaderHeader;
+ 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();
-
- 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));
- }
+ // Opens the input file using the run loader
+ // Not used anymore
}
//____________________________________________________________________________
-void AliJetKineReader::FillMomentumArray(Int_t event)
+Bool_t AliJetKineReader::FillMomentumArray()
{
-//
-// Fill momentum array for event
-//
+ // Fill momentum array for event
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();
- // Get event from runloader
- fRunLoader->GetEvent(event);
// Get the stack
- AliStack* stack = fRunLoader->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();
- fAliHeader = fRunLoader->GetHeader();
-
+ Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
+ Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
+ fAliHeader = fMCEvent->Header();
+
+
+ TLorentzVector p4;
// Loop over particles
- Int_t* flag = new Int_t[nt];
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;
- // Skip non-final state particles, neutrinos and particles with pt < pt_min
-
- if (
- (status != 1)
- || (pdg == 12 || pdg == 14 || pdg == 16)
- || (pt < ptMin)
- ) 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)->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();
+ 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);
- } // Fast TPC
-
- // Fast simulation of EMCAL if requested
+ } // End fast TPC
+
+ // Fast simulation of EMCAL if requested
if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimEMCAL()) {
- Float_t charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
- // Charged particles
- if (charge != 0) {
+ 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);
- } // charged
- // Neutral particles
- // Exclude K0L, n, nbar
+ } // end "if" charged particles
+ // Neutral particles (exclude K0L, n, nbar)
if (pdg == kNeutron || pdg == kK0Long) continue;
- } // Fast EMCAL
+ } // End fast EMCAL
- fMass = part->GetCalcMass();
- fPdgC = part->GetPdgCode();
- // Fill momentum array
+ // Fill momentum array
Float_t r = p/p0;
-// printf("smeared %13.3f %13.3f %13.3f\n", p0, p, r);
-
Float_t px = r * part->Px();
Float_t py = r * part->Py();
Float_t pz = r * part->Pz();
- Float_t m = part->GetMass();
+ 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);
- new ((*fMomentumArray)[goodTrack])
- TLorentzVector(px, py, pz, e);
+ // all tracks are signal ... ?
+ sflag[goodTrack]=1;
+ cflag[goodTrack]=0;
+ if (pt > ptMin) cflag[goodTrack] = 1; // track surviving pt cut
goodTrack++;
- }
- fSignalFlag.Set(goodTrack,flag);
- printf("\nNumber of good tracks in event %d= %d \n",event,goodTrack);
+ } // track loop
+
+// 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;
}
Float_t AliJetKineReader::SmearMomentum(Int_t ind, Float_t p)
{
-//
-// The relative momentum error, i.e.
-// (delta p)/p = sqrt (a**2 + (b*p)**2) * 10**-2,
-// where typically a = 0.75 and b = 0.16 - 0.24 depending on multiplicity
-// (the lower value is for dn/d(eta) about 2000, and the higher one for 8000)
-//
-// If we include information from TRD b will be by a factor 2/3 smaller.
-//
-// ind = 1: high multiplicity
-// ind = 2: low multiplicity
-// ind = 3: high multiplicity + TRD
-// ind = 4: low multiplicity + TRD
-
- Float_t pSmeared;
- Float_t a = 0.75;
- Float_t b = 0.12;
-
- if (ind == 1) b = 0.12;
- if (ind == 2) b = 0.08;
- if (ind == 3) b = 0.12;
- if (ind == 4) b = 0.08;
-
- Float_t sigma = p * TMath::Sqrt(a * a + b * b * p * p)*0.01;
- pSmeared = p + gRandom->Gaus(0., sigma);
- return pSmeared;
+ // The relative momentum error, i.e.
+ // (delta p)/p = sqrt (a**2 + (b*p)**2) * 10**-2,
+ // where typically a = 0.75 and b = 0.16 - 0.24 depending on multiplicity
+ // (the lower value is for dn/d(eta) about 2000, and the higher one for 8000)
+ //
+ // If we include information from TRD, b will be a factor 2/3 smaller.
+ //
+ // ind = 1: high multiplicity
+ // ind = 2: low multiplicity
+ // ind = 3: high multiplicity + TRD
+ // ind = 4: low multiplicity + TRD
+
+ Float_t pSmeared;
+ Float_t a = 0.75;
+ Float_t b = 0.12;
+
+ if (ind == 1) b = 0.12;
+ if (ind == 2) b = 0.08;
+ if (ind == 3) b = 0.12;
+ if (ind == 4) b = 0.08;
+
+ Float_t sigma = p * TMath::Sqrt(a * a + b * b * p * p)*0.01;
+ pSmeared = p + gRandom->Gaus(0., sigma);
+ return pSmeared;
}
+
+
Bool_t AliJetKineReader::Efficiency(Float_t p, Float_t /*eta*/, Float_t phi)
{
-//
-// Fast simulation of geometrical acceptance and tracking efficiency
-//
-// Tracking
- Float_t eff = 0.99;
- if (p < 0.5) eff -= (0.5-p)*0.2/0.3;
-// Geometry
- if (p > 0.5) {
- phi *= 180. / TMath::Pi();
- // Sector number 0 - 17
- Int_t isec = Int_t(phi / 20.);
- // Sector centre
- Float_t phi0 = isec * 20. + 10.;
- Float_t phir = TMath::Abs(phi-phi0);
- // 2 deg of dead space
- if (phir > 9.) eff = 0.;
- }
- if (gRandom->Rndm() > eff) {
- return kFALSE;
+ // Fast simulation of geometrical acceptance and tracking efficiency
+
+ // Tracking
+ Float_t eff = 0.99;
+ if (p < 0.5) eff -= (0.5-p)*0.2/0.3;
+ // Geometry
+ if (p > 0.5) {
+ phi *= 180. / TMath::Pi();
+ // Sector number 0 - 17
+ Int_t isec = Int_t(phi / 20.);
+ // Sector centre
+ Float_t phi0 = isec * 20. + 10.;
+ Float_t phir = TMath::Abs(phi-phi0);
+ // 2 deg of dead space
+ if (phir > 9.) eff = 0.;
+ }
+
+ if (gRandom->Rndm() > eff) {
+ return kFALSE;
+ } else {
+ return kTRUE;
+ }
+
+}
+
+TClonesArray* AliJetKineReader::GetGeneratedJets()
+{
+ //
+ // Get generated jets from mc header
+ //
+ if (!fGenJets) {
+ fGenJets = new TClonesArray("AliAODJets", 100);
} else {
- return kTRUE;
+ fGenJets->Clear();
}
+
+
+ AliHeader* alih = GetAliHeader();
+ if (alih == 0) return 0;
+ AliGenEventHeader * genh = alih->GenEventHeader();
+ if (genh == 0) return 0;
+
+ Int_t nj =((AliGenPythiaEventHeader*)genh)->NTriggerJets();
+ for (Int_t i = 0; i < nj; i++) {
+ Float_t p[4];
+ ((AliGenPythiaEventHeader*)genh)->TriggerJet(i,p);
+ new ((*fGenJets)[i]) AliAODJet(p[0], p[1], p[2], p[3]);
+ }
+
+ 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