AliPhysics  6cf2591 (6cf2591)
 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.), fJetEmbeddingArrayName(), fJetEmbeddingArray(0), fRandomPsi3(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  Int_t accTracks = 0;
108 
109  // Add events in input array
110  if(fInputArray)
111  for(Int_t iPart=0; iPart<fInputArray->GetEntries(); iPart++)
112  {
113  // Remove particles contained in jet array (on demand)
114  if(fJetRemovalArray && IsParticleInJet(iPart))
115  continue;
116 
117  // Take only particles from the randomization acceptance
118  AliAODTrack* inputParticle = static_cast<AliAODTrack*>(fInputArray->At(iPart));
119  if(fRandomizeInPhi && (inputParticle->Phi() < fMinPhi || inputParticle->Phi() >= fMaxPhi) )
120  continue;
121  if( (fRandomizeInTheta || fRandomizeInEta) && (inputParticle->Eta() < fMinEta || inputParticle->Eta() >= fMaxEta) )
122  continue;
123 
124  new ((*fOutputArray)[accTracks]) AliAODTrack(*((AliAODTrack*)fInputArray->At(iPart)));
125 
126  // Randomize on demand
127  AliAODTrack* particle = static_cast<AliAODTrack*>(fOutputArray->At(accTracks));
128  RandomizeTrack(particle);
129 
130  accTracks++;
131  }
132 
133  // Add particles for embedding (on demand)
135  for(Int_t iPart=0; iPart<fJetEmbeddingArray->GetEntries(); iPart++)
136  {
137  // Take only particles from the randomization acceptance
138  AliPicoTrack* inputParticle = static_cast<AliPicoTrack*>(fJetEmbeddingArray->At(iPart));
139  if(fRandomizeInPhi && (inputParticle->Phi() < fMinPhi || inputParticle->Phi() >= fMaxPhi) )
140  continue;
141  if( (fRandomizeInTheta || fRandomizeInEta) && (inputParticle->Eta() < fMinEta || inputParticle->Eta() >= fMaxEta) )
142  continue;
143 
144  new ((*fOutputArray)[accTracks]) AliAODTrack(*(GetAODTrack(inputParticle)));
145 
146  // Randomize on demand
147  AliAODTrack* particle = static_cast<AliAODTrack*>(fOutputArray->At(accTracks));
148  RandomizeTrack(particle);
149 
150  accTracks++;
151  }
152 
153  // std::cout << Form("%i particles from jets removed out of %i tracks. ", fInputArray->GetEntries()-accTracks, fInputArray->GetEntries()) << std::endl;
154  return kTRUE;
155 }
156 
157 //_____________________________________________________________________________________________________
159 {
160  if(fRandomizeInPhi)
161  particle->SetPhi(fMinPhi + fRandom->Rndm()*(fMaxPhi-fMinPhi));
163  {
164  Double_t minTheta = 2.*atan(exp(-fMinEta));
165  Double_t maxTheta = 2.*atan(exp(-fMaxEta));
166  particle->SetTheta(minTheta + fRandom->Rndm()*(maxTheta-minTheta));
167  }
168  if(fRandomizeInEta)
169  {
170  Double_t randomEta = fMinEta + fRandom->Rndm()*(fMaxEta-fMinEta);
171  Double_t randomTheta = 2.*atan(exp(-randomEta));
172  particle->SetTheta(randomTheta);
173  }
174 
175  if(fRandomizeInPt)
176  particle->SetPt(fMinPt + fRandom->Rndm()*(fMaxPt-fMinPt));
177 
179  particle->SetPhi(AddFlow(particle->Phi(), particle->Pt()));
180 }
181 
182 //_____________________________________________________________________________________________________
184 {
185  AliAODTrack* newTrack = new AliAODTrack();
186  newTrack->SetPt(track->Pt());
187  newTrack->SetTheta(2.*atan(exp(-track->Eta()))); // there is no setter for eta
188  newTrack->SetPhi(track->Phi());
189  newTrack->SetCharge(track->Charge());
190  newTrack->SetLabel(track->GetLabel());
191 
192  // Hybrid tracks (compatible with LHC11h)
193  UInt_t filterMap = BIT(8) | BIT(9);
194  newTrack->SetIsHybridGlobalConstrainedGlobal();
195  newTrack->SetFilterMap(filterMap);
196  return newTrack;
197 }
198 
199 //_____________________________________________________________________________________________________
201 {
202  for(Int_t i=0; i<fJetRemovalArray->GetEntries(); i++)
203  {
204  AliEmcalJet* tmpJet = static_cast<AliEmcalJet*>(fJetRemovalArray->At(i));
205  Double_t tmpPt = tmpJet->Pt() - tmpJet->Area()*GetExternalRho();
206 
207  if(tmpPt >= fJetRemovalPtThreshold)
208  if(tmpJet->ContainsTrack(part)>=0)
209  return kTRUE;
210  }
211 
212  return kFALSE;
213 }
214 
215 //_____________________________________________________________________________________________________
217 {
218  // Get rho from event.
219  AliRhoParameter* rho = 0;
220  if (!fJetRemovalRhoObj.IsNull())
221  {
222  rho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fJetRemovalRhoObj.Data()));
223  if (!rho) {
224  AliError(Form("%s: Could not retrieve rho with name %s!", GetName(), fJetRemovalRhoObj.Data()));
225  return 0;
226  }
227  }
228  else
229  return 0;
230 
231  return (rho->GetVal());
232 }
233 
234 
235 //_____________________________________________________________________________________________________
237 {
238  // adapted from AliFlowTrackSimple
239  Double_t precisionPhi = 1e-10;
240  Int_t maxNumberOfIterations = 200;
241 
242  Double_t phi0=phi;
243  Double_t f=0.;
244  Double_t fp=0.;
245  Double_t phiprev=0.;
246  Int_t ptBin = 0;
247 
248  // Evaluate V2 for track pt/centrality
249  Double_t v2 = 0;
250  if(fDistributionV2)
251  {
252  ptBin = fDistributionV2->GetXaxis()->FindBin(pt);
253  if(ptBin>fDistributionV2->GetNbinsX())
254  v2 = fDistributionV2->GetBinContent(fDistributionV2->GetNbinsX(), fDistributionV2->GetYaxis()->FindBin(fCent));
255  else if(ptBin>0)
256  v2 = fDistributionV2->GetBinContent(ptBin, fDistributionV2->GetYaxis()->FindBin(fCent));
257  }
258 
259  // Evaluate V3 for track pt/centrality
260  Double_t v3 = 0;
261  if(fDistributionV3)
262  {
263  ptBin = fDistributionV3->GetXaxis()->FindBin(pt);
264  if(ptBin>fDistributionV3->GetNbinsX())
265  v3 = fDistributionV3->GetBinContent(fDistributionV3->GetNbinsX(), fDistributionV3->GetYaxis()->FindBin(fCent));
266  else if(ptBin>0)
267  v3 = fDistributionV3->GetBinContent(ptBin, fDistributionV3->GetYaxis()->FindBin(fCent));
268  }
269 
270  // Evaluate V4 for track pt/centrality
271  Double_t v4 = 0;
272  if(fDistributionV4)
273  {
274  ptBin = fDistributionV4->GetXaxis()->FindBin(pt);
275  if(ptBin>fDistributionV4->GetNbinsX())
276  v4 = fDistributionV4->GetBinContent(fDistributionV4->GetNbinsX(), fDistributionV4->GetYaxis()->FindBin(fCent));
277  else if(ptBin>0)
278  v4 = fDistributionV4->GetBinContent(ptBin, fDistributionV4->GetYaxis()->FindBin(fCent));
279  }
280 
281  // Evaluate V5 for track pt/centrality
282  Double_t v5 = 0;
283  if(fDistributionV5)
284  {
285  ptBin = fDistributionV5->GetXaxis()->FindBin(pt);
286  if(ptBin>fDistributionV5->GetNbinsX())
287  v5 = fDistributionV5->GetBinContent(fDistributionV5->GetNbinsX(), fDistributionV5->GetYaxis()->FindBin(fCent));
288  else if(ptBin>0)
289  v5 = fDistributionV5->GetBinContent(ptBin, fDistributionV5->GetYaxis()->FindBin(fCent));
290  }
291 
292  // Add all v's
293  for (Int_t i=0; i<maxNumberOfIterations; i++)
294  {
295  phiprev=phi; //store last value for comparison
296  f = phi-phi0
297  + v2*TMath::Sin(2.*(phi-(fEPV0+(TMath::Pi()/2.))))
298  +2./3.*v3*TMath::Sin(3.*(phi-fRandomPsi3))
299  +0.5 *v4*TMath::Sin(4.*(phi-fRandomPsi4))
300  +0.4 *v5*TMath::Sin(5.*(phi-fRandomPsi5));
301 
302  fp = 1.0+2.0*(
303  +v2*TMath::Cos(2.*(phi-(fEPV0+(TMath::Pi()/2.))))
304  +v3*TMath::Cos(3.*(phi-fRandomPsi3))
305  +v4*TMath::Cos(4.*(phi-fRandomPsi4))
306  +v5*TMath::Cos(5.*(phi-fRandomPsi5))); //first derivative
307 
308  phi -= f/fp;
309  if (TMath::AreEqualAbs(phiprev,phi,precisionPhi)) break;
310  }
311 
312  return phi;
313 }
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:123
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
eventwise calculated psi 5
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
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:40
Double_t Eta() const
Definition: AliPicoTrack.h:37
TString fJetRemovalArrayName
Name of array to rho object.
TString fJetEmbeddingArrayName
threshold at which jets given in fInputJetArray 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
TClonesArray * fJetEmbeddingArray
Name of the TClonesArray containing tracks for embedding.
Double_t Pt() const
Definition: AliEmcalJet.h:102
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
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
AliAODTrack * GetAODTrack(AliPicoTrack *track)
bool Bool_t
Definition: External.C:53
Double_t fMaxEta
range for eta for randomization
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.