AliPhysics  41af4b0 (41af4b0)
AliJetEmbeddingTask.cxx
Go to the documentation of this file.
1 // $Id$
2 //
3 // Jet embedding task.
4 //
5 // Author: S.Aiola, C.Loizides
6 
7 #include <TRandom3.h>
8 #include <TFile.h>
9 #include <AliLog.h>
10 #include <TGrid.h>
11 #include <THnSparse.h>
12 #include <AliPicoTrack.h>
13 
14 #include "AliJetEmbeddingTask.h"
15 
16 ClassImp(AliJetEmbeddingTask)
17 
18 //________________________________________________________________________
21  fMassless(kFALSE),
22  fNeutralFraction(0),
23  fNeutralMass(0.135),
24  fMass(0.1396),
25  fPathMinputFile(""),
26  fPathpTinputFile(""),
27  fMinputName(""),
28  fpTinputName(""),
29  fFromTree(0),
30  fPathTreeinputFile(""),
31  fTreeinputName("fTreeJet"),
32  fBranchJDetName("fJetDet."),
33  fBranchJParName("fJetPart."),
34  fTreeJet4Vect(0),
35  fCurrentEntry(0),
36  fInput(0),
37  fRandomEntry(0),
38  fNBins(10),
39  fDownscale(0),
40  fPtRanges(0),
41  fNevPerBin(0),
42  fCount(0),
43  fCurrentBin(0),
44  fGoBack(0),
45  fhPartJet(0),
46  fhEtaPart(0),
47  fhPhiPart(0),
48  fhTreeEntriesUsed(0),
49  fNTreeEntries(0),
50  fOldCutOnDetPt(0)
51 
52 {
53  // Default constructor.
54  SetSuffix("Embedded");
55  DefineOutput(1, TList::Class());
56 }
57 
58 //________________________________________________________________________
60  AliJetModelBaseTask(name, kTRUE),
61  fMassless(kFALSE),
62  fNeutralFraction(0),
63  fNeutralMass(0.135),
64  fMass(0.1396),
65  fPathMinputFile(""),
66  fPathpTinputFile(""),
67  fMinputName(""),
68  fpTinputName(""),
69  fFromTree(0),
70  fPathTreeinputFile(""),
71  fTreeinputName("fTreeJet"),
72  fBranchJDetName("fJetDet."),
73  fBranchJParName("fJetPart."),
74  fTreeJet4Vect(0),
75  fCurrentEntry(0),
76  fInput(0),
77  fRandomEntry(0),
78  fNBins(10),
79  fDownscale(0),
80  fPtRanges(0),
81  fNevPerBin(0),
82  fCount(0),
83  fCurrentBin(0),
84  fGoBack(0),
85  fhPartJet(0),
86  fhEtaPart(0),
87  fhPhiPart(0),
88  fhTreeEntriesUsed(0),
89  fNTreeEntries(0),
90  fOldCutOnDetPt(0)
91 {
92  // Standard constructor.
93  SetSuffix("Embedded");
94  DefineOutput(2, TList::Class());
95 }
96 
97 //________________________________________________________________________
99 {
100  // Destructor
101  delete[] fDownscale; fDownscale = 0;
102  delete[] fPtRanges; fPtRanges = 0;
103 }
104 
105 //________________________________________________________________________
106 
108 
110 
111  OpenFile(1);
112  fInput = new TList();
113  fInput->SetOwner();
114 
115  delete gRandom;
116  gRandom = new TRandom3(0); //random seed
117 
118  if(!fPathTreeinputFile.IsNull()){
120  if(!fTreeJet4Vect) AliFatal("Something went wrong in setting the tree");
121  Printf("Emb task");
122  TLorentzVector *detjet = 0;
123  fTreeJet4Vect->SetBranchAddress(fBranchJDetName, &detjet);
124  fTreeJet4Vect->GetEntry(0);
125  fTreeJet4Vect->Show();
126  fNTreeEntries = fTreeJet4Vect->GetEntries();
127  fhTreeEntriesUsed = new TH1F("fhTreeEntriesUsed", "Entries;Entry in TTree", fNTreeEntries, 0, fNTreeEntries-1);
129 
130  fCurrentEntry = gRandom->Integer(fNTreeEntries); //in each worker it starts from a different entry
131  Printf("Entries %lld, start from %d", fNTreeEntries, fCurrentEntry);
132  }
133 
134  if(!fPathMinputFile.IsNull() && fPathpTinputFile.IsNull()){
136  if(!fHMassDistrib) AliFatal("Something went wrong in setting the M distribution");
137  fInput->Add(fHMassDistrib);
138  }
139 
140  if(fPathMinputFile.IsNull() && !fPathpTinputFile.IsNull()){
142  if(!fPtSpectrum) AliFatal("Something went wrong in setting the pT distribution");
143  fInput->Add(fPtSpectrum);
144  }
145 
146  if(!fPathMinputFile.IsNull() && !fPathpTinputFile.IsNull()){
148  if(!fHMassDistrib) AliFatal("Something went wrong in setting the M distribution");
149  if(!fPtSpectrum) AliFatal("Something went wrong in setting the pT distribution");
150  fInput->Add(fHMassDistrib);
151  fInput->Add(fPtSpectrum);
152  }
153  PostData(2, fInput);
154 
155  const Int_t ndim = 4;
156  Int_t nbins = 80, mind = -30, maxd = 30, minm = 0, maxm = 80, minpt = 0, maxpt = 120;
157  const Int_t nbinshnsp[ndim] = {nbins, nbins, nbins, nbins};
158  const Double_t minhnsp[ndim] = {(Double_t)mind, (Double_t)mind, (Double_t)minm, (Double_t)minpt};
159  const Double_t maxhnsp[ndim] = {(Double_t)maxd, (Double_t)maxd, (Double_t)maxm, (Double_t)maxpt};
160  TString title = "Check part level;#it{M}_{det} - #it{M}_{par} (GeV); #it{p}_{T, det} - #it{p}_{T, par} (GeV/#it{c}); #it{M} (GeV); #it{p}_{T} (GeV/#it{c})";
161 
162  fhPartJet = new THnSparseF("fhPartJet", title.Data(), ndim, nbinshnsp, minhnsp, maxhnsp);
163  fhPartJet->Sumw2();
164  fOutput->Add(fhPartJet);
165 
166  fhEtaPart = new TH1F("fhEtaPart" , "#eta distribution part; #eta", 100, -1, 1);
167  fhEtaPart->Sumw2();
168  fOutput->Add(fhEtaPart);
169  fhPhiPart = new TH1F("fhPhiPart" , "#varphi distribution; #varphi", 100, (-1)*TMath::Pi(), TMath::Pi());
170  fhPhiPart->Sumw2();
171  fOutput->Add(fhPhiPart);
172 
173 }
174 
175 //________________________________________________________________________
177 {
178  // Embed particles.
179  if (fNClusters > 0 && fOutClusters) {
180  if (fCopyArray)
181  CopyClusters();
182  for (Int_t i = 0; i < fNClusters; ++i) {
183  AddCluster();
184  }
185  }
186 
187  if (fNTracks > 0 && fOutTracks) {
188  if (fCopyArray)
189  CopyTracks();
190  for (Int_t i = 0; i < fNTracks; ++i) {
191 
192  Short_t charge = 1;
193  if(fNeutralFraction>0.) {
194  Double_t rnd = gRandom->Rndm();
195  if(rnd<fNeutralFraction) {
196  charge = 0;
197  //mass = fNeutralMass;
198  }
199  }
200  // Add track from tree of 4-vectors (jet reco) and save the particle level somewhere
201  if(fFromTree){
202 
203  if(!fTreeJet4Vect || fBranchJDetName.IsNull()) {
204  AliFatal(Form("Tree or branch name not found"));
205  }
206  TLorentzVector *jetDet = 0;
207  TLorentzVector *jetPar = 0;
208  fTreeJet4Vect->ResetBranchAddresses();
209  //fTreeJet4Vect->SetBranchAddress(fBranchJDetName.Data(), &jetDet, &bDet);
210  fTreeJet4Vect->SetBranchAddress(fBranchJDetName.Data(), &jetDet);
211  fTreeJet4Vect->SetBranchAddress(fBranchJParName.Data(), &jetPar);
212 
213  Double_t pTemb = -1, pTpar = -1;
215  else {
216  fCurrentEntry = 0;
217  AliWarning("Starting from first entry again");
218  fTreeJet4Vect->GetEntry(fCurrentEntry);
219  }
220 
221  pTpar = jetPar->Pt(); //cut on the particle level pT
222  if(fOldCutOnDetPt) pTpar = jetDet->Pt(); //cut on the detector level pT (as it was before, keeping the functionality)
223  // selected pT range
224  if((fPtMin != 0 && fPtMax != 0) && !fRandomEntry) {
225  Int_t countWhile = 0;
226  while(!(pTpar > fPtMin && pTpar < fPtMax) && countWhile < fNTreeEntries){
227  fCurrentEntry++;
229  else {
230  fCurrentEntry = 0;
231  AliWarning("Starting from first entry again");
232  fTreeJet4Vect->GetEntry(fCurrentEntry);
233  }
234  pTpar = jetPar->Pt(); //cut on the particle level pT
235  if(fOldCutOnDetPt) pTpar = jetDet->Pt(); //cut on the detector level pT (as it was before, keeping the functionality)
236  countWhile++;
237  }
238  }
239 
241  // Add the track that complies with the settings
242  AddTrack(jetDet->Pt(), jetDet->Eta(), jetDet->Phi(),0,0,0,0,kFALSE, fCurrentEntry, charge, jetDet->M());
243  //Printf("Embedded det %.2f, part %.2f", jetDet->Pt(), jetPar->Pt());
244 
245  Double_t x[4] = {jetDet->M() - jetPar->M(), jetDet->Pt() - jetPar->Pt(), jetPar->M(), jetPar->Pt()};
246  fhPartJet->Fill(x);
247 
248  fhEtaPart->Fill(jetPar->Eta());
249  fhPhiPart->Fill(jetPar->Phi());
250  fCount++; // count the number of track embedded in the current pT range
251  fCurrentEntry++; //increase for next iteration
252 
253  } else { //other inputs
254 
255  Double_t mass = fMass;
256  Short_t charge = 1;
257  if(fNeutralFraction>0.) {
258  Double_t rnd = gRandom->Rndm();
259  if(rnd<fNeutralFraction) {
260  charge = 0;
261  mass = fNeutralMass;
262  }
263  }
264  if(fMassless) mass = 0.;
265  if(fMassFromDistr) mass = -999;
266 
267  AddTrack(-999,-999,-999,0,0,0,0,kFALSE,0,charge,mass);
268  }
269  }
270  }
271 
273 }
274 
275 //________________________________________________________________________
276 
278 
279  if(fCount > fNevPerBin) {
280  //Printf("%d = fCount, Increasing fCurrentBin %d -> %d", fCount, fCurrentBin, fCurrentBin+1);
281  fCurrentBin++;
282  fCount = 0;
283 
284  }
285 
286  if (fCurrentBin >= fNBins) {
287  AliError(Form("Bin %d out of bound %d, set to fNBins - 1 = %d", fCurrentBin, fNBins, fNBins - 1));
289  if (fCurrentBin < 0) fCurrentBin = fNBins - 1;
290  fGoBack++;
291  }
292  return fDownscale[fCurrentBin];
293 }
294 
295 
296 //________________________________________________________________________
298  if(!tree){
299  AliError("Null tree");
300  return;
301  }
302  fFromTree = kTRUE;
303  fTreeJet4Vect = (TTree*)tree->Clone(Form("%sCpEmb", tree->GetName()));
304  AliInfo(Form("Input tree set %d (%p -> %p)", fTreeJet4Vect->GetNbranches(), tree, fTreeJet4Vect));
305 
306  return;
307 }
308 
309 //________________________________________________________________________
310 
312 
313  if(filename.Contains("alien")) {
314  TGrid::Connect("alien://");
315  }
316  TFile *f = TFile::Open(filename);
317  if(!f){
318  Printf("File %s not found, cannot SetTree", filename.Data());
319  return;
320  }
321 
322  TTree *tree = dynamic_cast<TTree*>(f->Get(treename));
323  if(!tree){
324  Printf("Tree %s not found!!!", treename.Data());
325  f->ls();
326  return;
327  }
328  SetTree(tree);
329 
330  //f->Close();
331  //delete f;
332 
333  return;
334 }
335 
336 //________________________________________________________________________
337 
339 
340  if(fNBins == 0){
341  AliError("Set number of bins first! SetNBinsEmbedding(nbins);");
342  return;
343  }
344  Printf("Creating array of factors with size %d", fNBins);
345  fDownscale = new Float_t[fNBins];
346  for(Int_t i = 0; i<fNBins; i++) {
347  fDownscale[i] = rej[i];
348  Printf("Bin %d -> Factor = %e", i, fDownscale[i]);
349  }
350  return;
351 }
352 
353 
355  Printf("fGoBack = %d", fGoBack);
356 }
357 
Int_t charge
Class for track embedding into an event.
Int_t fCurrentEntry
Current TTree entry.
const char * filename
Definition: TestFCM.C:1
TString fBranchJDetName
name of the detector level jet branch in the TTree
double Double_t
Definition: External.C:58
Int_t fCount
counts number of embedded tracks in the current pT bin
void SetpTDistributionFromFile(TString filename, TString histoname)
TString fPathpTinputFile
path to the file where the external input pT distribution is (can be from alien)
const char * title
Definition: MakeQAPdf.C:27
THnSparse * fhPartJet
! control histogram particle level correponsding to embedded from tree
void SetTree(TTree *tree)
Float_t fPtMin
pt minimum value
Int_t fNevPerBin
number of embedded tracks per pT bins
Double_t mass
TList * fOutput
! output list for QA histograms
TClonesArray * fOutClusters
! output cluster collection
TTree * fTreeJet4Vect
! tree containing the jet 4-vectors (input for embed.)
Float_t * fDownscale
factor to exclude randomly entries from the embedding
TH1F * fhPhiPart
! Phi particle corresponding to embedded from tree
Base class for embedding into an event.
Int_t fCurrentBin
the current pT bin
Long64_t fNTreeEntries
Tree entry number.
TRandom * gRandom
TString fPathTreeinputFile
path to the file where the external input Tree is (can be from alien)
Int_t fGoBack
how many times fCurrentBin is set back to 0
TClonesArray * fOutTracks
! output track collection
void Terminate(Option_t *option="")
Double_t fMass
assign this mass to particle
void SetTreeFromFile(TString filenameM, TString treename)
Int_t fNTracks
how many tracks are being processed
Float_t * fPtRanges
low lims of the pT ranges corresponding to the fDownscale bin
int Int_t
Definition: External.C:63
Bool_t fMassless
make particles massless
float Float_t
Definition: External.C:68
Bool_t fOldCutOnDetPt
Apply cut on the embedded track at detector level (this was default before, not we apply the cut at p...
TString fMinputName
name of the external input Mass histogram
Double_t M() const
Definition: AliPicoTrack.h:36
Double_t fNeutralMass
assign this mass to neutral particles
Int_t fNClusters
how many clusters are being processed
TH1F * fhEtaPart
! Eta particle corresponding to embedded from tree
Bool_t fFromTree
if true use TTree of jet 4-vectors as input single tracks
TList * fInput
! Input histograms saved in this list
Double_t fNeutralFraction
assign charge==0 to fraction of particles
TH1F * fhTreeEntriesUsed
! Entries of the TTree used
Bool_t fCopyArray
whether or not the array will be copied to a new one before modelling
TString fpTinputName
name of the external input pT histogram
short Short_t
Definition: External.C:23
Bool_t fRandomEntry
draw random number to extract the entry number in the tree
TString fTreeinputName
name of the external input Tree
void FillHistograms()
do jet model action
AliVCluster * AddCluster(Double_t e=-1, Double_t eta=-999, Double_t phi=-1, Int_t label=0)
add a cell with given energy, position and times
Bool_t fMassFromDistr
draw the particle mass from fHMassDistrib
TH1F * fPtSpectrum
pt spectrum to extract random pt values
TH1F * fHMassDistrib
shape of mass distribution of embedded tracks
const char Option_t
Definition: External.C:48
void SetSuffix(const char *s)
const Int_t nbins
AliPicoTrack * AddTrack(Double_t pt=-999, Double_t eta=-999, Double_t phi=-999, Byte_t type=0, Double_t etaemc=0, Double_t phiemc=0, Double_t ptemc=0, Bool_t ise=kFALSE, Int_t label=0, Short_t charge=1, Double_t mass=0.1396)
add a cluster (copy)
TString fBranchJParName
name of the particle level jet branch in the TTree
void SetRejection(Float_t *rej)
void SetMassDistributionFromFile(TString filename, TString histoname)
void SetMassAndPtDistributionFromFile(TString filenameM, TString filenamepT, TString histonameM, TString histonamepT)
TList * OpenFile(const char *fname)
Definition: DrawAnaELoss.C:65
TString fPathMinputFile
path to the file where the external input Mass distribution is (can be from alien) ...
void Run()
intialize task
Float_t fPtMax
pt maximum value