]> 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
fast global merger installed. To use the previous version of the global merger, set...
[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"
3a9d4bcf 40#include "AliESDEvent.h"
1f329128 41#include "AliAODHandler.h"
42#include "AliAODInputHandler.h"
43#include "AliAODJet.h"
44#include "AliAODMCParticle.h"
45#include "AliAODTrack.h"
3a9d4bcf 46#include "AliESDtrack.h"
1f329128 47#include "AliKFVertex.h"
48#include "AliMCEvent.h"
49#include "AliMCEventHandler.h"
50#include "AliStack.h"
51
52#include "AliAnalysisHelperJetTasks.h"
53#include "AliGenPythiaEventHeader.h"
54#include "AliInputEventHandler.h"
55#include "AliLog.h"
56#include "AliStack.h"
57
58////////////////////////////////////////////////
59//---------------------------------------------
60// Class for transverse regions analysis
61//---------------------------------------------
62////////////////////////////////////////////////
63
64
65using namespace std;
66
67ClassImp(AliAnalyseUE)
68
69//-------------------------------------------------------------------
70AliAnalyseUE::AliAnalyseUE() :
71 TObject(),
72 //fTaskUE(0),
73 fkAOD(0x0),
3a9d4bcf 74 fkMC(0x0),
75 fkESD(0x0),
1f329128 76 fDebug(0),
77 fSimulateChJetPt(kFALSE),
3a9d4bcf 78 fStack(0x0),
1f329128 79 fAnaType(1),
80 fAreaReg(1.5393), // Pi*0.7*0.7
81 fConeRadius(0.7),
82 fFilterBit(0xFF),
83 fRegionType(1),
84 fUseChargeHadrons(kFALSE),
85 fUseChPartJet(kFALSE),
86 fUsePositiveCharge(kTRUE),
87 fUseSingleCharge(kFALSE),
88 fOrdering(1),
89 fJet1EtaCut(0.2),
90 fJet2DeltaPhiCut(2.616), // 150 degrees
91 fJet2RatioPtCut(0.8),
92 fJet3PtCut(15.),
93 fTrackEtaCut(0.9),
94 fTrackPtCut(0.),
95 fHistos(0x0),
96 fSumPtRegionPosit(0.),
97 fSumPtRegionNegat(0.),
98 fSumPtRegionForward(0.),
99 fSumPtRegionBackward(0.),
100 fMaxPartPtRegion(0.),
101 fNTrackRegionPosit(0),
102 fNTrackRegionNegat(0),
103 fNTrackRegionForward(0),
104 fNTrackRegionBackward(0),
3a9d4bcf 105 fSettingsTree(0x0),
106 fLtLabel(-999),
107 fLtMCLabel(-999)
1f329128 108{
109 // constructor
110}
111
112
113//-------------------------------------------------------------------
114AliAnalyseUE::AliAnalyseUE(const AliAnalyseUE & original) :
115 TObject(original),
116 //fTaskUE(original.fTaskUE),
117 fkAOD(original.fkAOD),
3a9d4bcf 118 fkMC(original.fkMC),
119 fkESD(original.fkESD),
1f329128 120 fDebug(original.fDebug),
121 fSimulateChJetPt(original.fSimulateChJetPt),
3a9d4bcf 122 fStack(original.fStack),
1f329128 123 fAnaType(original.fAnaType),
124 fAreaReg(original.fAreaReg),
125 fConeRadius(original.fConeRadius),
126 fFilterBit(original.fFilterBit),
127 fRegionType(original.fRegionType),
128 fUseChargeHadrons(original.fUseChargeHadrons),
129 fUseChPartJet(original.fUseChPartJet),
130 fUsePositiveCharge(original.fUsePositiveCharge),
131 fUseSingleCharge(original.fUseSingleCharge),
132 fOrdering(original.fOrdering),
133 fJet1EtaCut(original.fJet1EtaCut),
134 fJet2DeltaPhiCut(original.fJet2DeltaPhiCut),
135 fJet2RatioPtCut(original.fJet2RatioPtCut),
136 fJet3PtCut(original.fJet3PtCut),
137 fTrackEtaCut(original.fTrackEtaCut),
138 fTrackPtCut(original.fTrackPtCut),
139 fHistos(original.fHistos),
140 fSumPtRegionPosit(original.fSumPtRegionPosit),
141 fSumPtRegionNegat(original.fSumPtRegionNegat),
142 fSumPtRegionForward(original.fSumPtRegionForward),
143 fSumPtRegionBackward(original.fSumPtRegionBackward),
144 fMaxPartPtRegion(original.fMaxPartPtRegion),
145 fNTrackRegionPosit(original.fNTrackRegionPosit),
146 fNTrackRegionNegat(original.fNTrackRegionNegat),
147 fNTrackRegionForward(original.fNTrackRegionForward),
148 fNTrackRegionBackward(original.fNTrackRegionBackward),
3a9d4bcf 149 fSettingsTree(original.fSettingsTree),
150 fLtLabel(original.fLtLabel),
151 fLtMCLabel(original.fLtMCLabel)
1f329128 152{
153 //copy constructor
154}
155
156//-------------------------------------------------------------------
157AliAnalyseUE & AliAnalyseUE::operator = (const AliAnalyseUE & /*source*/)
158{
159 // assignment operator
160 return *this;
161}
162
163
164//-------------------------------------------------------------------
165AliAnalyseUE::~AliAnalyseUE(){
166
167 //clear memory
1f329128 168
169
170
171}
172
173
174//-------------------------------------------------------------------
175void AliAnalyseUE::AnalyseMC(TVector3 *jetVect,AliMCEvent *mcEvent, AliGenPythiaEventHeader *pythiaGenHeader,Int_t conePosition, Bool_t useAliStack, Bool_t constrainDistance, Double_t minDistance){
176
177 // Execute the analysis in case of MC input
178 fSumPtRegionPosit = 0.;
179 fSumPtRegionNegat = 0.;
180 fSumPtRegionForward = 0.;
181 fSumPtRegionBackward = 0.;
182 fMaxPartPtRegion = 0.;
183 fNTrackRegionPosit = 0;
184 fNTrackRegionNegat = 0;
185 fNTrackRegionForward = 0;
186 fNTrackRegionBackward = 0;
187
188 static Double_t const kPI = TMath::Pi();
189 static Double_t const k270rad = 270.*kPI/180.;
190
191 //Get Jets from MC header
192 Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
193 AliAODJet pythiaGenJets[4];
194 TVector3 jetVectnew[4];
195 Int_t iCount = 0;
196 for(int ip = 0;ip < nPythiaGenJets;++ip){
197 if (iCount>3) break;
198 Float_t p[4];
199 pythiaGenHeader->TriggerJet(ip,p);
200 TVector3 tempVect(p[0],p[1],p[2]);
201 if ( TMath::Abs(tempVect.Eta())>fJet1EtaCut ) continue;
202 pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
203 jetVectnew[iCount].SetXYZ(pythiaGenJets[iCount].Px(), pythiaGenJets[iCount].Py(), pythiaGenJets[iCount].Pz());
204 iCount++;
205 }
206
207 if (!iCount) return;// no jet in eta acceptance
208
209 //Search the index of the nearest MC jet to the leading jet reconstructed from the input data
210 Int_t index = 0;
211 if (constrainDistance){
212 Float_t deltaR = 0.;
213 Float_t dRTemp = 0.;
214 for (Int_t i=0; i<iCount; i++){
215 if (!i) {
216 dRTemp = jetVectnew[i].DeltaR(jetVect[0]);
217 index = i;
218 }
219
220 deltaR = jetVectnew[i].DeltaR(jetVect[0]);
221 if (deltaR < dRTemp){
222 index = i;
223 dRTemp = deltaR;
224 }
225 }
226
227 if (jetVectnew[index].DeltaR(jetVect[0]) > minDistance) return;
228 }
229
230 //Let's add some taste to jet and simulate pt of charged alone
231 //eta and phi are kept as original
232 //Play a Normal Distribution
233 Float_t random = 1.;
234 if (fSimulateChJetPt){
235 while(1){
236 random = gRandom->Gaus(0.6,0.25);
237 if (random > 0. && random < 1. &&
238 (random * jetVectnew[index].Pt()>6.)) break;
239 }
240 }
241
242 //Set new Pt & Fill histogram accordingly
243 Double_t maxPtJet1 = random * jetVectnew[index].Pt();
244
245 fHistos->FillHistogram("hEleadingPt", maxPtJet1 );
246
247 if (useAliStack){//Try Stack Information to perform UE analysis
248
249 AliStack* mcStack = mcEvent->Stack();//Load Stack
250 Int_t nTracksMC = mcStack->GetNtrack();
251 for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
252 //Cuts
253 if(!(mcStack->IsPhysicalPrimary(iTracks))) continue;
254
255 TParticle* mctrk = mcStack->Particle(iTracks);
256
257 Double_t charge = mctrk->GetPDG()->Charge();
258 Double_t pT = mctrk->Pt();
259 Double_t eta = mctrk->Eta();
260 Int_t pdgCode = mctrk->GetPdgCode();
261
262 if (!TrackMCSelected(charge, pT, eta, pdgCode))continue;
263
264 TVector3 partVect(mctrk->Px(), mctrk->Py(), mctrk->Pz());
265 Double_t deltaPhi = jetVectnew[index].DeltaPhi(partVect)+k270rad;
266 if( deltaPhi > 2.*TMath::Pi() ) deltaPhi-= 2.*TMath::Pi();
267 fHistos->FillHistogram("hdNdEtaPhiDist",deltaPhi, maxPtJet1 );
268
269 fHistos->FillHistogram("hFullRegPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
270
271 //We are not interested on stack organization but don't loose track of info
272
273 TVector3 tempVector = jetVectnew[0];
274 jetVectnew[0] = jetVectnew[index];
275 jetVectnew[index] = tempVector;
276
277 Int_t region = IsTrackInsideRegion( jetVectnew, &partVect, conePosition );
278
279 if (region == 1) {
280 if( fMaxPartPtRegion < mctrk->Pt() ) fMaxPartPtRegion = mctrk->Pt();
281 fSumPtRegionPosit += mctrk->Pt();
282 fNTrackRegionPosit++;
283 fHistos->FillHistogram("hTransRegPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
284 }
285 if (region == -1) {
286 if( fMaxPartPtRegion < mctrk->Pt() ) fMaxPartPtRegion = mctrk->Pt();
287 fSumPtRegionNegat += mctrk->Pt();
288 fNTrackRegionNegat++;
289 fHistos->FillHistogram("hTransRegPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
290 }
291 if (region == 2){ //forward
292 fSumPtRegionForward += mctrk->Pt();
293 fNTrackRegionForward++;
294 fHistos->FillHistogram("hRegForwardPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
295 }
296 if (region == -2){ //backward
297 fSumPtRegionBackward += mctrk->Pt();
298 fNTrackRegionBackward++;
299 fHistos->FillHistogram("hRegBackwardPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
300 }
301 } //end loop on stack particles
302 }else{//Try mc Particle
303
304 TClonesArray* farray = (TClonesArray*)fkAOD->FindListObject("mcparticles");
305
306 Int_t ntrks = farray->GetEntries();
307 if (fDebug>1) AliInfo(Form("In UE MC analysis tracks %d \n",ntrks));
308 for (Int_t i =0 ; i < ntrks; i++){
309 AliAODMCParticle* mctrk = (AliAODMCParticle*)farray->At(i);
310 //Cuts
311 if (!(mctrk->IsPhysicalPrimary())) continue;
312 //if (!(mctrk->IsPrimary())) continue;
313
314 Double_t charge = mctrk->Charge();
315 Double_t pT = mctrk->Pt();
316 Double_t eta = mctrk->Eta();
317 Int_t pdgCode = mctrk->GetPdgCode();
318
319 if (!TrackMCSelected(charge, pT, eta, pdgCode))continue;
320
321 TVector3 partVect(mctrk->Px(), mctrk->Py(), mctrk->Pz());
322
323 Double_t deltaPhi = jetVectnew[index].DeltaPhi(partVect)+k270rad;
324 if( deltaPhi > 2.*TMath::Pi() ) deltaPhi-= 2.*TMath::Pi();
325 fHistos->FillHistogram("hdNdEtaPhiDist", deltaPhi, maxPtJet1 );
326
327 fHistos->FillHistogram("hFullRegPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
328
329 //We are not interested on stack organization but don't loose track of info
330 TVector3 tempVector = jetVectnew[0];
331 jetVectnew[0] = jetVectnew[index];
332 jetVectnew[index] = tempVector;
333
334 Int_t region = IsTrackInsideRegion( jetVectnew, &partVect, conePosition );
335
336 if (region == 1) { //right
337 if( fMaxPartPtRegion < mctrk->Pt() ) fMaxPartPtRegion = mctrk->Pt();
338 fSumPtRegionPosit += mctrk->Pt();
339 fNTrackRegionPosit++;
340 fHistos->FillHistogram("hTransRegPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
341 }
342 if (region == -1) { //left
343 if( fMaxPartPtRegion < mctrk->Pt() ) fMaxPartPtRegion = mctrk->Pt();
344 fSumPtRegionNegat += mctrk->Pt();
345 fNTrackRegionNegat++;
346 fHistos->FillHistogram("hTransRegPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
347 }
348 if (region == 2){ //forward
349 fSumPtRegionForward += mctrk->Pt();
350 fNTrackRegionForward++;
351 fHistos->FillHistogram("hRegForwardPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
352 }
353 if (region == -2){ //backward
354 fSumPtRegionBackward += mctrk->Pt();
355 fNTrackRegionBackward++;
356 fHistos->FillHistogram("hRegBackwardPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
357 }
358
359 }//end loop AliAODMCParticle tracks
360 }
361}
362
363
364
365//-------------------------------------------------------------------
366Bool_t AliAnalyseUE::AnaTypeSelection(TVector3 *jetVect ){
367
368 // Cut events by jets topology
369 // anaType:
370 // 1 = inclusive,
371 // - Jet1 |eta| < jet1EtaCut
372 // 2 = back to back inclusive
373 // - fulfill case 1
374 // - |Jet1.Phi - Jet2.Phi| > jet2DeltaPhiCut
375 // - Jet2.Pt/Jet1Pt > jet2RatioPtCut
376 // 3 = back to back exclusive
377 // - fulfill case 2
378 // - Jet3.Pt < jet3PtCut
379
380 Double_t eta=jetVect[0].Eta();
381 if( TMath::Abs(eta) > fJet1EtaCut) {
382 if( fDebug > 1 ) AliInfo("\n Skipping Event...Jet1 |eta| > fJet1EtaCut");
383 return kFALSE;
384 }
385 // back to back inclusive
386 if( fAnaType > 1 && fAnaType < 4 && jetVect[1].Pt() < 0. ) {
387 if( fDebug > 1 ) AliInfo("\n Skipping Event... no second Jet found");
388 return kFALSE;
389 }
390 if( fAnaType > 1 && fAnaType < 4 && jetVect[1].Pt() > 0. ) {
391 if( TMath::Abs(jetVect[0].DeltaPhi(jetVect[1])) < fJet2DeltaPhiCut ||
392 jetVect[1].Pt()/jetVect[0].Pt() < fJet2RatioPtCut ) {
393 if( fDebug > 1 ) AliInfo("\n Skipping Event... |Jet1.Phi - Jet2.Phi| < fJet2DeltaPhiCut");
394 return kFALSE;
395 }
396 }
397 // back to back exclusive
398 if( fAnaType > 2 && fAnaType < 4 && jetVect[2].Pt() > 0. ) {
399 if( jetVect[2].Pt() > fJet3PtCut ) {
400 if( fDebug > 1 ) AliInfo("\n Skipping Event... Jet3.Pt > fJet3PtCut ");
401 return kFALSE;
402 }
403 }
404 return kTRUE;
405}
406
407
1f329128 408//-------------------------------------------------------------------
409void AliAnalyseUE::FillRegions(Bool_t isNorm2Area, TVector3 *jetVect){
410
411 // Fill the different topological regions
412 Double_t maxPtJet1 = jetVect[0].Pt();
413 static Double_t const kPI = TMath::Pi();
414 static Double_t const k120rad = 120.*kPI/180.;
415 Double_t const kMyTolerance = 0.0000001;
416
417 //Area for Normalization
418 // Forward and backward
419 Double_t normArea = 1.;
420 // Transverse
421 if (isNorm2Area) {
422 SetRegionArea(jetVect);
423 normArea = 2.*fTrackEtaCut*k120rad ;
424 } else fAreaReg = 1.;
425
426 Double_t avePosRegion = (fNTrackRegionPosit) ? fSumPtRegionPosit/fNTrackRegionPosit : 0.;
427 Double_t aveNegRegion = (fNTrackRegionNegat) ? fSumPtRegionNegat/fNTrackRegionNegat : 0.;
428 if( avePosRegion > aveNegRegion ) {
429 FillAvePartPtRegion( maxPtJet1, avePosRegion/fAreaReg, aveNegRegion/fAreaReg );
430 } else {
431 FillAvePartPtRegion( maxPtJet1, aveNegRegion/fAreaReg, avePosRegion/fAreaReg );
432 }
433
434 //How quantities will be sorted before Fill Min and Max Histogram
435 // 1=Plots will be CDF-like
436 // 2=Plots will be Marchesini-like
437 // 3=Minimum zone is selected as the one having lowest pt per track
438 if( fOrdering == 1 ) {
439 if( fSumPtRegionPosit > fSumPtRegionNegat ) {
440 FillSumPtRegion( maxPtJet1, fSumPtRegionPosit/fAreaReg, fSumPtRegionNegat/fAreaReg );
441 } else {
442 FillSumPtRegion( maxPtJet1, fSumPtRegionNegat/fAreaReg, fSumPtRegionPosit/fAreaReg );
443 }
444 if (fNTrackRegionPosit > fNTrackRegionNegat ) {
445 FillMultRegion( maxPtJet1, fNTrackRegionPosit/fAreaReg, fNTrackRegionNegat/fAreaReg, fSumPtRegionNegat/fAreaReg );
446 } else {
447 FillMultRegion( maxPtJet1, fNTrackRegionNegat/fAreaReg, fNTrackRegionPosit/fAreaReg, fSumPtRegionPosit/fAreaReg );
448 }
449 } else if( fOrdering == 2 ) {
450 if (fSumPtRegionPosit > fSumPtRegionNegat) {
451 FillSumPtRegion( maxPtJet1, fSumPtRegionPosit/fAreaReg, fSumPtRegionNegat/fAreaReg );
452 FillMultRegion( maxPtJet1, fNTrackRegionPosit/fAreaReg, fNTrackRegionNegat/fAreaReg, fSumPtRegionNegat/fAreaReg );
453 } else {
454 FillSumPtRegion( maxPtJet1, fSumPtRegionNegat/fAreaReg, fSumPtRegionPosit/fAreaReg );
455 FillMultRegion( maxPtJet1, fNTrackRegionNegat/fAreaReg, fNTrackRegionPosit/fAreaReg, fSumPtRegionPosit/fAreaReg );
456 }
457 } else if( fOrdering == 3 ){
458 if (avePosRegion > aveNegRegion) {
459 FillSumPtRegion( maxPtJet1, fSumPtRegionPosit/fAreaReg, fSumPtRegionNegat/fAreaReg );
460 FillMultRegion( maxPtJet1, fNTrackRegionPosit/fAreaReg, fNTrackRegionNegat/fAreaReg, fSumPtRegionNegat/fAreaReg );
461 }else{
462 FillSumPtRegion( maxPtJet1, fSumPtRegionNegat/fAreaReg, fSumPtRegionPosit/fAreaReg );
463 FillMultRegion( maxPtJet1, fNTrackRegionNegat/fAreaReg, fNTrackRegionPosit/fAreaReg, fSumPtRegionPosit/fAreaReg );
464 }
465 }
466 fHistos->FillHistogram("hRegionMaxPartPtMaxVsEt",maxPtJet1, fMaxPartPtRegion);
467
468 // Compute pedestal like magnitudes
469 fHistos->FillHistogram("hRegionDiffSumPtVsEt",maxPtJet1, (TMath::Abs(fSumPtRegionPosit-fSumPtRegionNegat)/(2.0*fAreaReg))+kMyTolerance);
470 fHistos->FillHistogram("hRegionAveSumPtVsEt", maxPtJet1, (fSumPtRegionPosit+fSumPtRegionNegat)/(2.0*fAreaReg));
471
472 // Transverse as a whole
473 fHistos->FillHistogram("hRegTransMult", maxPtJet1, fNTrackRegionPosit + fNTrackRegionNegat, (fNTrackRegionPosit + fNTrackRegionNegat)/(2.0*fAreaReg));
474 fHistos->FillHistogram("hRegTransSumPtVsMult",maxPtJet1, fNTrackRegionPosit + fNTrackRegionNegat , (fSumPtRegionNegat + fSumPtRegionPosit)/(2.0 *fAreaReg));
475
476 // Fill Histograms for Forward and away side w.r.t. leading jet direction
477 // Pt dependence
478 //fHistos->FillHistogram("hRegForwardSumPtVsEt",maxPtJet1, fSumPtRegionForward/normArea );
479 //fHistos->FillHistogram("hRegForwardMultVsEt",maxPtJet1, fNTrackRegionForward/normArea );
480 //fHistos->FillHistogram("hRegBackwardSumPtVsEt",maxPtJet1, fSumPtRegionBackward/normArea );
481 //fHistos->FillHistogram("hRegBackwardMultVsEt",maxPtJet1, fNTrackRegionBackward/normArea);
482
483 // Multiplicity dependence
484 fHistos->FillHistogram("hRegForwardMult", maxPtJet1, fNTrackRegionForward, fNTrackRegionForward/normArea);
485 fHistos->FillHistogram("hRegForwardSumPtvsMult", maxPtJet1, fNTrackRegionForward,fSumPtRegionForward/normArea);
486 fHistos->FillHistogram("hRegBackwardMult", maxPtJet1, fNTrackRegionBackward, fNTrackRegionBackward/normArea );
487 fHistos->FillHistogram("hRegBackwardSumPtvsMult", maxPtJet1, fNTrackRegionBackward,fSumPtRegionBackward/normArea);
488}
489
1f329128 490
491//-------------------------------------------------------------------
3a9d4bcf 492void AliAnalyseUE::FindMaxMinRegions(TVector3 *jetVect, Int_t conePosition, Int_t mctrue=0, Int_t eff=0){
493
494 // If mctrue = 1 consider branch AliAODMCParticles
495 // If eff = 1 track cuts for efficiency studies
496
1f329128 497 // Identify the different topological zones
498 fSumPtRegionPosit = 0.;
499 fSumPtRegionNegat = 0.;
500 fSumPtRegionForward = 0.;
501 fSumPtRegionBackward = 0.;
502 fMaxPartPtRegion = 0.;
503 fNTrackRegionPosit = 0;
504 fNTrackRegionNegat = 0;
505 fNTrackRegionForward = 0;
506 fNTrackRegionBackward = 0;
507 static Double_t const kPI = TMath::Pi();
508 static Double_t const kTWOPI = 2.*kPI;
509 static Double_t const k270rad = 270.*kPI/180.;
510 Double_t const kMyTolerance = 0.0000001;
511
3a9d4bcf 512 Int_t nTracks=0;
513 TClonesArray *tca = 0;
514 if (!mctrue) {
515 nTracks = fkAOD->GetNTracks();
516 if (fDebug > 1) AliInfo(Form(" ==== AOD tracks = %d \n ",nTracks));
517 }else{
518 tca = dynamic_cast<TClonesArray*>(fkAOD->FindListObject(AliAODMCParticle::StdBranchName()));
519 if(!tca){
e4cb1c80 520 Printf("%s:%d No AODMC Branch found !!!",(char*)__FILE__,__LINE__);
521 return;
522 }
3a9d4bcf 523 nTracks = tca->GetEntriesFast();
524 if (fDebug > 1) AliInfo(Form(" ==== AOD MC particles = %d \n ",nTracks));
525 }
1f329128 526
3a9d4bcf 527 //If UE task d0 distribution is not filled
528 Int_t flag_tmp=0;
529 if (fHistos->GetHistograms()->FindObject("hDCAxy")) flag_tmp = 1;
530
1f329128 531 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
3a9d4bcf 532 AliVParticle *part;
533 AliAODTrack *partRECO=0;
534 AliAODMCParticle *partMC=0;
535 Double_t charge=-999.;
536 Double_t pt=-999.;
537 Double_t eta=-999.;
538 Int_t pdgcode=-999;
539 TString suffix("");
540
541 // get track reco or true
542 if (!mctrue){
543 partRECO = dynamic_cast<AliAODTrack*>(fkAOD->GetTrack(ipart));
544 part = partRECO;
545 }
546 else {
547 partMC = dynamic_cast<AliAODMCParticle*>(tca->At(ipart));
548 part = partMC;
549 charge = partMC->Charge();
550 pt = partMC->Pt();
551 eta = partMC->Eta();
552 pdgcode = partMC->GetPdgCode();
553 suffix.Append("MC");
554 }
555
1f329128 556
1f329128 557 if (fDebug > 1) AliInfo(Form(" ==== AOD track = %d pt = %f charge = %d \n ",ipart,part->Pt(),part->Charge()));
3a9d4bcf 558
1f329128 559 // track selection
3a9d4bcf 560 if (!mctrue && !eff ){
561 if( !TrackSelected(partRECO)) continue; //track selection for data and MC reconstructed
562 if (flag_tmp){
563 if (fkESD && fkESD->GetNumberOfTracks() ){
564 AliInfo("READING ESD *************************************************");
565 Int_t id = partRECO->GetID();
566 AliESDtrack *trackESD;
567 trackESD = (AliESDtrack*)fkESD->GetTrack(id);
568 Float_t d0;
569 Float_t z;
570 trackESD->GetImpactParameters(d0,z);
571 fHistos->FillHistogram("hDCAxy", d0, jetVect[0].Pt());
572 }else AliInfo("NO TRACKS ************************************************")
573 }
574 }
575
576 if (!mctrue && eff){
577 if (!TrackSelectedEfficiency(partRECO )) continue; //track selection for MC reconstructed for efficiency studies
578 if(fkESD && fkESD->GetNumberOfTracks()){
579 Int_t id = partRECO->GetID();
580 AliESDtrack * trackESD = (AliESDtrack*) fkESD->GetTrack(id);
581 Float_t d0;
582 Float_t z;
583 trackESD->GetImpactParameters(d0,z);
584 fHistos->FillHistogram("hDCAxyPrimary", d0, jetVect[0].Pt());
585 }
586 }
587
588 if (mctrue){
589 if ( !(TrackMCSelected(charge, pt, eta, pdgcode) && partMC->IsPhysicalPrimary())) continue; //track selection for MC true
590 }
1f329128 591
592 TVector3 partVect(part->Px(), part->Py(), part->Pz());
593 Bool_t isFlagPart = kTRUE;
594 Double_t deltaPhi = jetVect[0].DeltaPhi(partVect)+k270rad;
595 if( deltaPhi > kTWOPI ) deltaPhi-= kTWOPI;
3a9d4bcf 596 if (fAnaType != 4 ) fHistos->FillHistogram(Form("hdNdEtaPhiDist%s",suffix.Data()),deltaPhi, jetVect[0].Pt());
1f329128 597 else if (TMath::Abs(deltaPhi-k270rad) >= kMyTolerance && TMath::Abs(jetVect[0].Eta()-partVect.Eta()) >= kMyTolerance){
3a9d4bcf 598 fHistos->FillHistogram(Form("hdNdEtaPhiDist%s",suffix.Data()),deltaPhi, jetVect[0].Pt());
1f329128 599 isFlagPart = kFALSE;
600 }
601
3a9d4bcf 602 fHistos->FillHistogram(Form("hFullRegPartPtDistVsEt%s",suffix.Data()), part->Pt(), jetVect[0].Pt());
1f329128 603
604 Int_t region = IsTrackInsideRegion( jetVect, &partVect, conePosition );
605 if (region == 1) {
606 if( fMaxPartPtRegion < part->Pt() ) fMaxPartPtRegion = part->Pt();
607 fSumPtRegionPosit += part->Pt();
608 fNTrackRegionPosit++;
3a9d4bcf 609 fHistos->FillHistogram(Form("hTransRegPartPtDistVsEt%s",suffix.Data()),part->Pt(), jetVect[0].Pt());
1f329128 610 }
611 if (region == -1) {
612 if( fMaxPartPtRegion < part->Pt() ) fMaxPartPtRegion = part->Pt();
613 fSumPtRegionNegat += part->Pt();
614 fNTrackRegionNegat++;
3a9d4bcf 615 fHistos->FillHistogram(Form("hTransRegPartPtDistVsEt%s",suffix.Data()),part->Pt(), jetVect[0].Pt());
1f329128 616 }
617 if (region == 2){ //forward
618 fSumPtRegionForward += part->Pt();
619 fNTrackRegionForward++;
3a9d4bcf 620 fHistos->FillHistogram(Form("hRegForwardPartPtDistVsEt%s",suffix.Data()),part->Pt(), jetVect[0].Pt());
1f329128 621 }
622 if (region == -2){ //backward
623 fSumPtRegionBackward += part->Pt();
624 fNTrackRegionBackward++;
3a9d4bcf 625 fHistos->FillHistogram(Form("hRegBackwardPartPtDistVsEt%s",suffix.Data()),part->Pt(), jetVect[0].Pt());
1f329128 626 }
627 }//end loop AOD tracks
628
629}
630
3a9d4bcf 631//-------------------------------------------------------------------
632/*TVector3 AliAnalyseUE::GetLeadingTracksMC(AliMCEvent *mcEvent){
633
634 fkMC = mcEvent;
635
636 Double_t maxPtJet1 = 0.;
637 Int_t index1 = -1;
638
639 TVector3 jetVect[3];
640
641 jetVect[0].SetPtEtaPhi(-1.,-1.,-1.);
642 jetVect[1].SetPtEtaPhi(-1.,-1.,-1.);
643 jetVect[2].SetPtEtaPhi(-1.,-1.,-1.);
644
645 Int_t nJets = 0;
646
647 TObjArray* tracks = SortChargedParticlesMC();
648 if( tracks ) {
649 nJets = tracks->GetEntriesFast();
650 if( nJets > 0 ) {
651 index1 = 0;
652 TParticle* jet = (TParticle*)tracks->At(0);
653 maxPtJet1 = jet->Pt();
654 jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
655 }
656 tracks->Clear();
657 delete tracks;
658 }
659 return *jetVect;
660
661}
662*/
663
664//-------------------------------------------------------------------
665TVector3 AliAnalyseUE::GetLeadingTracksMC(AliMCEvent *mcEvent){
666
667 fkMC = mcEvent;
668
669 Double_t maxPtJet1 = 0.;
670 Int_t index1 = -1;
671
672 TVector3 jetVect[3];
673
674 jetVect[0].SetPtEtaPhi(-1.,-1.,-1.);
675 jetVect[1].SetPtEtaPhi(-1.,-1.,-1.);
676 jetVect[2].SetPtEtaPhi(-1.,-1.,-1.);
677
678 Int_t nJets = 0;
679
680 TObjArray* tracks = SortChargedParticlesMC();
681 if( tracks ) {
682 nJets = tracks->GetEntriesFast();
683 if( nJets > 0 ) {
684 index1 = 0;
685 AliAODMCParticle* jet = (AliAODMCParticle*)tracks->At(0);
686 fLtMCLabel = jet->GetLabel();
687 maxPtJet1 = jet->Pt();
688 jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
689 }
690 tracks->Clear();
691 delete tracks;
692 }
693 return *jetVect;
694
695}
696
1f329128 697//-------------------------------------------------------------------
698TVector3 AliAnalyseUE::GetOrderedClusters(TString aodBranch, Bool_t chargedJets, Double_t chJetPtMin){
699
700 // jets from AOD, on-the-fly or leading particle
701 Double_t maxPtJet1 = 0.;
702 Int_t index1 = -1;
703 Double_t maxPtJet2 = 0.; // jet 2 need for back to back inclusive
704 Int_t index2 = -1;
705 Double_t maxPtJet3 = 0.; // jet 3 need for back to back exclusive
706 Int_t index3 = -1;
707 TVector3 jetVect[3];
708
709 jetVect[0].SetPtEtaPhi(-1.,-1.,-1.);
710 jetVect[1].SetPtEtaPhi(-1.,-1.,-1.);
711 jetVect[2].SetPtEtaPhi(-1.,-1.,-1.);
712
713 Int_t nJets = 0;
714 //TClonesArray* fArrayJets;
715 TObjArray* arrayJets;
716 // 1) JETS FROM AOD BRANCH (standard, non-standard or delta)
717 if (!chargedJets && fAnaType != 4 ) {
718 AliInfo(" ==== Read AODs !");
719 AliInfo(Form(" ==== Reading Branch: %s ", aodBranch.Data()));
720 arrayJets = (TObjArray*)fkAOD->GetList()->FindObject(aodBranch.Data());
721 if (!arrayJets){
722 AliFatal(" No jet-array! ");
723 return *jetVect;
724 }
725
726 // Find Leading Jets 1,2,3
727 // (could be skipped if Jets are sort by Pt...)
728 nJets=arrayJets->GetEntries();
729 for( Int_t i=0; i<nJets; ++i ) {
730 AliAODJet* jet = (AliAODJet*)arrayJets->At(i);
731 Double_t jetPt = jet->Pt();//*1.666; // FIXME Jet Pt Correction ?????!!!
732
733 if( jetPt > maxPtJet1 ) {
734 maxPtJet3 = maxPtJet2; index3 = index2;
735 maxPtJet2 = maxPtJet1; index2 = index1;
736 maxPtJet1 = jetPt; index1 = i;
737 } else if( jetPt > maxPtJet2 ) {
738 maxPtJet3 = maxPtJet2; index3 = index2;
739 maxPtJet2 = jetPt; index2 = i;
740 } else if( jetPt > maxPtJet3 ) {
741 maxPtJet3 = jetPt; index3 = i;
742 }
743 }
744
745 if( index1 != -1 ) {
746 AliAODJet *jet =(AliAODJet*) arrayJets->At(index1);
747 if(jet)jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
748 }
749 if( index2 != -1 ) {
750 AliAODJet* jet= (AliAODJet*) arrayJets->At(index2);
751 if(jet)jetVect[1].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
752 }
753 if( index3 != -1 ) {
754 AliAODJet* jet = (AliAODJet*) arrayJets->At(index3);
755 if(jet)jetVect[2].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
756 }
757
758 }
759
760
761 // 2) ON-THE-FLY CDF ALGORITHM
762 if (chargedJets){
1f329128 763 arrayJets = FindChargedParticleJets(chJetPtMin);
764 if( arrayJets ) {
765 nJets = arrayJets->GetEntriesFast();
766 if( nJets > 0 ) {
767 index1 = 0;
768 AliAODJet* jet = (AliAODJet*)arrayJets->At(0);
769 maxPtJet1 = jet->Pt();
770 jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
771 }
772 if( nJets > 1 ) {
773 index2 = 1;
774 AliAODJet* jet = (AliAODJet*)arrayJets->At(1);
775 maxPtJet2 = jet->Pt();
776 jetVect[1].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
777 }
778 if( nJets > 2 ) {
779 index3 = 2;
780 AliAODJet* jet = (AliAODJet*)arrayJets->At(2);
781 maxPtJet3 = jet->Pt();
782 jetVect[2].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
783 }
784
785 arrayJets->Delete();
786 delete arrayJets;
787 }
788 }
789
790
791 // 3) LEADING PARTICLE
792 if( fAnaType == 4 ){
793 TObjArray* tracks = SortChargedParticles();
794 if( tracks ) {
795 nJets = tracks->GetEntriesFast();
796 if( nJets > 0 ) {
797 index1 = 0;
798 AliAODTrack* jet = (AliAODTrack*)tracks->At(0);
3a9d4bcf 799 fLtLabel = jet->GetLabel();
800 maxPtJet1 = jet->Pt();
1f329128 801 jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
802 }
803 tracks->Clear();
804 delete tracks;
805 }
806
807 }
ddbad1f5 808 fHistos->FillHistogram("hNJets",nJets);
1f329128 809
810 return *jetVect;
811
812}
813
814
815//-------------------------------------------------------------------
816void AliAnalyseUE::Initialize(AliAnalysisTaskUE& taskUE){
3a9d4bcf 817/*void AliAnalyseUE::Initialize(AliAnalysisTask& task){// when correction task is in trunk
818 if (task->InheritsFrom("AliAnalysisTaskUE")){
819 AliAnalysisTaskUE *taskUE = dynamic_cast<AliAnalysisTaskUE*> task;
820 }
821 else if (task->InheritsFrom("AliAnalysisTaskCorrectionsUE")){
822 AliAnalysisTaskCorrectionsUE *taskUE = dynamic_cast<AliAnalysisTaskCorrectionsUE*> task;
823 }
824
825*/
1f329128 826 //Get principal settings from current instance of UE analysis-task
827 fAnaType = taskUE.GetAnaTopology();
828 fkAOD = taskUE.GetAOD();
829 fConeRadius = taskUE.GetConeRadius();
830 fDebug = taskUE.GetDebugLevel();
831 fFilterBit = taskUE.GetFilterBit();
832 fJet1EtaCut = taskUE.GetJet1EtaCut();
833 fJet2DeltaPhiCut = taskUE.GetJet2DeltaPhiCut();
834 fJet2RatioPtCut = taskUE.GetJet2RatioPtCut();
835 fJet3PtCut = taskUE.GetJet3PtCut();
836 fOrdering = taskUE.GetPtSumOrdering() ;
837 fRegionType = taskUE.GetRegionType();
838 fSimulateChJetPt = taskUE.GetSimulateChJetPt();
839 fTrackEtaCut = taskUE.GetTrackEtaCut();
840 fTrackPtCut = taskUE.GetTrackPtCut();
ddbad1f5 841 fHistos = taskUE.fHistosUE;
1f329128 842 fUseChargeHadrons = taskUE.GetUseChargeHadrons();
843 fUseChPartJet = taskUE.GetUseChPartJet();
844 fUsePositiveCharge = taskUE.GetUseNegativeChargeType();
845 fUseSingleCharge = taskUE.GetUseSingleCharge();
846
1f329128 847}
848
3a9d4bcf 849//-------------------------------------------------------------------
850void AliAnalyseUE::Initialize(Int_t anaType, AliAODEvent* aod,Double_t coneRadius, Int_t debug, Int_t filterBit, Double_t jet1EtaCut, Double_t jet2DeltaPhiCut, Double_t jet2RatioPtCut, Double_t jet3PtCut, Int_t ordering, Int_t regionType,Bool_t simulateChJetPt, Double_t trackEtaCut, Double_t trackPtCut, Bool_t useChargeHadrons, Bool_t useChPartJet, Bool_t usePositiveCharge, Bool_t useSingleCharge, AliHistogramsUE* histos){
851
852 fAnaType = anaType;
853 fkAOD = aod;
854 fConeRadius = coneRadius;
855 fDebug = debug;
856 fFilterBit = filterBit;
857 fJet1EtaCut = jet1EtaCut;
858 fJet2DeltaPhiCut = jet2DeltaPhiCut;
859 fJet2RatioPtCut = jet2RatioPtCut;
860 fJet3PtCut = jet3PtCut;
861 fOrdering = ordering;
862 fRegionType = regionType;
863 fSimulateChJetPt = simulateChJetPt;
864 fTrackEtaCut = trackEtaCut;
865 fTrackPtCut = trackPtCut;
866 fUseChargeHadrons = useChargeHadrons;
867 fUseChPartJet = useChPartJet;
868 fUsePositiveCharge = usePositiveCharge;
869 fUseSingleCharge = useSingleCharge;
870 fHistos = histos;
871}
872
873
874//-------------------------------------------------------------------
875Bool_t AliAnalyseUE::TriggerSelection(AliInputEventHandler* inputHandler){
876
877 //Use AliPhysicsSelection to select good events
878 if( !inputHandler->InheritsFrom("AliAODInputHandler") ) { // input AOD
879 if (inputHandler->IsEventSelected()) {
880 if (fDebug > 1) AliInfo(" Trigger Selection: event ACCEPTED ... ");
881 } else {
882 if (fDebug > 1) AliInfo(" Trigger Selection: event REJECTED ... ");
883 return kFALSE;
884 }
885 }
886
887 return kTRUE;
888
889}
890
891
1f329128 892//-------------------------------------------------------------------
1f329128 893Bool_t AliAnalyseUE::VertexSelection(AliAODEvent *aod, Int_t tracks, Double_t zed ){
894
895 //Require 1 vertex (no TPC stand-alone) with a minimum number of tracks and z-coordinate in a limited range
896 Int_t nVertex = aod->GetNumberOfVertices();
897 if( nVertex > 0 ) { // Only one vertex (reject pileup)
898 AliAODVertex* vertex = (AliAODVertex*)aod->GetPrimaryVertex();
899 Int_t nTracksPrim = vertex->GetNContributors();
900 Double_t zVertex = vertex->GetZ();
901 if (fDebug > 1) AliInfo(Form(" Vertex in = %f with %d particles by %s data ...",zVertex,nTracksPrim,vertex->GetName()));
902 // Select a quality vertex by number of tracks?
903 if( nTracksPrim < tracks || TMath::Abs(zVertex) > zed ) {
904 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
905 return kFALSE;
906 }
907 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
908 } else {
909 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
910 return kFALSE;
911 }
912
913 return kTRUE;
914}
915
ddbad1f5 916//-------------------------------------------------------------------
917Bool_t AliAnalyseUE::VertexSelectionOld(AliAODEvent *aod ){
918
919 AliKFVertex primVtx(*(aod->GetPrimaryVertex()));
920 Int_t nTracksPrim=primVtx.GetNContributors();
921 if (fDebug > 1) AliInfo(Form(" Primary-vertex Selection: %d",nTracksPrim));
922 if(!nTracksPrim){
923 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
924 return kFALSE;
925 }
926 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED ...");
927
928 return kTRUE;
929}
930
1f329128 931// PRIVATE METHODS **************************************************
932
933TObjArray* AliAnalyseUE::FindChargedParticleJets( Double_t chJetPtMin )
934{
935 // Return a TObjArray of "charged particle jets"
936
937 // Charged particle jet definition from reference:
938 // "Charged jet evolution and the underlying event
939 // in proton-antiproton collisions at 1.8 TeV"
940 // PHYSICAL REVIEW D 65 092002, CDF Collaboration
941
942 // We defined "jets" as circular regions in eta-phi space with
943 // radius defined by R = sqrt( (eta-eta0)^2 +(phi-phi0)^2 ).
944 // Our jet algorithm is as follows:
945 // 1- Order all charged particles according to their pT .
946 // 2- Start with the highest pT particle and include in the jet all
947 // particles within the radius R=0.7 considering each particle
948 // in the order of decreasing pT and recalculating the centroid
949 // of the jet after each new particle is added to the jet .
950 // 3- Go to the next highest pT particle not already included in
951 // a jet and add to the jet all particles not already included in
952 // a jet within R=0.7.
953 // 4- Continue until all particles are in a jet.
954 // We defined the transverse momentum of the jet to be
955 // the scalar pT sum of all the particles within the jet, where pT
956 // is measured with respect to the beam axis
957
958 // 1 - Order all charged particles according to their pT .
959 Int_t nTracks = fkAOD->GetNTracks();
960 if( !nTracks ) return 0;
961 TObjArray tracks(nTracks);
962
963 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
964 AliAODTrack* part = fkAOD->GetTrack( ipart );
965 if( !part->TestFilterBit(fFilterBit) ) continue; // track cut selection
966 if( !part->Charge() ) continue;
967 if( part->Pt() < fTrackPtCut ) continue;
968 tracks.AddLast(part);
969 }
970 QSortTracks( tracks, 0, tracks.GetEntriesFast() );
971
972 nTracks = tracks.GetEntriesFast();
973 if( !nTracks ) return 0;
974
975 TObjArray *jets = new TObjArray(nTracks);
976 TIter itrack(&tracks);
977 while( nTracks ) {
978 // 2- Start with the highest pT particle ...
979 Float_t px,py,pz,pt;
980 AliAODTrack* track = (AliAODTrack*)itrack.Next();
981 if( !track ) continue;
982 px = track->Px();
983 py = track->Py();
984 pz = track->Pz();
985 pt = track->Pt(); // Use the energy member to store Pt
986 jets->AddLast( new TLorentzVector(px, py, pz, pt) );
987 tracks.Remove( track );
988 TLorentzVector* jet = (TLorentzVector*)jets->Last();
989 jet->SetPtEtaPhiE( 1., jet->Eta(), jet->Phi(), pt );
990 // 3- Go to the next highest pT particle not already included...
991 AliAODTrack* track1;
992 Double_t fPt = jet->E();
993 while ( (track1 = (AliAODTrack*)(itrack.Next())) ) {
994 Double_t tphi = track1->Phi(); // here Phi is from 0 <-> 2Pi
995 if (tphi > TMath::Pi()) tphi -= 2. * TMath::Pi(); // convert to -Pi <-> Pi
996 Double_t dphi = TVector2::Phi_mpi_pi(jet->Phi()-tphi);
997 Double_t r = TMath::Sqrt( (jet->Eta()-track1->Eta())*(jet->Eta()-track1->Eta()) +dphi*dphi );
998 if( r < fConeRadius ) {
999 fPt = jet->E()+track1->Pt(); // Scalar sum of Pt
1000 // recalculating the centroid
1001 Double_t eta = jet->Eta()*jet->E()/fPt + track1->Eta()*track1->Pt()/fPt;
1002 Double_t phi = jet->Phi()*jet->E()/fPt + tphi*track1->Pt()/fPt;
1003 jet->SetPtEtaPhiE( 1., eta, phi, fPt );
1004 tracks.Remove( track1 );
1005 }
1006 }
1007
1008 tracks.Compress();
1009 nTracks = tracks.GetEntries();
1010 // 4- Continue until all particles are in a jet.
1011 itrack.Reset();
1012 } // end while nTracks
1013
1014 // Convert to AODjets....
1015 Int_t njets = jets->GetEntriesFast();
1016 TObjArray* aodjets = new TObjArray(njets);
1017 aodjets->SetOwner(kTRUE);
1018 for(Int_t ijet=0; ijet<njets; ++ijet) {
1019 TLorentzVector* jet = (TLorentzVector*)jets->At(ijet);
1020 if (jet->E() < chJetPtMin) continue;
1021 Float_t px, py,pz,en; // convert to 4-vector
1022 px = jet->E() * TMath::Cos(jet->Phi()); // Pt * cos(phi)
1023 py = jet->E() * TMath::Sin(jet->Phi()); // Pt * sin(phi)
1024 pz = jet->E() / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-jet->Eta())));
1025 en = TMath::Sqrt(px * px + py * py + pz * pz);
1026
1027 aodjets->AddLast( new AliAODJet(px, py, pz, en) );
1028 }
1029 jets->Delete();
1030 delete jets;
1031
1032 // Order jets according to their pT .
1033 QSortTracks( *aodjets, 0, aodjets->GetEntriesFast() );
1034
1035 // debug
1036 if (fDebug>3) AliInfo(Form(" %d Charged jets found\n",njets));
1037
1038 return aodjets;
1039}
1040
1041
1042//____________________________________________________________________
1043void AliAnalyseUE::FillAvePartPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin )
1044{
1045
1046 // Fill average particle Pt of control regions
1047
1048 // Max cone
1049 fHistos->FillHistogram("hRegionAvePartPtMaxVsEt", leadingE, ptMax);
1050 // Min cone
1051 fHistos->FillHistogram("hRegionAvePartPtMinVsEt", leadingE, ptMin);
1052 // MAke distributions for UE comparison with MB data
1053 fHistos->FillHistogram("hMinRegAvePt", ptMin);
1054
1055}
1056
1057//____________________________________________________________________
1058void AliAnalyseUE::FillMultRegion( Double_t leadingE, Double_t nTrackPtmax, Double_t nTrackPtmin, Double_t ptMin )
1059{
1060
1061 // Fill Nch multiplicity of control regions
1062
1063 // Max cone
1064 fHistos->FillHistogram("hRegionMultMaxVsEt", leadingE, nTrackPtmax);
1065 fHistos->FillHistogram("hRegionMultMax", nTrackPtmax);
1066
1067 // Min cone
1068 fHistos->FillHistogram("hRegionMultMinVsEt", leadingE, nTrackPtmin );
1069 fHistos->FillHistogram("hRegionMultMin", nTrackPtmin);
1070
1071 // MAke distributions for UE comparison with MB data
1072 fHistos->FillHistogram("hMinRegSumPtvsMult", nTrackPtmin,ptMin);
1073
1074}
1075
1076//____________________________________________________________________
1077void AliAnalyseUE::FillSumPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin )
1078{
1079 // Fill sumPt of control regions
1080
1081 // Max cone
1082 fHistos->FillHistogram("hRegionSumPtMaxVsEt", leadingE, ptMax);
1083
1084 // Min cone
1085 fHistos->FillHistogram("hRegionSumPtMinVsEt", leadingE, ptMin);
1086
1087 // MAke distributions for UE comparison with MB data
1088 fHistos->FillHistogram("hMinRegSumPt", ptMin);
1089 fHistos->FillHistogram("hMinRegSumPtJetPtBin", leadingE, ptMin);
1090 fHistos->FillHistogram("hMaxRegSumPtJetPtBin", leadingE, ptMax);
1091
1092}
1093
1094//-------------------------------------------------------------------
1095Int_t AliAnalyseUE::IsTrackInsideRegion(TVector3 *jetVect, TVector3 *partVect, Int_t conePosition)
1096{
1097 // return de region in delta phi
1098 // -1 negative delta phi
1099 // 1 positive delta phi
1100 // 0 outside region
1101 static const Double_t k60rad = 60.*TMath::Pi()/180.;
1102 static const Double_t k120rad = 120.*TMath::Pi()/180.;
1103
1104 Int_t region = 0;
1105 if( fRegionType == 1 ) {
1106 if( TMath::Abs(partVect->Eta()) > fTrackEtaCut ) return 0;
1107 // transverse regions
1108 if (jetVect[0].DeltaPhi(*partVect) < -k60rad && jetVect[0].DeltaPhi(*partVect) > -k120rad ) region = -1; //left
1109 if (jetVect[0].DeltaPhi(*partVect) > k60rad && jetVect[0].DeltaPhi(*partVect) < k120rad ) region = 1; //right
1110 if (TMath::Abs(jetVect[0].DeltaPhi(*partVect)) < k60rad ) region = 2; //forward
1111 if (TMath::Abs(jetVect[0].DeltaPhi(*partVect)) > k120rad ) region = -2; //backward
1112
1113 } else if( fRegionType == 2 ) {
1114 // Cone regions
1115 Double_t deltaR = 0.;
1116
1117 TVector3 positVect,negatVect;
1118 if (conePosition==1){
1119 positVect.SetMagThetaPhi(1, 2.*atan(exp(-jetVect[0].Eta())), jetVect[0].Phi()+TMath::PiOver2());
1120 negatVect.SetMagThetaPhi(1, 2.*atan(exp(-jetVect[0].Eta())), jetVect[0].Phi()-TMath::PiOver2());
1121 }else if (conePosition==2){
1122 if(fAnaType<2) AliFatal("Prevent error in Analysis type there might be only 1 jet. To avoid overflow better Correct UE config");
1123 positVect.SetMagThetaPhi(1, 2.*atan(exp(-(jetVect[0].Eta()+jetVect[1].Eta())/2.)), jetVect[0].Phi()+TMath::PiOver2());
1124 negatVect.SetMagThetaPhi(1, 2.*atan(exp(-(jetVect[0].Eta()+jetVect[1].Eta())/2.)), jetVect[0].Phi()-TMath::PiOver2());
1125 }else if (conePosition==3){
1126 if(fAnaType<2) AliFatal("Prevent error in Analysis type there might be only 1 jet. To avoid overflow better Correct UE config");
1127 Double_t weightEta = jetVect[0].Eta() * jetVect[0].Pt()/(jetVect[0].Pt() + jetVect[1].Pt()) +
1128 jetVect[1].Eta() * jetVect[1].Pt()/(jetVect[0].Pt() + jetVect[1].Pt());
1129 //Double_t weightEta = jetVect[0].Eta() * jetVect[0].Mag()/(jetVect[0].Mag() + jetVect[1].Mag()) +
1130 // jetVect[1].Eta() * jetVect[1].Mag()/(jetVect[0].Mag() + jetVect[1].Mag());
1131 positVect.SetMagThetaPhi(1, 2.*atan(exp(-weightEta)), jetVect[0].Phi()+TMath::PiOver2());
1132 negatVect.SetMagThetaPhi(1, 2.*atan(exp(-weightEta)), jetVect[0].Phi()-TMath::PiOver2());
1133 }
1134 if (TMath::Abs(positVect.DeltaPhi(*partVect)) < fConeRadius ) {
1135 region = 1;
1136 deltaR = positVect.DrEtaPhi(*partVect);
1137 } else if (TMath::Abs(negatVect.DeltaPhi(*partVect)) < fConeRadius) {
1138 region = -1;
1139 deltaR = negatVect.DrEtaPhi(*partVect);
1140 }
1141
1142 if (deltaR > fConeRadius) region = 0;
1143
1144 } else
1145 AliError("Unknow region type");
1146
1147 return region;
1148 }
1149
1150
1151
1152//-------------------------------------------------------------------
1153void AliAnalyseUE::QSortTracks(TObjArray &a, Int_t first, Int_t last)
1154{
1155 // Sort array of TObjArray of tracks by Pt using a quicksort algorithm.
1156
1157 static TObject *tmp;
1158 static int i; // "static" to save stack space
1159 int j;
1160
1161 while (last - first > 1) {
1162 i = first;
1163 j = last;
1164 for (;;) {
1165 while (++i < last && ((AliVParticle*)a[i])->Pt() > ((AliVParticle*)a[first])->Pt() )
1166 ;
1167 while (--j > first && ((AliVParticle*)a[j])->Pt() < ((AliVParticle*)a[first])->Pt() )
1168 ;
1169 if (i >= j)
1170 break;
1171
1172 tmp = a[i];
1173 a[i] = a[j];
1174 a[j] = tmp;
1175 }
1176 if (j == first) {
1177 ++first;
1178 continue;
1179 }
1180 tmp = a[first];
1181 a[first] = a[j];
1182 a[j] = tmp;
1183 if (j - first < last - (j + 1)) {
1184 QSortTracks(a, first, j);
1185 first = j + 1; // QSortTracks(j + 1, last);
1186 } else {
1187 QSortTracks(a, j + 1, last);
1188 last = j; // QSortTracks(first, j);
1189 }
1190 }
1191}
1192
1193
3a9d4bcf 1194//-------------------------------------------------------------------
1195void AliAnalyseUE::QSortTracksMC(TObjArray &a, Int_t first, Int_t last)
1196{
1197 // Sort array of TObjArray of tracks by Pt using a quicksort algorithm.
1198
1199 static TObject *tmp;
1200 static int i; // "static" to save stack space
1201 int j;
1202
1203 while (last - first > 1) {
1204 i = first;
1205 j = last;
1206 for (;;) {
1207 while (++i < last && ((AliAODMCParticle*)a[i])->Pt() > ((AliAODMCParticle*)a[first])->Pt() )
1208 ;
1209 while (--j > first && ((AliAODMCParticle*)a[j])->Pt() < ((AliAODMCParticle*)a[first])->Pt() )
1210 ;
1211 if (i >= j)
1212 break;
1213
1214 tmp = a[i];
1215 a[i] = a[j];
1216 a[j] = tmp;
1217 }
1218 if (j == first) {
1219 ++first;
1220 continue;
1221 }
1222 tmp = a[first];
1223 a[first] = a[j];
1224 a[j] = tmp;
1225 if (j - first < last - (j + 1)) {
1226 QSortTracksMC(a, first, j);
1227 first = j + 1; // QSortTracks(j + 1, last);
1228 } else {
1229 QSortTracksMC(a, j + 1, last);
1230 last = j; // QSortTracks(first, j);
1231 }
1232 }
1233}
1234
1235
1236
1f329128 1237
1238
1239//-------------------------------------------------------------------
1240void AliAnalyseUE::SetRegionArea(TVector3 *jetVect)
1241{
1242 // Set region area
1243 Double_t areaCorrFactor=0.;
1244 Double_t deltaEta = 0.;
1245 if (fRegionType==1) fAreaReg = 2.*fTrackEtaCut*60.*TMath::Pi()/180.;
1246 else if (fRegionType==2){
1247 deltaEta = 0.9-TMath::Abs(jetVect[0].Eta());
1248 if (deltaEta>fConeRadius) fAreaReg = TMath::Pi()*fConeRadius*fConeRadius;
1249 else{
1250 areaCorrFactor = fConeRadius*fConeRadius*TMath::ACos(deltaEta/fConeRadius) -
1251 (deltaEta*fConeRadius)*TMath::Sqrt( 1. - deltaEta*deltaEta/(fConeRadius*fConeRadius));
1252 fAreaReg=TMath::Pi()*fConeRadius*fConeRadius-areaCorrFactor;
1253 }
1254 }else AliWarning("Unknown Region Type");
1255 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));
1256}
1257
1258
1259//____________________________________________________________________
1260TObjArray* AliAnalyseUE::SortChargedParticles()
1261{
1262 // return an array with all charged particles ordered according to their pT .
1263 Int_t nTracks = fkAOD->GetNTracks();
1264 if( !nTracks ) return 0;
1265 TObjArray* tracks = new TObjArray(nTracks);
1266
1267 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
1268 AliAODTrack* part = fkAOD->GetTrack( ipart );
1269 if( !part->TestFilterBit( fFilterBit ) ) continue; // track cut selection
1270 if( !part->Charge() ) continue;
1271 if( part->Pt() < fTrackPtCut ) continue;
1272 tracks->AddLast( part );
1273 }
1274 QSortTracks( *tracks, 0, tracks->GetEntriesFast() );
1275
1276 nTracks = tracks->GetEntriesFast();
1277 if( !nTracks ) return 0;
1278
1279 return tracks;
1280 }
1281
3a9d4bcf 1282//____________________________________________________________________
1283/*TObjArray* AliAnalyseUE::SortChargedParticlesMC()
1284{
1285 // return an array with all charged particles ordered according to their pT .
1286 AliStack *mcStack = fkMC->Stack();
1287 Int_t nTracks = mcStack->GetNtrack();
1288 if( !nTracks ) return 0;
1289 TObjArray* tracks = new TObjArray(nTracks);
1290
1291 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
1292 if (!(mcStack->IsPhysicalPrimary(ipart))) continue;
1293 TParticle *part = mcStack->Particle(ipart);
1294
1295 if( !part->GetPDG()->Charge() ) continue;
1296 if( part->Pt() < fTrackPtCut ) continue;
1297 tracks->AddLast( part );
1298 }
1299 QSortTracksMC( *tracks, 0, tracks->GetEntriesFast() );
1300
1301 nTracks = tracks->GetEntriesFast();
1302 if( !nTracks ) return 0;
1303
1304 return tracks;
1305 }
1306*/
1307
1308//____________________________________________________________________
1309TObjArray* AliAnalyseUE::SortChargedParticlesMC()
1310{
1311 // return an array with all charged particles ordered according to their pT .
1312 TClonesArray *tca = dynamic_cast<TClonesArray*>(fkAOD->FindListObject(AliAODMCParticle::StdBranchName()));
1313 if(!tca){
1314 Printf("No branch!!!");
1315 return 0;
1316 }
1317 Int_t nTracks = tca->GetEntriesFast();
1318 if( !nTracks ) return 0;
1319 TObjArray* tracks = new TObjArray(nTracks);
1320
1321 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
1322 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(ipart));
1323 if(!part->IsPhysicalPrimary())continue;
1324 if( part->Pt() < fTrackPtCut ) continue;
1325 if(!part->Charge()) continue;
1326 tracks->AddLast( part );
1327 }
1328 QSortTracksMC( *tracks, 0, tracks->GetEntriesFast() );
1329
1330 nTracks = tracks->GetEntriesFast();
1331 if( !nTracks ) return 0;
1332
1333 return tracks;
1334 }
1335
1f329128 1336
1337//____________________________________________________________________
19699b9c 1338Bool_t AliAnalyseUE::TrackMCSelected(Double_t charge, Double_t pT, Double_t eta, Int_t pdgCode)const{
1f329128 1339
1340 // MC track selection
1341 Double_t const kMyTolerance = 0.0000001;
1342 //if (charge == 0. || charge == -99.) return kFALSE;
1343 if (charge < kMyTolerance || charge + 99 < kMyTolerance) return kFALSE;
1344
1345 if ( fUseSingleCharge ) { // Charge selection
1346 if ( fUsePositiveCharge && charge < 0.) return kFALSE; // keep Positives
1347 if ( !fUsePositiveCharge && charge > 0.) return kFALSE; // keep Negatives
1348 }
1349
1350 //Kinematics cuts on particle
1351 if ((pT < fTrackPtCut) || (TMath::Abs(eta) > fTrackEtaCut )) return kFALSE;
1352
1353 Bool_t isHadron = TMath::Abs(pdgCode)==211 ||
1354 TMath::Abs(pdgCode)==2212 ||
1355 TMath::Abs(pdgCode)==321;
1356
1357 if ( fUseChargeHadrons && !isHadron ) return kFALSE;
1358
1359 return kTRUE;
1360
1361}
1362
1363//____________________________________________________________________
19699b9c 1364Bool_t AliAnalyseUE::TrackSelected(AliAODTrack* part)const{
1f329128 1365
1366 // Real track selection
1367 if ( !part->TestFilterBit(fFilterBit) ) return kFALSE; // track cut selection
1368 if ( !part->IsPrimaryCandidate()) return kFALSE; // reject whatever is not linked to collision point
1369 // PID Selection: Reject everything but hadrons
1370 Bool_t isHadron = part->GetMostProbablePID()==AliAODTrack::kPion ||
1371 part->GetMostProbablePID()==AliAODTrack::kKaon ||
1372 part->GetMostProbablePID()==AliAODTrack::kProton;
1373 if ( fUseChargeHadrons && !isHadron ) return kFALSE;
1374
1375 if ( !part->Charge() ) return kFALSE; //Only charged
1376 if ( fUseSingleCharge ) { // Charge selection
1377 if ( fUsePositiveCharge && part->Charge() < 0.) return kFALSE; // keep Positives
1378 if ( !fUsePositiveCharge && part->Charge() > 0.) return kFALSE; // keep Negatives
1379 }
1380
1381 if ( part->Pt() < fTrackPtCut ) return kFALSE;
1382 if( TMath::Abs(part->Eta()) > fTrackEtaCut ) return kFALSE;
1383
1384 return kTRUE;
1385}
1386
3a9d4bcf 1387//____________________________________________________________________
1388Bool_t AliAnalyseUE::TrackSelectedEfficiency(AliAODTrack* part)const{
1389
1390 if (!fStack) AliWarning("Attention: stack not set from mother task");
1391 // Real track selection
1392 if ( !part->TestFilterBit(fFilterBit) ) return kFALSE; // track cut selection
1393 if ( !part->IsPrimaryCandidate()) return kFALSE; // reject whatever is not linked to collision point
1394 // PID Selection: Reject everything but hadrons
1395 Bool_t isHadron = part->GetMostProbablePID()==AliAODTrack::kPion ||
1396 part->GetMostProbablePID()==AliAODTrack::kKaon ||
1397 part->GetMostProbablePID()==AliAODTrack::kProton;
1398 if ( fUseChargeHadrons && !isHadron ) return kFALSE;
1399
1400 if ( !part->Charge() ) return kFALSE; //Only charged
1401 if ( fUseSingleCharge ) { // Charge selection
1402 if ( fUsePositiveCharge && part->Charge() < 0.) return kFALSE; // keep Positives
1403 if ( !fUsePositiveCharge && part->Charge() > 0.) return kFALSE; // keep Negatives
1404 }
1405
1406 /*if ( part->Pt() < fTrackPtCut ) return kFALSE;
1407 if( TMath::Abs(part->Eta()) > fTrackEtaCut ) return kFALSE;*/
1408
1409 //Check if there was a primary fulfilling the required cuts
1410 Double_t charge=-999.;
1411 Double_t pt=-999.;
1412 Double_t eta=-999.;
1413 Int_t pdgcode=-999;
1414
1415 Int_t label = part->GetLabel();
1416 TClonesArray *tca=0;
1417 tca = dynamic_cast<TClonesArray*>(fkAOD->FindListObject(AliAODMCParticle::StdBranchName()));
1418 //TParticle *partMC = fStack->ParticleFromTreeK(label);
1419 AliAODMCParticle *partMC=dynamic_cast<AliAODMCParticle*>(tca->At(TMath::Abs(label)));
1420 if (!partMC->IsPhysicalPrimary())return kFALSE;
1421 //charge = ((TParticlePDG*)partMC->GetPDG())->Charge();
1422 charge = partMC->Charge();
1423 pt = partMC->Pt();
1424 eta = partMC->Eta();
1425 pdgcode = partMC->GetPdgCode();
1426
1427 if (!TrackMCSelected(charge, pt, eta, pdgcode)) return kFALSE;
1428 return kTRUE;
1429}
1430
1f329128 1431