AliPhysics  fe039ad (fe039ad)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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), 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  new ((*fOutputArray)[accTracks]) AliAODTrack(*((AliAODTrack*)fInputArray->At(iPart)));
129 
130  // Randomize on demand
131  AliAODTrack* particle = static_cast<AliAODTrack*>(fOutputArray->At(accTracks));
132  RandomizeTrack(particle);
133 
134  accTracks++;
135  }
136 
137  // Add particles for embedding (on demand)
139  for(Int_t iPart=0; iPart<fJetEmbeddingArray->GetEntries(); iPart++)
140  {
141  // Take only particles from the randomization acceptance
142  AliPicoTrack* inputParticle = static_cast<AliPicoTrack*>(fJetEmbeddingArray->At(iPart));
143  if(fRandomizeInPhi && (inputParticle->Phi() < fMinPhi || inputParticle->Phi() >= fMaxPhi) )
144  continue;
145  if( (fRandomizeInTheta || fRandomizeInEta) && (inputParticle->Eta() < fMinEta || inputParticle->Eta() >= fMaxEta) )
146  continue;
147 
148  new ((*fOutputArray)[accTracks]) AliAODTrack(*(GetAODTrack(inputParticle)));
149 
150  // Randomize on demand
151  AliAODTrack* particle = static_cast<AliAODTrack*>(fOutputArray->At(accTracks));
152  RandomizeTrack(particle);
153 
154  accTracks++;
155  }
156 
157  // std::cout << Form("%i particles from jets removed out of %i tracks. ", fInputArray->GetEntries()-accTracks, fInputArray->GetEntries()) << std::endl;
158  return kTRUE;
159 }
160 
161 //_____________________________________________________________________________________________________
163 {
164  if(fRandomizeInPhi)
165  particle->SetPhi(fMinPhi + fRandom->Rndm()*(fMaxPhi-fMinPhi));
167  {
168  Double_t minTheta = 2.*atan(exp(-fMinEta));
169  Double_t maxTheta = 2.*atan(exp(-fMaxEta));
170  particle->SetTheta(minTheta + fRandom->Rndm()*(maxTheta-minTheta));
171  }
172  if(fRandomizeInEta)
173  {
174  Double_t randomEta = fMinEta + fRandom->Rndm()*(fMaxEta-fMinEta);
175  Double_t randomTheta = 2.*atan(exp(-randomEta));
176  particle->SetTheta(randomTheta);
177  }
178 
179  if(fRandomizeInPt)
180  particle->SetPt(fMinPt + fRandom->Rndm()*(fMaxPt-fMinPt));
181 
183  particle->SetPhi(AddFlow(particle->Phi(), particle->Pt()));
184 }
185 
186 //_____________________________________________________________________________________________________
188 {
189  AliAODTrack* newTrack = new AliAODTrack();
190  newTrack->SetPt(track->Pt());
191  newTrack->SetTheta(2.*atan(exp(-track->Eta()))); // there is no setter for eta
192  newTrack->SetPhi(track->Phi());
193  newTrack->SetCharge(track->Charge());
194  newTrack->SetLabel(track->GetLabel());
195 
196  // Hybrid tracks (compatible with LHC11h)
197  UInt_t filterMap = BIT(8) | BIT(9);
198  newTrack->SetIsHybridGlobalConstrainedGlobal();
199  newTrack->SetFilterMap(filterMap);
200  return newTrack;
201 }
202 
203 //_____________________________________________________________________________________________________
205 {
206  // Check if particle w/ index 'part' is contained in one of the removal jets
207  // The removal jets can be jets with a pt theshold or the leading jets
208  for(Int_t i=0; i<fJetRemovalArray->GetEntries(); i++)
209  {
210  AliEmcalJet* tmpJet = static_cast<AliEmcalJet*>(fJetRemovalArray->At(i));
211 
212  // Check if to remove leading jets
214  {
215  // leading jet removal mode
216  if( (fJetRemovalNLeadingJets == 1) && (tmpJet == fLeadingJet))
217  {
218  if(tmpJet->ContainsTrack(part)>=0)
219  return kTRUE;
220  else
221  return kFALSE; // we know we are done here: the track is not contained in leading jet
222  }
223  // leading or leading removal mode
224  else if ( (fJetRemovalNLeadingJets == 2) && ((tmpJet == fLeadingJet) || (tmpJet == fSubleadingJet)))
225  {
226  if(tmpJet->ContainsTrack(part)>=0)
227  return kTRUE;
228  }
229  }
230  // Check if to remove jets above threshold
231  else
232  {
233  Double_t tmpPt = tmpJet->Pt() - tmpJet->Area()*GetExternalRho();
234  if(tmpPt >= fJetRemovalPtThreshold)
235  if(tmpJet->ContainsTrack(part)>=0)
236  return kTRUE;
237  }
238  }
239  return kFALSE;
240 }
241 
242 //_____________________________________________________________________________________________________
244 {
245  // Get rho from event.
246  AliRhoParameter* rho = 0;
247  if (!fJetRemovalRhoObj.IsNull())
248  {
249  rho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fJetRemovalRhoObj.Data()));
250  if (!rho) {
251  AliError(Form("%s: Could not retrieve rho with name %s!", GetName(), fJetRemovalRhoObj.Data()));
252  return 0;
253  }
254  }
255  else
256  return 0;
257 
258  return (rho->GetVal());
259 }
260 
261 
262 //_____________________________________________________________________________________________________
264 {
265  // adapted from AliFlowTrackSimple
266  Double_t precisionPhi = 1e-10;
267  Int_t maxNumberOfIterations = 200;
268 
269  Double_t phi0=phi;
270  Double_t f=0.;
271  Double_t fp=0.;
272  Double_t phiprev=0.;
273  Int_t ptBin = 0;
274 
275  // Evaluate V2 for track pt/centrality
276  Double_t v2 = 0;
277  if(fDistributionV2)
278  {
279  ptBin = fDistributionV2->GetXaxis()->FindBin(pt);
280  if(ptBin>fDistributionV2->GetNbinsX())
281  v2 = fDistributionV2->GetBinContent(fDistributionV2->GetNbinsX(), fDistributionV2->GetYaxis()->FindBin(fCent));
282  else if(ptBin>0)
283  v2 = fDistributionV2->GetBinContent(ptBin, fDistributionV2->GetYaxis()->FindBin(fCent));
284  }
285 
286  // Evaluate V3 for track pt/centrality
287  Double_t v3 = 0;
288  if(fDistributionV3)
289  {
290  ptBin = fDistributionV3->GetXaxis()->FindBin(pt);
291  if(ptBin>fDistributionV3->GetNbinsX())
292  v3 = fDistributionV3->GetBinContent(fDistributionV3->GetNbinsX(), fDistributionV3->GetYaxis()->FindBin(fCent));
293  else if(ptBin>0)
294  v3 = fDistributionV3->GetBinContent(ptBin, fDistributionV3->GetYaxis()->FindBin(fCent));
295  }
296 
297  // Evaluate V4 for track pt/centrality
298  Double_t v4 = 0;
299  if(fDistributionV4)
300  {
301  ptBin = fDistributionV4->GetXaxis()->FindBin(pt);
302  if(ptBin>fDistributionV4->GetNbinsX())
303  v4 = fDistributionV4->GetBinContent(fDistributionV4->GetNbinsX(), fDistributionV4->GetYaxis()->FindBin(fCent));
304  else if(ptBin>0)
305  v4 = fDistributionV4->GetBinContent(ptBin, fDistributionV4->GetYaxis()->FindBin(fCent));
306  }
307 
308  // Evaluate V5 for track pt/centrality
309  Double_t v5 = 0;
310  if(fDistributionV5)
311  {
312  ptBin = fDistributionV5->GetXaxis()->FindBin(pt);
313  if(ptBin>fDistributionV5->GetNbinsX())
314  v5 = fDistributionV5->GetBinContent(fDistributionV5->GetNbinsX(), fDistributionV5->GetYaxis()->FindBin(fCent));
315  else if(ptBin>0)
316  v5 = fDistributionV5->GetBinContent(ptBin, fDistributionV5->GetYaxis()->FindBin(fCent));
317  }
318 
319  // Add all v's
320  for (Int_t i=0; i<maxNumberOfIterations; i++)
321  {
322  phiprev=phi; //store last value for comparison
323  f = phi-phi0
324  + v2*TMath::Sin(2.*(phi-(fEPV0+(TMath::Pi()/2.))))
325  +2./3.*v3*TMath::Sin(3.*(phi-fRandomPsi3))
326  +0.5 *v4*TMath::Sin(4.*(phi-fRandomPsi4))
327  +0.4 *v5*TMath::Sin(5.*(phi-fRandomPsi5));
328 
329  fp = 1.0+2.0*(
330  +v2*TMath::Cos(2.*(phi-(fEPV0+(TMath::Pi()/2.))))
331  +v3*TMath::Cos(3.*(phi-fRandomPsi3))
332  +v4*TMath::Cos(4.*(phi-fRandomPsi4))
333  +v5*TMath::Cos(5.*(phi-fRandomPsi5))); //first derivative
334 
335  phi -= f/fp;
336  if (TMath::AreEqualAbs(phiprev,phi,precisionPhi)) break;
337  }
338 
339  return phi;
340 }
341 
342 //_____________________________________________________________________________________________________
344 {
345  // Customized from AliJetContainer::GetLeadingJet()
346  // Get the leading+subleading jet; the sorting is according to pt-A*rho
347 
348  jetLeading = 0;
349  jetSubLeading = 0;
350 
351  Double_t tmpLeadingPt = 0;
352  Double_t tmpSubleadingPt = 0;
353 
354  for(Int_t i=0; i<fJetRemovalArray->GetEntries(); i++)
355  {
356  AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJetRemovalArray->At(i));
357  if ( (jet->Pt()-jet->Area()*GetExternalRho()) > tmpLeadingPt )
358  {
359  jetSubLeading = jetLeading;
360  jetLeading = jet;
361  tmpSubleadingPt = tmpLeadingPt;
362  tmpLeadingPt = jet->Pt()-jet->Area()*GetExternalRho();
363  }
364  else if ( (jet->Pt()-jet->Area()*GetExternalRho()) > tmpSubleadingPt )
365  {
366  jetSubLeading = jet;
367  tmpSubleadingPt = jet->Pt()-jet->Area()*GetExternalRho();
368  }
369  }
370 }
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'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's position in theta
Double_t fMinPhi
randomize the particle's position in Pt
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 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
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's position in azimuth
Short_t Charge() const
Definition: AliPicoTrack.h:39
TH2D * fDistributionV3
Distribution for v2 in bins of pt and centrality.