]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/JetTasks/AliAnalyseUE.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Corrected memory leaks and go back to old vertex selection
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalyseUE.cxx
CommitLineData
1f329128 1/*************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: A.Abrahantes, E.Lopez, S.Vallero *
5 * Version 1.0 *
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#include <TROOT.h>
16#include <TBranch.h>
17#include <TCanvas.h>
18#include <TChain.h>
19#include <TFile.h>
20#include <TH1F.h>
21#include <TH1I.h>
22#include <TH2F.h>
23#include <TList.h>
24#include <TLorentzVector.h>
25#include <TMath.h>
26#include <TObjArray.h>
27#include <TProfile.h>
28#include <TRandom.h>
29#include <TSystem.h>
30#include <TTree.h>
31#include <TVector3.h>
32
33#include "AliAnalyseUE.h"
34#include "AliAnalysisTaskUE.h"
35#include "AliAnalysisTask.h"
36#include "AliHistogramsUE.h"
37
38#include "AliAnalysisManager.h"
39#include "AliAODEvent.h"
40#include "AliAODHandler.h"
41#include "AliAODInputHandler.h"
42#include "AliAODJet.h"
43#include "AliAODMCParticle.h"
44#include "AliAODTrack.h"
45#include "AliKFVertex.h"
46#include "AliMCEvent.h"
47#include "AliMCEventHandler.h"
48#include "AliStack.h"
49
50#include "AliAnalysisHelperJetTasks.h"
51#include "AliGenPythiaEventHeader.h"
52#include "AliInputEventHandler.h"
53#include "AliLog.h"
54#include "AliStack.h"
55
56////////////////////////////////////////////////
57//---------------------------------------------
58// Class for transverse regions analysis
59//---------------------------------------------
60////////////////////////////////////////////////
61
62
63using namespace std;
64
65ClassImp(AliAnalyseUE)
66
67//-------------------------------------------------------------------
68AliAnalyseUE::AliAnalyseUE() :
69 TObject(),
70 //fTaskUE(0),
71 fkAOD(0x0),
72 fDebug(0),
73 fSimulateChJetPt(kFALSE),
74 fAnaType(1),
75 fAreaReg(1.5393), // Pi*0.7*0.7
76 fConeRadius(0.7),
77 fFilterBit(0xFF),
78 fRegionType(1),
79 fUseChargeHadrons(kFALSE),
80 fUseChPartJet(kFALSE),
81 fUsePositiveCharge(kTRUE),
82 fUseSingleCharge(kFALSE),
83 fOrdering(1),
84 fJet1EtaCut(0.2),
85 fJet2DeltaPhiCut(2.616), // 150 degrees
86 fJet2RatioPtCut(0.8),
87 fJet3PtCut(15.),
88 fTrackEtaCut(0.9),
89 fTrackPtCut(0.),
90 fHistos(0x0),
91 fSumPtRegionPosit(0.),
92 fSumPtRegionNegat(0.),
93 fSumPtRegionForward(0.),
94 fSumPtRegionBackward(0.),
95 fMaxPartPtRegion(0.),
96 fNTrackRegionPosit(0),
97 fNTrackRegionNegat(0),
98 fNTrackRegionForward(0),
99 fNTrackRegionBackward(0),
100 fSettingsTree(0x0)
101
102{
103 // constructor
104}
105
106
107//-------------------------------------------------------------------
108AliAnalyseUE::AliAnalyseUE(const AliAnalyseUE & original) :
109 TObject(original),
110 //fTaskUE(original.fTaskUE),
111 fkAOD(original.fkAOD),
112 fDebug(original.fDebug),
113 fSimulateChJetPt(original.fSimulateChJetPt),
114 fAnaType(original.fAnaType),
115 fAreaReg(original.fAreaReg),
116 fConeRadius(original.fConeRadius),
117 fFilterBit(original.fFilterBit),
118 fRegionType(original.fRegionType),
119 fUseChargeHadrons(original.fUseChargeHadrons),
120 fUseChPartJet(original.fUseChPartJet),
121 fUsePositiveCharge(original.fUsePositiveCharge),
122 fUseSingleCharge(original.fUseSingleCharge),
123 fOrdering(original.fOrdering),
124 fJet1EtaCut(original.fJet1EtaCut),
125 fJet2DeltaPhiCut(original.fJet2DeltaPhiCut),
126 fJet2RatioPtCut(original.fJet2RatioPtCut),
127 fJet3PtCut(original.fJet3PtCut),
128 fTrackEtaCut(original.fTrackEtaCut),
129 fTrackPtCut(original.fTrackPtCut),
130 fHistos(original.fHistos),
131 fSumPtRegionPosit(original.fSumPtRegionPosit),
132 fSumPtRegionNegat(original.fSumPtRegionNegat),
133 fSumPtRegionForward(original.fSumPtRegionForward),
134 fSumPtRegionBackward(original.fSumPtRegionBackward),
135 fMaxPartPtRegion(original.fMaxPartPtRegion),
136 fNTrackRegionPosit(original.fNTrackRegionPosit),
137 fNTrackRegionNegat(original.fNTrackRegionNegat),
138 fNTrackRegionForward(original.fNTrackRegionForward),
139 fNTrackRegionBackward(original.fNTrackRegionBackward),
140 fSettingsTree(original.fSettingsTree)
141{
142 //copy constructor
143}
144
145//-------------------------------------------------------------------
146AliAnalyseUE & AliAnalyseUE::operator = (const AliAnalyseUE & /*source*/)
147{
148 // assignment operator
149 return *this;
150}
151
152
153//-------------------------------------------------------------------
154AliAnalyseUE::~AliAnalyseUE(){
155
156 //clear memory
157 delete[] fkAOD;
158 fkAOD = NULL;
159
160
161
162}
163
164
165//-------------------------------------------------------------------
166void AliAnalyseUE::AnalyseMC(TVector3 *jetVect,AliMCEvent *mcEvent, AliGenPythiaEventHeader *pythiaGenHeader,Int_t conePosition, Bool_t useAliStack, Bool_t constrainDistance, Double_t minDistance){
167
168 // Execute the analysis in case of MC input
169 fSumPtRegionPosit = 0.;
170 fSumPtRegionNegat = 0.;
171 fSumPtRegionForward = 0.;
172 fSumPtRegionBackward = 0.;
173 fMaxPartPtRegion = 0.;
174 fNTrackRegionPosit = 0;
175 fNTrackRegionNegat = 0;
176 fNTrackRegionForward = 0;
177 fNTrackRegionBackward = 0;
178
179 static Double_t const kPI = TMath::Pi();
180 static Double_t const k270rad = 270.*kPI/180.;
181
182 //Get Jets from MC header
183 Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
184 AliAODJet pythiaGenJets[4];
185 TVector3 jetVectnew[4];
186 Int_t iCount = 0;
187 for(int ip = 0;ip < nPythiaGenJets;++ip){
188 if (iCount>3) break;
189 Float_t p[4];
190 pythiaGenHeader->TriggerJet(ip,p);
191 TVector3 tempVect(p[0],p[1],p[2]);
192 if ( TMath::Abs(tempVect.Eta())>fJet1EtaCut ) continue;
193 pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
194 jetVectnew[iCount].SetXYZ(pythiaGenJets[iCount].Px(), pythiaGenJets[iCount].Py(), pythiaGenJets[iCount].Pz());
195 iCount++;
196 }
197
198 if (!iCount) return;// no jet in eta acceptance
199
200 //Search the index of the nearest MC jet to the leading jet reconstructed from the input data
201 Int_t index = 0;
202 if (constrainDistance){
203 Float_t deltaR = 0.;
204 Float_t dRTemp = 0.;
205 for (Int_t i=0; i<iCount; i++){
206 if (!i) {
207 dRTemp = jetVectnew[i].DeltaR(jetVect[0]);
208 index = i;
209 }
210
211 deltaR = jetVectnew[i].DeltaR(jetVect[0]);
212 if (deltaR < dRTemp){
213 index = i;
214 dRTemp = deltaR;
215 }
216 }
217
218 if (jetVectnew[index].DeltaR(jetVect[0]) > minDistance) return;
219 }
220
221 //Let's add some taste to jet and simulate pt of charged alone
222 //eta and phi are kept as original
223 //Play a Normal Distribution
224 Float_t random = 1.;
225 if (fSimulateChJetPt){
226 while(1){
227 random = gRandom->Gaus(0.6,0.25);
228 if (random > 0. && random < 1. &&
229 (random * jetVectnew[index].Pt()>6.)) break;
230 }
231 }
232
233 //Set new Pt & Fill histogram accordingly
234 Double_t maxPtJet1 = random * jetVectnew[index].Pt();
235
236 fHistos->FillHistogram("hEleadingPt", maxPtJet1 );
237
238 if (useAliStack){//Try Stack Information to perform UE analysis
239
240 AliStack* mcStack = mcEvent->Stack();//Load Stack
241 Int_t nTracksMC = mcStack->GetNtrack();
242 for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
243 //Cuts
244 if(!(mcStack->IsPhysicalPrimary(iTracks))) continue;
245
246 TParticle* mctrk = mcStack->Particle(iTracks);
247
248 Double_t charge = mctrk->GetPDG()->Charge();
249 Double_t pT = mctrk->Pt();
250 Double_t eta = mctrk->Eta();
251 Int_t pdgCode = mctrk->GetPdgCode();
252
253 if (!TrackMCSelected(charge, pT, eta, pdgCode))continue;
254
255 TVector3 partVect(mctrk->Px(), mctrk->Py(), mctrk->Pz());
256 Double_t deltaPhi = jetVectnew[index].DeltaPhi(partVect)+k270rad;
257 if( deltaPhi > 2.*TMath::Pi() ) deltaPhi-= 2.*TMath::Pi();
258 fHistos->FillHistogram("hdNdEtaPhiDist",deltaPhi, maxPtJet1 );
259
260 fHistos->FillHistogram("hFullRegPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
261
262 //We are not interested on stack organization but don't loose track of info
263
264 TVector3 tempVector = jetVectnew[0];
265 jetVectnew[0] = jetVectnew[index];
266 jetVectnew[index] = tempVector;
267
268 Int_t region = IsTrackInsideRegion( jetVectnew, &partVect, conePosition );
269
270 if (region == 1) {
271 if( fMaxPartPtRegion < mctrk->Pt() ) fMaxPartPtRegion = mctrk->Pt();
272 fSumPtRegionPosit += mctrk->Pt();
273 fNTrackRegionPosit++;
274 fHistos->FillHistogram("hTransRegPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
275 }
276 if (region == -1) {
277 if( fMaxPartPtRegion < mctrk->Pt() ) fMaxPartPtRegion = mctrk->Pt();
278 fSumPtRegionNegat += mctrk->Pt();
279 fNTrackRegionNegat++;
280 fHistos->FillHistogram("hTransRegPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
281 }
282 if (region == 2){ //forward
283 fSumPtRegionForward += mctrk->Pt();
284 fNTrackRegionForward++;
285 fHistos->FillHistogram("hRegForwardPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
286 }
287 if (region == -2){ //backward
288 fSumPtRegionBackward += mctrk->Pt();
289 fNTrackRegionBackward++;
290 fHistos->FillHistogram("hRegBackwardPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
291 }
292 } //end loop on stack particles
293 }else{//Try mc Particle
294
295 TClonesArray* farray = (TClonesArray*)fkAOD->FindListObject("mcparticles");
296
297 Int_t ntrks = farray->GetEntries();
298 if (fDebug>1) AliInfo(Form("In UE MC analysis tracks %d \n",ntrks));
299 for (Int_t i =0 ; i < ntrks; i++){
300 AliAODMCParticle* mctrk = (AliAODMCParticle*)farray->At(i);
301 //Cuts
302 if (!(mctrk->IsPhysicalPrimary())) continue;
303 //if (!(mctrk->IsPrimary())) continue;
304
305 Double_t charge = mctrk->Charge();
306 Double_t pT = mctrk->Pt();
307 Double_t eta = mctrk->Eta();
308 Int_t pdgCode = mctrk->GetPdgCode();
309
310 if (!TrackMCSelected(charge, pT, eta, pdgCode))continue;
311
312 TVector3 partVect(mctrk->Px(), mctrk->Py(), mctrk->Pz());
313
314 Double_t deltaPhi = jetVectnew[index].DeltaPhi(partVect)+k270rad;
315 if( deltaPhi > 2.*TMath::Pi() ) deltaPhi-= 2.*TMath::Pi();
316 fHistos->FillHistogram("hdNdEtaPhiDist", deltaPhi, maxPtJet1 );
317
318 fHistos->FillHistogram("hFullRegPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
319
320 //We are not interested on stack organization but don't loose track of info
321 TVector3 tempVector = jetVectnew[0];
322 jetVectnew[0] = jetVectnew[index];
323 jetVectnew[index] = tempVector;
324
325 Int_t region = IsTrackInsideRegion( jetVectnew, &partVect, conePosition );
326
327 if (region == 1) { //right
328 if( fMaxPartPtRegion < mctrk->Pt() ) fMaxPartPtRegion = mctrk->Pt();
329 fSumPtRegionPosit += mctrk->Pt();
330 fNTrackRegionPosit++;
331 fHistos->FillHistogram("hTransRegPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
332 }
333 if (region == -1) { //left
334 if( fMaxPartPtRegion < mctrk->Pt() ) fMaxPartPtRegion = mctrk->Pt();
335 fSumPtRegionNegat += mctrk->Pt();
336 fNTrackRegionNegat++;
337 fHistos->FillHistogram("hTransRegPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
338 }
339 if (region == 2){ //forward
340 fSumPtRegionForward += mctrk->Pt();
341 fNTrackRegionForward++;
342 fHistos->FillHistogram("hRegForwardPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
343 }
344 if (region == -2){ //backward
345 fSumPtRegionBackward += mctrk->Pt();
346 fNTrackRegionBackward++;
347 fHistos->FillHistogram("hRegBackwardPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
348 }
349
350 }//end loop AliAODMCParticle tracks
351 }
352}
353
354
355
356//-------------------------------------------------------------------
357Bool_t AliAnalyseUE::AnaTypeSelection(TVector3 *jetVect ){
358
359 // Cut events by jets topology
360 // anaType:
361 // 1 = inclusive,
362 // - Jet1 |eta| < jet1EtaCut
363 // 2 = back to back inclusive
364 // - fulfill case 1
365 // - |Jet1.Phi - Jet2.Phi| > jet2DeltaPhiCut
366 // - Jet2.Pt/Jet1Pt > jet2RatioPtCut
367 // 3 = back to back exclusive
368 // - fulfill case 2
369 // - Jet3.Pt < jet3PtCut
370
371 Double_t eta=jetVect[0].Eta();
372 if( TMath::Abs(eta) > fJet1EtaCut) {
373 if( fDebug > 1 ) AliInfo("\n Skipping Event...Jet1 |eta| > fJet1EtaCut");
374 return kFALSE;
375 }
376 // back to back inclusive
377 if( fAnaType > 1 && fAnaType < 4 && jetVect[1].Pt() < 0. ) {
378 if( fDebug > 1 ) AliInfo("\n Skipping Event... no second Jet found");
379 return kFALSE;
380 }
381 if( fAnaType > 1 && fAnaType < 4 && jetVect[1].Pt() > 0. ) {
382 if( TMath::Abs(jetVect[0].DeltaPhi(jetVect[1])) < fJet2DeltaPhiCut ||
383 jetVect[1].Pt()/jetVect[0].Pt() < fJet2RatioPtCut ) {
384 if( fDebug > 1 ) AliInfo("\n Skipping Event... |Jet1.Phi - Jet2.Phi| < fJet2DeltaPhiCut");
385 return kFALSE;
386 }
387 }
388 // back to back exclusive
389 if( fAnaType > 2 && fAnaType < 4 && jetVect[2].Pt() > 0. ) {
390 if( jetVect[2].Pt() > fJet3PtCut ) {
391 if( fDebug > 1 ) AliInfo("\n Skipping Event... Jet3.Pt > fJet3PtCut ");
392 return kFALSE;
393 }
394 }
395 return kTRUE;
396}
397
398
1f329128 399//-------------------------------------------------------------------
400void AliAnalyseUE::FillRegions(Bool_t isNorm2Area, TVector3 *jetVect){
401
402 // Fill the different topological regions
403 Double_t maxPtJet1 = jetVect[0].Pt();
404 static Double_t const kPI = TMath::Pi();
405 static Double_t const k120rad = 120.*kPI/180.;
406 Double_t const kMyTolerance = 0.0000001;
407
408 //Area for Normalization
409 // Forward and backward
410 Double_t normArea = 1.;
411 // Transverse
412 if (isNorm2Area) {
413 SetRegionArea(jetVect);
414 normArea = 2.*fTrackEtaCut*k120rad ;
415 } else fAreaReg = 1.;
416
417 Double_t avePosRegion = (fNTrackRegionPosit) ? fSumPtRegionPosit/fNTrackRegionPosit : 0.;
418 Double_t aveNegRegion = (fNTrackRegionNegat) ? fSumPtRegionNegat/fNTrackRegionNegat : 0.;
419 if( avePosRegion > aveNegRegion ) {
420 FillAvePartPtRegion( maxPtJet1, avePosRegion/fAreaReg, aveNegRegion/fAreaReg );
421 } else {
422 FillAvePartPtRegion( maxPtJet1, aveNegRegion/fAreaReg, avePosRegion/fAreaReg );
423 }
424
425 //How quantities will be sorted before Fill Min and Max Histogram
426 // 1=Plots will be CDF-like
427 // 2=Plots will be Marchesini-like
428 // 3=Minimum zone is selected as the one having lowest pt per track
429 if( fOrdering == 1 ) {
430 if( fSumPtRegionPosit > fSumPtRegionNegat ) {
431 FillSumPtRegion( maxPtJet1, fSumPtRegionPosit/fAreaReg, fSumPtRegionNegat/fAreaReg );
432 } else {
433 FillSumPtRegion( maxPtJet1, fSumPtRegionNegat/fAreaReg, fSumPtRegionPosit/fAreaReg );
434 }
435 if (fNTrackRegionPosit > fNTrackRegionNegat ) {
436 FillMultRegion( maxPtJet1, fNTrackRegionPosit/fAreaReg, fNTrackRegionNegat/fAreaReg, fSumPtRegionNegat/fAreaReg );
437 } else {
438 FillMultRegion( maxPtJet1, fNTrackRegionNegat/fAreaReg, fNTrackRegionPosit/fAreaReg, fSumPtRegionPosit/fAreaReg );
439 }
440 } else if( fOrdering == 2 ) {
441 if (fSumPtRegionPosit > fSumPtRegionNegat) {
442 FillSumPtRegion( maxPtJet1, fSumPtRegionPosit/fAreaReg, fSumPtRegionNegat/fAreaReg );
443 FillMultRegion( maxPtJet1, fNTrackRegionPosit/fAreaReg, fNTrackRegionNegat/fAreaReg, fSumPtRegionNegat/fAreaReg );
444 } else {
445 FillSumPtRegion( maxPtJet1, fSumPtRegionNegat/fAreaReg, fSumPtRegionPosit/fAreaReg );
446 FillMultRegion( maxPtJet1, fNTrackRegionNegat/fAreaReg, fNTrackRegionPosit/fAreaReg, fSumPtRegionPosit/fAreaReg );
447 }
448 } else if( fOrdering == 3 ){
449 if (avePosRegion > aveNegRegion) {
450 FillSumPtRegion( maxPtJet1, fSumPtRegionPosit/fAreaReg, fSumPtRegionNegat/fAreaReg );
451 FillMultRegion( maxPtJet1, fNTrackRegionPosit/fAreaReg, fNTrackRegionNegat/fAreaReg, fSumPtRegionNegat/fAreaReg );
452 }else{
453 FillSumPtRegion( maxPtJet1, fSumPtRegionNegat/fAreaReg, fSumPtRegionPosit/fAreaReg );
454 FillMultRegion( maxPtJet1, fNTrackRegionNegat/fAreaReg, fNTrackRegionPosit/fAreaReg, fSumPtRegionPosit/fAreaReg );
455 }
456 }
457 fHistos->FillHistogram("hRegionMaxPartPtMaxVsEt",maxPtJet1, fMaxPartPtRegion);
458
459 // Compute pedestal like magnitudes
460 fHistos->FillHistogram("hRegionDiffSumPtVsEt",maxPtJet1, (TMath::Abs(fSumPtRegionPosit-fSumPtRegionNegat)/(2.0*fAreaReg))+kMyTolerance);
461 fHistos->FillHistogram("hRegionAveSumPtVsEt", maxPtJet1, (fSumPtRegionPosit+fSumPtRegionNegat)/(2.0*fAreaReg));
462
463 // Transverse as a whole
464 fHistos->FillHistogram("hRegTransMult", maxPtJet1, fNTrackRegionPosit + fNTrackRegionNegat, (fNTrackRegionPosit + fNTrackRegionNegat)/(2.0*fAreaReg));
465 fHistos->FillHistogram("hRegTransSumPtVsMult",maxPtJet1, fNTrackRegionPosit + fNTrackRegionNegat , (fSumPtRegionNegat + fSumPtRegionPosit)/(2.0 *fAreaReg));
466
467 // Fill Histograms for Forward and away side w.r.t. leading jet direction
468 // Pt dependence
469 //fHistos->FillHistogram("hRegForwardSumPtVsEt",maxPtJet1, fSumPtRegionForward/normArea );
470 //fHistos->FillHistogram("hRegForwardMultVsEt",maxPtJet1, fNTrackRegionForward/normArea );
471 //fHistos->FillHistogram("hRegBackwardSumPtVsEt",maxPtJet1, fSumPtRegionBackward/normArea );
472 //fHistos->FillHistogram("hRegBackwardMultVsEt",maxPtJet1, fNTrackRegionBackward/normArea);
473
474 // Multiplicity dependence
475 fHistos->FillHistogram("hRegForwardMult", maxPtJet1, fNTrackRegionForward, fNTrackRegionForward/normArea);
476 fHistos->FillHistogram("hRegForwardSumPtvsMult", maxPtJet1, fNTrackRegionForward,fSumPtRegionForward/normArea);
477 fHistos->FillHistogram("hRegBackwardMult", maxPtJet1, fNTrackRegionBackward, fNTrackRegionBackward/normArea );
478 fHistos->FillHistogram("hRegBackwardSumPtvsMult", maxPtJet1, fNTrackRegionBackward,fSumPtRegionBackward/normArea);
479}
480
1f329128 481
482//-------------------------------------------------------------------
483void AliAnalyseUE::FindMaxMinRegions(TVector3 *jetVect, Int_t conePosition){
484
485 // Identify the different topological zones
486 fSumPtRegionPosit = 0.;
487 fSumPtRegionNegat = 0.;
488 fSumPtRegionForward = 0.;
489 fSumPtRegionBackward = 0.;
490 fMaxPartPtRegion = 0.;
491 fNTrackRegionPosit = 0;
492 fNTrackRegionNegat = 0;
493 fNTrackRegionForward = 0;
494 fNTrackRegionBackward = 0;
495 static Double_t const kPI = TMath::Pi();
496 static Double_t const kTWOPI = 2.*kPI;
497 static Double_t const k270rad = 270.*kPI/180.;
498 Double_t const kMyTolerance = 0.0000001;
499
500 Int_t nTracks = fkAOD->GetNTracks();
501 if (fDebug > 1) AliInfo(Form(" ==== AOD tracks = %d \n ",nTracks));
502
503 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
504
505 AliAODTrack* part = fkAOD->GetTrack( ipart );
506 if (fDebug > 1) AliInfo(Form(" ==== AOD track = %d pt = %f charge = %d \n ",ipart,part->Pt(),part->Charge()));
507 // track selection
508 if (! TrackSelected(part)) continue;
509
510
511 TVector3 partVect(part->Px(), part->Py(), part->Pz());
512 Bool_t isFlagPart = kTRUE;
513 Double_t deltaPhi = jetVect[0].DeltaPhi(partVect)+k270rad;
514 if( deltaPhi > kTWOPI ) deltaPhi-= kTWOPI;
515 if (fAnaType != 4 ) fHistos->FillHistogram("hdNdEtaPhiDist",deltaPhi, jetVect[0].Pt());
516 else if (TMath::Abs(deltaPhi-k270rad) >= kMyTolerance && TMath::Abs(jetVect[0].Eta()-partVect.Eta()) >= kMyTolerance){
517 fHistos->FillHistogram("hdNdEtaPhiDist",deltaPhi, jetVect[0].Pt());
518 isFlagPart = kFALSE;
519 }
520
521 fHistos->FillHistogram("hFullRegPartPtDistVsEt", part->Pt(), jetVect[0].Pt());
522
523 Int_t region = IsTrackInsideRegion( jetVect, &partVect, conePosition );
524 if (region == 1) {
525 if( fMaxPartPtRegion < part->Pt() ) fMaxPartPtRegion = part->Pt();
526 fSumPtRegionPosit += part->Pt();
527 fNTrackRegionPosit++;
528 fHistos->FillHistogram("hTransRegPartPtDistVsEt",part->Pt(), jetVect[0].Pt());
529 }
530 if (region == -1) {
531 if( fMaxPartPtRegion < part->Pt() ) fMaxPartPtRegion = part->Pt();
532 fSumPtRegionNegat += part->Pt();
533 fNTrackRegionNegat++;
534 fHistos->FillHistogram("hTransRegPartPtDistVsEt",part->Pt(), jetVect[0].Pt());
535 }
536 if (region == 2){ //forward
537 fSumPtRegionForward += part->Pt();
538 fNTrackRegionForward++;
539 fHistos->FillHistogram("hRegForwardPartPtDistVsEt",part->Pt(), jetVect[0].Pt());
540 }
541 if (region == -2){ //backward
542 fSumPtRegionBackward += part->Pt();
543 fNTrackRegionBackward++;
544 fHistos->FillHistogram("hRegBackwardPartPtDistVsEt",part->Pt(), jetVect[0].Pt());
545 }
546 }//end loop AOD tracks
547
548}
549
1f329128 550//-------------------------------------------------------------------
551TVector3 AliAnalyseUE::GetOrderedClusters(TString aodBranch, Bool_t chargedJets, Double_t chJetPtMin){
552
553 // jets from AOD, on-the-fly or leading particle
554 Double_t maxPtJet1 = 0.;
555 Int_t index1 = -1;
556 Double_t maxPtJet2 = 0.; // jet 2 need for back to back inclusive
557 Int_t index2 = -1;
558 Double_t maxPtJet3 = 0.; // jet 3 need for back to back exclusive
559 Int_t index3 = -1;
560 TVector3 jetVect[3];
561
562 jetVect[0].SetPtEtaPhi(-1.,-1.,-1.);
563 jetVect[1].SetPtEtaPhi(-1.,-1.,-1.);
564 jetVect[2].SetPtEtaPhi(-1.,-1.,-1.);
565
566 Int_t nJets = 0;
567 //TClonesArray* fArrayJets;
568 TObjArray* arrayJets;
569 // 1) JETS FROM AOD BRANCH (standard, non-standard or delta)
570 if (!chargedJets && fAnaType != 4 ) {
571 AliInfo(" ==== Read AODs !");
572 AliInfo(Form(" ==== Reading Branch: %s ", aodBranch.Data()));
573 arrayJets = (TObjArray*)fkAOD->GetList()->FindObject(aodBranch.Data());
574 if (!arrayJets){
575 AliFatal(" No jet-array! ");
576 return *jetVect;
577 }
578
579 // Find Leading Jets 1,2,3
580 // (could be skipped if Jets are sort by Pt...)
581 nJets=arrayJets->GetEntries();
582 for( Int_t i=0; i<nJets; ++i ) {
583 AliAODJet* jet = (AliAODJet*)arrayJets->At(i);
584 Double_t jetPt = jet->Pt();//*1.666; // FIXME Jet Pt Correction ?????!!!
585
586 if( jetPt > maxPtJet1 ) {
587 maxPtJet3 = maxPtJet2; index3 = index2;
588 maxPtJet2 = maxPtJet1; index2 = index1;
589 maxPtJet1 = jetPt; index1 = i;
590 } else if( jetPt > maxPtJet2 ) {
591 maxPtJet3 = maxPtJet2; index3 = index2;
592 maxPtJet2 = jetPt; index2 = i;
593 } else if( jetPt > maxPtJet3 ) {
594 maxPtJet3 = jetPt; index3 = i;
595 }
596 }
597
598 if( index1 != -1 ) {
599 AliAODJet *jet =(AliAODJet*) arrayJets->At(index1);
600 if(jet)jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
601 }
602 if( index2 != -1 ) {
603 AliAODJet* jet= (AliAODJet*) arrayJets->At(index2);
604 if(jet)jetVect[1].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
605 }
606 if( index3 != -1 ) {
607 AliAODJet* jet = (AliAODJet*) arrayJets->At(index3);
608 if(jet)jetVect[2].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
609 }
610
611 }
612
613
614 // 2) ON-THE-FLY CDF ALGORITHM
615 if (chargedJets){
1f329128 616 arrayJets = FindChargedParticleJets(chJetPtMin);
617 if( arrayJets ) {
618 nJets = arrayJets->GetEntriesFast();
619 if( nJets > 0 ) {
620 index1 = 0;
621 AliAODJet* jet = (AliAODJet*)arrayJets->At(0);
622 maxPtJet1 = jet->Pt();
623 jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
624 }
625 if( nJets > 1 ) {
626 index2 = 1;
627 AliAODJet* jet = (AliAODJet*)arrayJets->At(1);
628 maxPtJet2 = jet->Pt();
629 jetVect[1].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
630 }
631 if( nJets > 2 ) {
632 index3 = 2;
633 AliAODJet* jet = (AliAODJet*)arrayJets->At(2);
634 maxPtJet3 = jet->Pt();
635 jetVect[2].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
636 }
637
638 arrayJets->Delete();
639 delete arrayJets;
640 }
641 }
642
643
644 // 3) LEADING PARTICLE
645 if( fAnaType == 4 ){
646 TObjArray* tracks = SortChargedParticles();
647 if( tracks ) {
648 nJets = tracks->GetEntriesFast();
649 if( nJets > 0 ) {
650 index1 = 0;
651 AliAODTrack* jet = (AliAODTrack*)tracks->At(0);
652 maxPtJet1 = jet->Pt();
653 jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
654 }
655 tracks->Clear();
656 delete tracks;
657 }
658
659 }
ddbad1f5 660 fHistos->FillHistogram("hNJets",nJets);
1f329128 661
662 return *jetVect;
663
664}
665
666
667//-------------------------------------------------------------------
668void AliAnalyseUE::Initialize(AliAnalysisTaskUE& taskUE){
669
670 //Get principal settings from current instance of UE analysis-task
671 fAnaType = taskUE.GetAnaTopology();
672 fkAOD = taskUE.GetAOD();
673 fConeRadius = taskUE.GetConeRadius();
674 fDebug = taskUE.GetDebugLevel();
675 fFilterBit = taskUE.GetFilterBit();
676 fJet1EtaCut = taskUE.GetJet1EtaCut();
677 fJet2DeltaPhiCut = taskUE.GetJet2DeltaPhiCut();
678 fJet2RatioPtCut = taskUE.GetJet2RatioPtCut();
679 fJet3PtCut = taskUE.GetJet3PtCut();
680 fOrdering = taskUE.GetPtSumOrdering() ;
681 fRegionType = taskUE.GetRegionType();
682 fSimulateChJetPt = taskUE.GetSimulateChJetPt();
683 fTrackEtaCut = taskUE.GetTrackEtaCut();
684 fTrackPtCut = taskUE.GetTrackPtCut();
ddbad1f5 685 fHistos = taskUE.fHistosUE;
1f329128 686 fUseChargeHadrons = taskUE.GetUseChargeHadrons();
687 fUseChPartJet = taskUE.GetUseChPartJet();
688 fUsePositiveCharge = taskUE.GetUseNegativeChargeType();
689 fUseSingleCharge = taskUE.GetUseSingleCharge();
690
1f329128 691}
692
693//-------------------------------------------------------------------
1f329128 694Bool_t AliAnalyseUE::VertexSelection(AliAODEvent *aod, Int_t tracks, Double_t zed ){
695
696 //Require 1 vertex (no TPC stand-alone) with a minimum number of tracks and z-coordinate in a limited range
697 Int_t nVertex = aod->GetNumberOfVertices();
698 if( nVertex > 0 ) { // Only one vertex (reject pileup)
699 AliAODVertex* vertex = (AliAODVertex*)aod->GetPrimaryVertex();
700 Int_t nTracksPrim = vertex->GetNContributors();
701 Double_t zVertex = vertex->GetZ();
702 if (fDebug > 1) AliInfo(Form(" Vertex in = %f with %d particles by %s data ...",zVertex,nTracksPrim,vertex->GetName()));
703 // Select a quality vertex by number of tracks?
704 if( nTracksPrim < tracks || TMath::Abs(zVertex) > zed ) {
705 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
706 return kFALSE;
707 }
708 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
709 } else {
710 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
711 return kFALSE;
712 }
713
714 return kTRUE;
715}
716
ddbad1f5 717//-------------------------------------------------------------------
718Bool_t AliAnalyseUE::VertexSelectionOld(AliAODEvent *aod ){
719
720 AliKFVertex primVtx(*(aod->GetPrimaryVertex()));
721 Int_t nTracksPrim=primVtx.GetNContributors();
722 if (fDebug > 1) AliInfo(Form(" Primary-vertex Selection: %d",nTracksPrim));
723 if(!nTracksPrim){
724 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
725 return kFALSE;
726 }
727 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED ...");
728
729 return kTRUE;
730}
731
1f329128 732// PRIVATE METHODS **************************************************
733
734TObjArray* AliAnalyseUE::FindChargedParticleJets( Double_t chJetPtMin )
735{
736 // Return a TObjArray of "charged particle jets"
737
738 // Charged particle jet definition from reference:
739 // "Charged jet evolution and the underlying event
740 // in proton-antiproton collisions at 1.8 TeV"
741 // PHYSICAL REVIEW D 65 092002, CDF Collaboration
742
743 // We defined "jets" as circular regions in eta-phi space with
744 // radius defined by R = sqrt( (eta-eta0)^2 +(phi-phi0)^2 ).
745 // Our jet algorithm is as follows:
746 // 1- Order all charged particles according to their pT .
747 // 2- Start with the highest pT particle and include in the jet all
748 // particles within the radius R=0.7 considering each particle
749 // in the order of decreasing pT and recalculating the centroid
750 // of the jet after each new particle is added to the jet .
751 // 3- Go to the next highest pT particle not already included in
752 // a jet and add to the jet all particles not already included in
753 // a jet within R=0.7.
754 // 4- Continue until all particles are in a jet.
755 // We defined the transverse momentum of the jet to be
756 // the scalar pT sum of all the particles within the jet, where pT
757 // is measured with respect to the beam axis
758
759 // 1 - Order all charged particles according to their pT .
760 Int_t nTracks = fkAOD->GetNTracks();
761 if( !nTracks ) return 0;
762 TObjArray tracks(nTracks);
763
764 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
765 AliAODTrack* part = fkAOD->GetTrack( ipart );
766 if( !part->TestFilterBit(fFilterBit) ) continue; // track cut selection
767 if( !part->Charge() ) continue;
768 if( part->Pt() < fTrackPtCut ) continue;
769 tracks.AddLast(part);
770 }
771 QSortTracks( tracks, 0, tracks.GetEntriesFast() );
772
773 nTracks = tracks.GetEntriesFast();
774 if( !nTracks ) return 0;
775
776 TObjArray *jets = new TObjArray(nTracks);
777 TIter itrack(&tracks);
778 while( nTracks ) {
779 // 2- Start with the highest pT particle ...
780 Float_t px,py,pz,pt;
781 AliAODTrack* track = (AliAODTrack*)itrack.Next();
782 if( !track ) continue;
783 px = track->Px();
784 py = track->Py();
785 pz = track->Pz();
786 pt = track->Pt(); // Use the energy member to store Pt
787 jets->AddLast( new TLorentzVector(px, py, pz, pt) );
788 tracks.Remove( track );
789 TLorentzVector* jet = (TLorentzVector*)jets->Last();
790 jet->SetPtEtaPhiE( 1., jet->Eta(), jet->Phi(), pt );
791 // 3- Go to the next highest pT particle not already included...
792 AliAODTrack* track1;
793 Double_t fPt = jet->E();
794 while ( (track1 = (AliAODTrack*)(itrack.Next())) ) {
795 Double_t tphi = track1->Phi(); // here Phi is from 0 <-> 2Pi
796 if (tphi > TMath::Pi()) tphi -= 2. * TMath::Pi(); // convert to -Pi <-> Pi
797 Double_t dphi = TVector2::Phi_mpi_pi(jet->Phi()-tphi);
798 Double_t r = TMath::Sqrt( (jet->Eta()-track1->Eta())*(jet->Eta()-track1->Eta()) +dphi*dphi );
799 if( r < fConeRadius ) {
800 fPt = jet->E()+track1->Pt(); // Scalar sum of Pt
801 // recalculating the centroid
802 Double_t eta = jet->Eta()*jet->E()/fPt + track1->Eta()*track1->Pt()/fPt;
803 Double_t phi = jet->Phi()*jet->E()/fPt + tphi*track1->Pt()/fPt;
804 jet->SetPtEtaPhiE( 1., eta, phi, fPt );
805 tracks.Remove( track1 );
806 }
807 }
808
809 tracks.Compress();
810 nTracks = tracks.GetEntries();
811 // 4- Continue until all particles are in a jet.
812 itrack.Reset();
813 } // end while nTracks
814
815 // Convert to AODjets....
816 Int_t njets = jets->GetEntriesFast();
817 TObjArray* aodjets = new TObjArray(njets);
818 aodjets->SetOwner(kTRUE);
819 for(Int_t ijet=0; ijet<njets; ++ijet) {
820 TLorentzVector* jet = (TLorentzVector*)jets->At(ijet);
821 if (jet->E() < chJetPtMin) continue;
822 Float_t px, py,pz,en; // convert to 4-vector
823 px = jet->E() * TMath::Cos(jet->Phi()); // Pt * cos(phi)
824 py = jet->E() * TMath::Sin(jet->Phi()); // Pt * sin(phi)
825 pz = jet->E() / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-jet->Eta())));
826 en = TMath::Sqrt(px * px + py * py + pz * pz);
827
828 aodjets->AddLast( new AliAODJet(px, py, pz, en) );
829 }
830 jets->Delete();
831 delete jets;
832
833 // Order jets according to their pT .
834 QSortTracks( *aodjets, 0, aodjets->GetEntriesFast() );
835
836 // debug
837 if (fDebug>3) AliInfo(Form(" %d Charged jets found\n",njets));
838
839 return aodjets;
840}
841
842
843//____________________________________________________________________
844void AliAnalyseUE::FillAvePartPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin )
845{
846
847 // Fill average particle Pt of control regions
848
849 // Max cone
850 fHistos->FillHistogram("hRegionAvePartPtMaxVsEt", leadingE, ptMax);
851 // Min cone
852 fHistos->FillHistogram("hRegionAvePartPtMinVsEt", leadingE, ptMin);
853 // MAke distributions for UE comparison with MB data
854 fHistos->FillHistogram("hMinRegAvePt", ptMin);
855
856}
857
858//____________________________________________________________________
859void AliAnalyseUE::FillMultRegion( Double_t leadingE, Double_t nTrackPtmax, Double_t nTrackPtmin, Double_t ptMin )
860{
861
862 // Fill Nch multiplicity of control regions
863
864 // Max cone
865 fHistos->FillHistogram("hRegionMultMaxVsEt", leadingE, nTrackPtmax);
866 fHistos->FillHistogram("hRegionMultMax", nTrackPtmax);
867
868 // Min cone
869 fHistos->FillHistogram("hRegionMultMinVsEt", leadingE, nTrackPtmin );
870 fHistos->FillHistogram("hRegionMultMin", nTrackPtmin);
871
872 // MAke distributions for UE comparison with MB data
873 fHistos->FillHistogram("hMinRegSumPtvsMult", nTrackPtmin,ptMin);
874
875}
876
877//____________________________________________________________________
878void AliAnalyseUE::FillSumPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin )
879{
880 // Fill sumPt of control regions
881
882 // Max cone
883 fHistos->FillHistogram("hRegionSumPtMaxVsEt", leadingE, ptMax);
884
885 // Min cone
886 fHistos->FillHistogram("hRegionSumPtMinVsEt", leadingE, ptMin);
887
888 // MAke distributions for UE comparison with MB data
889 fHistos->FillHistogram("hMinRegSumPt", ptMin);
890 fHistos->FillHistogram("hMinRegSumPtJetPtBin", leadingE, ptMin);
891 fHistos->FillHistogram("hMaxRegSumPtJetPtBin", leadingE, ptMax);
892
893}
894
895//-------------------------------------------------------------------
896Int_t AliAnalyseUE::IsTrackInsideRegion(TVector3 *jetVect, TVector3 *partVect, Int_t conePosition)
897{
898 // return de region in delta phi
899 // -1 negative delta phi
900 // 1 positive delta phi
901 // 0 outside region
902 static const Double_t k60rad = 60.*TMath::Pi()/180.;
903 static const Double_t k120rad = 120.*TMath::Pi()/180.;
904
905 Int_t region = 0;
906 if( fRegionType == 1 ) {
907 if( TMath::Abs(partVect->Eta()) > fTrackEtaCut ) return 0;
908 // transverse regions
909 if (jetVect[0].DeltaPhi(*partVect) < -k60rad && jetVect[0].DeltaPhi(*partVect) > -k120rad ) region = -1; //left
910 if (jetVect[0].DeltaPhi(*partVect) > k60rad && jetVect[0].DeltaPhi(*partVect) < k120rad ) region = 1; //right
911 if (TMath::Abs(jetVect[0].DeltaPhi(*partVect)) < k60rad ) region = 2; //forward
912 if (TMath::Abs(jetVect[0].DeltaPhi(*partVect)) > k120rad ) region = -2; //backward
913
914 } else if( fRegionType == 2 ) {
915 // Cone regions
916 Double_t deltaR = 0.;
917
918 TVector3 positVect,negatVect;
919 if (conePosition==1){
920 positVect.SetMagThetaPhi(1, 2.*atan(exp(-jetVect[0].Eta())), jetVect[0].Phi()+TMath::PiOver2());
921 negatVect.SetMagThetaPhi(1, 2.*atan(exp(-jetVect[0].Eta())), jetVect[0].Phi()-TMath::PiOver2());
922 }else if (conePosition==2){
923 if(fAnaType<2) AliFatal("Prevent error in Analysis type there might be only 1 jet. To avoid overflow better Correct UE config");
924 positVect.SetMagThetaPhi(1, 2.*atan(exp(-(jetVect[0].Eta()+jetVect[1].Eta())/2.)), jetVect[0].Phi()+TMath::PiOver2());
925 negatVect.SetMagThetaPhi(1, 2.*atan(exp(-(jetVect[0].Eta()+jetVect[1].Eta())/2.)), jetVect[0].Phi()-TMath::PiOver2());
926 }else if (conePosition==3){
927 if(fAnaType<2) AliFatal("Prevent error in Analysis type there might be only 1 jet. To avoid overflow better Correct UE config");
928 Double_t weightEta = jetVect[0].Eta() * jetVect[0].Pt()/(jetVect[0].Pt() + jetVect[1].Pt()) +
929 jetVect[1].Eta() * jetVect[1].Pt()/(jetVect[0].Pt() + jetVect[1].Pt());
930 //Double_t weightEta = jetVect[0].Eta() * jetVect[0].Mag()/(jetVect[0].Mag() + jetVect[1].Mag()) +
931 // jetVect[1].Eta() * jetVect[1].Mag()/(jetVect[0].Mag() + jetVect[1].Mag());
932 positVect.SetMagThetaPhi(1, 2.*atan(exp(-weightEta)), jetVect[0].Phi()+TMath::PiOver2());
933 negatVect.SetMagThetaPhi(1, 2.*atan(exp(-weightEta)), jetVect[0].Phi()-TMath::PiOver2());
934 }
935 if (TMath::Abs(positVect.DeltaPhi(*partVect)) < fConeRadius ) {
936 region = 1;
937 deltaR = positVect.DrEtaPhi(*partVect);
938 } else if (TMath::Abs(negatVect.DeltaPhi(*partVect)) < fConeRadius) {
939 region = -1;
940 deltaR = negatVect.DrEtaPhi(*partVect);
941 }
942
943 if (deltaR > fConeRadius) region = 0;
944
945 } else
946 AliError("Unknow region type");
947
948 return region;
949 }
950
951
952
953//-------------------------------------------------------------------
954void AliAnalyseUE::QSortTracks(TObjArray &a, Int_t first, Int_t last)
955{
956 // Sort array of TObjArray of tracks by Pt using a quicksort algorithm.
957
958 static TObject *tmp;
959 static int i; // "static" to save stack space
960 int j;
961
962 while (last - first > 1) {
963 i = first;
964 j = last;
965 for (;;) {
966 while (++i < last && ((AliVParticle*)a[i])->Pt() > ((AliVParticle*)a[first])->Pt() )
967 ;
968 while (--j > first && ((AliVParticle*)a[j])->Pt() < ((AliVParticle*)a[first])->Pt() )
969 ;
970 if (i >= j)
971 break;
972
973 tmp = a[i];
974 a[i] = a[j];
975 a[j] = tmp;
976 }
977 if (j == first) {
978 ++first;
979 continue;
980 }
981 tmp = a[first];
982 a[first] = a[j];
983 a[j] = tmp;
984 if (j - first < last - (j + 1)) {
985 QSortTracks(a, first, j);
986 first = j + 1; // QSortTracks(j + 1, last);
987 } else {
988 QSortTracks(a, j + 1, last);
989 last = j; // QSortTracks(first, j);
990 }
991 }
992}
993
994
995
996
997//-------------------------------------------------------------------
998void AliAnalyseUE::SetRegionArea(TVector3 *jetVect)
999{
1000 // Set region area
1001 Double_t areaCorrFactor=0.;
1002 Double_t deltaEta = 0.;
1003 if (fRegionType==1) fAreaReg = 2.*fTrackEtaCut*60.*TMath::Pi()/180.;
1004 else if (fRegionType==2){
1005 deltaEta = 0.9-TMath::Abs(jetVect[0].Eta());
1006 if (deltaEta>fConeRadius) fAreaReg = TMath::Pi()*fConeRadius*fConeRadius;
1007 else{
1008 areaCorrFactor = fConeRadius*fConeRadius*TMath::ACos(deltaEta/fConeRadius) -
1009 (deltaEta*fConeRadius)*TMath::Sqrt( 1. - deltaEta*deltaEta/(fConeRadius*fConeRadius));
1010 fAreaReg=TMath::Pi()*fConeRadius*fConeRadius-areaCorrFactor;
1011 }
1012 }else AliWarning("Unknown Region Type");
1013 if (fDebug>10) AliInfo(Form("\n dEta=%5.3f Angle =%5.3f Region Area = %5.3f Corr Factor=%5.4f \n",deltaEta,TMath::ACos(deltaEta/fConeRadius),fAreaReg,areaCorrFactor));
1014}
1015
1016
1017//____________________________________________________________________
1018TObjArray* AliAnalyseUE::SortChargedParticles()
1019{
1020 // return an array with all charged particles ordered according to their pT .
1021 Int_t nTracks = fkAOD->GetNTracks();
1022 if( !nTracks ) return 0;
1023 TObjArray* tracks = new TObjArray(nTracks);
1024
1025 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
1026 AliAODTrack* part = fkAOD->GetTrack( ipart );
1027 if( !part->TestFilterBit( fFilterBit ) ) continue; // track cut selection
1028 if( !part->Charge() ) continue;
1029 if( part->Pt() < fTrackPtCut ) continue;
1030 tracks->AddLast( part );
1031 }
1032 QSortTracks( *tracks, 0, tracks->GetEntriesFast() );
1033
1034 nTracks = tracks->GetEntriesFast();
1035 if( !nTracks ) return 0;
1036
1037 return tracks;
1038 }
1039
1040
1041//____________________________________________________________________
1042const Bool_t AliAnalyseUE::TrackMCSelected(Double_t charge, Double_t pT, Double_t eta, Int_t pdgCode){
1043
1044 // MC track selection
1045 Double_t const kMyTolerance = 0.0000001;
1046 //if (charge == 0. || charge == -99.) return kFALSE;
1047 if (charge < kMyTolerance || charge + 99 < kMyTolerance) return kFALSE;
1048
1049 if ( fUseSingleCharge ) { // Charge selection
1050 if ( fUsePositiveCharge && charge < 0.) return kFALSE; // keep Positives
1051 if ( !fUsePositiveCharge && charge > 0.) return kFALSE; // keep Negatives
1052 }
1053
1054 //Kinematics cuts on particle
1055 if ((pT < fTrackPtCut) || (TMath::Abs(eta) > fTrackEtaCut )) return kFALSE;
1056
1057 Bool_t isHadron = TMath::Abs(pdgCode)==211 ||
1058 TMath::Abs(pdgCode)==2212 ||
1059 TMath::Abs(pdgCode)==321;
1060
1061 if ( fUseChargeHadrons && !isHadron ) return kFALSE;
1062
1063 return kTRUE;
1064
1065}
1066
1067//____________________________________________________________________
1068const Bool_t AliAnalyseUE::TrackSelected(AliAODTrack* part){
1069
1070 // Real track selection
1071 if ( !part->TestFilterBit(fFilterBit) ) return kFALSE; // track cut selection
1072 if ( !part->IsPrimaryCandidate()) return kFALSE; // reject whatever is not linked to collision point
1073 // PID Selection: Reject everything but hadrons
1074 Bool_t isHadron = part->GetMostProbablePID()==AliAODTrack::kPion ||
1075 part->GetMostProbablePID()==AliAODTrack::kKaon ||
1076 part->GetMostProbablePID()==AliAODTrack::kProton;
1077 if ( fUseChargeHadrons && !isHadron ) return kFALSE;
1078
1079 if ( !part->Charge() ) return kFALSE; //Only charged
1080 if ( fUseSingleCharge ) { // Charge selection
1081 if ( fUsePositiveCharge && part->Charge() < 0.) return kFALSE; // keep Positives
1082 if ( !fUsePositiveCharge && part->Charge() > 0.) return kFALSE; // keep Negatives
1083 }
1084
1085 if ( part->Pt() < fTrackPtCut ) return kFALSE;
1086 if( TMath::Abs(part->Eta()) > fTrackEtaCut ) return kFALSE;
1087
1088 return kTRUE;
1089}
1090
1091