X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=JETAN%2FAliJetKineReader.cxx;h=4fd19fb1b806e1e4b63286d8f6c3c3f4c08a159a;hb=2e0a9a3b302da2ede89c1ad6ab85ac9253679479;hp=06406a8062b691564952bb1d10ef8371f0773341;hpb=48bb52d83fb00b554792532c4df12572bfbbd49d;p=u%2Fmrichter%2FAliRoot.git diff --git a/JETAN/AliJetKineReader.cxx b/JETAN/AliJetKineReader.cxx index 06406a8062b..4fd19fb1b80 100644 --- a/JETAN/AliJetKineReader.cxx +++ b/JETAN/AliJetKineReader.cxx @@ -13,38 +13,42 @@ * 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 #include #include #include -#include #include #include #include #include // 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 } //____________________________________________________________________________ @@ -52,143 +56,229 @@ AliJetKineReader::AliJetKineReader() 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(); - - fRunLoader->GetEvent(event); - AliStack* stack = fRunLoader->Stack(); + // Get the stack + AliStack* stack = fMCEvent->Stack(); + // Number of primaries Int_t nt = stack->GetNprimary(); - // Get cuts set by user + + // 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; - if ( - (status != 1) - || (pdg == 12 || pdg == 14) - || (pt < ptMin) - ) 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; + Float_t charge = + TDatabasePDG::Instance()->GetParticle(pdg)->Charge(); + if (charge == 0) continue; // Simulate efficiency - if (!Efficiency(p, eta, phi)) continue; + if (!Efficiency(p0, eta, phi)) continue; // Simulate resolution - SmearMomentum(3, p); - } + 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 + // Neutral particles (exclude K0L, n, nbar) + if (pdg == kNeutron || pdg == kK0Long) continue; + } // End fast EMCAL - fMass = part->GetCalcMass(); - fPdgC = part->GetPdgCode(); - // Fill momentum array + // 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(); + 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()Px(),part->Py(),part->Pz(), part->Energy()); + // 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 %d \n", 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; +}