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