AliPhysics  8d00e07 (8d00e07)
AliAnalysisTaskParticleRandomizer.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2016, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: R. Haake. *
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 #include <iostream>
17 #include <TRandom3.h>
18 #include <AliLog.h>
19 #include <TString.h>
20 #include <TMath.h>
21 #include <TClonesArray.h>
22 #include <AliAODTrack.h>
23 #include <AliPicoTrack.h>
24 #include <AliEmcalJet.h>
25 #include <AliRhoParameter.h>
26 #include "TH1D.h"
27 #include "TH2D.h"
28 #include "AliAnalysisTaskEmcal.h"
30 
34 
35 //_____________________________________________________________________________________________________
37  AliAnalysisTaskEmcal("AliAnalysisTaskParticleRandomizer", kFALSE), fRandomizeInPhi(1), fRandomizeInEta(0), fRandomizeInTheta(0), fRandomizeInPt(0), fTrackEfficiency(1.0), fMinPhi(0), fMaxPhi(TMath::TwoPi()), fMinEta(-0.9), fMaxEta(+0.9), fMinPt(0), fMaxPt(120), fDistributionV2(0), fDistributionV3(0), fDistributionV4(0), fDistributionV5(0), fInputArrayName(), fOutputArrayName(), fInputArray(0), fOutputArray(0), fJetRemovalRhoObj(), fJetRemovalArrayName(), fJetRemovalArray(0), fJetRemovalPtThreshold(999.), fJetRemovalNLeadingJets(0), fJetEmbeddingArrayName(), fJetEmbeddingArray(0), fRandomPsi3(0), fLeadingJet(0), fSubleadingJet(0), fRandom()
38 {
39 // constructor
40 }
41 
42 //_____________________________________________________________________________________________________
44 {
45 // destructor
46 }
47 
48 
49 //_____________________________________________________________________________________________________
51 {
52  // Check user input
53  if(fInputArrayName.IsNull())
54  AliWarning(Form("Name of input array not given!"));
55  if(fOutputArrayName.IsNull())
56  AliFatal(Form("Name of output array not given!"));
57 
58  fRandom = new TRandom3(0);
59 }
60 
61 //_____________________________________________________________________________________________________
63 {
64  // Check if arrays are OK
65  fInputArray = static_cast<TClonesArray*>(InputEvent()->FindListObject(Form("%s", fInputArrayName.Data())));
66  if(!fInputArrayName.IsNull() && !fInputArray)
67  AliFatal(Form("Input array '%s' not found!", fInputArrayName.Data()));
68 
69  // On demand, load also jets
70  if(!fJetRemovalArrayName.IsNull())
71  {
72  fJetRemovalArray = static_cast<TClonesArray*>(InputEvent()->FindListObject(Form("%s", fJetRemovalArrayName.Data())));
73  if(!fJetRemovalArray)
74  AliError(Form("Jet array '%s' demanded but not found in event!", fJetRemovalArrayName.Data()));
75  }
76 
77  // On demand, load array for embedding
78  if(!fJetEmbeddingArrayName.IsNull())
79  {
80  fJetEmbeddingArray = static_cast<TClonesArray*>(InputEvent()->FindListObject(Form("%s", fJetEmbeddingArrayName.Data())));
82  AliError(Form("Embedding array '%s' demanded but not found in event!", fJetEmbeddingArrayName.Data()));
83  }
84 
85  if((InputEvent()->FindListObject(Form("%s", fOutputArrayName.Data()))))
86  AliFatal(Form("Output array '%s' already exists in the event! Rename it.", fInputArrayName.Data()));
87 
88  if(fInputArray)
89  if(strcmp(fInputArray->GetClass()->GetName(), "AliAODTrack"))
90  AliError(Form("Track type %s not yet supported. Use AliAODTrack", fInputArray->GetClass()->GetName()));
91 
92  // Copy the input array to the output array
93  fOutputArray = new TClonesArray("AliAODTrack");
94  fOutputArray->SetName(fOutputArrayName.Data());
95  InputEvent()->AddObject(fOutputArray);
96 
98 }
99 
100 //_____________________________________________________________________________________________________
102 {
103  fRandomPsi3 = fRandom->Rndm()*TMath::Pi(); // once per event, create a random value dedicated for Psi3
104  fRandomPsi4 = fRandom->Rndm()*TMath::Pi(); // once per event, create a random value dedicated for Psi4
105  fRandomPsi5 = fRandom->Rndm()*TMath::Pi(); // once per event, create a random value dedicated for Psi5
106 
107  // Get leading jets on demand
110 
111  Int_t accTracks = 0;
112 
113  // Add events in input array
114  if(fInputArray)
115  for(Int_t iPart=0; iPart<fInputArray->GetEntries(); iPart++)
116  {
117  // Remove particles contained in jet array (on demand)
118  if(fJetRemovalArray && IsParticleInJet(iPart))
119  continue;
120 
121  // Take only particles from the randomization acceptance
122  AliAODTrack* inputParticle = static_cast<AliAODTrack*>(fInputArray->At(iPart));
123  if(fRandomizeInPhi && (inputParticle->Phi() < fMinPhi || inputParticle->Phi() >= fMaxPhi) )
124  continue;
125  if( (fRandomizeInTheta || fRandomizeInEta) && (inputParticle->Eta() < fMinEta || inputParticle->Eta() >= fMaxEta) )
126  continue;
127 
128  // Discard tracks due to lowered tracking efficiency
129  if (fTrackEfficiency < 1.0)
130  if (fTrackEfficiency < fRandom->Rndm())
131  continue;
132 
133  new ((*fOutputArray)[accTracks]) AliAODTrack(*((AliAODTrack*)fInputArray->At(iPart)));
134 
135  // Randomize on demand
136  AliAODTrack* particle = static_cast<AliAODTrack*>(fOutputArray->At(accTracks));
137  RandomizeTrack(particle);
138 
139  accTracks++;
140  }
141 
142  // Add particles for embedding (on demand)
144  for(Int_t iPart=0; iPart<fJetEmbeddingArray->GetEntries(); iPart++)
145  {
146  // Take only particles from the randomization acceptance
147  AliPicoTrack* inputParticle = static_cast<AliPicoTrack*>(fJetEmbeddingArray->At(iPart));
148  if(fRandomizeInPhi && (inputParticle->Phi() < fMinPhi || inputParticle->Phi() >= fMaxPhi) )
149  continue;
150  if( (fRandomizeInTheta || fRandomizeInEta) && (inputParticle->Eta() < fMinEta || inputParticle->Eta() >= fMaxEta) )
151  continue;
152 
153  // Discard tracks due to lowered tracking efficiency
154  if (fTrackEfficiency < 1.0)
155  if (fTrackEfficiency < fRandom->Rndm())
156  continue;
157 
158  new ((*fOutputArray)[accTracks]) AliAODTrack(*(GetAODTrack(inputParticle)));
159 
160  // Randomize on demand
161  AliAODTrack* particle = static_cast<AliAODTrack*>(fOutputArray->At(accTracks));
162  RandomizeTrack(particle);
163 
164  accTracks++;
165  }
166 
167  // std::cout << Form("%i particles from jets removed out of %i tracks. ", fInputArray->GetEntries()-accTracks, fInputArray->GetEntries()) << std::endl;
168  return kTRUE;
169 }
170 
171 //_____________________________________________________________________________________________________
173 {
174  if(fRandomizeInPhi)
175  particle->SetPhi(fMinPhi + fRandom->Rndm()*(fMaxPhi-fMinPhi));
177  {
178  Double_t minTheta = 2.*atan(exp(-fMinEta));
179  Double_t maxTheta = 2.*atan(exp(-fMaxEta));
180  particle->SetTheta(minTheta + fRandom->Rndm()*(maxTheta-minTheta));
181  }
182  if(fRandomizeInEta)
183  {
184  Double_t randomEta = fMinEta + fRandom->Rndm()*(fMaxEta-fMinEta);
185  Double_t randomTheta = 2.*atan(exp(-randomEta));
186  particle->SetTheta(randomTheta);
187  }
188 
189  if(fRandomizeInPt)
190  particle->SetPt(fMinPt + fRandom->Rndm()*(fMaxPt-fMinPt));
191 
193  particle->SetPhi(AddFlow(particle->Phi(), particle->Pt()));
194 }
195 
196 //_____________________________________________________________________________________________________
198 {
199  AliAODTrack* newTrack = new AliAODTrack();
200  newTrack->SetPt(track->Pt());
201  newTrack->SetTheta(2.*atan(exp(-track->Eta()))); // there is no setter for eta
202  newTrack->SetPhi(track->Phi());
203  newTrack->SetCharge(track->Charge());
204  newTrack->SetLabel(track->GetLabel());
205 
206  // Hybrid tracks (compatible with LHC11h)
207  UInt_t filterMap = BIT(8) | BIT(9);
208  newTrack->SetIsHybridGlobalConstrainedGlobal();
209  newTrack->SetFilterMap(filterMap);
210  return newTrack;
211 }
212 
213 //_____________________________________________________________________________________________________
215 {
216  // Check if particle w/ index 'part' is contained in one of the removal jets
217  // The removal jets can be jets with a pt theshold or the leading jets
218  for(Int_t i=0; i<fJetRemovalArray->GetEntries(); i++)
219  {
220  AliEmcalJet* tmpJet = static_cast<AliEmcalJet*>(fJetRemovalArray->At(i));
221 
222  // Check if to remove leading jets
224  {
225  // leading jet removal mode
226  if( (fJetRemovalNLeadingJets == 1) && (tmpJet == fLeadingJet))
227  {
228  if(tmpJet->ContainsTrack(part)>=0)
229  return kTRUE;
230  else
231  return kFALSE; // we know we are done here: the track is not contained in leading jet
232  }
233  // leading or leading removal mode
234  else if ( (fJetRemovalNLeadingJets == 2) && ((tmpJet == fLeadingJet) || (tmpJet == fSubleadingJet)))
235  {
236  if(tmpJet->ContainsTrack(part)>=0)
237  return kTRUE;
238  }
239  }
240  // Check if to remove jets above threshold
241  else
242  {
243  Double_t tmpPt = tmpJet->Pt() - tmpJet->Area()*GetExternalRho();
244  if(tmpPt >= fJetRemovalPtThreshold)
245  if(tmpJet->ContainsTrack(part)>=0)
246  return kTRUE;
247  }
248  }
249  return kFALSE;
250 }
251 
252 //_____________________________________________________________________________________________________
254 {
255  // Get rho from event.
256  AliRhoParameter* rho = 0;
257  if (!fJetRemovalRhoObj.IsNull())
258  {
259  rho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fJetRemovalRhoObj.Data()));
260  if (!rho) {
261  AliError(Form("%s: Could not retrieve rho with name %s!", GetName(), fJetRemovalRhoObj.Data()));
262  return 0;
263  }
264  }
265  else
266  return 0;
267 
268  return (rho->GetVal());
269 }
270 
271 
272 //_____________________________________________________________________________________________________
274 {
275  // adapted from AliFlowTrackSimple
276  Double_t precisionPhi = 1e-10;
277  Int_t maxNumberOfIterations = 200;
278 
279  Double_t phi0=phi;
280  Double_t f=0.;
281  Double_t fp=0.;
282  Double_t phiprev=0.;
283  Int_t ptBin = 0;
284 
285  // Evaluate V2 for track pt/centrality
286  Double_t v2 = 0;
287  if(fDistributionV2)
288  {
289  ptBin = fDistributionV2->GetXaxis()->FindBin(pt);
290  if(ptBin>fDistributionV2->GetNbinsX())
291  v2 = fDistributionV2->GetBinContent(fDistributionV2->GetNbinsX(), fDistributionV2->GetYaxis()->FindBin(fCent));
292  else if(ptBin>0)
293  v2 = fDistributionV2->GetBinContent(ptBin, fDistributionV2->GetYaxis()->FindBin(fCent));
294  }
295 
296  // Evaluate V3 for track pt/centrality
297  Double_t v3 = 0;
298  if(fDistributionV3)
299  {
300  ptBin = fDistributionV3->GetXaxis()->FindBin(pt);
301  if(ptBin>fDistributionV3->GetNbinsX())
302  v3 = fDistributionV3->GetBinContent(fDistributionV3->GetNbinsX(), fDistributionV3->GetYaxis()->FindBin(fCent));
303  else if(ptBin>0)
304  v3 = fDistributionV3->GetBinContent(ptBin, fDistributionV3->GetYaxis()->FindBin(fCent));
305  }
306 
307  // Evaluate V4 for track pt/centrality
308  Double_t v4 = 0;
309  if(fDistributionV4)
310  {
311  ptBin = fDistributionV4->GetXaxis()->FindBin(pt);
312  if(ptBin>fDistributionV4->GetNbinsX())
313  v4 = fDistributionV4->GetBinContent(fDistributionV4->GetNbinsX(), fDistributionV4->GetYaxis()->FindBin(fCent));
314  else if(ptBin>0)
315  v4 = fDistributionV4->GetBinContent(ptBin, fDistributionV4->GetYaxis()->FindBin(fCent));
316  }
317 
318  // Evaluate V5 for track pt/centrality
319  Double_t v5 = 0;
320  if(fDistributionV5)
321  {
322  ptBin = fDistributionV5->GetXaxis()->FindBin(pt);
323  if(ptBin>fDistributionV5->GetNbinsX())
324  v5 = fDistributionV5->GetBinContent(fDistributionV5->GetNbinsX(), fDistributionV5->GetYaxis()->FindBin(fCent));
325  else if(ptBin>0)
326  v5 = fDistributionV5->GetBinContent(ptBin, fDistributionV5->GetYaxis()->FindBin(fCent));
327  }
328 
329  // Add all v's
330  for (Int_t i=0; i<maxNumberOfIterations; i++)
331  {
332  phiprev=phi; //store last value for comparison
333  f = phi-phi0
334  + v2*TMath::Sin(2.*(phi-(fEPV0+(TMath::Pi()/2.))))
335  +2./3.*v3*TMath::Sin(3.*(phi-fRandomPsi3))
336  +0.5 *v4*TMath::Sin(4.*(phi-fRandomPsi4))
337  +0.4 *v5*TMath::Sin(5.*(phi-fRandomPsi5));
338 
339  fp = 1.0+2.0*(
340  +v2*TMath::Cos(2.*(phi-(fEPV0+(TMath::Pi()/2.))))
341  +v3*TMath::Cos(3.*(phi-fRandomPsi3))
342  +v4*TMath::Cos(4.*(phi-fRandomPsi4))
343  +v5*TMath::Cos(5.*(phi-fRandomPsi5))); //first derivative
344 
345  phi -= f/fp;
346  if (TMath::AreEqualAbs(phiprev,phi,precisionPhi)) break;
347  }
348 
349  return phi;
350 }
351 
352 //_____________________________________________________________________________________________________
354 {
355  // Customized from AliJetContainer::GetLeadingJet()
356  // Get the leading+subleading jet; the sorting is according to pt-A*rho
357 
358  jetLeading = 0;
359  jetSubLeading = 0;
360 
361  Double_t tmpLeadingPt = 0;
362  Double_t tmpSubleadingPt = 0;
363 
364  for(Int_t i=0; i<fJetRemovalArray->GetEntries(); i++)
365  {
366  AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJetRemovalArray->At(i));
367  if ( (jet->Pt()-jet->Area()*GetExternalRho()) > tmpLeadingPt )
368  {
369  jetSubLeading = jetLeading;
370  jetLeading = jet;
371  tmpSubleadingPt = tmpLeadingPt;
372  tmpLeadingPt = jet->Pt()-jet->Area()*GetExternalRho();
373  }
374  else if ( (jet->Pt()-jet->Area()*GetExternalRho()) > tmpSubleadingPt )
375  {
376  jetSubLeading = jet;
377  tmpSubleadingPt = jet->Pt()-jet->Area()*GetExternalRho();
378  }
379  }
380 }
Double_t fMaxPt
range for Pt for randomization
TH2D * fDistributionV4
Distribution for v3 in bins of pt and centrality.
Double_t fMinEta
range for phi for randomization
Double_t Area() const
Definition: AliEmcalJet.h:130
Bool_t fRandomizeInTheta
randomize the particle&#39;s position in pseudorap
double Double_t
Definition: External.C:58
Double_t fRandomPsi4
eventwise calculated psi 3
TRandom3 * fRandom
! random number generator
Int_t fJetRemovalNLeadingJets
threshold at which jets given in fInputJetArray will be removed
TClonesArray * fInputArray
Name of the destination TClonesArray.
TH2D * fDistributionV5
Distribution for v4 in bins of pt and centrality.
Base task in the EMCAL framework.
Bool_t fRandomizeInPt
randomize the particle&#39;s position in theta
Double_t fMinPhi
Artificial tracking efficiency factor.
void GetLeadingJets(AliEmcalJet *&jetLeading, AliEmcalJet *&jetSubLeading)
Double_t fEPV0
!event plane V0
TClonesArray * fJetRemovalArray
Name of the TClonesArray containing jets for removal that will be loaded.
TString part
use mixed event to constrain combinatorial background
Definition: InvMassFit.C:52
TString fOutputArrayName
Name of the TClonesArray that will be loaded.
TClonesArray * fOutputArray
! Destination TClonesArray
Int_t GetLabel() const
Definition: AliPicoTrack.h:41
Double_t Eta() const
Definition: AliPicoTrack.h:37
TString fJetRemovalArrayName
Name of array to rho object.
TString fJetEmbeddingArrayName
if this set via ActivateLeadingJetRemoval, the first n leading jets will be removed ...
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
Int_t ContainsTrack(AliVParticle *track, TClonesArray *tracks) const
Double_t Phi() const
Definition: AliPicoTrack.h:33
Double_t fCent
!event centrality
Double_t fRandomPsi5
eventwise calculated psi 4
Double_t Pt() const
Definition: AliPicoTrack.h:23
Double_t fMinPt
range for eta for randomization
AliEmcalJet * fSubleadingJet
! subleading jet (calculated event-by-event)
TClonesArray * fJetEmbeddingArray
Name of the TClonesArray containing tracks for embedding.
Double_t Pt() const
Definition: AliEmcalJet.h:109
TString fInputArrayName
Distribution for v5 in bins of pt and centrality.
Double_t fMaxPhi
range for phi for randomization
TH2D * fDistributionV2
range for Pt for randomization
Double_t fTrackEfficiency
randomize the particle&#39;s position in Pt
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:51
virtual void ExecOnce()
Perform steps needed to initialize the analysis.
void ExecOnce()
Perform steps needed to initialize the analysis.
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.
AliAODTrack * GetAODTrack(AliPicoTrack *track)
bool Bool_t
Definition: External.C:53
Double_t fMaxEta
range for eta for randomization
AliEmcalJet * fLeadingJet
eventwise calculated psi 5
Bool_t fRandomizeInEta
randomize the particle&#39;s position in azimuth
Short_t Charge() const
Definition: AliPicoTrack.h:39
TH2D * fDistributionV3
Distribution for v2 in bins of pt and centrality.