1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //-------------------------------------------------------------------------
17 // Jet kinematics reader
18 // MC reader for jet analysis
19 // Author: Andreas Morsch (andreas.morsch@cern.ch)
20 //-------------------------------------------------------------------------
23 #include <TClonesArray.h>
25 #include <TParticle.h>
26 #include <TParticlePDG.h>
27 #include <TLorentzVector.h>
29 #include <TDatabasePDG.h>
32 #include "AliAODJet.h"
34 #include "AliJetKineReaderHeader.h"
35 #include "AliJetKineReader.h"
36 #include "AliMCEventHandler.h"
37 #include "AliMCEvent.h"
39 #include "AliHeader.h"
40 #include "AliGenPythiaEventHeader.h"
43 ClassImp(AliJetKineReader)
45 AliJetKineReader::AliJetKineReader():
51 // Default constructor
54 //____________________________________________________________________________
56 AliJetKineReader::~AliJetKineReader()
62 //____________________________________________________________________________
64 void AliJetKineReader::OpenInputFiles()
66 // Opens the input file using the run loader
70 //____________________________________________________________________________
72 Bool_t AliJetKineReader::FillMomentumArray()
74 // Fill momentum array for event
77 // temporary storage of signal and cut flags
78 Int_t* sflag = new Int_t[50000];
79 Int_t* cflag = new Int_t[50000];
83 AliStack* stack = fMCEvent->Stack();
84 // Number of primaries
85 Int_t nt = stack->GetNprimary();
87 // Get cuts set by user and header
88 Double_t ptMin = ((AliJetKineReaderHeader*) fReaderHeader)->GetPtCut();
89 Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
90 Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
91 fAliHeader = fMCEvent->Header();
95 // Loop over particles
96 for (Int_t it = 0; it < nt; it++) {
97 TParticle *part = stack->Particle(it);
98 Int_t status = part->GetStatusCode();
99 Int_t pdg = TMath::Abs(part->GetPdgCode());
100 Float_t pt = part->Pt();
102 // Skip non-final state particles, neutrinos
103 // and particles with pt < pt_min
105 || (pdg == 12 || pdg == 14 || pdg == 16)) continue;
107 Float_t p = part->P();
109 Float_t eta = part->Eta();
110 Float_t phi = part->Phi();
112 // Fast simulation of TPC if requested
113 if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimTPC()) {
114 // Charged particles only
116 TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
117 if (charge == 0) continue;
118 // Simulate efficiency
119 if (!Efficiency(p0, eta, phi)) continue;
120 // Simulate resolution
121 p = SmearMomentum(4, p0);
124 // Fast simulation of EMCAL if requested
125 if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimEMCAL()) {
127 TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
128 // Charged particles only
130 // Simulate efficiency
131 if (!Efficiency(p0, eta, phi)) continue;
132 // Simulate resolution
133 p = SmearMomentum(4, p0);
134 } // end "if" charged particles
135 // Neutral particles (exclude K0L, n, nbar)
136 if (pdg == kNeutron || pdg == kK0Long) continue;
139 // Fill momentum array
141 Float_t px = r * part->Px();
142 Float_t py = r * part->Py();
143 Float_t pz = r * part->Pz();
144 TParticlePDG *pPDG = part->GetPDG();
147 // this is very rare...
148 // Is avoided by AliPDG::AddParticlesToPdgDataBase();
149 // but this should be called only once... (how in proof?)
150 // Calucluate mass with unsmeared momentum values
151 m = part->Energy()*part->Energy() -
152 (px * px + py * py + pz * pz)/r/r;
153 if(m>0)m = TMath::Sqrt(m);
155 AliInfo(Form("Unknown Particle using %d calculated mass m = %3.3f",part->GetPdgCode(),m));
161 Float_t e = TMath::Sqrt(px * px + py * py + pz * pz + m * m);
162 p4 = TLorentzVector(px, py, pz, e);
163 if ( (p4.Eta()>etaMax) || (p4.Eta()<etaMin)) continue;
164 new ((*fMomentumArray)[goodTrack]) TLorentzVector(p4);
166 // all tracks are signal ... ?
169 if (pt > ptMin) cflag[goodTrack] = 1; // track surviving pt cut
173 // set the signal flags
174 fSignalFlag.Set(goodTrack, sflag);
175 fCutFlag.Set(goodTrack, cflag);
176 AliInfo(Form(" Number of good tracks %d \n", goodTrack));
185 Float_t AliJetKineReader::SmearMomentum(Int_t ind, Float_t p)
187 // The relative momentum error, i.e.
188 // (delta p)/p = sqrt (a**2 + (b*p)**2) * 10**-2,
189 // where typically a = 0.75 and b = 0.16 - 0.24 depending on multiplicity
190 // (the lower value is for dn/d(eta) about 2000, and the higher one for 8000)
192 // If we include information from TRD, b will be a factor 2/3 smaller.
194 // ind = 1: high multiplicity
195 // ind = 2: low multiplicity
196 // ind = 3: high multiplicity + TRD
197 // ind = 4: low multiplicity + TRD
203 if (ind == 1) b = 0.12;
204 if (ind == 2) b = 0.08;
205 if (ind == 3) b = 0.12;
206 if (ind == 4) b = 0.08;
208 Float_t sigma = p * TMath::Sqrt(a * a + b * b * p * p)*0.01;
209 pSmeared = p + gRandom->Gaus(0., sigma);
214 Bool_t AliJetKineReader::Efficiency(Float_t p, Float_t /*eta*/, Float_t phi)
216 // Fast simulation of geometrical acceptance and tracking efficiency
220 if (p < 0.5) eff -= (0.5-p)*0.2/0.3;
223 phi *= 180. / TMath::Pi();
224 // Sector number 0 - 17
225 Int_t isec = Int_t(phi / 20.);
227 Float_t phi0 = isec * 20. + 10.;
228 Float_t phir = TMath::Abs(phi-phi0);
229 // 2 deg of dead space
230 if (phir > 9.) eff = 0.;
233 if (gRandom->Rndm() > eff) {
241 TClonesArray* AliJetKineReader::GetGeneratedJets()
244 // Get generated jets from mc header
247 fGenJets = new TClonesArray("AliAODJets", 100);
253 AliHeader* alih = GetAliHeader();
254 if (alih == 0) return 0;
255 AliGenEventHeader * genh = alih->GenEventHeader();
256 if (genh == 0) return 0;
258 Int_t nj =((AliGenPythiaEventHeader*)genh)->NTriggerJets();
259 for (Int_t i = 0; i < nj; i++) {
261 ((AliGenPythiaEventHeader*)genh)->TriggerJet(i,p);
262 new ((*fGenJets)[i]) AliAODJet(p[0], p[1], p[2], p[3]);
268 void AliJetKineReader::SetInputEvent(TObject* /*esd*/, TObject* /*aod*/, TObject* mc)
270 // Connect the input event
271 fMCEvent = (AliMCEvent*) mc;