1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 //--------------------------------------------------
19 // Filling of CalTrkTrack objects in the MC reader task
21 // Author: magali.estienne@subatech.in2p3.fr
22 // alexandre.shabetai@cern.ch
23 //-------------------------------------------------
25 // --- ROOT system ---
26 #include <TDatabasePDG.h>
30 // --- AliRoot system ---
31 #include "AliJetFillCalTrkTrackMC.h"
32 #include "AliAODMCParticle.h"
33 #include "AliJetKineReaderHeader.h"
34 #include "AliVEvent.h"
35 #include "AliMCEvent.h"
39 ClassImp(AliJetFillCalTrkTrackMC)
41 /////////////////////////////////////////////////////////////////////////
43 AliJetFillCalTrkTrackMC::AliJetFillCalTrkTrackMC():
44 AliJetFillCalTrkEvent(),
46 fApplyMIPCorrection(kTRUE),
53 //-----------------------------------------------------------------------
54 AliJetFillCalTrkTrackMC::AliJetFillCalTrkTrackMC(AliVEvent* evt):
55 AliJetFillCalTrkEvent(),
57 fApplyMIPCorrection(kTRUE),
64 //-----------------------------------------------------------------------
65 AliJetFillCalTrkTrackMC::AliJetFillCalTrkTrackMC(const AliJetFillCalTrkTrackMC &det):
66 AliJetFillCalTrkEvent(det),
67 fHadCorr(det.fHadCorr),
68 fApplyMIPCorrection(det.fApplyMIPCorrection),
70 fMCEvent(det.fMCEvent)
75 //-----------------------------------------------------------------------
76 AliJetFillCalTrkTrackMC& AliJetFillCalTrkTrackMC::operator=(const AliJetFillCalTrkTrackMC& other)
80 fHadCorr = other.fHadCorr;
81 fApplyMIPCorrection = other.fApplyMIPCorrection;
83 fMCEvent = other.fMCEvent;
90 //-----------------------------------------------------------------------
91 AliJetFillCalTrkTrackMC::~AliJetFillCalTrkTrackMC()
96 //-----------------------------------------------------------------------
97 void AliJetFillCalTrkTrackMC::Exec(Option_t const * /*option*/)
100 // Fill AliJetFillCalTrkTrackMC the with the charged particle information
102 fDebug = fReaderHeader->GetDebug();
103 fOpt = fReaderHeader->GetDetector();
104 TString datatype = fReaderHeader->GetDataType();
109 if (datatype.Contains("AODMC2B")) { type=kTrackAODMCChargedAcceptance; }
110 else if (datatype.Contains("AODMC2")) { type=kTrackAODMCCharged; }
111 else if (datatype.Contains("AODMC")) { type=kTrackAODMCAll; }
113 else if (!datatype.Contains("AOD") && datatype.Contains("MC2")) {type=kTrackKineCharged;}
114 else if (!datatype.Contains("AOD") && datatype.Contains("MC")) {type=kTrackKineAll;}
116 else { cout<< "Unknown Data type !" << endl; }
118 // temporary storage of signal and pt cut flag
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();
129 if(fDebug>2) Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
131 Int_t goodTrackStd = 0;
132 Int_t goodTrackNonStd = 0;
134 if (type == kTrackKineAll || type == kTrackKineCharged){
139 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
141 TClonesArray *mcArray = dynamic_cast<TClonesArray*>(fVEvt->FindListObject(AliAODMCParticle::StdBranchName()));
143 Printf("%s:%d No MC particle branch found",(char*)__FILE__,__LINE__);
146 for(int it = 0;it < mcArray->GetEntriesFast();++it){
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);
159 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
160 if(part->Charge()==0)continue;
161 if(kTrackAODMCCharged){
162 fCalTrkEvent->AddCalTrkTrack(part,cflag,sflag);
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);
174 if(fDebug>0) printf("Number of good tracks selected: %5d \n", goodTrackStd+goodTrackNonStd);
179 //-----------------------------------------------------------------------
180 Bool_t AliJetFillCalTrkTrackMC::FillKine()
186 if(!fMCEvent) {cout<<"could not get MCEvent!"<<endl; return kFALSE;}
188 Int_t nt = fMCEvent->GetNumberOfTracks();
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();
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();
203 if( (pdg == 12 || pdg == 14 || pdg == 16)) continue;
205 Float_t p = part->P();
207 Float_t eta = part->Eta();
208 Float_t phi = part->Phi();
209 Float_t charge = part->Charge();
211 if (((AliJetKineReaderHeader*)fReaderHeader)->ChargedOnly()) {
212 if (charge == 0) continue;
213 } // End charged only
215 // Fast simulation of EMCAL if requested
216 if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimEMCAL()) {
217 // Charged particles only
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;
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);
238 // Fill momentum array
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();
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;
249 // Flag used to store the r factor
251 Bool_t cflag = kFALSE, sflag = kTRUE;
252 if (pt > ptMin) cflag = kTRUE; // track surviving pt cut
254 fCalTrkEvent->AddCalTrkTrackKine(part,cflag,sflag,ptReso);
260 if(fDebug>0) printf(" Number of good tracks %d \n", goodTrack);
266 //-----------------------------------------------------------------------
267 Float_t AliJetFillCalTrkTrackMC::SmearMomentum(Int_t ind, Float_t p)
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)
274 // If we include information from TRD, b will be a factor 2/3 smaller.
276 // ind = 1: high multiplicity
277 // ind = 2: low multiplicity
278 // ind = 3: high multiplicity + TRD
279 // ind = 4: low multiplicity + TRD
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;
290 Float_t sigma = p * TMath::Sqrt(a * a + b * b * p * p)*0.01;
291 pSmeared = p + gRandom->Gaus(0., sigma);
296 //-----------------------------------------------------------------------
297 Bool_t AliJetFillCalTrkTrackMC::Efficiency(Float_t p, Float_t /*eta*/, Float_t phi)
299 // Fast simulation of geometrical acceptance and tracking efficiency
303 if (p < 0.5) eff -= (0.5-p)*0.2/0.3;
306 phi *= 180. / TMath::Pi();
307 // Sector number 0 - 17
308 Int_t isec = Int_t(phi / 20.);
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.;
316 if (gRandom->Rndm() > eff) {