]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/JetTasks/AliAnalysisTaskUE.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Possibility to copy the number of the TPC clusters from an AOD track to an ESD track...
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskUE.cxx
CommitLineData
f3050824 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id:$ */
6ef5bfa4 17
f3050824 18#include <TROOT.h>
19#include <TSystem.h>
20#include <TChain.h>
21#include <TFile.h>
22#include <TList.h>
23#include <TH1I.h>
24#include <TH1F.h>
25#include <TH2F.h>
623ed0a0 26#include <TProfile.h>
f3050824 27#include <TCanvas.h>
28#include <TVector3.h>
29#include <TLorentzVector.h>
30#include <TMath.h>
6ef5bfa4 31#include <TTree.h>
32#include <TBranch.h>
623ed0a0 33#include <TRandom.h>
f3050824 34
35#include "AliAnalysisTaskUE.h"
36#include "AliAnalysisManager.h"
37#include "AliMCEventHandler.h"
623ed0a0 38#include "AliMCEvent.h"
f3050824 39#include "AliAODEvent.h"
40#include "AliAODInputHandler.h"
41#include "AliAODHandler.h"
42#include "AliStack.h"
43#include "AliAODJet.h"
44#include "AliAODTrack.h"
623ed0a0 45#include "AliAODMCParticle.h"
1c796df8 46#include "AliKFVertex.h"
f3050824 47
623ed0a0 48#include "AliGenPythiaEventHeader.h"
49#include "AliAnalysisHelperJetTasks.h"
216601f0 50#include "AliInputEventHandler.h"
623ed0a0 51#include "AliStack.h"
f3050824 52#include "AliLog.h"
53
02480db2 54//__________________________________________________________________
f3050824 55// Analysis class for Underlying Event studies
56//
57// Look for correlations on the tranverse regions to
58// the leading charged jet
59//
02480db2 60// This class needs as input AOD with track and Jets.
61// The output is a list of histograms
f3050824 62//
63// AOD can be either connected to the InputEventHandler
64// for a chain of AOD files
65// or
66// to the OutputEventHandler
67// for a chain of ESD files, so this case class should be
68// in the train after the Jet finder
69//
70// Arian.Abrahantes.Quintana@cern.ch
71// Ernesto.Lopez.Torres@cern.ch
72//
73
74
75ClassImp( AliAnalysisTaskUE)
76
77////////////////////////////////////////////////////////////////////////
78
79
80//____________________________________________________________________
81AliAnalysisTaskUE:: AliAnalysisTaskUE(const char* name):
94ee88d4 82AliAnalysisTask(name,""),
7fa8b2da 83fTrigger(0),
6b7d2b7e 84fDebug(0),
c4396fd4 85fDeltaAOD(kFALSE),
86fDeltaAODBranch(""),
e1060f2b 87fAODBranch("jets"),
6b7d2b7e 88fArrayJets(0x0),
6ef5bfa4 89fAOD(0x0),
90fAODjets(0x0),
91fListOfHistos(0x0),
92fBinsPtInHist(30),
93fMinJetPtInHist(0.),
baf36668 94fMaxJetPtInHist(300.),
95fIsNorm2Area(kTRUE),
623ed0a0 96fUseMCParticleBranch(kFALSE),
97fConstrainDistance(kTRUE),
98fMinDistance(0.2),
99fSimulateChJetPt(kFALSE),
100fUseAliStack(kTRUE),
94ee88d4 101fnTracksVertex(3), // QA tracks pointing to principal vertex (= 3 default)
102fZVertex(5.),
6ef5bfa4 103fAnaType(1),
104fRegionType(1),
105fConeRadius(0.7),
623ed0a0 106fConePosition(1),
6ef5bfa4 107fAreaReg(1.5393), // Pi*0.7*0.7
108fUseChPartJet(kFALSE),
02480db2 109fUseChPartMaxPt(kFALSE),
9d9a3c85 110fUseChargeHadrons(kFALSE),
6ef5bfa4 111fUseSingleCharge(kFALSE),
112fUsePositiveCharge(kTRUE),
113fOrdering(1),
114fFilterBit(0xFF),
115fJetsOnFly(kFALSE),
623ed0a0 116fChJetPtMin(5.0),
6ef5bfa4 117fJet1EtaCut(0.2),
118fJet2DeltaPhiCut(2.616), // 150 degrees
119fJet2RatioPtCut(0.8),
120fJet3PtCut(15.),
121fTrackPtCut(0.),
122fTrackEtaCut(0.9),
c4396fd4 123fAvgTrials(1),
6ef5bfa4 124fhNJets(0x0),
125fhEleadingPt(0x0),
126fhMinRegPtDist(0x0),
127fhRegionMultMin(0x0),
128fhMinRegAvePt(0x0),
129fhMinRegSumPt(0x0),
130fhMinRegMaxPtPart(0x0),
131fhMinRegSumPtvsMult(0x0),
de6f8090 132fhdNdEtaPhiDist(0x0),
6ef5bfa4 133fhFullRegPartPtDistVsEt(0x0),
134fhTransRegPartPtDistVsEt(0x0),
135fhRegionSumPtMaxVsEt(0x0),
136fhRegionMultMax(0x0),
137fhRegionMultMaxVsEt(0x0),
138fhRegionSumPtMinVsEt(0x0), //fhRegionMultMin(0x0),
e1060f2b 139fhRegionMultMinVsEt(0x0),
140fhRegionAveSumPtVsEt(0x0),
141fhRegionDiffSumPtVsEt(0x0),
6ef5bfa4 142fhRegionAvePartPtMaxVsEt(0x0),
143fhRegionAvePartPtMinVsEt(0x0),
144fhRegionMaxPartPtMaxVsEt(0x0),
94ee88d4 145//fhRegForwardSumPtVsEt(0x0),
146//fhRegForwardMultVsEt(0x0),
147//fhRegBackwardSumPtVsEt(0x0),
148//fhRegBackwardMultVsEt(0x0),
baf36668 149fhRegForwardMult(0x0),
150fhRegForwardSumPtvsMult(0x0),
151fhRegBackwardMult(0x0),
152fhRegBackwardSumPtvsMult(0x0),
153fhRegForwardPartPtDistVsEt(0x0),
154fhRegBackwardPartPtDistVsEt(0x0),
155fhRegTransMult(0x0),
156fhRegTransSumPtVsMult(0x0),
5713f021 157fhMinRegSumPtJetPtBin(0x0),
158fhMaxRegSumPtJetPtBin(0x0),
94ee88d4 159fhVertexMult(0x0),
623ed0a0 160fh1Xsec(0x0),
161fh1Trials(0x0),
6ef5bfa4 162fSettingsTree(0x0)//, fhValidRegion(0x0)
f3050824 163{
164 // Default constructor
165 // Define input and output slots here
166 // Input slot #0 works with a TChain
167 DefineInput(0, TChain::Class());
168 // Output slot #0 writes into a TList container
169 DefineOutput(0, TList::Class());
94ee88d4 170
f3050824 171}
172
623ed0a0 173//______________________________________________________________
174Bool_t AliAnalysisTaskUE::Notify()
175{
176 //
177 // Implemented Notify() to read the cross sections
178 // and number of trials from pyxsec.root
519378fb 179 // Copy from AliAnalysisTaskJFSystematics
856f9829 180 fAvgTrials = 1;
623ed0a0 181 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
519378fb 182 Float_t xsection = 0;
183 Float_t ftrials = 1;
623ed0a0 184 if(tree){
185 TFile *curfile = tree->GetCurrentFile();
186 if (!curfile) {
187 Error("Notify","No current file");
188 return kFALSE;
189 }
190 if(!fh1Xsec||!fh1Trials){
191 Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
192 return kFALSE;
193 }
c4396fd4 194 AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
623ed0a0 195 fh1Xsec->Fill("<#sigma>",xsection);
9d9a3c85 196
c4396fd4 197 // construct average trials
198 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
a221560c 199 if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
c4396fd4 200 }
201
623ed0a0 202 return kTRUE;
203}
204
f3050824 205//____________________________________________________________________
206void AliAnalysisTaskUE::ConnectInputData(Option_t* /*option*/)
207{
6ef5bfa4 208 // Connect the input data
f3050824 209
6ef5bfa4 210 // We need AOD with tracks and jets.
211 // Since AOD can be either connected to the InputEventHandler (input chain fron AOD files)
1e996f55 212 // or to the OutputEventHandler ( AOD is create by a previus task in the train )
6ef5bfa4 213 // we need to check where it is and get the pointer to AODEvent in the right way
f3050824 214
c4396fd4 215 // Delta AODs are also implemented
216
217
218 if (fDebug > 1) AliInfo("ConnectInputData() ");
f3050824 219
220 TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
221
c4396fd4 222 if( handler && handler->InheritsFrom("AliAODInputHandler") ) { // input AOD
223 fAOD = ((AliAODInputHandler*)handler)->GetEvent();
224 if (fDebug > 1) AliInfo(" ==== Tracks from AliAODInputHandler");
02480db2 225 // Case when jets are reconstructed on the fly from AOD tracks
226 // (the Jet Finder is using the AliJetAODReader) of InputEventHandler
227 // and put in the OutputEventHandler AOD. Useful whe you want to reconstruct jets with
228 // different parameters to default ones stored in the AOD or to use a different algorithm
229 if( fJetsOnFly ) {
230 handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
231 if( handler && handler->InheritsFrom("AliAODHandler") ) {
232 fAODjets = ((AliAODHandler*)handler)->GetAOD();
233 if (fDebug > 1) AliInfo(" ==== Jets from AliAODHandler (on the fly)");
234 }
235 } else {
236 fAODjets = fAOD;
237 if (fDebug > 1) AliInfo(" ==== Jets from AliAODInputHandler");
238 }
c4396fd4 239 } else { //output AOD
6ef5bfa4 240 handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
241 if( handler && handler->InheritsFrom("AliAODHandler") ) {
c4396fd4 242 fAOD = ((AliAODHandler*)handler)->GetAOD();
6ef5bfa4 243 fAODjets = fAOD;
c4396fd4 244 if (fDebug > 1) AliInfo(" ==== Tracks and Jets from AliAODHandler");
6ef5bfa4 245 } else {
246 AliFatal("I can't get any AOD Event Handler");
247 return;
248 }
249 }
f3050824 250}
251
252//____________________________________________________________________
253void AliAnalysisTaskUE::CreateOutputObjects()
254{
6ef5bfa4 255 // Create the output container
256 //
f3050824 257 if (fDebug > 1) AliInfo("CreateOutPutData()");
258 //
259 // Histograms
1e996f55 260
7fa8b2da 261 // OpenFile(0);
f3050824 262 CreateHistos();
94ee88d4 263 fListOfHistos->SetOwner(kTRUE);
264 PostData(0,fListOfHistos);
1e996f55 265
f3050824 266}
267
f3050824 268//____________________________________________________________________
269void AliAnalysisTaskUE::Exec(Option_t */*option*/)
270{
1c796df8 271 //Trigger selection ************************************************
94ee88d4 272 /*AliInputEventHandler* inputHandler = (AliInputEventHandler*)
216601f0 273 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
02480db2 274 if (inputHandler->IsEventSelected()) {
275 if (fDebug > 1) AliInfo(" Trigger Selection: event ACCEPTED ... ");
276 } else {
277 if (fDebug > 1) AliInfo(" Trigger Selection: event REJECTED ... ");
278 return;
279 }
94ee88d4 280 */
1c796df8 281 //Event selection (vertex) *****************************************
94ee88d4 282
283 Int_t nVertex = fAOD->GetNumberOfVertices();
284 if( nVertex > 0 ) { // Only one vertex (reject pileup??)
285 AliAODVertex* vertex = (AliAODVertex*)fAOD->GetPrimaryVertex();
286 TString tpcvertex("TPCVertex");
287 if (vertex->GetName() == tpcvertex){
288 if (fDebug > 1) AliInfo(Form("Primary vertex selection: %s event REJECTED ...",vertex->GetName()));
289 return;
290 }
291 Int_t nTracksPrim = vertex->GetNContributors();
292 Double_t zVertex = vertex->GetZ();
293 if (fDebug > 1) AliInfo(Form(" Vertex in = %f with %d particles by %s data ...",zVertex,nTracksPrim,vertex->GetName()));
294 // Select a quality vertex by number of tracks?
295 if( nTracksPrim < fnTracksVertex || TMath::Abs(zVertex) > fZVertex ) {
296 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
297 return;
298 }
299 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
300 } else {
301 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
302 return;
303 }
304 /*AliKFVertex primVtx(*(fAOD->GetPrimaryVertex()));
1c796df8 305 Int_t nTracksPrim=primVtx.GetNContributors();
306 if (fDebug > 1) AliInfo(Form(" Primary-vertex Selection: %d",nTracksPrim));
307 if(!nTracksPrim){
02480db2 308 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
309 return;
310 }
1c796df8 311 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED ...");
94ee88d4 312 */
6ef5bfa4 313 // Execute analysis for current event
314 //
315 if ( fDebug > 3 ) AliInfo( " Processing event..." );
519378fb 316 // fetch the pythia header info and get the trials
317 AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
318 Float_t nTrials = 1;
c4396fd4 319 if (mcHandler) {
320 AliMCEvent* mcEvent = mcHandler->MCEvent();
321 if (mcEvent) {
322 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
323 if(pythiaGenHeader){
324 nTrials = pythiaGenHeader->Trials();
325 }
519378fb 326 }
327 }
02480db2 328 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
329 AnalyseUE();
330
f3050824 331 // Post the data
332 PostData(0, fListOfHistos);
333}
f3050824 334
335//____________________________________________________________________
336void AliAnalysisTaskUE::AnalyseUE()
337{
02480db2 338 Double_t const kMyTolerance = 0.0000001;
f3050824 339 //
340 // Look for correlations on the tranverse regions to
341 // the leading charged jet
342 //
6ef5bfa4 343
f3050824 344 // ------------------------------------------------
345 // Find Leading Jets 1,2,3
346 // (could be skipped if Jets are sort by Pt...)
347 Double_t maxPtJet1 = 0.;
348 Int_t index1 = -1;
349 Double_t maxPtJet2 = 0.; // jet 2 need for back to back inclusive
350 Int_t index2 = -1;
351 Double_t maxPtJet3 = 0.; // jet 3 need for back to back exclusive
352 Int_t index3 = -1;
6ef5bfa4 353 TVector3 jetVect[3];
354 Int_t nJets = 0;
355
c4396fd4 356
02480db2 357 if( !fUseChPartJet && fAnaType != 4 ) {
6ef5bfa4 358
359 // Use AOD Jets
c4396fd4 360 if(fDeltaAOD){
361 if (fDebug > 1) AliInfo(" ==== Jets From Delta-AODs !");
362 if (fDebug > 1) AliInfo(Form(" ==== Reading Branch: %s ",fDeltaAODBranch.Data()));
02480db2 363 fArrayJets = (TClonesArray*)fAODjets->GetList()->FindObject(fDeltaAODBranch.Data());
02480db2 364 } else {
e1060f2b 365 if (fDebug > 1) AliInfo(" ==== Read Standard-AODs !");
366 if (fDebug > 1) AliInfo(Form(" ==== Reading Branch: %s ",fAODBranch.Data()));
94ee88d4 367 fArrayJets = (TClonesArray*)fAODjets->GetList()->FindObject(fAODBranch.Data());
368 }
369 if (!fArrayJets){
370 AliFatal(" No jet-array! ");
371 return;
c4396fd4 372 }
94ee88d4 373 nJets=fArrayJets->GetEntries();
c4396fd4 374 //printf("AOD %d jets \n", nJets);
375
6ef5bfa4 376 for( Int_t i=0; i<nJets; ++i ) {
94ee88d4 377 AliAODJet* jet = (AliAODJet*)fArrayJets->At(i);
c4396fd4 378 Double_t jetPt = jet->Pt();//*1.666; // FIXME Jet Pt Correction ?????!!!
9d9a3c85 379
6ef5bfa4 380 if( jetPt > maxPtJet1 ) {
02480db2 381 maxPtJet3 = maxPtJet2; index3 = index2;
382 maxPtJet2 = maxPtJet1; index2 = index1;
383 maxPtJet1 = jetPt; index1 = i;
6ef5bfa4 384 } else if( jetPt > maxPtJet2 ) {
02480db2 385 maxPtJet3 = maxPtJet2; index3 = index2;
386 maxPtJet2 = jetPt; index2 = i;
6ef5bfa4 387 } else if( jetPt > maxPtJet3 ) {
02480db2 388 maxPtJet3 = jetPt; index3 = i;
6ef5bfa4 389 }
390 }
9d9a3c85 391
6ef5bfa4 392 if( index1 != -1 ) {
94ee88d4 393 AliAODJet *jet =(AliAODJet*) fArrayJets->At(index1);
7f80aa7b 394 if(jet)jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
6ef5bfa4 395 }
396 if( index2 != -1 ) {
94ee88d4 397 AliAODJet* jet= (AliAODJet*) fArrayJets->At(index2);
7f80aa7b 398 if(jet)jetVect[1].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
6ef5bfa4 399 }
400 if( index3 != -1 ) {
94ee88d4 401 AliAODJet* jet = (AliAODJet*) fArrayJets->At(index3);
7f80aa7b 402 if(jet)jetVect[2].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
6ef5bfa4 403 }
404
02480db2 405 }
406
407 if( fUseChPartJet ) {
408
6ef5bfa4 409 // Use "Charged Particle Jets"
410 TObjArray* jets = FindChargedParticleJets();
411 if( jets ) {
412 nJets = jets->GetEntriesFast();
413 if( nJets > 0 ) {
02480db2 414 index1 = 0;
415 AliAODJet* jet = (AliAODJet*)jets->At(0);
416 maxPtJet1 = jet->Pt();
417 jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
6ef5bfa4 418 }
419 if( nJets > 1 ) {
02480db2 420 index2 = 1;
421 AliAODJet* jet = (AliAODJet*)jets->At(1);
9d9a3c85 422 maxPtJet2 = jet->Pt();
02480db2 423 jetVect[1].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
6ef5bfa4 424 }
425 if( nJets > 2 ) {
02480db2 426 index3 = 2;
427 AliAODJet* jet = (AliAODJet*)jets->At(2);
9d9a3c85 428 maxPtJet3 = jet->Pt();
02480db2 429 jetVect[2].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
6ef5bfa4 430 }
431
432 jets->Delete();
433 delete jets;
02480db2 434 if (fDebug > 1) AliInfo(" ==== Jets Created in OUR class !: CDF algorithm ");
f3050824 435 }
436 }
c4396fd4 437
02480db2 438
439 if( fAnaType == 4 ) {
440
441 // Use "Max Pt Charge Particle"
442 TObjArray* tracks = SortChargedParticles();
443 if( tracks ) {
444 nJets = tracks->GetEntriesFast();
445 if( nJets > 0 ) {
446 index1 = 0;
447 AliAODTrack* jet = (AliAODTrack*)tracks->At(0);
448 maxPtJet1 = jet->Pt();
449 jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
450 }
451 tracks->Clear();
452 delete tracks;
453 }
454 }
455
456
457 fhNJets->Fill(nJets);
f3050824 458
6ef5bfa4 459 if( fDebug > 1 ) {
460 if( index1 < 0 ) {
461 AliInfo("\n Skipping Event, not jet found...");
462 return;
463 } else
464 AliInfo(Form("\n Pt Leading Jet = %6.1f eta=%5.3f ", maxPtJet1, jetVect[0].Eta() ));
465 }
466
f3050824 467
468 // ----------------------------------------------
6ef5bfa4 469 // Cut events by jets topology
f3050824 470 // fAnaType:
6ef5bfa4 471 // 1 = inclusive,
f3050824 472 // - Jet1 |eta| < fJet1EtaCut
473 // 2 = back to back inclusive
474 // - fulfill case 1
475 // - |Jet1.Phi - Jet2.Phi| > fJet2DeltaPhiCut
476 // - Jet2.Pt/Jet1Pt > fJet2RatioPtCut
6ef5bfa4 477 // 3 = back to back exclusive
f3050824 478 // - fulfill case 2
479 // - Jet3.Pt < fJet3PtCut
480
6ef5bfa4 481 if( index1 < 0 || TMath::Abs(jetVect[0].Eta()) > fJet1EtaCut) {
482 if( fDebug > 1 ) AliInfo("\n Skipping Event...Jet1 |eta| > fJet1EtaCut");
f3050824 483 return;
484 }
f3050824 485 // back to back inclusive
02480db2 486 if( fAnaType > 1 && fAnaType < 4 && index2 == -1 ) {
6ef5bfa4 487 if( fDebug > 1 ) AliInfo("\n Skipping Event... no second Jet found");
488 return;
489 }
02480db2 490 if( fAnaType > 1 && fAnaType < 4 && index2 > -1 ) {
f3050824 491 if( TMath::Abs(jetVect[0].DeltaPhi(jetVect[1])) < fJet2DeltaPhiCut ||
6ef5bfa4 492 maxPtJet2/maxPtJet1 < fJet2RatioPtCut ) {
493 if( fDebug > 1 ) AliInfo("\n Skipping Event... |Jet1.Phi - Jet2.Phi| < fJet2DeltaPhiCut");
494 return;
f3050824 495 }
f3050824 496 }
f3050824 497 // back to back exclusive
02480db2 498 if( fAnaType > 2 && fAnaType < 4 && index3 > -1 ) {
6ef5bfa4 499 if( maxPtJet3 > fJet3PtCut ) {
500 if( fDebug > 1 ) AliInfo("\n Skipping Event... Jet3.Pt > fJet3PtCut ");
f3050824 501 return;
502 }
f3050824 503 }
f3050824 504
623ed0a0 505 //fhEleadingPt->Fill( maxPtJet1 );
baf36668 506
f3050824 507 // ----------------------------------------------
508 // Find max and min regions
509 Double_t sumPtRegionPosit = 0.;
510 Double_t sumPtRegionNegat = 0.;
baf36668 511 Double_t sumPtRegionForward = 0;
512 Double_t sumPtRegionBackward = 0;
f3050824 513 Double_t maxPartPtRegion = 0.;
514 Int_t nTrackRegionPosit = 0;
515 Int_t nTrackRegionNegat = 0;
baf36668 516 Int_t nTrackRegionForward = 0;
517 Int_t nTrackRegionBackward = 0;
9d9a3c85 518 static Double_t const kPI = TMath::Pi();
519 static Double_t const kTWOPI = 2.*kPI;
520 static Double_t const k270rad = 270.*kPI/180.;
baf36668 521 static Double_t const k120rad = 120.*kPI/180.;
522
523 //Area for Normalization Purpose at Display histos
524 // Forward and backward
525 Double_t normArea = 1.;
526 // Transverse
527 if (fIsNorm2Area) {
528 SetRegionArea(jetVect);
529 normArea = 2.*fTrackEtaCut*k120rad ;
530 } else fAreaReg = 1.;
f3050824 531
623ed0a0 532 if (!fUseMCParticleBranch){
94ee88d4 533 AliAODVertex* vertex = (AliAODVertex*)fAOD->GetPrimaryVertex();
534 fhVertexMult->Fill( vertex->GetNContributors() );
535
623ed0a0 536 fhEleadingPt->Fill( maxPtJet1 );
537 Int_t nTracks = fAOD->GetNTracks();
02480db2 538 if (fDebug > 1) AliInfo(Form(" ==== AOD tracks = %d \n ",nTracks));
6ef5bfa4 539
623ed0a0 540 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
541
542 AliAODTrack* part = fAOD->GetTrack( ipart );
02480db2 543 if (fDebug > 1) AliInfo(Form(" ==== AOD track = %d pt = %f charge = %d \n ",ipart,part->Pt(),part->Charge()));
623ed0a0 544 if ( !part->TestFilterBit(fFilterBit) ) continue; // track cut selection
545 if (!part->IsPrimaryCandidate()) continue; // reject whatever is not linked to collision point
546 // PID Selection: Reject everything but hadrons
547 Bool_t isHadron = part->GetMostProbablePID()==AliAODTrack::kPion ||
548 part->GetMostProbablePID()==AliAODTrack::kKaon ||
549 part->GetMostProbablePID()==AliAODTrack::kProton;
9d9a3c85 550 if ( fUseChargeHadrons && !isHadron ) continue;
623ed0a0 551
552 if ( !part->Charge() ) continue; //Only charged
553 if ( fUseSingleCharge ) { // Charge selection
554 if ( fUsePositiveCharge && part->Charge() < 0.) continue; // keep Positives
555 if ( !fUsePositiveCharge && part->Charge() > 0.) continue; // keep Negatives
556 }
557
558 if ( part->Pt() < fTrackPtCut ) continue;
559 if( TMath::Abs(part->Eta()) > fTrackEtaCut ) continue;
560
561 TVector3 partVect(part->Px(), part->Py(), part->Pz());
02480db2 562 Bool_t isFlagPart = kTRUE;
623ed0a0 563 Double_t deltaPhi = jetVect[0].DeltaPhi(partVect)+k270rad;
9d9a3c85 564 if( deltaPhi > kTWOPI ) deltaPhi-= kTWOPI;
02480db2 565 if (fAnaType != 4 ) fhdNdEtaPhiDist->Fill( deltaPhi, maxPtJet1 );
566 else if (TMath::Abs(deltaPhi-k270rad) >= kMyTolerance && TMath::Abs(jetVect[0].Eta()-partVect.Eta()) >= kMyTolerance){
567 fhdNdEtaPhiDist->Fill( deltaPhi, maxPtJet1 );
568 isFlagPart = kFALSE;
569 }
623ed0a0 570
571 fhFullRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 );
572
573 Int_t region = IsTrackInsideRegion( jetVect, &partVect );
574
baf36668 575 if (region == 1) {
623ed0a0 576 if( maxPartPtRegion < part->Pt() ) maxPartPtRegion = part->Pt();
577 sumPtRegionPosit += part->Pt();
578 nTrackRegionPosit++;
579 fhTransRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 );
580 }
baf36668 581 if (region == -1) {
623ed0a0 582 if( maxPartPtRegion < part->Pt() ) maxPartPtRegion = part->Pt();
583 sumPtRegionNegat += part->Pt();
584 nTrackRegionNegat++;
585 fhTransRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 );
586 }
baf36668 587 if (region == 2){ //forward
588 sumPtRegionForward += part->Pt();
589 nTrackRegionForward++;
590 fhRegForwardPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 );
591 }
592 if (region == -2){ //backward
593 sumPtRegionBackward += part->Pt();
594 nTrackRegionBackward++;
595 fhRegBackwardPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 );
596 }
623ed0a0 597 }//end loop AOD tracks
598 }
599 else {
6ef5bfa4 600
623ed0a0 601 // this is the part we only use when we have MC information
602 // More than a test for values of it also resumes the reconstruction efficiency of jets
603 // As commented bellow if available for the data, we try to pair reconstructed jets with simulated ones
604 // afterwards we kept angular variables of MC jet to perform UE analysis over MC particles
605 // TODO: Handle Multiple jet environment. 06/2009 just suited for inclusive jet condition ( fAnaType = 1 )
6ef5bfa4 606
623ed0a0 607 AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
608 if (!mcHandler) {
609 Printf("ERROR: Could not retrieve MC event handler");
610 return;
611 }
612
613 AliMCEvent* mcEvent = mcHandler->MCEvent();
614 if (!mcEvent) {
615 Printf("ERROR: Could not retrieve MC event");
616 return;
617 }
618 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
619 if(!pythiaGenHeader){
620 return;
621 }
f3050824 622
623ed0a0 623 //Get Jets from MC header
624 Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
625 AliAODJet pythiaGenJets[4];
626 TVector3 jetVectnew[4];
627 Int_t iCount = 0;
628 for(int ip = 0;ip < nPythiaGenJets;++ip){
629 if (iCount>3) break;
630 Float_t p[4];
631 pythiaGenHeader->TriggerJet(ip,p);
632 TVector3 tempVect(p[0],p[1],p[2]);
633 if ( TMath::Abs(tempVect.Eta())>fJet1EtaCut ) continue;
634 pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
635 jetVectnew[iCount].SetXYZ(pythiaGenJets[iCount].Px(), pythiaGenJets[iCount].Py(), pythiaGenJets[iCount].Pz());
636 iCount++;
637 }
6ef5bfa4 638
623ed0a0 639 if (!iCount) return;// no jet in eta acceptance
f3050824 640
623ed0a0 641 //Search the index of the nearest MC jet to the leading jet reconstructed from the input data
642 Int_t index = 0;
643 if (fConstrainDistance){
644 Float_t deltaR = 0.;
645 Float_t dRTemp = 0.;
646 for (Int_t i=0; i<iCount; i++){
647 if (!i) {
648 dRTemp = jetVectnew[i].DeltaR(jetVect[0]);
649 index = i;
650 }
651 deltaR = jetVectnew[i].DeltaR(jetVect[0]);
652 if (deltaR < dRTemp){
653 index = i;
654 dRTemp = deltaR;
655 }
656 }
657
658 if (jetVectnew[index].DeltaR(jetVect[0]) > fMinDistance) return;
659 }
660 //Let's add some taste to jet and simulate pt of charged alone
661 //eta and phi are kept as original
662 //Play a Normal Distribution
663 Float_t random = 1.;
664 if (fSimulateChJetPt){
665 while(1){
666 random = gRandom->Gaus(0.6,0.25);
667 if (random > 0. && random < 1. &&
668 (random * jetVectnew[index].Pt()>6.)) break;
669 }
f3050824 670 }
623ed0a0 671
672 //Set new Pt & Fill histogram accordingly
673 maxPtJet1 = random * jetVectnew[index].Pt();
674
675
676 fhEleadingPt->Fill( maxPtJet1 );
677
678 if (fUseAliStack){//Try Stack Information to perform UE analysis
679
680 AliStack* mcStack = mcEvent->Stack();//Load Stack
681 Int_t nTracksMC = mcStack->GetNtrack();
682 for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
683 //Cuts
684 if(!(mcStack->IsPhysicalPrimary(iTracks))) continue;
685
686 TParticle* mctrk = mcStack->Particle(iTracks);
687
688 Double_t charge = mctrk->GetPDG()->Charge();
689 if (charge == 0) continue;
690
691 if ( fUseSingleCharge ) { // Charge selection
692 if ( fUsePositiveCharge && charge < 0.) continue; // keep Positives
693 if ( !fUsePositiveCharge && charge > 0.) continue; // keep Negatives
694 }
695
696 //Kinematics cuts on particle
697 if ((mctrk->Pt() < fTrackPtCut) || (TMath::Abs(mctrk->Eta()) > fTrackEtaCut )) continue;
698
699 Bool_t isHadron = TMath::Abs(mctrk->GetPdgCode())==211 ||
700 TMath::Abs(mctrk->GetPdgCode())==2212 ||
701 TMath::Abs(mctrk->GetPdgCode())==321;
702
9d9a3c85 703 if ( fUseChargeHadrons && !isHadron ) continue;
623ed0a0 704
705 TVector3 partVect(mctrk->Px(), mctrk->Py(), mctrk->Pz());
706
707 Double_t deltaPhi = jetVectnew[index].DeltaPhi(partVect)+k270rad;
708 if( deltaPhi > 2.*TMath::Pi() ) deltaPhi-= 2.*TMath::Pi();
9d9a3c85 709 fhdNdEtaPhiDist->Fill( deltaPhi, maxPtJet1 );
623ed0a0 710
711 fhFullRegPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 );
712
713 //We are not interested on stack organization but don't loose track of info
714 TVector3 tempVector = jetVectnew[0];
715 jetVectnew[0] = jetVectnew[index];
716 jetVectnew[index] = tempVector;
717
718 Int_t region = IsTrackInsideRegion( jetVectnew, &partVect );
719
baf36668 720 if (region == 1) {
623ed0a0 721 if( maxPartPtRegion < mctrk->Pt() ) maxPartPtRegion = mctrk->Pt();
722 sumPtRegionPosit += mctrk->Pt();
723 nTrackRegionPosit++;
724 fhTransRegPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 );
725 }
baf36668 726 if (region == -1) {
623ed0a0 727 if( maxPartPtRegion < mctrk->Pt() ) maxPartPtRegion = mctrk->Pt();
728 sumPtRegionNegat += mctrk->Pt();
729 nTrackRegionNegat++;
730 fhTransRegPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 );
731 }
baf36668 732 if (region == 2){ //forward
733 sumPtRegionForward += mctrk->Pt();
734 nTrackRegionForward++;
735 fhRegForwardPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 );
736 }
737 if (region == -2){ //backward
738 sumPtRegionBackward += mctrk->Pt();
739 nTrackRegionBackward++;
740 fhRegBackwardPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 );
741 }
742
623ed0a0 743 } // end loop stack Particles
744
745 }else{//Try mc Particle
746
747 TClonesArray* farray = (TClonesArray*)fAOD->FindListObject("mcparticles");
748
749 Int_t ntrks = farray->GetEntries();
750 if (fDebug>1) AliInfo(Form("In UE MC analysis tracks %d \n",ntrks));
751 for(Int_t i =0 ; i < ntrks; i++){
752 AliAODMCParticle* mctrk = (AliAODMCParticle*)farray->At(i);
753 //Cuts
754 if (!(mctrk->IsPhysicalPrimary())) continue;
755 //if (!(mctrk->IsPrimary())) continue;
756
757 if (mctrk->Charge() == 0 || mctrk->Charge()==-99) continue;
758
759 if (mctrk->Pt() < fTrackPtCut ) continue;
760 if( TMath::Abs(mctrk->Eta()) > fTrackEtaCut ) continue;
761
762 Bool_t isHadron = TMath::Abs(mctrk->GetPdgCode())==211 ||
763 TMath::Abs(mctrk->GetPdgCode())==2212 ||
764 TMath::Abs(mctrk->GetPdgCode())==321;
765
9d9a3c85 766 if ( fUseChargeHadrons && !isHadron ) continue;
623ed0a0 767
768 TVector3 partVect(mctrk->Px(), mctrk->Py(), mctrk->Pz());
769
770 Double_t deltaPhi = jetVectnew[index].DeltaPhi(partVect)+k270rad;
771 if( deltaPhi > 2.*TMath::Pi() ) deltaPhi-= 2.*TMath::Pi();
9d9a3c85 772 fhdNdEtaPhiDist->Fill( deltaPhi, maxPtJet1 );
623ed0a0 773
774 fhFullRegPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 );
775
776 //We are not interested on stack organization but don't loose track of info
777 TVector3 tempVector = jetVectnew[0];
778 jetVectnew[0] = jetVectnew[index];
779 jetVectnew[index] = tempVector;
780
781 Int_t region = IsTrackInsideRegion( jetVectnew, &partVect );
782
baf36668 783 if (region == 1) { //right
623ed0a0 784 if( maxPartPtRegion < mctrk->Pt() ) maxPartPtRegion = mctrk->Pt();
785 sumPtRegionPosit += mctrk->Pt();
786 nTrackRegionPosit++;
787 fhTransRegPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 );
788 }
baf36668 789 if (region == -1) { //left
623ed0a0 790 if( maxPartPtRegion < mctrk->Pt() ) maxPartPtRegion = mctrk->Pt();
791 sumPtRegionNegat += mctrk->Pt();
792 nTrackRegionNegat++;
793 fhTransRegPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 );
794 }
baf36668 795 if (region == 2){ //forward
796 sumPtRegionForward += mctrk->Pt();
797 nTrackRegionForward++;
798 fhRegForwardPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 );
799 }
800 if (region == -2){ //backward
801 sumPtRegionBackward += mctrk->Pt();
802 nTrackRegionBackward++;
803 fhRegBackwardPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 );
804 }
623ed0a0 805
806 }//end loop AliAODMCParticle tracks
f3050824 807 }
808 }
02480db2 809
810 Double_t avePosRegion = (nTrackRegionPosit) ? sumPtRegionPosit/nTrackRegionPosit : 0.;
811 Double_t aveNegRegion = (nTrackRegionNegat) ? sumPtRegionNegat/nTrackRegionNegat : 0.;
812 if( avePosRegion > aveNegRegion ) {
813 FillAvePartPtRegion( maxPtJet1, avePosRegion/fAreaReg, aveNegRegion/fAreaReg );
814 } else {
815 FillAvePartPtRegion( maxPtJet1, aveNegRegion/fAreaReg, avePosRegion/fAreaReg );
816 }
817
f3050824 818 //How quantities will be sorted before Fill Min and Max Histogram
819 // 1=Plots will be CDF-like
820 // 2=Plots will be Marchesini-like
02480db2 821 // 3=Minimum zone is selected as the one having lowest pt per track
6ef5bfa4 822 if( fOrdering == 1 ) {
823 if( sumPtRegionPosit > sumPtRegionNegat ) {
824 FillSumPtRegion( maxPtJet1, sumPtRegionPosit/fAreaReg, sumPtRegionNegat/fAreaReg );
f3050824 825 } else {
6ef5bfa4 826 FillSumPtRegion( maxPtJet1, sumPtRegionNegat/fAreaReg, sumPtRegionPosit/fAreaReg );
f3050824 827 }
828 if (nTrackRegionPosit > nTrackRegionNegat ) {
6ef5bfa4 829 FillMultRegion( maxPtJet1, nTrackRegionPosit/fAreaReg, nTrackRegionNegat/fAreaReg, sumPtRegionNegat/fAreaReg );
f3050824 830 } else {
6ef5bfa4 831 FillMultRegion( maxPtJet1, nTrackRegionNegat/fAreaReg, nTrackRegionPosit/fAreaReg, sumPtRegionPosit/fAreaReg );
f3050824 832 }
6ef5bfa4 833 } else if( fOrdering == 2 ) {
f3050824 834 if (sumPtRegionPosit > sumPtRegionNegat) {
6ef5bfa4 835 FillSumPtRegion( maxPtJet1, sumPtRegionPosit/fAreaReg, sumPtRegionNegat/fAreaReg );
836 FillMultRegion( maxPtJet1, nTrackRegionPosit/fAreaReg, nTrackRegionNegat/fAreaReg, sumPtRegionNegat/fAreaReg );
f3050824 837 } else {
6ef5bfa4 838 FillSumPtRegion( maxPtJet1, sumPtRegionNegat/fAreaReg, sumPtRegionPosit/fAreaReg );
839 FillMultRegion( maxPtJet1, nTrackRegionNegat/fAreaReg, nTrackRegionPosit/fAreaReg, sumPtRegionPosit/fAreaReg );
f3050824 840 }
02480db2 841 } else if( fOrdering == 3 ){
842 if (avePosRegion > aveNegRegion) {
843 FillSumPtRegion( maxPtJet1, sumPtRegionPosit/fAreaReg, sumPtRegionNegat/fAreaReg );
844 FillMultRegion( maxPtJet1, nTrackRegionPosit/fAreaReg, nTrackRegionNegat/fAreaReg, sumPtRegionNegat/fAreaReg );
845 }else{
846 FillSumPtRegion( maxPtJet1, sumPtRegionNegat/fAreaReg, sumPtRegionPosit/fAreaReg );
847 FillMultRegion( maxPtJet1, nTrackRegionNegat/fAreaReg, nTrackRegionPosit/fAreaReg, sumPtRegionPosit/fAreaReg );
848 }
f3050824 849 }
f3050824 850 fhRegionMaxPartPtMaxVsEt->Fill(maxPtJet1, maxPartPtRegion );
6ef5bfa4 851
852 // Compute pedestal like magnitudes
853 fhRegionDiffSumPtVsEt->Fill(maxPtJet1, TMath::Abs(sumPtRegionPosit-sumPtRegionNegat)/(2.0*fAreaReg));
854 fhRegionAveSumPtVsEt->Fill(maxPtJet1, (sumPtRegionPosit+sumPtRegionNegat)/(2.0*fAreaReg));
baf36668 855 // Transverse as a whole
94ee88d4 856 fhRegTransMult->Fill( maxPtJet1, nTrackRegionPosit + nTrackRegionNegat, (nTrackRegionPosit + nTrackRegionNegat)/(2.0*fAreaReg));
857 fhRegTransSumPtVsMult->Fill(maxPtJet1, nTrackRegionPosit + nTrackRegionNegat , (sumPtRegionNegat + sumPtRegionPosit)/(2.0*fAreaReg) );
baf36668 858
859 // Fill Histograms for Forward and away side w.r.t. leading jet direction
860 // Pt dependence
94ee88d4 861 //fhRegForwardSumPtVsEt->Fill( maxPtJet1, sumPtRegionForward/normArea );
862 //fhRegForwardMultVsEt->Fill( maxPtJet1, nTrackRegionForward/normArea );
863 //fhRegBackwardSumPtVsEt->Fill( maxPtJet1, sumPtRegionBackward/normArea );
864 //fhRegBackwardMultVsEt->Fill( maxPtJet1, nTrackRegionBackward/normArea );
baf36668 865 // Multiplicity dependence
94ee88d4 866 fhRegForwardMult->Fill(maxPtJet1, nTrackRegionForward, nTrackRegionForward/normArea);
867 fhRegForwardSumPtvsMult->Fill(maxPtJet1, nTrackRegionForward,sumPtRegionForward/normArea);
868 fhRegBackwardMult->Fill(maxPtJet1, nTrackRegionBackward, nTrackRegionBackward/normArea );
869 fhRegBackwardSumPtvsMult->Fill(maxPtJet1, nTrackRegionBackward,sumPtRegionBackward/normArea);
6ef5bfa4 870
f3050824 871}
872
f3050824 873//____________________________________________________________________
874void AliAnalysisTaskUE::FillSumPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin )
875{
6ef5bfa4 876 // Fill sumPt of control regions
877
878 // Max cone
879 fhRegionSumPtMaxVsEt->Fill( leadingE, ptMax );
880 // Min cone
881 fhRegionSumPtMinVsEt->Fill( leadingE, ptMin );
882 // MAke distributions for UE comparison with MB data
883 fhMinRegSumPt->Fill(ptMin);
5713f021 884 fhMinRegSumPtJetPtBin->Fill(leadingE, ptMin);
885 fhMaxRegSumPtJetPtBin->Fill(leadingE, ptMax);
f3050824 886}
887
888//____________________________________________________________________
889void AliAnalysisTaskUE::FillAvePartPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin )
890{
6ef5bfa4 891 // Fill average particle Pt of control regions
892
893 // Max cone
894 fhRegionAvePartPtMaxVsEt->Fill( leadingE, ptMax );
895 // Min cone
896 fhRegionAvePartPtMinVsEt->Fill( leadingE, ptMin );
897 // MAke distributions for UE comparison with MB data
898 fhMinRegAvePt->Fill(ptMin);
f3050824 899}
900
901//____________________________________________________________________
902void AliAnalysisTaskUE::FillMultRegion( Double_t leadingE, Double_t nTrackPtmax, Double_t nTrackPtmin, Double_t ptMin )
903{
6ef5bfa4 904 // Fill Nch multiplicity of control regions
905
906 // Max cone
907 fhRegionMultMaxVsEt->Fill( leadingE, nTrackPtmax );
908 fhRegionMultMax->Fill( nTrackPtmax );
909 // Min cone
910 fhRegionMultMinVsEt->Fill( leadingE, nTrackPtmin );
911 fhRegionMultMin->Fill( nTrackPtmin );
912 // MAke distributions for UE comparison with MB data
913 fhMinRegSumPtvsMult->Fill(nTrackPtmin,ptMin);
f3050824 914}
915
916//____________________________________________________________________
917Int_t AliAnalysisTaskUE::IsTrackInsideRegion(TVector3 *jetVect, TVector3 *partVect)
918{
919 // return de region in delta phi
920 // -1 negative delta phi
921 // 1 positive delta phi
922 // 0 outside region
923 static const Double_t k60rad = 60.*TMath::Pi()/180.;
924 static const Double_t k120rad = 120.*TMath::Pi()/180.;
925
926 Int_t region = 0;
927 if( fRegionType == 1 ) {
928 if( TMath::Abs(partVect->Eta()) > fTrackEtaCut ) return 0;
929 // transverse regions
baf36668 930 if (jetVect[0].DeltaPhi(*partVect) < -k60rad && jetVect[0].DeltaPhi(*partVect) > -k120rad ) region = -1; //left
931 if (jetVect[0].DeltaPhi(*partVect) > k60rad && jetVect[0].DeltaPhi(*partVect) < k120rad ) region = 1; //right
932 if (TMath::Abs(jetVect[0].DeltaPhi(*partVect)) < k60rad ) region = 2; //forward
933 if (TMath::Abs(jetVect[0].DeltaPhi(*partVect)) > k120rad ) region = -2; //backward
6ef5bfa4 934
f3050824 935 } else if( fRegionType == 2 ) {
936 // Cone regions
937 Double_t deltaR = 0.;
938
939 TVector3 positVect,negatVect;
623ed0a0 940 if (fConePosition==1){
941 positVect.SetMagThetaPhi(1, 2.*atan(exp(-jetVect[0].Eta())), jetVect[0].Phi()+TMath::PiOver2());
942 negatVect.SetMagThetaPhi(1, 2.*atan(exp(-jetVect[0].Eta())), jetVect[0].Phi()-TMath::PiOver2());
943 }else if (fConePosition==2){
944 if(fAnaType<2) AliFatal("Prevent error in Analysis type there might be only 1 jet. To avoid overflow better Correct UE config");
945 positVect.SetMagThetaPhi(1, 2.*atan(exp(-(jetVect[0].Eta()+jetVect[1].Eta())/2.)), jetVect[0].Phi()+TMath::PiOver2());
946 negatVect.SetMagThetaPhi(1, 2.*atan(exp(-(jetVect[0].Eta()+jetVect[1].Eta())/2.)), jetVect[0].Phi()-TMath::PiOver2());
947 }else if (fConePosition==3){
948 if(fAnaType<2) AliFatal("Prevent error in Analysis type there might be only 1 jet. To avoid overflow better Correct UE config");
949 Double_t weightEta = jetVect[0].Eta() * jetVect[0].Pt()/(jetVect[0].Pt() + jetVect[1].Pt()) +
950 jetVect[1].Eta() * jetVect[1].Pt()/(jetVect[0].Pt() + jetVect[1].Pt());
951 //Double_t weightEta = jetVect[0].Eta() * jetVect[0].Mag()/(jetVect[0].Mag() + jetVect[1].Mag()) +
952 // jetVect[1].Eta() * jetVect[1].Mag()/(jetVect[0].Mag() + jetVect[1].Mag());
953 positVect.SetMagThetaPhi(1, 2.*atan(exp(-weightEta)), jetVect[0].Phi()+TMath::PiOver2());
954 negatVect.SetMagThetaPhi(1, 2.*atan(exp(-weightEta)), jetVect[0].Phi()-TMath::PiOver2());
955 }
f3050824 956 if (TMath::Abs(positVect.DeltaPhi(*partVect)) < fConeRadius ) {
6ef5bfa4 957 region = 1;
958 deltaR = positVect.DrEtaPhi(*partVect);
f3050824 959 } else if (TMath::Abs(negatVect.DeltaPhi(*partVect)) < fConeRadius) {
6ef5bfa4 960 region = -1;
961 deltaR = negatVect.DrEtaPhi(*partVect);
f3050824 962 }
963
964 if (deltaR > fConeRadius) region = 0;
6ef5bfa4 965
f3050824 966 } else
967 AliError("Unknow region type");
6ef5bfa4 968
f3050824 969 // For debug (to be removed)
baf36668 970 if( fDebug > 5 && region != 0 ) AliInfo(Form("Delta Phi = %3.2f region = %d \n", jetVect[0].DeltaPhi(*partVect), region)); //fhValidRegion->Fill( partVect->Eta()-jetVect[0].Eta(), jetVect[0].DeltaPhi(*partVect) );
f3050824 971
972 return region;
973}
974
02480db2 975//____________________________________________________________________
976TObjArray* AliAnalysisTaskUE::SortChargedParticles()
977{
978 // return an array with all charged particles ordered according to their pT .
979 Int_t nTracks = fAOD->GetNTracks();
980 if( !nTracks ) return 0;
981 TObjArray* tracks = new TObjArray(nTracks);
982
983 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
984 AliAODTrack* part = fAOD->GetTrack( ipart );
985 if( !part->TestFilterBit( fFilterBit ) ) continue; // track cut selection
986 if( !part->Charge() ) continue;
987 if( part->Pt() < fTrackPtCut ) continue;
988 tracks->AddLast( part );
989 }
990 QSortTracks( *tracks, 0, tracks->GetEntriesFast() );
991
992 nTracks = tracks->GetEntriesFast();
993 if( !nTracks ) return 0;
994
995 return tracks;
996}
6ef5bfa4 997
998//____________________________________________________________________
999TObjArray* AliAnalysisTaskUE::FindChargedParticleJets()
1000{
1001 // Return a TObjArray of "charged particle jets"
1002 //
02480db2 1003 // Charged particle jet definition from reference:
6ef5bfa4 1004 //
1005 // "Charged jet evolution and the underlying event
1006 // in proton-antiproton collisions at 1.8 TeV"
1007 // PHYSICAL REVIEW D 65 092002, CDF Collaboration
1008 //
02480db2 1009 // We defined "jets" as circular regions in eta-phi space with
1010 // radius defined by R = sqrt( (eta-eta0)^2 +(phi-phi0)^2 ).
6ef5bfa4 1011 // Our jet algorithm is as follows:
1012 // 1- Order all charged particles according to their pT .
1013 // 2- Start with the highest pT particle and include in the jet all
1014 // particles within the radius R=0.7 considering each particle
1015 // in the order of decreasing pT and recalculating the centroid
1016 // of the jet after each new particle is added to the jet .
1017 // 3- Go to the next highest pT particle not already included in
1018 // a jet and add to the jet all particles not already included in
1019 // a jet within R=0.7.
1020 // 4- Continue until all particles are in a jet.
02480db2 1021 // We defined the transverse momentum of the jet to be
6ef5bfa4 1022 // the scalar pT sum of all the particles within the jet, where pT
1023 // is measured with respect to the beam axis
1024
1025 // 1 - Order all charged particles according to their pT .
1026 Int_t nTracks = fAOD->GetNTracks();
1027 if( !nTracks ) return 0;
1028 TObjArray tracks(nTracks);
1029
1030 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
1031 AliAODTrack* part = fAOD->GetTrack( ipart );
1032 if( !part->TestFilterBit(fFilterBit) ) continue; // track cut selection
1033 if( !part->Charge() ) continue;
1034 if( part->Pt() < fTrackPtCut ) continue;
1035 tracks.AddLast(part);
1036 }
1037 QSortTracks( tracks, 0, tracks.GetEntriesFast() );
1038
1039 nTracks = tracks.GetEntriesFast();
1040 if( !nTracks ) return 0;
9d9a3c85 1041
6ef5bfa4 1042 TObjArray *jets = new TObjArray(nTracks);
1043 TIter itrack(&tracks);
1044 while( nTracks ) {
1045 // 2- Start with the highest pT particle ...
1046 Float_t px,py,pz,pt;
1047 AliAODTrack* track = (AliAODTrack*)itrack.Next();
1048 if( !track ) continue;
1049 px = track->Px();
1050 py = track->Py();
1051 pz = track->Pz();
1052 pt = track->Pt(); // Use the energy member to store Pt
1053 jets->AddLast( new TLorentzVector(px, py, pz, pt) );
1054 tracks.Remove( track );
1055 TLorentzVector* jet = (TLorentzVector*)jets->Last();
9d9a3c85 1056 jet->SetPtEtaPhiE( 1., jet->Eta(), jet->Phi(), pt );
6ef5bfa4 1057 // 3- Go to the next highest pT particle not already included...
1058 AliAODTrack* track1;
9d9a3c85 1059 Double_t fPt = jet->E();
6ef5bfa4 1060 while ( (track1 = (AliAODTrack*)(itrack.Next())) ) {
9d9a3c85 1061 Double_t tphi = track1->Phi(); // here Phi is from 0 <-> 2Pi
1062 if (tphi > TMath::Pi()) tphi -= 2. * TMath::Pi(); // convert to -Pi <-> Pi
1063 Double_t dphi = TVector2::Phi_mpi_pi(jet->Phi()-tphi);
6ef5bfa4 1064 Double_t r = TMath::Sqrt( (jet->Eta()-track1->Eta())*(jet->Eta()-track1->Eta()) +
1065 dphi*dphi );
1066 if( r < fConeRadius ) {
9d9a3c85 1067 fPt = jet->E()+track1->Pt(); // Scalar sum of Pt
6ef5bfa4 1068 // recalculating the centroid
623ed0a0 1069 Double_t eta = jet->Eta()*jet->E()/fPt + track1->Eta()*track1->Pt()/fPt;
9d9a3c85 1070 Double_t phi = jet->Phi()*jet->E()/fPt + tphi*track1->Pt()/fPt;
de6f8090 1071 jet->SetPtEtaPhiE( 1., eta, phi, fPt );
6ef5bfa4 1072 tracks.Remove( track1 );
1073 }
1074 }
94ee88d4 1075
6ef5bfa4 1076 tracks.Compress();
1077 nTracks = tracks.GetEntries();
1078 // 4- Continue until all particles are in a jet.
1079 itrack.Reset();
1080 } // end while nTracks
94ee88d4 1081
6ef5bfa4 1082 // Convert to AODjets....
1083 Int_t njets = jets->GetEntriesFast();
1084 TObjArray* aodjets = new TObjArray(njets);
1085 aodjets->SetOwner(kTRUE);
1086 for(Int_t ijet=0; ijet<njets; ++ijet) {
1087 TLorentzVector* jet = (TLorentzVector*)jets->At(ijet);
623ed0a0 1088 if (jet->E() < fChJetPtMin) continue;
6ef5bfa4 1089 Float_t px, py,pz,en; // convert to 4-vector
1090 px = jet->E() * TMath::Cos(jet->Phi()); // Pt * cos(phi)
1091 py = jet->E() * TMath::Sin(jet->Phi()); // Pt * sin(phi)
1092 pz = jet->E() / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-jet->Eta())));
1093 en = TMath::Sqrt(px * px + py * py + pz * pz);
9d9a3c85 1094
6ef5bfa4 1095 aodjets->AddLast( new AliAODJet(px, py, pz, en) );
1096 }
5cb0f438 1097 jets->Delete();
1098 delete jets;
94ee88d4 1099
6ef5bfa4 1100 // Order jets according to their pT .
1101 QSortTracks( *aodjets, 0, aodjets->GetEntriesFast() );
1102
1103 // debug
1104 if (fDebug>3) AliInfo(Form(" %d Charged jets found\n",njets));
1105
1106 return aodjets;
1107}
1108
1109//____________________________________________________________________
1110void AliAnalysisTaskUE::QSortTracks(TObjArray &a, Int_t first, Int_t last)
1111{
1112 // Sort array of TObjArray of tracks by Pt using a quicksort algorithm.
1113
1114 static TObject *tmp;
1115 static int i; // "static" to save stack space
1116 int j;
1117
1118 while (last - first > 1) {
1119 i = first;
1120 j = last;
1121 for (;;) {
1122 while (++i < last && ((AliVParticle*)a[i])->Pt() > ((AliVParticle*)a[first])->Pt() )
1123 ;
1124 while (--j > first && ((AliVParticle*)a[j])->Pt() < ((AliVParticle*)a[first])->Pt() )
1125 ;
1126 if (i >= j)
1127 break;
1128
1129 tmp = a[i];
1130 a[i] = a[j];
1131 a[j] = tmp;
1132 }
1133 if (j == first) {
1134 ++first;
1135 continue;
1136 }
1137 tmp = a[first];
1138 a[first] = a[j];
1139 a[j] = tmp;
1140 if (j - first < last - (j + 1)) {
1141 QSortTracks(a, first, j);
1142 first = j + 1; // QSortTracks(j + 1, last);
1143 } else {
1144 QSortTracks(a, j + 1, last);
1145 last = j; // QSortTracks(first, j);
1146 }
1147 }
1148}
1149
1150//____________________________________________________________________
1151void AliAnalysisTaskUE::SetRegionArea(TVector3 *jetVect)
1152{
1153 Double_t fAreaCorrFactor=0.;
1154 Double_t deltaEta = 0.;
1155 if (fRegionType==1) fAreaReg = 2.*fTrackEtaCut*60.*TMath::Pi()/180.;
1156 else if (fRegionType==2){
1157 deltaEta = 0.9-TMath::Abs(jetVect[0].Eta());
1158 if (deltaEta>fConeRadius) fAreaReg = TMath::Pi()*fConeRadius*fConeRadius;
1159 else{
1160 fAreaCorrFactor = fConeRadius*fConeRadius*TMath::ACos(deltaEta/fConeRadius) -
1161 (deltaEta*fConeRadius)*TMath::Sqrt( 1. - deltaEta*deltaEta/(fConeRadius*fConeRadius));
1162 fAreaReg=TMath::Pi()*fConeRadius*fConeRadius-fAreaCorrFactor;
1163 }
1164 }else AliWarning("Unknown Rgion Type");
623ed0a0 1165 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,fAreaCorrFactor));
6ef5bfa4 1166}
1167
f3050824 1168//____________________________________________________________________
1169void AliAnalysisTaskUE::CreateHistos()
1170{
1171 fListOfHistos = new TList();
1172
1173
9d9a3c85 1174 fhNJets = new TH1F("fhNJets", "n Jet", 20, 0, 20);
f3050824 1175 fhNJets->SetXTitle("# of jets");
1176 fhNJets->Sumw2();
f3050824 1177 fListOfHistos->Add( fhNJets ); // At(0)
1178
1179 fhEleadingPt = new TH1F("hEleadingPt", "Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
1180 fhEleadingPt->SetXTitle("P_{T} (GeV/c)");
1181 fhEleadingPt->SetYTitle("dN/dP_{T}");
1182 fhEleadingPt->Sumw2();
1183 fListOfHistos->Add( fhEleadingPt ); // At(1)
1184
1185 fhMinRegPtDist = new TH1F("hMinRegPtDist", "P_{T} distribution in Min zone", 50,0.,20.);
1186 fhMinRegPtDist->SetXTitle("P_{T} (GeV/c)");
1187 fhMinRegPtDist->SetYTitle("dN/dP_{T}");
1188 fhMinRegPtDist->Sumw2();
1189 fListOfHistos->Add( fhMinRegPtDist ); // At(2)
1190
1191 fhRegionMultMin = new TH1F("hRegionMultMin", "N_{ch}^{90, min}", 21, -0.5, 20.5);
1192 fhRegionMultMin->SetXTitle("N_{ch tracks}");
1193 fhRegionMultMin->Sumw2();
1194 fListOfHistos->Add( fhRegionMultMin ); // At(3)
1195
1196 fhMinRegAvePt = new TH1F("hMinRegAvePt", "#LTp_{T}#GT", 50, 0., 20.);
1197 fhMinRegAvePt->SetXTitle("P_{T} (GeV/c)");
1198 fhMinRegAvePt->Sumw2();
1199 fListOfHistos->Add( fhMinRegAvePt ); // At(4)
1200
1201 fhMinRegSumPt = new TH1F("hMinRegSumPt", "#Sigma p_{T} ", 50, 0., 20.);
1202 fhMinRegSumPt->SetYTitle("Ed^{3}N_{tracks}/dp^{3} (c^{3}/GeV^{2})");
1203 fhMinRegSumPt->SetXTitle("#Sigma p_{T} (GeV/c)");
1204 fhMinRegSumPt->Sumw2();
1205 fListOfHistos->Add( fhMinRegSumPt ); // At(5)
6ef5bfa4 1206
f3050824 1207 fhMinRegMaxPtPart = new TH1F("hMinRegMaxPtPart", "max(p_{T})|_{event} ", 50, 0., 20.);
1208 fhMinRegMaxPtPart->SetYTitle("Ed^{3}N_{tracks}/dp^{3} (c^{3}/GeV^{2})");
1209 fhMinRegMaxPtPart->SetXTitle("p_{T} (GeV/c)");
1210 fhMinRegMaxPtPart->Sumw2();
1211 fListOfHistos->Add( fhMinRegMaxPtPart ); // At(6)
1212
1213 fhMinRegSumPtvsMult = new TH1F("hMinRegSumPtvsMult", "#Sigma p_{T} vs. Multiplicity ", 21, -0.5, 20.5);
1214 fhMinRegSumPtvsMult->SetYTitle("#Sigma p_{T} (GeV/c)");
1215 fhMinRegSumPtvsMult->SetXTitle("N_{charge}");
1216 fhMinRegSumPtvsMult->Sumw2();
1217 fListOfHistos->Add( fhMinRegSumPtvsMult ); // At(7);
6ef5bfa4 1218
9d9a3c85 1219 fhdNdEtaPhiDist = new TH2F("hdNdEtaPhiDist", Form("Charge particle density |#eta|<%3.1f vs #Delta#phi", fTrackEtaCut),
1220 62, 0., 2.*TMath::Pi(), fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
de6f8090 1221 fhdNdEtaPhiDist->SetXTitle("#Delta#phi");
9d9a3c85 1222 fhdNdEtaPhiDist->SetYTitle("Leading Jet P_{T}");
de6f8090 1223 fhdNdEtaPhiDist->Sumw2();
1224 fListOfHistos->Add( fhdNdEtaPhiDist ); // At(8)
6ef5bfa4 1225
f3050824 1226 // Can be use to get part pt distribution for differente Jet Pt bins
9d9a3c85 1227 fhFullRegPartPtDistVsEt = new TH2F("hFullRegPartPtDistVsEt", Form( "dN/dP_{T} |#eta|<%3.1f vs Leading Jet P_{T}", fTrackEtaCut),
27d185bd 1228 200,0.,100., fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
f3050824 1229 fhFullRegPartPtDistVsEt->SetYTitle("Leading Jet P_{T}");
1230 fhFullRegPartPtDistVsEt->SetXTitle("p_{T}");
1231 fhFullRegPartPtDistVsEt->Sumw2();
1232 fListOfHistos->Add( fhFullRegPartPtDistVsEt ); // At(9)
6ef5bfa4 1233
1234 // Can be use to get part pt distribution for differente Jet Pt bins
9d9a3c85 1235 fhTransRegPartPtDistVsEt = new TH2F("hTransRegPartPtDistVsEt", Form( "dN/dP_{T} in tranvese regions |#eta|<%3.1f vs Leading Jet P_{T}", fTrackEtaCut),
27d185bd 1236 200,0.,100., fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
f3050824 1237 fhTransRegPartPtDistVsEt->SetYTitle("Leading Jet P_{T}");
1238 fhTransRegPartPtDistVsEt->SetXTitle("p_{T}");
1239 fhTransRegPartPtDistVsEt->Sumw2();
1240 fListOfHistos->Add( fhTransRegPartPtDistVsEt ); // At(10)
1241
1242
1243 fhRegionSumPtMaxVsEt = new TH1F("hRegionSumPtMaxVsEt", "P_{T}^{90, max} vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
1244 fhRegionSumPtMaxVsEt->SetXTitle("P_{T} (GeV/c)");
1245 fhRegionSumPtMaxVsEt->Sumw2();
1246 fListOfHistos->Add( fhRegionSumPtMaxVsEt ); // At(11)
6ef5bfa4 1247
f3050824 1248 fhRegionSumPtMinVsEt = new TH1F("hRegionSumPtMinVsEt", "P_{T}^{90, min} vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
1249 fhRegionSumPtMinVsEt->SetXTitle("P_{T} (GeV/c)");
1250 fhRegionSumPtMinVsEt->Sumw2();
1251 fListOfHistos->Add( fhRegionSumPtMinVsEt ); // At(12)
1252
1253 fhRegionMultMax = new TH1I("hRegionMultMax", "N_{ch}^{90, max}", 21, -0.5, 20.5);
1254 fhRegionMultMax->SetXTitle("N_{ch tracks}");
1255 fhRegionMultMax->Sumw2();
1256 fListOfHistos->Add( fhRegionMultMax ); // At(13)
6ef5bfa4 1257
f3050824 1258 fhRegionMultMaxVsEt = new TH1F("hRegionMultMaxVsEt", "N_{ch}^{90, max} vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
1259 fhRegionMultMaxVsEt->SetXTitle("E (GeV hRegionAveSumPtVsEt/c)");
1260 fhRegionMultMaxVsEt->Sumw2();
1261 fListOfHistos->Add( fhRegionMultMaxVsEt ); // At(14)
6ef5bfa4 1262
f3050824 1263 fhRegionMultMinVsEt = new TH1F("hRegionMultMinVsEt", "N_{ch}^{90, min} vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
1264 fhRegionMultMinVsEt->SetXTitle("E (GeV/c)");
1265 fhRegionMultMinVsEt->Sumw2();
1266 fListOfHistos->Add( fhRegionMultMinVsEt ); // At(15)
6ef5bfa4 1267
f3050824 1268 fhRegionAveSumPtVsEt = new TH1F("hRegionAveSumPtVsEt", "(P_{T}^{90, max} + P_{T}^{90, min})/2 vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
1269 fhRegionAveSumPtVsEt->SetXTitle("P_{T} (GeV/c)");
1270 fhRegionAveSumPtVsEt->Sumw2();
1271 fListOfHistos->Add( fhRegionAveSumPtVsEt ); // At(16)
6ef5bfa4 1272
f3050824 1273 fhRegionDiffSumPtVsEt= new TH1F("hRegionPtDiffVsEt", "(P_{T}^{90, max} - P_{T}^{90, min}) vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
1274 fhRegionDiffSumPtVsEt->SetXTitle("P_{T} (GeV/c)");
1275 fhRegionDiffSumPtVsEt->Sumw2();
1276 fListOfHistos->Add( fhRegionDiffSumPtVsEt ); // At(17)
1277
1278 fhRegionAvePartPtMaxVsEt = new TH1F("hRegionAvePartPtMaxVsEt", "#LTp_{T}#GT^{90, max} vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
1279 fhRegionAvePartPtMaxVsEt->SetXTitle("P_{T} (GeV/c)");
1280 fhRegionAvePartPtMaxVsEt->Sumw2();
1281 fListOfHistos->Add( fhRegionAvePartPtMaxVsEt ); // At(18)
6ef5bfa4 1282
f3050824 1283 fhRegionAvePartPtMinVsEt = new TH1F("hRegionAvePartPtMinVsEt", "#LTp_{T}#GT^{90, min} vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
1284 fhRegionAvePartPtMinVsEt->SetXTitle("P_{T} (GeV/c)");
1285 fhRegionAvePartPtMinVsEt->Sumw2();
1286 fListOfHistos->Add( fhRegionAvePartPtMinVsEt ); // At(19)
6ef5bfa4 1287
f3050824 1288 fhRegionMaxPartPtMaxVsEt = new TH1F("hRegionMaxPartPtMaxVsEt", "max(p_{T})^{90} vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
1289 fhRegionMaxPartPtMaxVsEt->SetXTitle("P_{T} (GeV/c)");
1290 fhRegionMaxPartPtMaxVsEt->Sumw2();
1291 fListOfHistos->Add( fhRegionMaxPartPtMaxVsEt ); // At(20)
94ee88d4 1292 /*
baf36668 1293 fhRegForwardSumPtVsEt = new TH1F("hRegForwardSumPtVsEt", "Forward #sum{p_{T}} vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
1294 fhRegForwardSumPtVsEt->SetXTitle("P_{T} (GeV/c)");
1295 fhRegForwardSumPtVsEt->Sumw2();
1296 fListOfHistos->Add( fhRegForwardSumPtVsEt ); // At(21)
1297
1298 fhRegForwardMultVsEt = new TH1F("hRegForwardMultVsEt", "Forward #sum{N_{ch}} vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
1299 fhRegForwardMultVsEt->SetXTitle("P_{T} (GeV/c)");
1300 fhRegForwardMultVsEt->Sumw2();
1301 fListOfHistos->Add( fhRegForwardMultVsEt ); // At(22)
1302
1303 fhRegBackwardSumPtVsEt = new TH1F("hRegBackwardSumPtVsEt", "Backward #sum{p_{T}} vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
1304 fhRegBackwardSumPtVsEt->SetXTitle("P_{T} (GeV/c)");
1305 fhRegBackwardSumPtVsEt->Sumw2();
1306 fListOfHistos->Add( fhRegBackwardSumPtVsEt ); // At(23)
1307
1308 fhRegBackwardMultVsEt = new TH1F("hRegBackwardMultVsEt", "Backward #sum{N_{ch}} vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
1309 fhRegBackwardMultVsEt->SetXTitle("P_{T} (GeV/c)");
1310 fhRegBackwardMultVsEt->Sumw2();
1311 fListOfHistos->Add( fhRegBackwardMultVsEt ); // At(24)
94ee88d4 1312 */
baf36668 1313 fhRegForwardMult = new TH2F("hRegForwardMult", "N_{ch}^{forward}",
1314 fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist, 21, -0.5, 20.5);
1315 fhRegForwardMult->SetXTitle("N_{ch tracks}");
1316 fhRegForwardMult->Sumw2();
1317 fListOfHistos->Add( fhRegForwardMult ); // At(25)
1318
1319 fhRegForwardSumPtvsMult = new TH2F("hRegForwardSumPtvsMult", "Forward #Sigma p_{T} vs. Multiplicity ",
1320 fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist, 21, -0.5, 20.5);
1321 fhRegForwardSumPtvsMult->SetYTitle("#Sigma p_{T} (GeV/c)");
1322 fhRegForwardSumPtvsMult->SetXTitle("N_{charge}");
1323 fhRegForwardSumPtvsMult->Sumw2();
1324 fListOfHistos->Add( fhRegForwardSumPtvsMult ); // At(26);
1325
1326 fhRegBackwardMult = new TH2F("hRegBackwardMult", "N_{ch}^{backward}",
1327 fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist, 21, -0.5, 20.5);
1328 fhRegBackwardMult->SetXTitle("N_{ch tracks}");
1329 fhRegBackwardMult->Sumw2();
1330 fListOfHistos->Add( fhRegBackwardMult ); // At(27)
1331
1332 fhRegBackwardSumPtvsMult = new TH2F("hRegBackwardSumPtvsMult", "Backward #Sigma p_{T} vs. Multiplicity ",
1333 fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist, 21, -0.5, 20.5);
1334 fhRegBackwardSumPtvsMult->SetYTitle("#Sigma p_{T} (GeV/c)");
1335 fhRegBackwardSumPtvsMult->SetXTitle("N_{charge}");
1336 fhRegBackwardSumPtvsMult->Sumw2();
1337 fListOfHistos->Add( fhRegBackwardSumPtvsMult ); // At(28);
1338
1339 fhRegForwardPartPtDistVsEt = new TH2F("hRegForwardPartPtDistVsEt", Form( "dN/dP_{T} |#eta|<%3.1f vs Leading Jet P_{T}", fTrackEtaCut),
27d185bd 1340 200,0.,100., fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
baf36668 1341 fhRegForwardPartPtDistVsEt->SetYTitle("Leading Jet P_{T}");
1342 fhRegForwardPartPtDistVsEt->SetXTitle("p_{T}");
1343 fhRegForwardPartPtDistVsEt->Sumw2();
1344 fListOfHistos->Add( fhRegForwardPartPtDistVsEt ); // At(29)
1345
1346 fhRegBackwardPartPtDistVsEt = new TH2F("hRegBackwardPartPtDistVsEt", Form( "dN/dP_{T} |#eta|<%3.1f vs Leading Jet P_{T}", fTrackEtaCut),
27d185bd 1347 200,0.,100., fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
baf36668 1348 fhRegBackwardPartPtDistVsEt->SetYTitle("Leading Jet P_{T}");
1349 fhRegBackwardPartPtDistVsEt->SetXTitle("p_{T}");
1350 fhRegBackwardPartPtDistVsEt->Sumw2();
1351 fListOfHistos->Add( fhRegBackwardPartPtDistVsEt ); // At(30)
1352
1353 fhRegTransMult = new TH2F("hRegTransMult", "N_{ch}^{transv}",
1354 fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist, 21, -0.5, 20.5);
1355 fhRegTransMult->SetXTitle("N_{ch tracks}");
1356 fhRegTransMult->Sumw2();
1357 fListOfHistos->Add( fhRegTransMult ); // At(31)
1358
1359 fhRegTransSumPtVsMult = new TH2F("hRegTransSumPtVsMult", "Transverse #Sigma p_{T} vs. Multiplicity ",
1360 fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist, 21, -0.5, 20.5);
1361 fhRegTransSumPtVsMult->SetYTitle("#Sigma p_{T} (GeV/c)");
1362 fhRegTransSumPtVsMult->SetXTitle("N_{charge}");
1363 fhRegTransSumPtVsMult->Sumw2();
1364 fListOfHistos->Add( fhRegTransSumPtVsMult ); // At(32);
1365
5713f021 1366 fhMinRegSumPtJetPtBin = new TH2F("hMinRegSumPtJetPtBin", "Transverse Min Reg #Sigma p_{T} per jet Pt Bin",
1367 fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist, 50, 0., 20.);
1368 fhMinRegSumPtJetPtBin->SetXTitle("Leading Jet P_{T}");
1369 fhMinRegSumPtJetPtBin->Sumw2();
1370 fListOfHistos->Add( fhMinRegSumPtJetPtBin ); // At(33)
1371
1372 fhMaxRegSumPtJetPtBin = new TH2F("hMaxRegSumPtJetPtBin", "Transverse Max Reg #Sigma p_{T} per jet Pt Bin",
1373 fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist, 50, 0., 20.);
1374 fhMaxRegSumPtJetPtBin->SetXTitle("Leading Jet P_{T}");
1375 fhMaxRegSumPtJetPtBin->Sumw2();
1376 fListOfHistos->Add( fhMaxRegSumPtJetPtBin ); // At(34)
1377
94ee88d4 1378 fhVertexMult = new TH1F("hVertexMult", "Multiplicity in Main Vertex", 81, -0.5 , 80.5);
1379 fhVertexMult->SetXTitle("Main Vertex Multiplicity");
1380 fhVertexMult->Sumw2();
1381 fListOfHistos->Add( fhVertexMult ); //At(35)
1382
623ed0a0 1383 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
1384 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
1385 fh1Xsec->Sumw2();
94ee88d4 1386 fListOfHistos->Add( fh1Xsec ); //At(36)
623ed0a0 1387
1388 fh1Trials = new TH1F("fh1Trials","trials from pyxsec.root",1,0,1);
1389 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
1390 fh1Trials->Sumw2();
94ee88d4 1391 fListOfHistos->Add( fh1Trials ); //At(37)
1e996f55 1392
6ef5bfa4 1393 fSettingsTree = new TTree("UEAnalysisSettings","Analysis Settings in UE estimation");
5c1827c4 1394 fSettingsTree->Branch("fFilterBit", &fFilterBit,"FilterBit/I");
1c796df8 1395 fSettingsTree->Branch("fTrigger", &fTrigger,"TriggerFlag/I");
1e996f55 1396 fSettingsTree->Branch("fConeRadius", &fConeRadius,"Rad/D");
1397 fSettingsTree->Branch("fJet1EtaCut", &fJet1EtaCut, "LeadJetEtaCut/D");
1398 fSettingsTree->Branch("fJet2DeltaPhiCut", &fJet2DeltaPhiCut, "DeltaPhi/D");
1399 fSettingsTree->Branch("fJet2RatioPtCut", &fJet2RatioPtCut, "Jet2Ratio/D");
1400 fSettingsTree->Branch("fJet3PtCut", &fJet3PtCut, "Jet3PtCut/D");
1401 fSettingsTree->Branch("fTrackPtCut", &fTrackPtCut, "TrackPtCut/D");
1402 fSettingsTree->Branch("fTrackEtaCut", &fTrackEtaCut, "TrackEtaCut/D");
1403 fSettingsTree->Branch("fAnaType", &fAnaType, "Ana/I");
1404 fSettingsTree->Branch("fRegionType", &fRegionType,"Reg/I");
1405 fSettingsTree->Branch("fOrdering", &fOrdering,"OrderMeth/I");
1406 fSettingsTree->Branch("fUseChPartJet", &fUseChPartJet,"UseChPart/O");
9d9a3c85 1407 fSettingsTree->Branch("fUseChargeHadrons", &fUseChargeHadrons,"UseChHadrons/O");
1e996f55 1408 fSettingsTree->Branch("fUseSingleCharge", &fUseSingleCharge,"UseSingleCh/O");
1409 fSettingsTree->Branch("fUsePositiveCharge", &fUsePositiveCharge,"UsePositiveCh/O");
1410 fSettingsTree->Fill();
1411
1412
94ee88d4 1413 fListOfHistos->Add( fSettingsTree ); // At(38)
6ef5bfa4 1414
1415 /*
1416 // For debug region selection
1417 fhValidRegion = new TH2F("hValidRegion", "dN/d#eta/d#phi",
1418 20, -2.,2., 62, -TMath::Pi(), TMath::Pi());
1419 fhValidRegion->SetYTitle("#Delta#phi");
1420 fhValidRegion->Sumw2();
1421 fListOfHistos->Add( fhValidRegion ); // At(15)
1422 */
f3050824 1423}
f3050824 1424
e1060f2b 1425
1426
6ef5bfa4 1427//____________________________________________________________________
1428void AliAnalysisTaskUE::Terminate(Option_t */*option*/)
1429{
1430 // Terminate analysis
1431 //
1432
1433 //Save Analysis Settings
1e996f55 1434 // WriteSettings();
6ef5bfa4 1435
1436 // Normalize histos to region area TODO:
1437 // Normalization done at Analysis, taking into account
1438 // area variations on per-event basis (cone case)
5713f021 1439
1440 //HIGH WARNING!!!!!: DO NOT SCALE ANY OF THE ORIGINAL HISTOGRAMS
1441 //MAKE A COPY, DRAW IT, And later sacale that copy. CAF Issue!!!!!
6ef5bfa4 1442
1443 // Draw some Test plot to the screen
1444 if (!gROOT->IsBatch()) {
1445
1446 // Update pointers reading them from the output slot
1447 fListOfHistos = dynamic_cast<TList*> (GetOutputData(0));
1448 if( !fListOfHistos ) {
1449 AliError("Histogram List is not available");
1450 return;
1451 }
9d9a3c85 1452 fhNJets = (TH1F*)fListOfHistos->At(0);
6ef5bfa4 1453 fhEleadingPt = (TH1F*)fListOfHistos->At(1);
9d9a3c85 1454 fhdNdEtaPhiDist = (TH2F*)fListOfHistos->At(8);
6ef5bfa4 1455 fhRegionSumPtMaxVsEt = (TH1F*)fListOfHistos->At(11);
1456 fhRegionSumPtMinVsEt = (TH1F*)fListOfHistos->At(12);
1457 fhRegionMultMaxVsEt = (TH1F*)fListOfHistos->At(14);
1458 fhRegionMultMinVsEt = (TH1F*)fListOfHistos->At(15);
1459 fhRegionAveSumPtVsEt = (TH1F*)fListOfHistos->At(16);
94ee88d4 1460 //fhRegForwardSumPtVsEt = (TH1F*)fListOfHistos->At(21);
1461 //fhRegBackwardSumPtVsEt = (TH1F*)fListOfHistos->At(23);
6ef5bfa4 1462
623ed0a0 1463 //fhValidRegion = (TH2F*)fListOfHistos->At(21);
1464
6ef5bfa4 1465 // Canvas
623ed0a0 1466 TCanvas* c1 = new TCanvas("c1",Form("sumPt dist (%s)", GetTitle()),60,60,1100,700);
6ef5bfa4 1467 c1->Divide(2,2);
1468 c1->cd(1);
1469 TH1F *h1r = new TH1F("hRegionEtvsSumPtMax" , "", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
1470 h1r->Divide(fhRegionSumPtMaxVsEt,fhEleadingPt,1,1);
1471 //h1r->Scale( areafactor );
1472 h1r->SetMarkerStyle(20);
1473 h1r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
1474 h1r->SetYTitle("P_{T}^{90, max}");
1475 h1r->DrawCopy("p");
1476
1477 c1->cd(2);
1478
1479 TH1F *h2r = new TH1F("hRegionEtvsSumPtMin" , "", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
1480 h2r->Divide(fhRegionSumPtMinVsEt,fhEleadingPt,1,1);
1481 //h2r->Scale( areafactor );
1482 h2r->SetMarkerStyle(20);
1483 h2r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
1484 h2r->SetYTitle("P_{T}^{90, min}");
1485 h2r->DrawCopy("p");
1486
1487 c1->cd(3);
1488 TH1F *h4r = new TH1F("hRegionEtvsDiffPt" , "", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
94ee88d4 1489 //TH1F *h41r = new TH1F("hRegForwvsDiffPt" , "", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
1490 //TH1F *h42r = new TH1F("hRegBackvsDiffPt" , "", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
1491 //h41r->Divide(fhRegForwardSumPtVsEt,fhEleadingPt,1,1);
1492 //h42r->Divide(fhRegBackwardSumPtVsEt,fhEleadingPt,1,1);
6ef5bfa4 1493 h4r->Divide(fhRegionAveSumPtVsEt,fhEleadingPt,1,1);
baf36668 1494 //h4r->Scale(2.); // make average
6ef5bfa4 1495 //h4r->Scale( areafactor );
1496 h4r->SetYTitle("#DeltaP_{T}^{90}");
1497 h4r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
1498 h4r->SetMarkerStyle(20);
1499 h4r->DrawCopy("p");
94ee88d4 1500 /*h41r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
baf36668 1501 h41r->SetMarkerStyle(22);
1502 h41r->DrawCopy("p same");
1503 h42r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
1504 h42r->SetMarkerStyle(23);
1505 h42r->SetMarkerColor(kRed);
1506 h42r->DrawCopy("p same");
94ee88d4 1507 */
6ef5bfa4 1508 c1->cd(4);
1509 TH1F *h5r = new TH1F("hRegionMultMaxVsEtleading", "", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
1510 TH1F *h6r = new TH1F("hRegionMultMinVsEtleading", "", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
1511
1512 h5r->Divide(fhRegionMultMaxVsEt,fhEleadingPt,1,1);
1513 h6r->Divide(fhRegionMultMinVsEt,fhEleadingPt,1,1);
1514 //h5r->Scale( areafactor );
1515 //h6r->Scale( areafactor );
1516 h5r->SetYTitle("N_{Tracks}^{90}");
1517 h5r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
1518 h5r->SetMarkerStyle(20);
1519 h5r->DrawCopy("p");
1520 h6r->SetMarkerStyle(21);
1521 h6r->SetMarkerColor(2);
1522 h6r->SetYTitle("N_{Tracks}^{90}");
1523 h6r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
1524 h6r->DrawCopy("p same");
1525 c1->Update();
1526
623ed0a0 1527 //Get Normalization
1528 fh1Xsec = (TProfile*)fListOfHistos->At(21);
1529 fh1Trials = (TH1F*)fListOfHistos->At(22);
1530
1531 Double_t xsec = fh1Xsec->GetBinContent(1);
1532 Double_t ntrials = fh1Trials->GetBinContent(1);
1533 Double_t normFactor = xsec/ntrials;
6b7d2b7e 1534 if(fDebug > 1)Printf("xSec %f nTrials %f Norm %f \n",xsec,ntrials,normFactor);
623ed0a0 1535
1536
5713f021 1537
6ef5bfa4 1538 TCanvas* c2 = new TCanvas("c2","Jet Pt dist",160,160,1200,800);
5713f021 1539 TH1 * copy = 0;
6ef5bfa4 1540 c2->Divide(2,2);
1541 c2->cd(1);
1542 fhEleadingPt->SetMarkerStyle(20);
1543 fhEleadingPt->SetMarkerColor(2);
5713f021 1544 //if( normFactor > 0.) fhEleadingPt->Scale(normFactor);
623ed0a0 1545 //fhEleadingPt->Draw("p");
5713f021 1546 copy = fhEleadingPt->DrawCopy("p");
1547 if( normFactor > 0.) copy->Scale(normFactor);
6ef5bfa4 1548 gPad->SetLogy();
1549
1550 c2->cd(2);
9d9a3c85 1551 Int_t xbin1 = fhdNdEtaPhiDist->GetYaxis()->FindFixBin(fMinJetPtInHist);
1552 Int_t xbin2 = fhdNdEtaPhiDist->GetYaxis()->FindFixBin(fMaxJetPtInHist);
1553 TH1D* dNdEtaPhiDistAllJets = fhdNdEtaPhiDist->ProjectionX("dNdEtaPhiDistAllJets",xbin1,xbin2);
1554 dNdEtaPhiDistAllJets->SetMarkerStyle(20);
1555 dNdEtaPhiDistAllJets->SetMarkerColor(2);
1556 dNdEtaPhiDistAllJets->DrawCopy("p");
6ef5bfa4 1557 gPad->SetLogy();
1558
1559 c2->cd(3);
1560 fhNJets->DrawCopy();
1561
623ed0a0 1562 //c2->cd(4);
1563 //fhValidRegion->DrawCopy("p");
1564
6ef5bfa4 1565 //fhTransRegPartPtDist = (TH1F*)fListOfHistos->At(2);
1566 fhRegionMultMin = (TH1F*)fListOfHistos->At(3);
1567 fhMinRegAvePt = (TH1F*)fListOfHistos->At(4);
1568 fhMinRegSumPt = (TH1F*)fListOfHistos->At(5);
1569 //fhMinRegMaxPtPart = (TH1F*)fListOfHistos->At(6);
1570 fhMinRegSumPtvsMult = (TH1F*)fListOfHistos->At(7);
1571
baf36668 1572
6ef5bfa4 1573 // Canvas
623ed0a0 1574 TCanvas* c3 = new TCanvas("c3"," pT dist",160,160,1200,800);
6ef5bfa4 1575 c3->Divide(2,2);
1576 c3->cd(1);
1577 /*fhTransRegPartPtDist->SetMarkerStyle(20);
1578 fhTransRegPartPtDist->SetMarkerColor(2);
1579 fhTransRegPartPtDist->Scale(areafactor/fhTransRegPartPtDist->GetEntries());
1580 fhTransRegPartPtDist->DrawCopy("p");
1581 gPad->SetLogy();
1582 */
1583 c3->cd(2);
1584 fhMinRegSumPt->SetMarkerStyle(20);
1585 fhMinRegSumPt->SetMarkerColor(2);
1586 //fhMinRegSumPt->Scale(areafactor);
1587 fhMinRegSumPt->DrawCopy("p");
1588 gPad->SetLogy();
1589
1590 c3->cd(3);
1591 fhMinRegAvePt->SetMarkerStyle(20);
1592 fhMinRegAvePt->SetMarkerColor(2);
1593 //fhMinRegAvePt->Scale(areafactor);
1594 fhMinRegAvePt->DrawCopy("p");
1595 gPad->SetLogy();
1596
1597 c3->cd(4);
1598
1599 TH1F *h7r = new TH1F("hRegionMultMinVsMult", "", 21, -0.5, 20.5);
1600
1601 h7r->Divide(fhMinRegSumPtvsMult,fhRegionMultMin,1,1);
1602
1603 h7r->SetMarkerStyle(20);
1604 h7r->SetMarkerColor(2);
1605 h7r->DrawCopy("p");
1606
1607 c3->Update();
1608
1609
1610 /* c2->cd(3);
1611 fhValidRegion->SetMarkerStyle(20);
1612 fhValidRegion->SetMarkerColor(2);
1613 fhValidRegion->DrawCopy("p");
1614 */
1615 c2->Update();
1616 } else {
1617 AliInfo(" Batch mode, not histograms will be shown...");
1618 }
1619
1620 if( fDebug > 1 )
1621 AliInfo("End analysis");
1622
1623}
f3050824 1624
6ef5bfa4 1625void AliAnalysisTaskUE::WriteSettings()
1e996f55 1626{
6ef5bfa4 1627 if (fDebug>5){
1628 AliInfo(" All Analysis Settings in Saved Tree");
1629 fSettingsTree->Scan();
1630 }
1631}