Compatibility with ROOT trunk
[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 <TLorentzVector.h>
28 #include <TSystem.h>
29 #include <TDatabasePDG.h>
30 #include <TRandom.h>
31 // From AliRoot ...
32 #include "AliAODJet.h"
33 #include "AliPDG.h"
34 #include "AliJetKineReaderHeader.h"
35 #include "AliJetKineReader.h"
36 #include "AliMCEventHandler.h"
37 #include "AliMCEvent.h"
38 #include "AliStack.h"
39 #include "AliHeader.h"
40 #include "AliGenPythiaEventHeader.h"
41 #include "AliLog.h"
42
43 ClassImp(AliJetKineReader)
44
45 AliJetKineReader::AliJetKineReader():
46     AliJetReader(),  
47     fAliHeader(0),
48     fMCEvent(0),
49     fGenJets(0)
50 {
51   // Default constructor
52 }
53
54 //____________________________________________________________________________
55
56 AliJetKineReader::~AliJetKineReader()
57 {
58   // Destructor
59     if (fAliHeader) {
60         delete fAliHeader;
61         fAliHeader = 0;
62     }
63 }
64
65 //____________________________________________________________________________
66
67 void AliJetKineReader::OpenInputFiles()
68 {
69   // Opens the input file using the run loader
70   // Not used anymore   
71 }
72
73 //____________________________________________________________________________
74
75 Bool_t AliJetKineReader::FillMomentumArray()
76 {
77   // Fill momentum array for event
78     Int_t goodTrack = 0;
79     // Clear array
80     // temporary storage of signal and cut flags
81     Int_t* sflag  = new Int_t[50000];
82     Int_t* cflag  = new Int_t[50000];
83   
84     ClearArray();
85     // Get the stack
86     AliStack* stack = fMCEvent->Stack();
87     // Number of primaries
88     Int_t nt = stack->GetNprimary();
89       
90     // Get cuts set by user and header
91     Double_t ptMin = ((AliJetKineReaderHeader*) fReaderHeader)->GetPtCut();
92     Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
93     Float_t etaMax = fReaderHeader->GetFiducialEtaMax();  
94     fAliHeader = fMCEvent->Header();
95       
96     
97     TLorentzVector p4;
98     // Loop over particles
99     for (Int_t it = 0; it < nt; it++) {
100         TParticle *part = stack->Particle(it);
101         Int_t   status  = part->GetStatusCode();
102         Int_t   pdg     = TMath::Abs(part->GetPdgCode());
103         Float_t pt      = part->Pt(); 
104           
105         // Skip non-final state particles, neutrinos 
106         // and particles with pt < pt_min 
107         if ((status != 1)            
108             || (pdg == 12 || pdg == 14 || pdg == 16)) continue; 
109         
110         Float_t p       = part->P();
111         Float_t p0      = p;
112         Float_t eta     = part->Eta();
113         Float_t phi     = part->Phi();
114
115
116         if (((AliJetKineReaderHeader*)fReaderHeader)->ChargedOnly()) {
117             // Charged particles only
118             Float_t charge = 
119                 TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
120             if (charge == 0)               continue;
121         } // End charged only
122
123         
124         // Fast simulation of TPC if requested
125         if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimTPC()) {
126             // Charged particles only
127             Float_t charge = 
128                 TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
129             if (charge == 0)               continue;
130             // Simulate efficiency
131             if (!Efficiency(p0, eta, phi)) continue;
132             // Simulate resolution
133             p = SmearMomentum(4, p0);
134         } // End fast TPC
135           
136           // Fast simulation of EMCAL if requested
137         if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimEMCAL()) {
138             Float_t charge = 
139                 TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
140             // Charged particles only
141             if (charge != 0){
142                 // Simulate efficiency
143                 if (!Efficiency(p0, eta, phi)) continue;
144                 // Simulate resolution
145                 p = SmearMomentum(4, p0);
146             } // end "if" charged particles
147               // Neutral particles (exclude K0L, n, nbar)
148             if (pdg == kNeutron || pdg == kK0Long) continue;
149         } // End fast EMCAL
150         
151           // Fill momentum array
152         Float_t r  = p/p0;
153         Float_t px = r * part->Px(); 
154         Float_t py = r * part->Py(); 
155         Float_t pz = r * part->Pz();
156         TParticlePDG *pPDG = part->GetPDG();
157         Float_t m = 0;
158         if(!pPDG){
159           // this is very rare...
160           // Is avoided by AliPDG::AddParticlesToPdgDataBase();
161           // but this should be called only once... (how in proof?)
162           // Calucluate mass with unsmeared momentum values
163           m  = part->Energy()*part->Energy() - 
164             (px * px + py * py + pz * pz)/r/r; 
165           if(m>0)m = TMath::Sqrt(m);
166           else m = 0;
167           AliInfo(Form("Unknown Particle using %d calculated mass m = %3.3f",part->GetPdgCode(),m));
168
169         }
170         else{
171           m  = pPDG->Mass();
172         }
173         Float_t e  = TMath::Sqrt(px * px + py * py + pz * pz + m * m);
174         p4 = TLorentzVector(px, py, pz, e);
175         if ( (p4.Eta()>etaMax) || (p4.Eta()<etaMin)) continue;
176         new ((*fMomentumArray)[goodTrack]) TLorentzVector(p4);
177         
178         // all tracks are signal ... ?
179         sflag[goodTrack]=1;
180         cflag[goodTrack]=0;
181         if (pt > ptMin) cflag[goodTrack] = 1; // track surviving pt cut
182         goodTrack++;
183     } // track loop
184
185 // set the signal flags
186     fSignalFlag.Set(goodTrack, sflag);
187     fCutFlag.Set(goodTrack, cflag);
188     AliInfo(Form(" Number of good tracks %d \n", goodTrack));
189     
190     delete[] sflag;
191     delete[] cflag;
192     
193     return kTRUE;
194 }
195
196
197 Float_t AliJetKineReader::SmearMomentum(Int_t ind, Float_t p)
198 {
199   //  The relative momentum error, i.e. 
200   //  (delta p)/p = sqrt (a**2 + (b*p)**2) * 10**-2,
201   //  where typically a = 0.75 and b = 0.16 - 0.24 depending on multiplicity
202   //  (the lower value is for dn/d(eta) about 2000, and the higher one for 8000)
203   //
204   //  If we include information from TRD, b will be a factor 2/3 smaller.
205   //
206   //  ind = 1: high multiplicity
207   //  ind = 2: low  multiplicity
208   //  ind = 3: high multiplicity + TRD
209   //  ind = 4: low  multiplicity + TRD
210   
211   Float_t pSmeared;
212   Float_t a = 0.75;
213   Float_t b = 0.12;
214   
215   if (ind == 1) b = 0.12;
216   if (ind == 2) b = 0.08;
217   if (ind == 3) b = 0.12;    
218   if (ind == 4) b = 0.08;    
219   
220   Float_t sigma =  p * TMath::Sqrt(a * a + b * b * p * p)*0.01;
221   pSmeared = p + gRandom->Gaus(0., sigma);
222   return pSmeared;
223 }
224
225
226 Bool_t AliJetKineReader::Efficiency(Float_t p, Float_t /*eta*/, Float_t phi)
227 {
228   // Fast simulation of geometrical acceptance and tracking efficiency
229  
230   //  Tracking
231   Float_t eff = 0.99;
232   if (p < 0.5) eff -= (0.5-p)*0.2/0.3;
233   // Geometry
234   if (p > 0.5) {
235     phi *= 180. / TMath::Pi();
236     // Sector number 0 - 17
237     Int_t isec  = Int_t(phi / 20.);
238     // Sector centre
239     Float_t phi0 = isec * 20. + 10.;
240     Float_t phir = TMath::Abs(phi-phi0);
241     // 2 deg of dead space
242     if (phir > 9.) eff = 0.;
243   }
244
245   if (gRandom->Rndm() > eff) {
246     return kFALSE;
247   } else {
248     return kTRUE;
249   }
250
251 }
252
253 TClonesArray*  AliJetKineReader::GetGeneratedJets()
254 {
255     //
256     // Get generated jets from mc header
257     //
258     if (!fGenJets) {
259         fGenJets = new TClonesArray("AliAODJets", 100);
260     } else {
261         fGenJets->Clear();
262     }
263     
264     
265     AliHeader* alih = GetAliHeader(); 
266     if (alih == 0) return 0;
267     AliGenEventHeader * genh = alih->GenEventHeader();
268     if (genh == 0) return 0;
269
270     Int_t nj =((AliGenPythiaEventHeader*)genh)->NTriggerJets(); 
271     for (Int_t i = 0; i < nj; i++) {
272         Float_t p[4];
273         ((AliGenPythiaEventHeader*)genh)->TriggerJet(i,p);
274         new ((*fGenJets)[i]) AliAODJet(p[0], p[1], p[2], p[3]);
275     }
276
277     return fGenJets;
278 }
279
280 void AliJetKineReader::SetInputEvent(const TObject* /*esd*/, const TObject* /*aod*/, const TObject* mc)
281 {
282     // Connect the input event
283     fMCEvent = (AliMCEvent*) mc;
284 }