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