AliPhysics  fe039ad (fe039ad)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliJetFastSimulation.cxx
Go to the documentation of this file.
1 // $Id$
2 //
3 // Jet model task to merge to existing branches
4 // only implemented for track branches
5 //
6 // Author: M. Verweij
7 
8 #include "AliJetFastSimulation.h"
9 
10 #include <TClonesArray.h>
11 #include <TFolder.h>
12 #include <TLorentzVector.h>
13 #include <TParticle.h>
14 #include <TParticlePDG.h>
15 #include <TRandom3.h>
16 #include <TProfile.h>
17 #include <TGrid.h>
18 #include <TFile.h>
19 #include <TF1.h>
20 #include "AliAnalysisManager.h"
21 #include "AliEMCALDigit.h"
22 #include "AliEMCALGeometry.h"
23 #include "AliEMCALRecPoint.h"
24 #include "AliGenerator.h"
25 #include "AliHeader.h"
26 #include "AliLog.h"
27 #include "AliPicoTrack.h"
28 #include "AliRun.h"
29 #include "AliRunLoader.h"
30 #include "AliStack.h"
31 #include "AliStack.h"
32 #include "AliVCluster.h"
33 #include "AliVEvent.h"
34 
35 ClassImp(AliJetFastSimulation)
36 
37 //________________________________________________________________________
40  fTracksOutName(""),
41  fTracksOut(0x0),
42  fNTrackClasses(2),
43  fRandom(0),
44  fEfficiencyFixed(1.1),
45  fMomResH1(0x0),
46  fMomResH2(0x0),
47  fMomResH3(0x0),
48  fMomResH1Fit(0x0),
49  fMomResH2Fit(0x0),
50  fMomResH3Fit(0x0),
51  fhEffH1(0x0),
52  fhEffH2(0x0),
53  fhEffH3(0x0),
54  fUseTrPtResolutionSmearing(kFALSE),
55  fUseDiceEfficiency(kFALSE),
56  fDiceEfficiencyMinPt(-1.),
57  fUncertEfficiency(1),
58  fUseTrPtResolutionFromOADB(kFALSE),
59  fUseTrEfficiencyFromOADB(kFALSE),
60  fPathTrPtResolution(""),
61  fPathTrEfficiency(""),
62  fHistPtDet(0),
63  fh2PtGenPtSmeared(0),
64  fp1Efficiency(0),
65  fp1PtResolution(0)
66 {
67  // Default constructor.
68  SetMakeGeneralHistograms(kTRUE);
69 }
70 
71 //________________________________________________________________________
73  AliAnalysisTaskEmcal(name,kTRUE),
74  fTracksOutName(""),
75  fTracksOut(0x0),
76  fNTrackClasses(2),
77  fRandom(0),
78  fEfficiencyFixed(1.1),
79  fMomResH1(0x0),
80  fMomResH2(0x0),
81  fMomResH3(0x0),
82  fMomResH1Fit(0x0),
83  fMomResH2Fit(0x0),
84  fMomResH3Fit(0x0),
85  fhEffH1(0x0),
86  fhEffH2(0x0),
87  fhEffH3(0x0),
88  fUseTrPtResolutionSmearing(kFALSE),
89  fUseDiceEfficiency(kFALSE),
90  fDiceEfficiencyMinPt(-1.),
91  fUncertEfficiency(1),
92  fUseTrPtResolutionFromOADB(kFALSE),
93  fUseTrEfficiencyFromOADB(kFALSE),
94  fPathTrPtResolution(""),
95  fPathTrEfficiency(""),
96  fHistPtDet(0),
97  fh2PtGenPtSmeared(0),
98  fp1Efficiency(0),
99  fp1PtResolution(0)
100 {
101  // Standard constructor.
103 }
104 
105 //________________________________________________________________________
107 {
108  // Destructor
109 
110  delete fRandom;
111 }
112 
113 //________________________________________________________________________
115 {
116  // Exec only once.
117 
119 
120  if(!fRandom) fRandom = new TRandom3(0);
121 
122  if (!fTracksOutName.IsNull()) {
123  fTracksOut = new TClonesArray("AliPicoTrack");
124  fTracksOut->SetName(fTracksOutName);
125  if (InputEvent()->FindListObject(fTracksOutName)) {
126  AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fTracksOutName.Data()));
127  return;
128  }
129  else {
130  InputEvent()->AddObject(fTracksOut);
131  }
132  }
133 }
134 //________________________________________________________________________
136  //initialize track response
140 }
141 
142 //________________________________________________________________________
144 {
146 
147  const Int_t nBinPt = 100;
148  Double_t binLimitsPt[nBinPt+1];
149  for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
150  if(iPt == 0){
151  binLimitsPt[iPt] = 0.0;
152  } else {// 1.0
153  binLimitsPt[iPt] = binLimitsPt[iPt-1] + 1.0;
154  }
155  }
156 
157  fHistPtDet = new TH1F("fHistpt","fHistPtDet;#it{p}_{T};N",nBinPt,binLimitsPt);
158  fOutput->Add(fHistPtDet);
159 
160  fh2PtGenPtSmeared = new TH2F("fh2PtGenPtSmeared","fh2PtGenPtSmeared",nBinPt,binLimitsPt,nBinPt,binLimitsPt);
162 
163  fp1Efficiency = new TProfile("fp1Efficiency","fp1Efficiency",nBinPt,binLimitsPt);
164  fOutput->Add(fp1Efficiency);
165 
166  fp1PtResolution = new TProfile("fp1PtResolution","fp1PtResolution",nBinPt,binLimitsPt);
167  fOutput->Add(fp1PtResolution);
168 
169  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
170 }
171 
172 
173 //________________________________________________________________________
175 {
176  //Check if information is provided detector level effects
177  if(!fMomResH1 && fNTrackClasses>0 )
179  if(!fMomResH2 && fNTrackClasses>1 )
181  if(!fMomResH3 && fNTrackClasses>2 )
183 
184  if(fEfficiencyFixed < 1.)
185  fUseDiceEfficiency = 1; // 1 is the default; 2 can be set by user but not implemented
186  else {
187  if(!fhEffH1 && fNTrackClasses>0 )
188  fUseDiceEfficiency = 0;
189  if(!fhEffH2 && fNTrackClasses>1 )
190  fUseDiceEfficiency = 0;
191  if(!fhEffH3 && fNTrackClasses>2 )
192  fUseDiceEfficiency = 0;
193  }
194 
195  fTracksOut->Delete();
196  SimulateTracks();
197  return kTRUE;
198 }
199 
200 //________________________________________________________________________
202 {
203  //Apply toy detector simulation to tracks
204  Int_t it = 0;
205  const Int_t nTracks = fTracks->GetEntriesFast();
206  for (Int_t i = 0; i < nTracks; ++i) {
207  AliPicoTrack *picotrack = static_cast<AliPicoTrack*>(fTracks->At(i));
208  if (!picotrack)
209  continue;
210 
211  Bool_t accept = kTRUE;
212  Double_t eff[3] = {0};
213  Double_t rnd = fRandom->Uniform(1.);
214  if(fUseDiceEfficiency) accept = DiceEfficiency(picotrack,eff,rnd);
215  if(!accept) continue;
216 
217  AliPicoTrack *track = NULL;
219  track = SmearPt(picotrack,eff,rnd);
220  (*fTracksOut)[it] = track;
221  } else
222  track = new ((*fTracksOut)[it]) AliPicoTrack(*picotrack);
223 
224  track->SetBit(TObject::kBitMask,1);
225  fHistPtDet->Fill(track->Pt());
226  it++;
227  }
228 }
229 
230 //________________________________________________________________________
232 {
233  // Dice to decide if particle is kept or not - toy model for efficiency
234  //
235  Double_t sumEff = 0.;
236  Double_t pT = 0.;
237 
238  if(fEfficiencyFixed<1.)
239  sumEff = fEfficiencyFixed;
240  else {
241  pT = vp->Pt();
242  Double_t pTtmp = pT;
243  if(pT>10.) pTtmp = 10.;
244  if(fhEffH1) eff[0] = fhEffH1->GetBinContent(fhEffH1->FindBin(pTtmp));
245  if(fhEffH2) eff[1] = fhEffH2->GetBinContent(fhEffH2->FindBin(pTtmp));
246  if(fhEffH3) eff[2] = fhEffH3->GetBinContent(fhEffH3->FindBin(pTtmp));
247 
248  sumEff = eff[0]+eff[1]+eff[2];
249  }
250  if(fUncertEfficiency!=1) sumEff=sumEff+fUncertEfficiency;
251  fp1Efficiency->Fill(vp->Pt(),sumEff);
252  if(rnd>sumEff && pT > fDiceEfficiencyMinPt) return kFALSE;
253  return kTRUE;
254 }
255 
256 //________________________________________________________________________
258 {
259  //Smear momentum of generated particle
260  Double_t smear = 1.;
261  Double_t pT = vp->Pt();
262  //Select hybrid track category
263 
264  //Sort efficiencies from large to small
265  Int_t cat[3] = {0};
266  TMath::Sort(3,eff,cat);
267  if(rnd<=eff[cat[2]])
268  smear = GetMomentumSmearing(cat[2],pT);
269  else if(rnd<=(eff[cat[2]]+eff[cat[1]]))
270  smear = GetMomentumSmearing(cat[1],pT);
271  else
272  smear = GetMomentumSmearing(cat[0],pT);
273 
274  fp1PtResolution->Fill(vp->Pt(),smear);
275 
276  Double_t sigma = vp->Pt()*smear;
277  Double_t pTrec = fRandom->Gaus(vp->Pt(),sigma);
278  fh2PtGenPtSmeared->Fill(vp->Pt(),pTrec);
279 
280  AliPicoTrack *picotrack = new AliPicoTrack(pTrec,
281  vp->Eta(),
282  vp->Phi(),
283  vp->Charge(),
284  vp->GetLabel(),
286  vp->GetTrackEtaOnEMCal(),
287  vp->GetTrackPhiOnEMCal(),
288  vp->GetTrackPtOnEMCal(),
289  vp->IsEMCAL(),
290  0.13957); //assume pion mass
291  return picotrack;
292 }
293 
294 //________________________________________________________________________
296 
297  //
298  // Get smearing on generated momentum
299  //
300 
301  TProfile *fMomRes = 0x0;
302  if(cat==1 && fMomResH1) fMomRes = (TProfile*)fMomResH1->Clone("fMomRes");
303  if(cat==2 && fMomResH2) fMomRes = (TProfile*)fMomResH2->Clone("fMomRes");
304  if(cat==3 && fMomResH3) fMomRes = (TProfile*)fMomResH3->Clone("fMomRes");
305 
306  if(!fMomRes)
307  return 0.;
308 
309  Double_t smear = 0.;
310  if(pt>20.) {
311  if(cat==1 && fMomResH1Fit) smear = fMomResH1Fit->Eval(pt);
312  if(cat==2 && fMomResH2Fit) smear = fMomResH2Fit->Eval(pt);
313  if(cat==3 && fMomResH3Fit) smear = fMomResH3Fit->Eval(pt);
314  }
315  else {
316  Int_t bin = fMomRes->FindBin(pt);
317  smear = fRandom->Gaus(fMomRes->GetBinContent(bin),fMomRes->GetBinError(bin));
318  }
319 
320  if(fMomRes) delete fMomRes;
321 
322  return smear;
323 }
324 
325 //________________________________________________________________________
327 
328  if (!gGrid && fPathTrPtResolution.Contains("alien://")) {
329  AliInfo("Trying to connect to AliEn ...");
330  TGrid::Connect("alien://");
331  }
332 
333  TFile *f = TFile::Open(fPathTrPtResolution.Data());
334  if(!f)return;
335  TProfile *fProfPtPtSigma1PtGlobSt = NULL;
336  TProfile *fProfPtPtSigma1PtGlobCnoSPD = NULL;
337  TProfile *fProfPtPtSigma1PtGlobCnoITS = NULL;
338  if(fNTrackClasses>0) fProfPtPtSigma1PtGlobSt = (TProfile*)f->Get("fProfPtPtSigma1PtGlobSt");
339  if(fNTrackClasses>1) fProfPtPtSigma1PtGlobCnoSPD = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoSPD");
340  if(fNTrackClasses>2) fProfPtPtSigma1PtGlobCnoITS = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoITS");
341 
342  SetSmearResolution(kTRUE);
343  SetMomentumResolutionHybrid(fProfPtPtSigma1PtGlobSt,fProfPtPtSigma1PtGlobCnoSPD,fProfPtPtSigma1PtGlobCnoITS);
344 }
345 
346 //________________________________________________________________________
348 
349  if (!gGrid && fPathTrPtResolution.Contains("alien://")) {
350  AliInfo("Trying to connect to AliEn ...");
351  TGrid::Connect("alien://");
352  }
353 
354  TFile *f = TFile::Open(fPathTrEfficiency.Data());
355  if(!f)return;
356  TH1D *hEffPosGlobSt = NULL;
357  TH1D *hEffPosGlobCnoSPD = NULL;
358  TH1D *hEffPosGlobCnoITS = NULL;
359  if(fNTrackClasses>0) hEffPosGlobSt = (TH1D*)f->Get("hEffPosGlobSt");
360  if(fNTrackClasses>1) hEffPosGlobCnoSPD = (TH1D*)f->Get("hEffPosGlobCnoSPD");
361  if(fNTrackClasses>2) hEffPosGlobCnoITS = (TH1D*)f->Get("hEffPosGlobCnoITS");
362 
363  SetDiceEfficiency(kTRUE);
364  SetEfficiencyHybrid(hEffPosGlobSt,hEffPosGlobCnoSPD,hEffPosGlobCnoITS);
365 }
366 
367 //________________________________________________________________________
368 void AliJetFastSimulation::SetMomentumResolutionHybrid(TProfile *p1, TProfile *p2, TProfile *p3) {
369  //
370  // set mom res profiles
371  //
372  if(fMomResH1) delete fMomResH1;
373  if(fMomResH2) delete fMomResH2;
374  if(fMomResH3) delete fMomResH3;
375 
376  if(p1) fMomResH1 = new TProfile(*p1);//(TProfile*)p1->Clone("fMomResH1");
377  if(p2) fMomResH2 = new TProfile(*p2);//(TProfile*)p2->Clone("fMomResH2");
378  if(p3) fMomResH3 = new TProfile(*p3);//(TProfile*)p3->Clone("fMomResH3");
379 }
380 
381 //________________________________________________________________________
383  //
384  // set tracking efficiency histos
385  //
386  if(h1) fhEffH1 = (TH1*)h1->Clone("fhEffH1");
387  if(h2) fhEffH2 = (TH1*)h2->Clone("fhEffH2");
388  if(h3) fhEffH3 = (TH1*)h3->Clone("fhEffH3");
389 }
390 
391 //________________________________________________________________________
393  //
394  // Fit linear function on momentum resolution at high pT
395  //
396 
397  if(!fMomResH1Fit && fMomResH1) {
398  fMomResH1Fit = new TF1("fMomResH1Fit","[0]+[1]*x",0.,200.);
399  fMomResH1->Fit(fMomResH1Fit,"LL V0","",5.,30.);
400  fMomResH1Fit ->SetRange(5.,100.);
401  }
402 
403  if(!fMomResH2Fit && fMomResH2) {
404  fMomResH2Fit = new TF1("fMomResH2Fit","[0]+[1]*x",0.,200.);
405  fMomResH2->Fit(fMomResH2Fit,"LL V0","",5.,30.);
406  fMomResH2Fit ->SetRange(5.,100.);
407  }
408 
409  if(!fMomResH3Fit && fMomResH3) {
410  fMomResH3Fit = new TF1("fMomResH3Fit","[0]+[1]*x",0.,200.);
411  fMomResH3->Fit(fMomResH3Fit,"LL V0","",5.,30.);
412  fMomResH3Fit ->SetRange(5.,100.);
413  }
414 
415 }
416 
417 
Double_t GetTrackPhiOnEMCal() const
Definition: AliPicoTrack.h:62
Bool_t DiceEfficiency(AliPicoTrack *vp, Double_t eff[3], Double_t rnd)
Int_t fNTrackClasses
output track collection
double Double_t
Definition: External.C:58
Definition: External.C:236
void SetDiceEfficiency(Int_t b)
Base task in the EMCAL framework.
virtual void UserCreateOutputObjects()
Double_t GetTrackPtOnEMCal() const
Definition: AliPicoTrack.h:64
Double_t fEfficiencyFixed
random number generator
void SetEfficiencyHybrid(TH1 *h1, TH1 *h2, TH1 *h3)
TProfile * fp1Efficiency
Control histo smeared momentum.
Bool_t IsEMCAL() const
Definition: AliPicoTrack.h:53
TH2F * fh2PtGenPtSmeared
pT spectrum of detector level particles
AliPicoTrack * SmearPt(AliPicoTrack *vp, Double_t eff[3], Double_t rnd)
Double_t * sigma
Int_t GetLabel() const
Definition: AliPicoTrack.h:41
Double_t Eta() const
Definition: AliPicoTrack.h:37
int Int_t
Definition: External.C:63
Definition: External.C:212
Double_t Phi() const
Definition: AliPicoTrack.h:33
void SetMomentumResolutionHybrid(TProfile *p1, TProfile *p2, TProfile *p3)
Byte_t GetTrackType() const
Definition: AliPicoTrack.h:43
Double_t Pt() const
Definition: AliPicoTrack.h:23
AliEmcalList * fOutput
!output list
TClonesArray * fTracks
!tracks
Bool_t Run()
Run function. This is the core function of the analysis and contains the user code. Therefore users have to implement this function.
void SetMakeGeneralHistograms(Bool_t g)
virtual void ExecOnce()
Perform steps needed to initialize the analysis.
Double_t GetMomentumSmearing(Int_t cat, Double_t pt)
void ExecOnce()
Perform steps needed to initialize the analysis.
void SetSmearResolution(Bool_t b)
void UserCreateOutputObjects()
Main initialization function on the worker.
Double_t GetTrackEtaOnEMCal() const
Definition: AliPicoTrack.h:63
bool Bool_t
Definition: External.C:53
TProfile * fp1PtResolution
Control profile efficiency.
Definition: External.C:196
Short_t Charge() const
Definition: AliPicoTrack.h:39