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