7797e8eb81b1f0e12857ab0dac3f6d5459e076a7
[u/mrichter/AliRoot.git] / JETAN / AliJetKineReader.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
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  **************************************************************************/
15
16 //-------------------------------------------------------------------------
17 // Jet kinematics reader 
18 // MC reader for jet analysis
19 // Author: Andreas Morsch (andreas.morsch@cern.ch)
20 //-------------------------------------------------------------------------
21
22 // From root ...
23 #include <TClonesArray.h>
24 #include <TPDGCode.h>
25 #include <TParticle.h>
26 #include <TParticlePDG.h>
27 //#include <TVector3.h>
28 #include <TLorentzVector.h>
29 #include <TSystem.h>
30 #include <TDatabasePDG.h>
31 #include <TRandom.h>
32 // From AliRoot ...
33 #include "AliJet.h"
34 #include "AliJetKineReaderHeader.h"
35 #include "AliJetKineReader.h"
36 #include "AliRunLoader.h"
37 #include "AliStack.h"
38 #include "AliHeader.h"
39 #include "AliGenPythiaEventHeader.h"
40
41 ClassImp(AliJetKineReader)
42
43 AliJetKineReader::AliJetKineReader():
44   fRunLoader(0x0),
45   fAliHeader(0x0),
46   fMass(0),
47   fPdgC(0)
48 {
49   // Default constructor
50   fReaderHeader = 0x0;
51 }
52
53 //____________________________________________________________________________
54
55 AliJetKineReader::~AliJetKineReader()
56 {
57   // Destructor
58   delete fReaderHeader;
59   delete fAliHeader;
60 }
61
62 //____________________________________________________________________________
63
64 void AliJetKineReader::OpenInputFiles()
65 {
66   // Opens the input file using the run loader
67   const char* dirName = fReaderHeader->GetDirectory();
68   char path[256];
69   sprintf(path, "%s/galice.root",dirName);
70   fRunLoader = AliRunLoader::Open(path);
71   fRunLoader->LoadKinematics();
72   fRunLoader->LoadHeader(); 
73   
74   // get number of events
75   Int_t nMax = fRunLoader->GetNumberOfEvents();
76   printf("\nTotal number of events = %d", nMax);
77   
78   // set number of events in header
79   if (fReaderHeader->GetLastEvent() == -1)
80     fReaderHeader->SetLastEvent(nMax);
81   else {
82     Int_t nUsr = fReaderHeader->GetLastEvent();
83     fReaderHeader->SetLastEvent(TMath::Min(nMax, nUsr));
84   }
85 }
86
87 //____________________________________________________________________________
88
89 Bool_t AliJetKineReader::FillMomentumArray(Int_t event)
90 {
91   // Fill momentum array for event
92  
93   Int_t goodTrack = 0;
94   // Clear array
95   ClearArray();
96   // Get event from runloader
97   fRunLoader->GetEvent(event);
98   // Get the stack
99   AliStack* stack = fRunLoader->Stack();
100   // Number of primaries
101   Int_t nt = stack->GetNprimary();
102
103   // Get cuts set by user and header
104   Double_t ptMin = ((AliJetKineReaderHeader*) fReaderHeader)->GetPtCut();
105   Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
106   Float_t etaMax = fReaderHeader->GetFiducialEtaMax();  
107   fAliHeader = fRunLoader->GetHeader();
108   
109   // temporary storage of signal and cut flags
110   Int_t* sflag  = new Int_t[nt];
111   Int_t* cflag  = new Int_t[nt];
112   
113   TLorentzVector p4;
114   // Loop over particles
115   for (Int_t it = 0; it < nt; it++) {
116     TParticle *part = stack->Particle(it);
117     Int_t   status  = part->GetStatusCode();
118     Int_t   pdg     = TMath::Abs(part->GetPdgCode());
119     Float_t pt      = part->Pt(); 
120     
121     // Skip non-final state particles, neutrinos 
122     // and particles with pt < pt_min 
123     if ((status != 1)            
124         || (pdg == 12 || pdg == 14 || pdg == 16)) continue; 
125     
126     Float_t p       = part->P();
127     Float_t p0      = p;
128     Float_t eta     = part->Eta();
129     Float_t phi     = part->Phi();
130     
131     // Fast simulation of TPC if requested
132     if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimTPC()) {
133       // Charged particles only
134       Float_t charge = 
135         TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
136       if (charge == 0)               continue;
137       // Simulate efficiency
138       if (!Efficiency(p0, eta, phi)) continue;
139       // Simulate resolution
140       p = SmearMomentum(4, p0);
141     } // End fast TPC
142     
143     // Fast simulation of EMCAL if requested
144     if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimEMCAL()) {
145       Float_t charge = 
146         TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
147       // Charged particles only
148       if (charge != 0){
149         // Simulate efficiency
150         if (!Efficiency(p0, eta, phi)) continue;
151         // Simulate resolution
152         p = SmearMomentum(4, p0);
153       } // end "if" charged particles
154       // Neutral particles (exclude K0L, n, nbar)
155       if (pdg == kNeutron || pdg == kK0Long) continue;
156     } // End fast EMCAL
157     
158     fMass = part->GetCalcMass();
159     fPdgC = part->GetPdgCode();
160
161     // Fill momentum array
162     Float_t r  = p/p0;
163     Float_t px = r * part->Px(); 
164     Float_t py = r * part->Py(); 
165     Float_t pz = r * part->Pz();
166     Float_t m  = part->GetMass();
167     Float_t e  = TMath::Sqrt(px * px + py * py + pz * pz + m * m);
168     p4 = TLorentzVector(px, py, pz, e);
169     if ( (p4.Eta()>etaMax) || (p4.Eta()<etaMin)) continue;
170     new ((*fMomentumArray)[goodTrack]) TLorentzVector(p4);
171     
172     // all tracks are signal ... ?
173     sflag[goodTrack]=1;
174     cflag[goodTrack]=0;
175     if (pt > ptMin) cflag[goodTrack]=1; // track surviving pt cut
176     goodTrack++;
177   }
178
179   // set the signal flags
180   fSignalFlag.Set(goodTrack,sflag);
181   fCutFlag.Set(goodTrack,cflag);
182   printf("\nIn event %d, number of good tracks %d \n", event, goodTrack);
183   return kTRUE;
184 }
185
186
187 Float_t AliJetKineReader::SmearMomentum(Int_t ind, Float_t p)
188 {
189   //  The relative momentum error, i.e. 
190   //  (delta p)/p = sqrt (a**2 + (b*p)**2) * 10**-2,
191   //  where typically a = 0.75 and b = 0.16 - 0.24 depending on multiplicity
192   //  (the lower value is for dn/d(eta) about 2000, and the higher one for 8000)
193   //
194   //  If we include information from TRD, b will be a factor 2/3 smaller.
195   //
196   //  ind = 1: high multiplicity
197   //  ind = 2: low  multiplicity
198   //  ind = 3: high multiplicity + TRD
199   //  ind = 4: low  multiplicity + TRD
200   
201   Float_t pSmeared;
202   Float_t a = 0.75;
203   Float_t b = 0.12;
204   
205   if (ind == 1) b = 0.12;
206   if (ind == 2) b = 0.08;
207   if (ind == 3) b = 0.12;    
208   if (ind == 4) b = 0.08;    
209   
210   Float_t sigma =  p * TMath::Sqrt(a * a + b * b * p * p)*0.01;
211   pSmeared = p + gRandom->Gaus(0., sigma);
212   return pSmeared;
213 }
214
215
216 Bool_t AliJetKineReader::Efficiency(Float_t p, Float_t /*eta*/, Float_t phi)
217 {
218   // Fast simulation of geometrical acceptance and tracking efficiency
219  
220   //  Tracking
221   Float_t eff = 0.99;
222   if (p < 0.5) eff -= (0.5-p)*0.2/0.3;
223   // Geometry
224   if (p > 0.5) {
225     phi *= 180. / TMath::Pi();
226     // Sector number 0 - 17
227     Int_t isec  = Int_t(phi / 20.);
228     // Sector centre
229     Float_t phi0 = isec * 20. + 10.;
230     Float_t phir = TMath::Abs(phi-phi0);
231     // 2 deg of dead space
232     if (phir > 9.) eff = 0.;
233   }
234
235   if (gRandom->Rndm() > eff) {
236     return kFALSE;
237   } else {
238     return kTRUE;
239   }
240
241 }
242
243 Bool_t AliJetKineReader::GetGenJets(AliJet* genJets)
244 {
245     // Get generated jets from mc header
246     AliHeader* alih = GetAliHeader(); 
247     if (alih == 0) return kFALSE;
248     AliGenEventHeader * genh = alih->GenEventHeader();
249     if (genh == 0) return kFALSE;
250     Int_t nj =((AliGenPythiaEventHeader*)genh)->NTriggerJets(); 
251     Int_t* m = new Int_t[nj];
252     Int_t* k = new Int_t[nj];
253     for (Int_t i=0; i< nj; i++) {
254         Float_t p[4];
255         ((AliGenPythiaEventHeader*)genh)->TriggerJet(i,p);
256         genJets->AddJet(p[0],p[1],p[2],p[3]);
257         m[i]=1;
258         k[i]=i;
259     }
260     genJets->SetNinput(nj);
261     genJets->SetMultiplicities(m);
262     genJets->SetInJet(k);
263     return kTRUE;
264 }