Adding include path to allow compilation of CleanGeom task
[u/mrichter/AliRoot.git] / JETAN / AliJetFillCalTrkTrackMC.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 /* $Id$ */
17
18 //--------------------------------------------------
19 // Filling of CalTrkTrack objects in the MC reader task
20 //
21 // Author: magali.estienne@subatech.in2p3.fr
22 //         alexandre.shabetai@cern.ch
23 //-------------------------------------------------
24
25 // --- ROOT system ---
26 #include <TDatabasePDG.h>
27 #include <TRandom.h>
28 #include <TPDGCode.h>
29
30 // --- AliRoot system ---
31 #include "AliJetFillCalTrkTrackMC.h"
32 #include "AliAODMCParticle.h"
33 #include "AliJetKineReaderHeader.h"
34 #include "AliVEvent.h"
35 #include "AliMCEvent.h"
36
37 using std::cout;
38 using std::endl;
39 ClassImp(AliJetFillCalTrkTrackMC)
40
41 /////////////////////////////////////////////////////////////////////////
42
43 AliJetFillCalTrkTrackMC::AliJetFillCalTrkTrackMC(): 
44   AliJetFillCalTrkEvent(),
45   fHadCorr(0),
46   fApplyMIPCorrection(kTRUE),
47   fVEvt(0x0),
48   fMCEvent(0x0)
49 {
50   // constructor
51 }
52
53 //-----------------------------------------------------------------------
54 AliJetFillCalTrkTrackMC::AliJetFillCalTrkTrackMC(AliVEvent* evt): 
55   AliJetFillCalTrkEvent(),
56   fHadCorr(0),
57   fApplyMIPCorrection(kTRUE),
58   fVEvt(evt),
59   fMCEvent(0x0)
60 {
61   // constructor
62 }
63
64 //-----------------------------------------------------------------------
65 AliJetFillCalTrkTrackMC::AliJetFillCalTrkTrackMC(const AliJetFillCalTrkTrackMC &det): 
66   AliJetFillCalTrkEvent(det),
67   fHadCorr(det.fHadCorr),
68   fApplyMIPCorrection(det.fApplyMIPCorrection),
69   fVEvt(det.fVEvt),
70   fMCEvent(det.fMCEvent)
71 {
72   // Copy constructor
73 }
74
75 //-----------------------------------------------------------------------
76 AliJetFillCalTrkTrackMC& AliJetFillCalTrkTrackMC::operator=(const AliJetFillCalTrkTrackMC& other)
77 {
78   // Assignment
79   if (this != &other) { 
80    fHadCorr = other.fHadCorr;
81    fApplyMIPCorrection = other.fApplyMIPCorrection;
82    fVEvt = other.fVEvt;
83    fMCEvent = other.fMCEvent;
84   }
85  
86   return (*this);
87
88 }
89
90 //-----------------------------------------------------------------------
91 AliJetFillCalTrkTrackMC::~AliJetFillCalTrkTrackMC()
92 {
93   // destructor
94 }
95
96 //-----------------------------------------------------------------------
97 void AliJetFillCalTrkTrackMC::Exec(Option_t const * /*option*/)
98 {
99   // Main method.
100   // Fill AliJetFillCalTrkTrackMC the with the charged particle information
101
102   fDebug = fReaderHeader->GetDebug();
103   fOpt = fReaderHeader->GetDetector();
104   TString datatype = fReaderHeader->GetDataType();
105   datatype.ToUpper();
106
107   Int_t type=0;
108
109   if      (datatype.Contains("AODMC2B")) { type=kTrackAODMCChargedAcceptance; }
110   else if (datatype.Contains("AODMC2"))  { type=kTrackAODMCCharged; }
111   else if (datatype.Contains("AODMC"))   { type=kTrackAODMCAll; }
112  
113   else if (!datatype.Contains("AOD") && datatype.Contains("MC2")) {type=kTrackKineCharged;} 
114   else if (!datatype.Contains("AOD") && datatype.Contains("MC"))  {type=kTrackKineAll;}
115   
116   else { cout<< "Unknown Data type !" << endl; }
117
118   // temporary storage of signal and pt cut flag
119   Bool_t sflag = 0;
120   Bool_t cflag = 0;
121
122   // get cuts set by user
123   Float_t ptMin  = fReaderHeader->GetPtCut();
124   Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
125   Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
126   Float_t phiMin = fReaderHeader->GetFiducialPhiMin();
127   Float_t phiMax = fReaderHeader->GetFiducialPhiMax();
128
129   if(fDebug>2) Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
130
131   Int_t goodTrackStd = 0;
132   Int_t goodTrackNonStd = 0;
133
134   if (type == kTrackKineAll || type == kTrackKineCharged){
135
136     FillKine();
137   } // Kine
138
139   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
140     if(!fVEvt) return; 
141     TClonesArray *mcArray = dynamic_cast<TClonesArray*>(fVEvt->FindListObject(AliAODMCParticle::StdBranchName()));
142     if(!mcArray) {
143       Printf("%s:%d No MC particle branch found",(char*)__FILE__,__LINE__);
144       return; 
145     }
146     for(int it = 0;it < mcArray->GetEntriesFast();++it){
147       cflag=sflag=0;
148       AliAODMCParticle *part = (AliAODMCParticle*)(mcArray->At(it));
149       if(!part->IsPhysicalPrimary())continue;
150       Int_t   pdg     = TMath::Abs(part->GetPdgCode());
151       // exclude neutrinos anyway
152       if((pdg == 12 || pdg == 14 || pdg == 16)) continue;
153       cflag = ( part->Pt()>ptMin ) ? 1 : 0;
154       sflag = ( TMath::Abs(part->GetLabel()) < 10000 ) ? 1 : 0;
155       if(type == kTrackAODMCAll){
156         fCalTrkEvent->AddCalTrkTrack(part,cflag,sflag);
157         goodTrackStd++;
158       }
159       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
160         if(part->Charge()==0)continue;
161         if(kTrackAODMCCharged){
162           fCalTrkEvent->AddCalTrkTrack(part,cflag,sflag);
163           goodTrackStd++;
164         }
165         else {
166           if((part->Eta()>etaMax) || (part->Eta()<etaMin)) continue;
167           Float_t phi = ( (part->Phi()) < 0) ? (part->Phi()) + 2. * TMath::Pi() : (part->Phi());
168           if((phi>phiMax) || (phi<phiMin)) continue;
169           fCalTrkEvent->AddCalTrkTrack(part,cflag,sflag);
170           goodTrackStd++;
171         }
172       }
173     }
174    if(fDebug>0) printf("Number of good tracks selected: %5d \n", goodTrackStd+goodTrackNonStd); 
175   } // AODMCparticle
176
177 }
178
179 //-----------------------------------------------------------------------
180 Bool_t AliJetFillCalTrkTrackMC::FillKine()
181 {
182   // Fill event
183   Int_t goodTrack = 0;
184   
185   // Get the stack
186   if(!fMCEvent) {cout<<"could not get MCEvent!"<<endl; return kFALSE;}
187   
188   Int_t nt = fMCEvent->GetNumberOfTracks();
189   
190   // Get cuts set by user and header
191   Double_t ptMin = ((AliJetKineReaderHeader*) fReaderHeader)->GetPtCut();
192   Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
193   Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
194   
195   TLorentzVector p4;
196   // Loop over particles
197   for (Int_t it = 0; it < nt; it++) {
198     if(!fMCEvent->IsPhysicalPrimary(it)) continue;
199     AliVParticle* part = fMCEvent->GetTrack(it);
200     Int_t   pdg    = TMath::Abs(part->PdgCode());
201     Float_t pt     = part->Pt();
202     
203     if( (pdg == 12 || pdg == 14 || pdg == 16)) continue;
204     
205     Float_t p      = part->P();
206     Float_t p0     = p;
207     Float_t eta    = part->Eta();
208     Float_t phi    = part->Phi();
209     Float_t charge = part->Charge();
210     
211     if (((AliJetKineReaderHeader*)fReaderHeader)->ChargedOnly()) {
212       if (charge == 0) continue;
213     } // End charged only
214     
215     // Fast simulation of EMCAL if requested
216     if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimEMCAL()) {
217       // Charged particles only
218       if (charge != 0){
219         // Simulate efficiency
220         if (!Efficiency(p0, eta, phi)) continue;
221         // Simulate resolution
222         p = SmearMomentum(4, p0);
223       } // end "if" charged particles
224       // Neutral particles (exclude K0L, n, nbar)
225       if (pdg == kNeutron || pdg == kK0Long) continue;
226     } // End fast EMCAL
227     
228     // Fast simulation of TPC if requested
229     if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimTPC()) {
230       // Charged particles only
231       if (charge == 0)               continue;
232       // Simulate efficiency
233       if (!Efficiency(p0, eta, phi)) continue;
234       // Simulate resolution
235       p = SmearMomentum(4, p0);
236     } // End fast TPC
237
238     // Fill momentum array
239     Float_t r  = p/p0;
240     Float_t px = r * part->Px();
241     Float_t py = r * part->Py();
242     Float_t pz = r * part->Pz();
243     Float_t m =  part->M();
244     
245     Float_t e  = TMath::Sqrt(px * px + py * py + pz * pz + m * m);
246     p4 = TLorentzVector(px, py, pz, e);
247     if ((p4.Eta()>etaMax) || (p4.Eta()<etaMin)) continue;
248
249     // Flag used to store the r factor
250     Float_t ptReso = r;
251     Bool_t cflag = kFALSE, sflag = kTRUE;
252     if (pt > ptMin) cflag = kTRUE; // track surviving pt cut
253
254     fCalTrkEvent->AddCalTrkTrackKine(part,cflag,sflag,ptReso);
255
256     goodTrack++;
257
258   } // track loop
259
260   if(fDebug>0) printf(" Number of good tracks %d \n", goodTrack);
261
262   return kTRUE;
263
264 }
265
266 //-----------------------------------------------------------------------
267 Float_t AliJetFillCalTrkTrackMC::SmearMomentum(Int_t ind, Float_t p)
268 {
269   //  The relative momentum error, i.e.
270   //  (delta p)/p = sqrt (a**2 + (b*p)**2) * 10**-2,
271   //  where typically a = 0.75 and b = 0.16 - 0.24 depending on multiplicity
272   //  (the lower value is for dn/d(eta) about 2000, and the higher one for 8000)
273   //
274   //  If we include information from TRD, b will be a factor 2/3 smaller.
275   //
276   //  ind = 1: high multiplicity
277   //  ind = 2: low  multiplicity
278   //  ind = 3: high multiplicity + TRD
279   //  ind = 4: low  multiplicity + TRD
280
281   Float_t pSmeared;
282   Float_t a = 0.75;
283   Float_t b = 0.12;
284
285   if (ind == 1) b = 0.12;
286   if (ind == 2) b = 0.08;
287   if (ind == 3) b = 0.12;
288   if (ind == 4) b = 0.08;
289
290   Float_t sigma =  p * TMath::Sqrt(a * a + b * b * p * p)*0.01;
291   pSmeared = p + gRandom->Gaus(0., sigma);
292   return pSmeared;
293
294 }
295
296 //-----------------------------------------------------------------------
297 Bool_t AliJetFillCalTrkTrackMC::Efficiency(Float_t p, Float_t /*eta*/, Float_t phi)
298 {
299   // Fast simulation of geometrical acceptance and tracking efficiency
300   
301   //  Tracking
302   Float_t eff = 0.99;
303   if (p < 0.5) eff -= (0.5-p)*0.2/0.3;
304   // Geometry
305   if (p > 0.5) {
306     phi *= 180. / TMath::Pi();
307     // Sector number 0 - 17
308     Int_t isec  = Int_t(phi / 20.);
309     // Sector centre
310     Float_t phi0 = isec * 20. + 10.;
311     Float_t phir = TMath::Abs(phi-phi0);
312     // 2 deg of dead space
313     if (phir > 9.) eff = 0.;
314   }
315   
316   if (gRandom->Rndm() > eff) {
317     return kFALSE;
318   } else {
319     return kTRUE;
320   }
321   
322 }
323
324