]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/AliAnalyseUE.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Transition PWG0 -> PWGUD
[u/mrichter/AliRoot.git] / PWGJE / 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;
b85b355c 549 if(!partMC)return;
550 charge = partMC->Charge();
3a9d4bcf 551 pt = partMC->Pt();
552 eta = partMC->Eta();
553 pdgcode = partMC->GetPdgCode();
554 suffix.Append("MC");
b85b355c 555 }
556
557 if(!part)return;
3a9d4bcf 558
1f329128 559 if (fDebug > 1) AliInfo(Form(" ==== AOD track = %d pt = %f charge = %d \n ",ipart,part->Pt(),part->Charge()));
3a9d4bcf 560
1f329128 561 // track selection
3a9d4bcf 562 if (!mctrue && !eff ){
563 if( !TrackSelected(partRECO)) continue; //track selection for data and MC reconstructed
564 if (flag_tmp){
565 if (fkESD && fkESD->GetNumberOfTracks() ){
566 AliInfo("READING ESD *************************************************");
567 Int_t id = partRECO->GetID();
568 AliESDtrack *trackESD;
569 trackESD = (AliESDtrack*)fkESD->GetTrack(id);
570 Float_t d0;
571 Float_t z;
572 trackESD->GetImpactParameters(d0,z);
573 fHistos->FillHistogram("hDCAxy", d0, jetVect[0].Pt());
13242232 574 }else AliInfo("NO TRACKS ************************************************") ;
575 }
3a9d4bcf 576 }
577
578 if (!mctrue && eff){
579 if (!TrackSelectedEfficiency(partRECO )) continue; //track selection for MC reconstructed for efficiency studies
580 if(fkESD && fkESD->GetNumberOfTracks()){
581 Int_t id = partRECO->GetID();
582 AliESDtrack * trackESD = (AliESDtrack*) fkESD->GetTrack(id);
583 Float_t d0;
584 Float_t z;
585 trackESD->GetImpactParameters(d0,z);
586 fHistos->FillHistogram("hDCAxyPrimary", d0, jetVect[0].Pt());
587 }
588 }
589
590 if (mctrue){
591 if ( !(TrackMCSelected(charge, pt, eta, pdgcode) && partMC->IsPhysicalPrimary())) continue; //track selection for MC true
592 }
1f329128 593
594 TVector3 partVect(part->Px(), part->Py(), part->Pz());
595 Bool_t isFlagPart = kTRUE;
596 Double_t deltaPhi = jetVect[0].DeltaPhi(partVect)+k270rad;
597 if( deltaPhi > kTWOPI ) deltaPhi-= kTWOPI;
3a9d4bcf 598 if (fAnaType != 4 ) fHistos->FillHistogram(Form("hdNdEtaPhiDist%s",suffix.Data()),deltaPhi, jetVect[0].Pt());
1f329128 599 else if (TMath::Abs(deltaPhi-k270rad) >= kMyTolerance && TMath::Abs(jetVect[0].Eta()-partVect.Eta()) >= kMyTolerance){
3a9d4bcf 600 fHistos->FillHistogram(Form("hdNdEtaPhiDist%s",suffix.Data()),deltaPhi, jetVect[0].Pt());
1f329128 601 isFlagPart = kFALSE;
602 }
603
3a9d4bcf 604 fHistos->FillHistogram(Form("hFullRegPartPtDistVsEt%s",suffix.Data()), part->Pt(), jetVect[0].Pt());
1f329128 605
606 Int_t region = IsTrackInsideRegion( jetVect, &partVect, conePosition );
607 if (region == 1) {
608 if( fMaxPartPtRegion < part->Pt() ) fMaxPartPtRegion = part->Pt();
609 fSumPtRegionPosit += part->Pt();
610 fNTrackRegionPosit++;
3a9d4bcf 611 fHistos->FillHistogram(Form("hTransRegPartPtDistVsEt%s",suffix.Data()),part->Pt(), jetVect[0].Pt());
1f329128 612 }
613 if (region == -1) {
614 if( fMaxPartPtRegion < part->Pt() ) fMaxPartPtRegion = part->Pt();
615 fSumPtRegionNegat += part->Pt();
616 fNTrackRegionNegat++;
3a9d4bcf 617 fHistos->FillHistogram(Form("hTransRegPartPtDistVsEt%s",suffix.Data()),part->Pt(), jetVect[0].Pt());
1f329128 618 }
619 if (region == 2){ //forward
620 fSumPtRegionForward += part->Pt();
621 fNTrackRegionForward++;
3a9d4bcf 622 fHistos->FillHistogram(Form("hRegForwardPartPtDistVsEt%s",suffix.Data()),part->Pt(), jetVect[0].Pt());
1f329128 623 }
624 if (region == -2){ //backward
625 fSumPtRegionBackward += part->Pt();
626 fNTrackRegionBackward++;
3a9d4bcf 627 fHistos->FillHistogram(Form("hRegBackwardPartPtDistVsEt%s",suffix.Data()),part->Pt(), jetVect[0].Pt());
1f329128 628 }
629 }//end loop AOD tracks
630
631}
632
3a9d4bcf 633//-------------------------------------------------------------------
634/*TVector3 AliAnalyseUE::GetLeadingTracksMC(AliMCEvent *mcEvent){
635
636 fkMC = mcEvent;
637
638 Double_t maxPtJet1 = 0.;
639 Int_t index1 = -1;
640
641 TVector3 jetVect[3];
642
643 jetVect[0].SetPtEtaPhi(-1.,-1.,-1.);
644 jetVect[1].SetPtEtaPhi(-1.,-1.,-1.);
645 jetVect[2].SetPtEtaPhi(-1.,-1.,-1.);
646
647 Int_t nJets = 0;
648
649 TObjArray* tracks = SortChargedParticlesMC();
650 if( tracks ) {
651 nJets = tracks->GetEntriesFast();
652 if( nJets > 0 ) {
653 index1 = 0;
654 TParticle* jet = (TParticle*)tracks->At(0);
655 maxPtJet1 = jet->Pt();
656 jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
657 }
658 tracks->Clear();
659 delete tracks;
660 }
661 return *jetVect;
662
663}
664*/
665
666//-------------------------------------------------------------------
667TVector3 AliAnalyseUE::GetLeadingTracksMC(AliMCEvent *mcEvent){
668
669 fkMC = mcEvent;
670
671 Double_t maxPtJet1 = 0.;
672 Int_t index1 = -1;
673
674 TVector3 jetVect[3];
675
676 jetVect[0].SetPtEtaPhi(-1.,-1.,-1.);
677 jetVect[1].SetPtEtaPhi(-1.,-1.,-1.);
678 jetVect[2].SetPtEtaPhi(-1.,-1.,-1.);
679
680 Int_t nJets = 0;
681
682 TObjArray* tracks = SortChargedParticlesMC();
683 if( tracks ) {
684 nJets = tracks->GetEntriesFast();
685 if( nJets > 0 ) {
686 index1 = 0;
687 AliAODMCParticle* jet = (AliAODMCParticle*)tracks->At(0);
688 fLtMCLabel = jet->GetLabel();
689 maxPtJet1 = jet->Pt();
690 jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
691 }
692 tracks->Clear();
693 delete tracks;
694 }
695 return *jetVect;
696
697}
698
1f329128 699//-------------------------------------------------------------------
700TVector3 AliAnalyseUE::GetOrderedClusters(TString aodBranch, Bool_t chargedJets, Double_t chJetPtMin){
701
702 // jets from AOD, on-the-fly or leading particle
703 Double_t maxPtJet1 = 0.;
704 Int_t index1 = -1;
705 Double_t maxPtJet2 = 0.; // jet 2 need for back to back inclusive
706 Int_t index2 = -1;
707 Double_t maxPtJet3 = 0.; // jet 3 need for back to back exclusive
708 Int_t index3 = -1;
709 TVector3 jetVect[3];
710
711 jetVect[0].SetPtEtaPhi(-1.,-1.,-1.);
712 jetVect[1].SetPtEtaPhi(-1.,-1.,-1.);
713 jetVect[2].SetPtEtaPhi(-1.,-1.,-1.);
714
715 Int_t nJets = 0;
716 //TClonesArray* fArrayJets;
717 TObjArray* arrayJets;
718 // 1) JETS FROM AOD BRANCH (standard, non-standard or delta)
719 if (!chargedJets && fAnaType != 4 ) {
720 AliInfo(" ==== Read AODs !");
721 AliInfo(Form(" ==== Reading Branch: %s ", aodBranch.Data()));
722 arrayJets = (TObjArray*)fkAOD->GetList()->FindObject(aodBranch.Data());
723 if (!arrayJets){
724 AliFatal(" No jet-array! ");
725 return *jetVect;
726 }
727
728 // Find Leading Jets 1,2,3
729 // (could be skipped if Jets are sort by Pt...)
730 nJets=arrayJets->GetEntries();
731 for( Int_t i=0; i<nJets; ++i ) {
732 AliAODJet* jet = (AliAODJet*)arrayJets->At(i);
733 Double_t jetPt = jet->Pt();//*1.666; // FIXME Jet Pt Correction ?????!!!
734
735 if( jetPt > maxPtJet1 ) {
736 maxPtJet3 = maxPtJet2; index3 = index2;
737 maxPtJet2 = maxPtJet1; index2 = index1;
738 maxPtJet1 = jetPt; index1 = i;
739 } else if( jetPt > maxPtJet2 ) {
740 maxPtJet3 = maxPtJet2; index3 = index2;
741 maxPtJet2 = jetPt; index2 = i;
742 } else if( jetPt > maxPtJet3 ) {
743 maxPtJet3 = jetPt; index3 = i;
744 }
745 }
746
747 if( index1 != -1 ) {
748 AliAODJet *jet =(AliAODJet*) arrayJets->At(index1);
749 if(jet)jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
750 }
751 if( index2 != -1 ) {
752 AliAODJet* jet= (AliAODJet*) arrayJets->At(index2);
753 if(jet)jetVect[1].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
754 }
755 if( index3 != -1 ) {
756 AliAODJet* jet = (AliAODJet*) arrayJets->At(index3);
757 if(jet)jetVect[2].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
758 }
759
760 }
761
762
763 // 2) ON-THE-FLY CDF ALGORITHM
764 if (chargedJets){
1f329128 765 arrayJets = FindChargedParticleJets(chJetPtMin);
766 if( arrayJets ) {
767 nJets = arrayJets->GetEntriesFast();
768 if( nJets > 0 ) {
769 index1 = 0;
770 AliAODJet* jet = (AliAODJet*)arrayJets->At(0);
771 maxPtJet1 = jet->Pt();
772 jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
773 }
774 if( nJets > 1 ) {
775 index2 = 1;
776 AliAODJet* jet = (AliAODJet*)arrayJets->At(1);
777 maxPtJet2 = jet->Pt();
778 jetVect[1].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
779 }
780 if( nJets > 2 ) {
781 index3 = 2;
782 AliAODJet* jet = (AliAODJet*)arrayJets->At(2);
783 maxPtJet3 = jet->Pt();
784 jetVect[2].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
785 }
786
787 arrayJets->Delete();
788 delete arrayJets;
789 }
790 }
791
792
793 // 3) LEADING PARTICLE
794 if( fAnaType == 4 ){
795 TObjArray* tracks = SortChargedParticles();
796 if( tracks ) {
797 nJets = tracks->GetEntriesFast();
798 if( nJets > 0 ) {
799 index1 = 0;
800 AliAODTrack* jet = (AliAODTrack*)tracks->At(0);
3a9d4bcf 801 fLtLabel = jet->GetLabel();
802 maxPtJet1 = jet->Pt();
1f329128 803 jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
804 }
805 tracks->Clear();
806 delete tracks;
807 }
808
809 }
ddbad1f5 810 fHistos->FillHistogram("hNJets",nJets);
1f329128 811
812 return *jetVect;
813
814}
815
816
817//-------------------------------------------------------------------
818void AliAnalyseUE::Initialize(AliAnalysisTaskUE& taskUE){
3a9d4bcf 819/*void AliAnalyseUE::Initialize(AliAnalysisTask& task){// when correction task is in trunk
820 if (task->InheritsFrom("AliAnalysisTaskUE")){
821 AliAnalysisTaskUE *taskUE = dynamic_cast<AliAnalysisTaskUE*> task;
822 }
823 else if (task->InheritsFrom("AliAnalysisTaskCorrectionsUE")){
824 AliAnalysisTaskCorrectionsUE *taskUE = dynamic_cast<AliAnalysisTaskCorrectionsUE*> task;
825 }
826
827*/
1f329128 828 //Get principal settings from current instance of UE analysis-task
829 fAnaType = taskUE.GetAnaTopology();
830 fkAOD = taskUE.GetAOD();
831 fConeRadius = taskUE.GetConeRadius();
832 fDebug = taskUE.GetDebugLevel();
833 fFilterBit = taskUE.GetFilterBit();
834 fJet1EtaCut = taskUE.GetJet1EtaCut();
835 fJet2DeltaPhiCut = taskUE.GetJet2DeltaPhiCut();
836 fJet2RatioPtCut = taskUE.GetJet2RatioPtCut();
837 fJet3PtCut = taskUE.GetJet3PtCut();
838 fOrdering = taskUE.GetPtSumOrdering() ;
839 fRegionType = taskUE.GetRegionType();
840 fSimulateChJetPt = taskUE.GetSimulateChJetPt();
841 fTrackEtaCut = taskUE.GetTrackEtaCut();
842 fTrackPtCut = taskUE.GetTrackPtCut();
ddbad1f5 843 fHistos = taskUE.fHistosUE;
1f329128 844 fUseChargeHadrons = taskUE.GetUseChargeHadrons();
845 fUseChPartJet = taskUE.GetUseChPartJet();
846 fUsePositiveCharge = taskUE.GetUseNegativeChargeType();
847 fUseSingleCharge = taskUE.GetUseSingleCharge();
848
1f329128 849}
850
3a9d4bcf 851//-------------------------------------------------------------------
852void 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){
853
854 fAnaType = anaType;
855 fkAOD = aod;
856 fConeRadius = coneRadius;
857 fDebug = debug;
858 fFilterBit = filterBit;
859 fJet1EtaCut = jet1EtaCut;
860 fJet2DeltaPhiCut = jet2DeltaPhiCut;
861 fJet2RatioPtCut = jet2RatioPtCut;
862 fJet3PtCut = jet3PtCut;
863 fOrdering = ordering;
864 fRegionType = regionType;
865 fSimulateChJetPt = simulateChJetPt;
866 fTrackEtaCut = trackEtaCut;
867 fTrackPtCut = trackPtCut;
868 fUseChargeHadrons = useChargeHadrons;
869 fUseChPartJet = useChPartJet;
870 fUsePositiveCharge = usePositiveCharge;
871 fUseSingleCharge = useSingleCharge;
872 fHistos = histos;
873}
874
875
876//-------------------------------------------------------------------
877Bool_t AliAnalyseUE::TriggerSelection(AliInputEventHandler* inputHandler){
878
879 //Use AliPhysicsSelection to select good events
880 if( !inputHandler->InheritsFrom("AliAODInputHandler") ) { // input AOD
881 if (inputHandler->IsEventSelected()) {
882 if (fDebug > 1) AliInfo(" Trigger Selection: event ACCEPTED ... ");
883 } else {
884 if (fDebug > 1) AliInfo(" Trigger Selection: event REJECTED ... ");
885 return kFALSE;
886 }
887 }
888
889 return kTRUE;
890
891}
892
893
1f329128 894//-------------------------------------------------------------------
1f329128 895Bool_t AliAnalyseUE::VertexSelection(AliAODEvent *aod, Int_t tracks, Double_t zed ){
896
897 //Require 1 vertex (no TPC stand-alone) with a minimum number of tracks and z-coordinate in a limited range
898 Int_t nVertex = aod->GetNumberOfVertices();
899 if( nVertex > 0 ) { // Only one vertex (reject pileup)
900 AliAODVertex* vertex = (AliAODVertex*)aod->GetPrimaryVertex();
901 Int_t nTracksPrim = vertex->GetNContributors();
902 Double_t zVertex = vertex->GetZ();
903 if (fDebug > 1) AliInfo(Form(" Vertex in = %f with %d particles by %s data ...",zVertex,nTracksPrim,vertex->GetName()));
904 // Select a quality vertex by number of tracks?
905 if( nTracksPrim < tracks || TMath::Abs(zVertex) > zed ) {
906 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
907 return kFALSE;
908 }
909 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
910 } else {
911 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
912 return kFALSE;
913 }
914
915 return kTRUE;
916}
917
ddbad1f5 918//-------------------------------------------------------------------
919Bool_t AliAnalyseUE::VertexSelectionOld(AliAODEvent *aod ){
920
921 AliKFVertex primVtx(*(aod->GetPrimaryVertex()));
922 Int_t nTracksPrim=primVtx.GetNContributors();
923 if (fDebug > 1) AliInfo(Form(" Primary-vertex Selection: %d",nTracksPrim));
924 if(!nTracksPrim){
925 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
926 return kFALSE;
927 }
928 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED ...");
929
930 return kTRUE;
931}
932
1f329128 933// PRIVATE METHODS **************************************************
934
935TObjArray* AliAnalyseUE::FindChargedParticleJets( Double_t chJetPtMin )
936{
937 // Return a TObjArray of "charged particle jets"
938
939 // Charged particle jet definition from reference:
940 // "Charged jet evolution and the underlying event
941 // in proton-antiproton collisions at 1.8 TeV"
942 // PHYSICAL REVIEW D 65 092002, CDF Collaboration
943
944 // We defined "jets" as circular regions in eta-phi space with
945 // radius defined by R = sqrt( (eta-eta0)^2 +(phi-phi0)^2 ).
946 // Our jet algorithm is as follows:
947 // 1- Order all charged particles according to their pT .
948 // 2- Start with the highest pT particle and include in the jet all
949 // particles within the radius R=0.7 considering each particle
950 // in the order of decreasing pT and recalculating the centroid
951 // of the jet after each new particle is added to the jet .
952 // 3- Go to the next highest pT particle not already included in
953 // a jet and add to the jet all particles not already included in
954 // a jet within R=0.7.
955 // 4- Continue until all particles are in a jet.
956 // We defined the transverse momentum of the jet to be
957 // the scalar pT sum of all the particles within the jet, where pT
958 // is measured with respect to the beam axis
959
960 // 1 - Order all charged particles according to their pT .
961 Int_t nTracks = fkAOD->GetNTracks();
962 if( !nTracks ) return 0;
963 TObjArray tracks(nTracks);
964
965 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
966 AliAODTrack* part = fkAOD->GetTrack( ipart );
967 if( !part->TestFilterBit(fFilterBit) ) continue; // track cut selection
968 if( !part->Charge() ) continue;
969 if( part->Pt() < fTrackPtCut ) continue;
970 tracks.AddLast(part);
971 }
972 QSortTracks( tracks, 0, tracks.GetEntriesFast() );
973
974 nTracks = tracks.GetEntriesFast();
975 if( !nTracks ) return 0;
976
977 TObjArray *jets = new TObjArray(nTracks);
978 TIter itrack(&tracks);
979 while( nTracks ) {
980 // 2- Start with the highest pT particle ...
981 Float_t px,py,pz,pt;
982 AliAODTrack* track = (AliAODTrack*)itrack.Next();
983 if( !track ) continue;
984 px = track->Px();
985 py = track->Py();
986 pz = track->Pz();
987 pt = track->Pt(); // Use the energy member to store Pt
988 jets->AddLast( new TLorentzVector(px, py, pz, pt) );
989 tracks.Remove( track );
990 TLorentzVector* jet = (TLorentzVector*)jets->Last();
991 jet->SetPtEtaPhiE( 1., jet->Eta(), jet->Phi(), pt );
992 // 3- Go to the next highest pT particle not already included...
993 AliAODTrack* track1;
994 Double_t fPt = jet->E();
995 while ( (track1 = (AliAODTrack*)(itrack.Next())) ) {
996 Double_t tphi = track1->Phi(); // here Phi is from 0 <-> 2Pi
997 if (tphi > TMath::Pi()) tphi -= 2. * TMath::Pi(); // convert to -Pi <-> Pi
998 Double_t dphi = TVector2::Phi_mpi_pi(jet->Phi()-tphi);
999 Double_t r = TMath::Sqrt( (jet->Eta()-track1->Eta())*(jet->Eta()-track1->Eta()) +dphi*dphi );
1000 if( r < fConeRadius ) {
1001 fPt = jet->E()+track1->Pt(); // Scalar sum of Pt
1002 // recalculating the centroid
1003 Double_t eta = jet->Eta()*jet->E()/fPt + track1->Eta()*track1->Pt()/fPt;
1004 Double_t phi = jet->Phi()*jet->E()/fPt + tphi*track1->Pt()/fPt;
1005 jet->SetPtEtaPhiE( 1., eta, phi, fPt );
1006 tracks.Remove( track1 );
1007 }
1008 }
1009
1010 tracks.Compress();
1011 nTracks = tracks.GetEntries();
1012 // 4- Continue until all particles are in a jet.
1013 itrack.Reset();
1014 } // end while nTracks
1015
1016 // Convert to AODjets....
1017 Int_t njets = jets->GetEntriesFast();
1018 TObjArray* aodjets = new TObjArray(njets);
1019 aodjets->SetOwner(kTRUE);
1020 for(Int_t ijet=0; ijet<njets; ++ijet) {
1021 TLorentzVector* jet = (TLorentzVector*)jets->At(ijet);
1022 if (jet->E() < chJetPtMin) continue;
1023 Float_t px, py,pz,en; // convert to 4-vector
1024 px = jet->E() * TMath::Cos(jet->Phi()); // Pt * cos(phi)
1025 py = jet->E() * TMath::Sin(jet->Phi()); // Pt * sin(phi)
1026 pz = jet->E() / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-jet->Eta())));
1027 en = TMath::Sqrt(px * px + py * py + pz * pz);
1028
1029 aodjets->AddLast( new AliAODJet(px, py, pz, en) );
1030 }
1031 jets->Delete();
1032 delete jets;
1033
1034 // Order jets according to their pT .
1035 QSortTracks( *aodjets, 0, aodjets->GetEntriesFast() );
1036
1037 // debug
1038 if (fDebug>3) AliInfo(Form(" %d Charged jets found\n",njets));
1039
1040 return aodjets;
1041}
1042
1043
1044//____________________________________________________________________
1045void AliAnalyseUE::FillAvePartPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin )
1046{
1047
1048 // Fill average particle Pt of control regions
1049
1050 // Max cone
1051 fHistos->FillHistogram("hRegionAvePartPtMaxVsEt", leadingE, ptMax);
1052 // Min cone
1053 fHistos->FillHistogram("hRegionAvePartPtMinVsEt", leadingE, ptMin);
1054 // MAke distributions for UE comparison with MB data
1055 fHistos->FillHistogram("hMinRegAvePt", ptMin);
1056
1057}
1058
1059//____________________________________________________________________
1060void AliAnalyseUE::FillMultRegion( Double_t leadingE, Double_t nTrackPtmax, Double_t nTrackPtmin, Double_t ptMin )
1061{
1062
1063 // Fill Nch multiplicity of control regions
1064
1065 // Max cone
1066 fHistos->FillHistogram("hRegionMultMaxVsEt", leadingE, nTrackPtmax);
1067 fHistos->FillHistogram("hRegionMultMax", nTrackPtmax);
1068
1069 // Min cone
1070 fHistos->FillHistogram("hRegionMultMinVsEt", leadingE, nTrackPtmin );
1071 fHistos->FillHistogram("hRegionMultMin", nTrackPtmin);
1072
1073 // MAke distributions for UE comparison with MB data
1074 fHistos->FillHistogram("hMinRegSumPtvsMult", nTrackPtmin,ptMin);
1075
1076}
1077
1078//____________________________________________________________________
1079void AliAnalyseUE::FillSumPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin )
1080{
1081 // Fill sumPt of control regions
1082
1083 // Max cone
1084 fHistos->FillHistogram("hRegionSumPtMaxVsEt", leadingE, ptMax);
1085
1086 // Min cone
1087 fHistos->FillHistogram("hRegionSumPtMinVsEt", leadingE, ptMin);
1088
1089 // MAke distributions for UE comparison with MB data
1090 fHistos->FillHistogram("hMinRegSumPt", ptMin);
1091 fHistos->FillHistogram("hMinRegSumPtJetPtBin", leadingE, ptMin);
1092 fHistos->FillHistogram("hMaxRegSumPtJetPtBin", leadingE, ptMax);
1093
1094}
1095
1096//-------------------------------------------------------------------
1097Int_t AliAnalyseUE::IsTrackInsideRegion(TVector3 *jetVect, TVector3 *partVect, Int_t conePosition)
1098{
1099 // return de region in delta phi
1100 // -1 negative delta phi
1101 // 1 positive delta phi
1102 // 0 outside region
1103 static const Double_t k60rad = 60.*TMath::Pi()/180.;
1104 static const Double_t k120rad = 120.*TMath::Pi()/180.;
1105
1106 Int_t region = 0;
1107 if( fRegionType == 1 ) {
1108 if( TMath::Abs(partVect->Eta()) > fTrackEtaCut ) return 0;
1109 // transverse regions
1110 if (jetVect[0].DeltaPhi(*partVect) < -k60rad && jetVect[0].DeltaPhi(*partVect) > -k120rad ) region = -1; //left
1111 if (jetVect[0].DeltaPhi(*partVect) > k60rad && jetVect[0].DeltaPhi(*partVect) < k120rad ) region = 1; //right
1112 if (TMath::Abs(jetVect[0].DeltaPhi(*partVect)) < k60rad ) region = 2; //forward
1113 if (TMath::Abs(jetVect[0].DeltaPhi(*partVect)) > k120rad ) region = -2; //backward
1114
1115 } else if( fRegionType == 2 ) {
1116 // Cone regions
1117 Double_t deltaR = 0.;
1118
1119 TVector3 positVect,negatVect;
1120 if (conePosition==1){
1121 positVect.SetMagThetaPhi(1, 2.*atan(exp(-jetVect[0].Eta())), jetVect[0].Phi()+TMath::PiOver2());
1122 negatVect.SetMagThetaPhi(1, 2.*atan(exp(-jetVect[0].Eta())), jetVect[0].Phi()-TMath::PiOver2());
1123 }else if (conePosition==2){
1124 if(fAnaType<2) AliFatal("Prevent error in Analysis type there might be only 1 jet. To avoid overflow better Correct UE config");
1125 positVect.SetMagThetaPhi(1, 2.*atan(exp(-(jetVect[0].Eta()+jetVect[1].Eta())/2.)), jetVect[0].Phi()+TMath::PiOver2());
1126 negatVect.SetMagThetaPhi(1, 2.*atan(exp(-(jetVect[0].Eta()+jetVect[1].Eta())/2.)), jetVect[0].Phi()-TMath::PiOver2());
1127 }else if (conePosition==3){
1128 if(fAnaType<2) AliFatal("Prevent error in Analysis type there might be only 1 jet. To avoid overflow better Correct UE config");
1129 Double_t weightEta = jetVect[0].Eta() * jetVect[0].Pt()/(jetVect[0].Pt() + jetVect[1].Pt()) +
1130 jetVect[1].Eta() * jetVect[1].Pt()/(jetVect[0].Pt() + jetVect[1].Pt());
1131 //Double_t weightEta = jetVect[0].Eta() * jetVect[0].Mag()/(jetVect[0].Mag() + jetVect[1].Mag()) +
1132 // jetVect[1].Eta() * jetVect[1].Mag()/(jetVect[0].Mag() + jetVect[1].Mag());
1133 positVect.SetMagThetaPhi(1, 2.*atan(exp(-weightEta)), jetVect[0].Phi()+TMath::PiOver2());
1134 negatVect.SetMagThetaPhi(1, 2.*atan(exp(-weightEta)), jetVect[0].Phi()-TMath::PiOver2());
1135 }
1136 if (TMath::Abs(positVect.DeltaPhi(*partVect)) < fConeRadius ) {
1137 region = 1;
1138 deltaR = positVect.DrEtaPhi(*partVect);
1139 } else if (TMath::Abs(negatVect.DeltaPhi(*partVect)) < fConeRadius) {
1140 region = -1;
1141 deltaR = negatVect.DrEtaPhi(*partVect);
1142 }
1143
1144 if (deltaR > fConeRadius) region = 0;
1145
1146 } else
1147 AliError("Unknow region type");
1148
1149 return region;
1150 }
1151
1152
1153
1154//-------------------------------------------------------------------
1155void AliAnalyseUE::QSortTracks(TObjArray &a, Int_t first, Int_t last)
1156{
1157 // Sort array of TObjArray of tracks by Pt using a quicksort algorithm.
1158
1159 static TObject *tmp;
1160 static int i; // "static" to save stack space
1161 int j;
1162
1163 while (last - first > 1) {
1164 i = first;
1165 j = last;
1166 for (;;) {
1167 while (++i < last && ((AliVParticle*)a[i])->Pt() > ((AliVParticle*)a[first])->Pt() )
1168 ;
1169 while (--j > first && ((AliVParticle*)a[j])->Pt() < ((AliVParticle*)a[first])->Pt() )
1170 ;
1171 if (i >= j)
1172 break;
1173
1174 tmp = a[i];
1175 a[i] = a[j];
1176 a[j] = tmp;
1177 }
1178 if (j == first) {
1179 ++first;
1180 continue;
1181 }
1182 tmp = a[first];
1183 a[first] = a[j];
1184 a[j] = tmp;
1185 if (j - first < last - (j + 1)) {
1186 QSortTracks(a, first, j);
1187 first = j + 1; // QSortTracks(j + 1, last);
1188 } else {
1189 QSortTracks(a, j + 1, last);
1190 last = j; // QSortTracks(first, j);
1191 }
1192 }
1193}
1194
1195
3a9d4bcf 1196//-------------------------------------------------------------------
1197void AliAnalyseUE::QSortTracksMC(TObjArray &a, Int_t first, Int_t last)
1198{
1199 // Sort array of TObjArray of tracks by Pt using a quicksort algorithm.
1200
1201 static TObject *tmp;
1202 static int i; // "static" to save stack space
1203 int j;
1204
1205 while (last - first > 1) {
1206 i = first;
1207 j = last;
1208 for (;;) {
1209 while (++i < last && ((AliAODMCParticle*)a[i])->Pt() > ((AliAODMCParticle*)a[first])->Pt() )
1210 ;
1211 while (--j > first && ((AliAODMCParticle*)a[j])->Pt() < ((AliAODMCParticle*)a[first])->Pt() )
1212 ;
1213 if (i >= j)
1214 break;
1215
1216 tmp = a[i];
1217 a[i] = a[j];
1218 a[j] = tmp;
1219 }
1220 if (j == first) {
1221 ++first;
1222 continue;
1223 }
1224 tmp = a[first];
1225 a[first] = a[j];
1226 a[j] = tmp;
1227 if (j - first < last - (j + 1)) {
1228 QSortTracksMC(a, first, j);
1229 first = j + 1; // QSortTracks(j + 1, last);
1230 } else {
1231 QSortTracksMC(a, j + 1, last);
1232 last = j; // QSortTracks(first, j);
1233 }
1234 }
1235}
1236
1237
1238
1f329128 1239
1240
1241//-------------------------------------------------------------------
1242void AliAnalyseUE::SetRegionArea(TVector3 *jetVect)
1243{
1244 // Set region area
1245 Double_t areaCorrFactor=0.;
1246 Double_t deltaEta = 0.;
1247 if (fRegionType==1) fAreaReg = 2.*fTrackEtaCut*60.*TMath::Pi()/180.;
1248 else if (fRegionType==2){
1249 deltaEta = 0.9-TMath::Abs(jetVect[0].Eta());
1250 if (deltaEta>fConeRadius) fAreaReg = TMath::Pi()*fConeRadius*fConeRadius;
1251 else{
1252 areaCorrFactor = fConeRadius*fConeRadius*TMath::ACos(deltaEta/fConeRadius) -
1253 (deltaEta*fConeRadius)*TMath::Sqrt( 1. - deltaEta*deltaEta/(fConeRadius*fConeRadius));
1254 fAreaReg=TMath::Pi()*fConeRadius*fConeRadius-areaCorrFactor;
1255 }
1256 }else AliWarning("Unknown Region Type");
1257 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));
1258}
1259
1260
1261//____________________________________________________________________
1262TObjArray* AliAnalyseUE::SortChargedParticles()
1263{
1264 // return an array with all charged particles ordered according to their pT .
1265 Int_t nTracks = fkAOD->GetNTracks();
1266 if( !nTracks ) return 0;
1267 TObjArray* tracks = new TObjArray(nTracks);
1268
1269 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
1270 AliAODTrack* part = fkAOD->GetTrack( ipart );
1271 if( !part->TestFilterBit( fFilterBit ) ) continue; // track cut selection
1272 if( !part->Charge() ) continue;
1273 if( part->Pt() < fTrackPtCut ) continue;
1274 tracks->AddLast( part );
1275 }
1276 QSortTracks( *tracks, 0, tracks->GetEntriesFast() );
1277
1278 nTracks = tracks->GetEntriesFast();
1279 if( !nTracks ) return 0;
1280
1281 return tracks;
1282 }
1283
3a9d4bcf 1284//____________________________________________________________________
1285/*TObjArray* AliAnalyseUE::SortChargedParticlesMC()
1286{
1287 // return an array with all charged particles ordered according to their pT .
1288 AliStack *mcStack = fkMC->Stack();
1289 Int_t nTracks = mcStack->GetNtrack();
1290 if( !nTracks ) return 0;
1291 TObjArray* tracks = new TObjArray(nTracks);
1292
1293 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
1294 if (!(mcStack->IsPhysicalPrimary(ipart))) continue;
1295 TParticle *part = mcStack->Particle(ipart);
1296
1297 if( !part->GetPDG()->Charge() ) continue;
1298 if( part->Pt() < fTrackPtCut ) continue;
1299 tracks->AddLast( part );
1300 }
1301 QSortTracksMC( *tracks, 0, tracks->GetEntriesFast() );
1302
1303 nTracks = tracks->GetEntriesFast();
1304 if( !nTracks ) return 0;
1305
1306 return tracks;
1307 }
1308*/
1309
1310//____________________________________________________________________
1311TObjArray* AliAnalyseUE::SortChargedParticlesMC()
1312{
1313 // return an array with all charged particles ordered according to their pT .
1314 TClonesArray *tca = dynamic_cast<TClonesArray*>(fkAOD->FindListObject(AliAODMCParticle::StdBranchName()));
1315 if(!tca){
1316 Printf("No branch!!!");
1317 return 0;
1318 }
1319 Int_t nTracks = tca->GetEntriesFast();
1320 if( !nTracks ) return 0;
1321 TObjArray* tracks = new TObjArray(nTracks);
1322
1323 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
1324 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(ipart));
b85b355c 1325 if(!part)continue;
3a9d4bcf 1326 if(!part->IsPhysicalPrimary())continue;
1327 if( part->Pt() < fTrackPtCut ) continue;
1328 if(!part->Charge()) continue;
1329 tracks->AddLast( part );
1330 }
1331 QSortTracksMC( *tracks, 0, tracks->GetEntriesFast() );
1332
1333 nTracks = tracks->GetEntriesFast();
1334 if( !nTracks ) return 0;
1335
1336 return tracks;
1337 }
1338
1f329128 1339
1340//____________________________________________________________________
19699b9c 1341Bool_t AliAnalyseUE::TrackMCSelected(Double_t charge, Double_t pT, Double_t eta, Int_t pdgCode)const{
1f329128 1342
1343 // MC track selection
1344 Double_t const kMyTolerance = 0.0000001;
1345 //if (charge == 0. || charge == -99.) return kFALSE;
1346 if (charge < kMyTolerance || charge + 99 < kMyTolerance) return kFALSE;
1347
1348 if ( fUseSingleCharge ) { // Charge selection
1349 if ( fUsePositiveCharge && charge < 0.) return kFALSE; // keep Positives
1350 if ( !fUsePositiveCharge && charge > 0.) return kFALSE; // keep Negatives
1351 }
1352
1353 //Kinematics cuts on particle
1354 if ((pT < fTrackPtCut) || (TMath::Abs(eta) > fTrackEtaCut )) return kFALSE;
1355
1356 Bool_t isHadron = TMath::Abs(pdgCode)==211 ||
1357 TMath::Abs(pdgCode)==2212 ||
1358 TMath::Abs(pdgCode)==321;
1359
1360 if ( fUseChargeHadrons && !isHadron ) return kFALSE;
1361
1362 return kTRUE;
1363
1364}
1365
1366//____________________________________________________________________
19699b9c 1367Bool_t AliAnalyseUE::TrackSelected(AliAODTrack* part)const{
1f329128 1368
1369 // Real track selection
1370 if ( !part->TestFilterBit(fFilterBit) ) return kFALSE; // track cut selection
1371 if ( !part->IsPrimaryCandidate()) return kFALSE; // reject whatever is not linked to collision point
1372 // PID Selection: Reject everything but hadrons
1373 Bool_t isHadron = part->GetMostProbablePID()==AliAODTrack::kPion ||
1374 part->GetMostProbablePID()==AliAODTrack::kKaon ||
1375 part->GetMostProbablePID()==AliAODTrack::kProton;
1376 if ( fUseChargeHadrons && !isHadron ) return kFALSE;
1377
1378 if ( !part->Charge() ) return kFALSE; //Only charged
1379 if ( fUseSingleCharge ) { // Charge selection
1380 if ( fUsePositiveCharge && part->Charge() < 0.) return kFALSE; // keep Positives
1381 if ( !fUsePositiveCharge && part->Charge() > 0.) return kFALSE; // keep Negatives
1382 }
1383
1384 if ( part->Pt() < fTrackPtCut ) return kFALSE;
1385 if( TMath::Abs(part->Eta()) > fTrackEtaCut ) return kFALSE;
1386
1387 return kTRUE;
1388}
1389
3a9d4bcf 1390//____________________________________________________________________
1391Bool_t AliAnalyseUE::TrackSelectedEfficiency(AliAODTrack* part)const{
1392
1393 if (!fStack) AliWarning("Attention: stack not set from mother task");
1394 // Real track selection
1395 if ( !part->TestFilterBit(fFilterBit) ) return kFALSE; // track cut selection
1396 if ( !part->IsPrimaryCandidate()) return kFALSE; // reject whatever is not linked to collision point
1397 // PID Selection: Reject everything but hadrons
1398 Bool_t isHadron = part->GetMostProbablePID()==AliAODTrack::kPion ||
1399 part->GetMostProbablePID()==AliAODTrack::kKaon ||
1400 part->GetMostProbablePID()==AliAODTrack::kProton;
1401 if ( fUseChargeHadrons && !isHadron ) return kFALSE;
1402
1403 if ( !part->Charge() ) return kFALSE; //Only charged
1404 if ( fUseSingleCharge ) { // Charge selection
1405 if ( fUsePositiveCharge && part->Charge() < 0.) return kFALSE; // keep Positives
1406 if ( !fUsePositiveCharge && part->Charge() > 0.) return kFALSE; // keep Negatives
1407 }
1408
1409 /*if ( part->Pt() < fTrackPtCut ) return kFALSE;
1410 if( TMath::Abs(part->Eta()) > fTrackEtaCut ) return kFALSE;*/
1411
1412 //Check if there was a primary fulfilling the required cuts
1413 Double_t charge=-999.;
1414 Double_t pt=-999.;
1415 Double_t eta=-999.;
1416 Int_t pdgcode=-999;
1417
1418 Int_t label = part->GetLabel();
1419 TClonesArray *tca=0;
1420 tca = dynamic_cast<TClonesArray*>(fkAOD->FindListObject(AliAODMCParticle::StdBranchName()));
b85b355c 1421 if(!tca)return kFALSE;
3a9d4bcf 1422 //TParticle *partMC = fStack->ParticleFromTreeK(label);
1423 AliAODMCParticle *partMC=dynamic_cast<AliAODMCParticle*>(tca->At(TMath::Abs(label)));
b85b355c 1424 if(!partMC)return kFALSE;
3a9d4bcf 1425 if (!partMC->IsPhysicalPrimary())return kFALSE;
1426 //charge = ((TParticlePDG*)partMC->GetPDG())->Charge();
1427 charge = partMC->Charge();
1428 pt = partMC->Pt();
1429 eta = partMC->Eta();
1430 pdgcode = partMC->GetPdgCode();
1431
1432 if (!TrackMCSelected(charge, pt, eta, pdgcode)) return kFALSE;
1433 return kTRUE;
1434}
1435
1f329128 1436