AliPhysics  a9863a5 (a9863a5)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
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  // exclude a fraction of the entries -- doesn't work very well
240  //if(fRandomEntry){
241  //
242  // if(fCurrentEntry < nentries) fTreeJet4Vect->GetEntry(fCurrentEntry);
243  // else {
244  // fCurrentEntry = 0;
245  // AliWarning("Starting from first entry again");
246  // fTreeJet4Vect->GetEntry(fCurrentEntry);
247  // }
248  // pTemb = jetDet->Pt();
249  //
250  // Float_t downscl = GetDownscalinigFactor();
251  // Float_t random = gRandom->Rndm();
252  //
253  // while (random > downscl){
254  // fCurrentEntry++;
255  // random = gRandom->Rndm();
256  // if(fCurrentEntry < nentries) fTreeJet4Vect->GetEntry(fCurrentEntry);
257  // else {
258  // fCurrentEntry = 0;
259  // AliWarning("Starting from first entry again");
260  // fTreeJet4Vect->GetEntry(fCurrentEntry);
261  // }
262  // pTemb = jetDet->Pt();
263  // if(pTemb < fPtRanges[fCurrentBin]) {
264  // random = gRandom->Rndm();
265  //
266  // }
267  //
268  // }
269  //
270  //}
271 
273  // Add the track that complies with the settings
274  AddTrack(jetDet->Pt(), jetDet->Eta(), jetDet->Phi(),0,0,0,0,kFALSE, fCurrentEntry, charge, jetDet->M());
275  //Printf("Embedded det %.2f, part %.2f", jetDet->Pt(), jetPar->Pt());
276 
277  Double_t x[4] = {jetDet->M() - jetPar->M(), jetDet->Pt() - jetPar->Pt(), jetPar->M(), jetPar->Pt()};
278  fhPartJet->Fill(x);
279 
280  fhEtaPart->Fill(jetPar->Eta());
281  fhPhiPart->Fill(jetPar->Phi());
282  fCount++; // count the number of track embedded in the current pT range
283  fCurrentEntry++; //increase for next iteration
284 
285  } else { //other inputs
286 
287  Double_t mass = fMass;
288  Short_t charge = 1;
289  if(fNeutralFraction>0.) {
290  Double_t rnd = gRandom->Rndm();
291  if(rnd<fNeutralFraction) {
292  charge = 0;
293  mass = fNeutralMass;
294  }
295  }
296  if(fMassless) mass = 0.;
297  if(fMassFromDistr) mass = -999;
298 
299  AddTrack(-999,-999,-999,0,0,0,0,kFALSE,0,charge,mass);
300  }
301  }
302  }
303 
305 }
306 
307 //________________________________________________________________________
308 
310 
311  if(fCount > fNevPerBin) {
312  //Printf("%d = fCount, Increasing fCurrentBin %d -> %d", fCount, fCurrentBin, fCurrentBin+1);
313  fCurrentBin++;
314  fCount = 0;
315 
316  }
317 
318  if (fCurrentBin >= fNBins) {
319  AliError(Form("Bin %d out of bound %d, set to fNBins - 1 = %d", fCurrentBin, fNBins, fNBins - 1));
321  if (fCurrentBin < 0) fCurrentBin = fNBins - 1;
322  fGoBack++;
323  }
324  return fDownscale[fCurrentBin];
325 }
326 
327 
328 //________________________________________________________________________
329 void AliJetEmbeddingTask::SetTree(TTree *tree) {
330  if(!tree){
331  AliError("Null tree");
332  return;
333  }
334  fFromTree = kTRUE;
335  fTreeJet4Vect = (TTree*)tree->Clone(Form("%sCpEmb", tree->GetName()));
336  AliInfo(Form("Input tree set %d (%p -> %p)", fTreeJet4Vect->GetNbranches(), tree, fTreeJet4Vect));
337 
338  return;
339 }
340 
341 //________________________________________________________________________
342 
343 void AliJetEmbeddingTask::SetTreeFromFile(TString filename, TString treename){
344 
345  if(filename.Contains("alien")) {
346  TGrid::Connect("alien://");
347  }
348  TFile *f = TFile::Open(filename);
349  if(!f){
350  Printf("File %s not found, cannot SetTree", filename.Data());
351  return;
352  }
353 
354  TTree *tree = dynamic_cast<TTree*>(f->Get(treename));
355  if(!tree){
356  Printf("Tree %s not found!!!", treename.Data());
357  f->ls();
358  return;
359  }
360  SetTree(tree);
361 
362  //f->Close();
363  //delete f;
364 
365  return;
366 }
367 
368 //________________________________________________________________________
369 
371 
372  if(fNBins == 0){
373  AliError("Set number of bins first! SetNBinsEmbedding(nbins);");
374  return;
375  }
376  Printf("Creating array of factors with size %d", fNBins);
377  fDownscale = new Float_t[fNBins];
378  for(Int_t i = 0; i<fNBins; i++) {
379  fDownscale[i] = rej[i];
380  Printf("Bin %d -> Factor = %e", i, fDownscale[i]);
381  }
382  return;
383 }
384 
385 
387  Printf("fGoBack = %d", fGoBack);
388 }
389 
Int_t charge
Class for track embedding into an event.
Int_t fCurrentEntry
Current TTree entry.
TString fBranchJDetName
name of the detector level jet branch in the TTree
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:26
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
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
Bool_t fMassless
make particles massless
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
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
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)
ClassImp(AliJetEmbeddingTask) AliJetEmbeddingTask
void SetMassAndPtDistributionFromFile(TString filenameM, TString filenamepT, TString histonameM, TString histonamepT)
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