]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PWG4/JetTasks/AliAnalysisTaskUE.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Missing packages produce a Fatal
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskUE.cxx
... / ...
CommitLineData
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:$ */
17
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>
26#include <TProfile.h>
27#include <TCanvas.h>
28#include <TVector3.h>
29#include <TLorentzVector.h>
30#include <TMath.h>
31#include <TTree.h>
32#include <TBranch.h>
33#include <TRandom.h>
34
35#include "AliAnalysisTaskUE.h"
36#include "AliAnalysisManager.h"
37#include "AliMCEventHandler.h"
38#include "AliMCEvent.h"
39#include "AliAODEvent.h"
40#include "AliAODInputHandler.h"
41#include "AliAODHandler.h"
42#include "AliStack.h"
43#include "AliAODJet.h"
44#include "AliAODTrack.h"
45#include "AliAODMCParticle.h"
46#include "AliKFVertex.h"
47
48#include "AliGenPythiaEventHeader.h"
49#include "AliAnalysisHelperJetTasks.h"
50#include "AliInputEventHandler.h"
51#include "AliStack.h"
52#include "AliLog.h"
53
54//__________________________________________________________________
55// Analysis class for Underlying Event studies
56//
57// Look for correlations on the tranverse regions to
58// the leading charged jet
59//
60// This class needs as input AOD with track and Jets.
61// The output is a list of histograms
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):
82AliAnalysisTask(name, ""),
83fTrigger(0),
84fDebug(0),
85fDeltaAOD(kFALSE),
86fDeltaAODBranch(""),
87fAODBranch("jets"),
88fArrayJets(0x0),
89fAOD(0x0),
90fAODjets(0x0),
91fListOfHistos(0x0),
92fBinsPtInHist(30),
93fMinJetPtInHist(0.),
94fMaxJetPtInHist(300.),
95fIsNorm2Area(kTRUE),
96fUseMCParticleBranch(kFALSE),
97fConstrainDistance(kTRUE),
98fMinDistance(0.2),
99fSimulateChJetPt(kFALSE),
100fUseAliStack(kTRUE),
101fAnaType(1),
102fRegionType(1),
103fConeRadius(0.7),
104fConePosition(1),
105fAreaReg(1.5393), // Pi*0.7*0.7
106fUseChPartJet(kFALSE),
107fUseChPartMaxPt(kFALSE),
108fUseChargeHadrons(kFALSE),
109fUseSingleCharge(kFALSE),
110fUsePositiveCharge(kTRUE),
111fOrdering(1),
112fFilterBit(0xFF),
113fJetsOnFly(kFALSE),
114fChJetPtMin(5.0),
115fJet1EtaCut(0.2),
116fJet2DeltaPhiCut(2.616), // 150 degrees
117fJet2RatioPtCut(0.8),
118fJet3PtCut(15.),
119fTrackPtCut(0.),
120fTrackEtaCut(0.9),
121fAvgTrials(1),
122fhNJets(0x0),
123fhEleadingPt(0x0),
124fhMinRegPtDist(0x0),
125fhRegionMultMin(0x0),
126fhMinRegAvePt(0x0),
127fhMinRegSumPt(0x0),
128fhMinRegMaxPtPart(0x0),
129fhMinRegSumPtvsMult(0x0),
130fhdNdEtaPhiDist(0x0),
131fhFullRegPartPtDistVsEt(0x0),
132fhTransRegPartPtDistVsEt(0x0),
133fhRegionSumPtMaxVsEt(0x0),
134fhRegionMultMax(0x0),
135fhRegionMultMaxVsEt(0x0),
136fhRegionSumPtMinVsEt(0x0), //fhRegionMultMin(0x0),
137fhRegionMultMinVsEt(0x0),
138fhRegionAveSumPtVsEt(0x0),
139fhRegionDiffSumPtVsEt(0x0),
140fhRegionAvePartPtMaxVsEt(0x0),
141fhRegionAvePartPtMinVsEt(0x0),
142fhRegionMaxPartPtMaxVsEt(0x0),
143fhRegForwardSumPtVsEt(0x0),
144fhRegForwardMultVsEt(0x0),
145fhRegBackwardSumPtVsEt(0x0),
146fhRegBackwardMultVsEt(0x0),
147fhRegForwardMult(0x0),
148fhRegForwardSumPtvsMult(0x0),
149fhRegBackwardMult(0x0),
150fhRegBackwardSumPtvsMult(0x0),
151fhRegForwardPartPtDistVsEt(0x0),
152fhRegBackwardPartPtDistVsEt(0x0),
153fhRegTransMult(0x0),
154fhRegTransSumPtVsMult(0x0),
155fhMinRegSumPtJetPtBin(0x0),
156fhMaxRegSumPtJetPtBin(0x0),
157fh1Xsec(0x0),
158fh1Trials(0x0),
159fSettingsTree(0x0)//, fhValidRegion(0x0)
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
169//______________________________________________________________
170Bool_t AliAnalysisTaskUE::Notify()
171{
172 //
173 // Implemented Notify() to read the cross sections
174 // and number of trials from pyxsec.root
175 // Copy from AliAnalysisTaskJFSystematics
176 fAvgTrials = 1;
177 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
178 Float_t xsection = 0;
179 Float_t ftrials = 1;
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 }
190 AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
191 fh1Xsec->Fill("<#sigma>",xsection);
192
193 // construct average trials
194 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
195 if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
196 }
197
198 return kTRUE;
199}
200
201//____________________________________________________________________
202void AliAnalysisTaskUE::ConnectInputData(Option_t* /*option*/)
203{
204 // Connect the input data
205
206 // We need AOD with tracks and jets.
207 // Since AOD can be either connected to the InputEventHandler (input chain fron AOD files)
208 // or to the OutputEventHandler ( AOD is create by a previus task in the train )
209 // we need to check where it is and get the pointer to AODEvent in the right way
210
211 // Delta AODs are also implemented
212
213
214 if (fDebug > 1) AliInfo("ConnectInputData() ");
215
216 TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
217
218 if( handler && handler->InheritsFrom("AliAODInputHandler") ) { // input AOD
219 fAOD = ((AliAODInputHandler*)handler)->GetEvent();
220 if (fDebug > 1) AliInfo(" ==== Tracks from AliAODInputHandler");
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 }
235 } else { //output AOD
236 handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
237 if( handler && handler->InheritsFrom("AliAODHandler") ) {
238 fAOD = ((AliAODHandler*)handler)->GetAOD();
239 fAODjets = fAOD;
240 if (fDebug > 1) AliInfo(" ==== Tracks and Jets from AliAODHandler");
241 } else {
242 AliFatal("I can't get any AOD Event Handler");
243 return;
244 }
245 }
246}
247
248//____________________________________________________________________
249void AliAnalysisTaskUE::CreateOutputObjects()
250{
251 // Create the output container
252 //
253 if (fDebug > 1) AliInfo("CreateOutPutData()");
254 //
255 // Histograms
256
257 // OpenFile(0);
258 CreateHistos();
259
260 fListOfHistos->SetOwner(kTRUE);
261 PostData(0, fListOfHistos);
262
263}
264
265
266//____________________________________________________________________
267void AliAnalysisTaskUE::Exec(Option_t */*option*/)
268{
269 //Trigger selection ************************************************
270 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
271 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
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 }
278
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){
284 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
285 return;
286 }
287 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED ...");
288
289 // Execute analysis for current event
290 //
291 if ( fDebug > 3 ) AliInfo( " Processing event..." );
292 // fetch the pythia header info and get the trials
293 AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
294 Float_t nTrials = 1;
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 }
302 }
303 }
304 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
305 AnalyseUE();
306
307 // Post the data
308 PostData(0, fListOfHistos);
309}
310
311//____________________________________________________________________
312void AliAnalysisTaskUE::AnalyseUE()
313{
314 Double_t const kMyTolerance = 0.0000001;
315 //
316 // Look for correlations on the tranverse regions to
317 // the leading charged jet
318 //
319
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;
329 TVector3 jetVect[3];
330 Int_t nJets = 0;
331
332
333 if( !fUseChPartJet && fAnaType != 4 ) {
334
335 // Use AOD Jets
336 if(fDeltaAOD){
337 if (fDebug > 1) AliInfo(" ==== Jets From Delta-AODs !");
338 if (fDebug > 1) AliInfo(Form(" ==== Reading Branch: %s ",fDeltaAODBranch.Data()));
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 {
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();
350 }
351 //printf("AOD %d jets \n", nJets);
352
353 for( Int_t i=0; i<nJets; ++i ) {
354 AliAODJet* jet;
355 if (fDeltaAOD){
356 jet =(AliAODJet*)fArrayJets->At(i);
357 }else{
358 jet = (AliAODJet*)((TClonesArray*)fAODjets->FindListObject(fAODBranch.Data()))->At(i);
359 }
360 Double_t jetPt = jet->Pt();//*1.666; // FIXME Jet Pt Correction ?????!!!
361
362 if( jetPt > maxPtJet1 ) {
363 maxPtJet3 = maxPtJet2; index3 = index2;
364 maxPtJet2 = maxPtJet1; index2 = index1;
365 maxPtJet1 = jetPt; index1 = i;
366 } else if( jetPt > maxPtJet2 ) {
367 maxPtJet3 = maxPtJet2; index3 = index2;
368 maxPtJet2 = jetPt; index2 = i;
369 } else if( jetPt > maxPtJet3 ) {
370 maxPtJet3 = jetPt; index3 = i;
371 }
372 }
373
374 if( index1 != -1 ) {
375 AliAODJet *jet = 0;
376 if (fDeltaAOD) {
377 jet =(AliAODJet*) fArrayJets->At(index1);
378 }else{
379 jet = (AliAODJet*)((TClonesArray*)fAODjets->FindListObject(fAODBranch.Data()))->At(index1);
380 }
381 if(jet)jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
382 }
383 if( index2 != -1 ) {
384 AliAODJet* jet = 0;
385 if (fDeltaAOD) {
386 jet= (AliAODJet*) fArrayJets->At(index2);
387 }else{
388 jet=(AliAODJet*) ((TClonesArray*)fAODjets->FindListObject(fAODBranch.Data()))->At(index2);
389 }
390 if(jet)jetVect[1].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
391 }
392 if( index3 != -1 ) {
393 AliAODJet* jet = 0;
394 if (fDeltaAOD) {
395 jet= (AliAODJet*) fArrayJets->At(index3);
396 }else{
397 ((TClonesArray*)fAODjets->FindListObject(fAODBranch.Data()))->At(index3);
398 }
399 if(jet)jetVect[2].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
400 }
401
402 }
403
404 if( fUseChPartJet ) {
405
406 // Use "Charged Particle Jets"
407 TObjArray* jets = FindChargedParticleJets();
408 if( jets ) {
409 nJets = jets->GetEntriesFast();
410 if( nJets > 0 ) {
411 index1 = 0;
412 AliAODJet* jet = (AliAODJet*)jets->At(0);
413 maxPtJet1 = jet->Pt();
414 jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
415 }
416 if( nJets > 1 ) {
417 index2 = 1;
418 AliAODJet* jet = (AliAODJet*)jets->At(1);
419 maxPtJet2 = jet->Pt();
420 jetVect[1].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
421 }
422 if( nJets > 2 ) {
423 index3 = 2;
424 AliAODJet* jet = (AliAODJet*)jets->At(2);
425 maxPtJet3 = jet->Pt();
426 jetVect[2].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
427 }
428
429 jets->Delete();
430 delete jets;
431 if (fDebug > 1) AliInfo(" ==== Jets Created in OUR class !: CDF algorithm ");
432 }
433 }
434
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);
455
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
464
465 // ----------------------------------------------
466 // Cut events by jets topology
467 // fAnaType:
468 // 1 = inclusive,
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
474 // 3 = back to back exclusive
475 // - fulfill case 2
476 // - Jet3.Pt < fJet3PtCut
477
478 if( index1 < 0 || TMath::Abs(jetVect[0].Eta()) > fJet1EtaCut) {
479 if( fDebug > 1 ) AliInfo("\n Skipping Event...Jet1 |eta| > fJet1EtaCut");
480 return;
481 }
482 // back to back inclusive
483 if( fAnaType > 1 && fAnaType < 4 && index2 == -1 ) {
484 if( fDebug > 1 ) AliInfo("\n Skipping Event... no second Jet found");
485 return;
486 }
487 if( fAnaType > 1 && fAnaType < 4 && index2 > -1 ) {
488 if( TMath::Abs(jetVect[0].DeltaPhi(jetVect[1])) < fJet2DeltaPhiCut ||
489 maxPtJet2/maxPtJet1 < fJet2RatioPtCut ) {
490 if( fDebug > 1 ) AliInfo("\n Skipping Event... |Jet1.Phi - Jet2.Phi| < fJet2DeltaPhiCut");
491 return;
492 }
493 }
494 // back to back exclusive
495 if( fAnaType > 2 && fAnaType < 4 && index3 > -1 ) {
496 if( maxPtJet3 > fJet3PtCut ) {
497 if( fDebug > 1 ) AliInfo("\n Skipping Event... Jet3.Pt > fJet3PtCut ");
498 return;
499 }
500 }
501
502 //fhEleadingPt->Fill( maxPtJet1 );
503
504 // ----------------------------------------------
505 // Find max and min regions
506 Double_t sumPtRegionPosit = 0.;
507 Double_t sumPtRegionNegat = 0.;
508 Double_t sumPtRegionForward = 0;
509 Double_t sumPtRegionBackward = 0;
510 Double_t maxPartPtRegion = 0.;
511 Int_t nTrackRegionPosit = 0;
512 Int_t nTrackRegionNegat = 0;
513 Int_t nTrackRegionForward = 0;
514 Int_t nTrackRegionBackward = 0;
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.;
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.;
528
529 if (!fUseMCParticleBranch){
530 fhEleadingPt->Fill( maxPtJet1 );
531 Int_t nTracks = fAOD->GetNTracks();
532 if (fDebug > 1) AliInfo(Form(" ==== AOD tracks = %d \n ",nTracks));
533
534 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
535
536 AliAODTrack* part = fAOD->GetTrack( ipart );
537 if (fDebug > 1) AliInfo(Form(" ==== AOD track = %d pt = %f charge = %d \n ",ipart,part->Pt(),part->Charge()));
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;
544 if ( fUseChargeHadrons && !isHadron ) continue;
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());
556 Bool_t isFlagPart = kTRUE;
557 Double_t deltaPhi = jetVect[0].DeltaPhi(partVect)+k270rad;
558 if( deltaPhi > kTWOPI ) deltaPhi-= kTWOPI;
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 }
564
565 fhFullRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 );
566
567 Int_t region = IsTrackInsideRegion( jetVect, &partVect );
568
569 if (region == 1) {
570 if( maxPartPtRegion < part->Pt() ) maxPartPtRegion = part->Pt();
571 sumPtRegionPosit += part->Pt();
572 nTrackRegionPosit++;
573 fhTransRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 );
574 }
575 if (region == -1) {
576 if( maxPartPtRegion < part->Pt() ) maxPartPtRegion = part->Pt();
577 sumPtRegionNegat += part->Pt();
578 nTrackRegionNegat++;
579 fhTransRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 );
580 }
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 }
591 }//end loop AOD tracks
592 }
593 else {
594
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 )
600
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 }
616
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 }
632
633 if (!iCount) return;// no jet in eta acceptance
634
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 }
664 }
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
697 if ( fUseChargeHadrons && !isHadron ) continue;
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();
703 fhdNdEtaPhiDist->Fill( deltaPhi, maxPtJet1 );
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
714 if (region == 1) {
715 if( maxPartPtRegion < mctrk->Pt() ) maxPartPtRegion = mctrk->Pt();
716 sumPtRegionPosit += mctrk->Pt();
717 nTrackRegionPosit++;
718 fhTransRegPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 );
719 }
720 if (region == -1) {
721 if( maxPartPtRegion < mctrk->Pt() ) maxPartPtRegion = mctrk->Pt();
722 sumPtRegionNegat += mctrk->Pt();
723 nTrackRegionNegat++;
724 fhTransRegPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 );
725 }
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
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
760 if ( fUseChargeHadrons && !isHadron ) continue;
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();
766 fhdNdEtaPhiDist->Fill( deltaPhi, maxPtJet1 );
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
777 if (region == 1) { //right
778 if( maxPartPtRegion < mctrk->Pt() ) maxPartPtRegion = mctrk->Pt();
779 sumPtRegionPosit += mctrk->Pt();
780 nTrackRegionPosit++;
781 fhTransRegPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 );
782 }
783 if (region == -1) { //left
784 if( maxPartPtRegion < mctrk->Pt() ) maxPartPtRegion = mctrk->Pt();
785 sumPtRegionNegat += mctrk->Pt();
786 nTrackRegionNegat++;
787 fhTransRegPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 );
788 }
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 }
799
800 }//end loop AliAODMCParticle tracks
801 }
802 }
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
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
815 // 3=Minimum zone is selected as the one having lowest pt per track
816 if( fOrdering == 1 ) {
817 if( sumPtRegionPosit > sumPtRegionNegat ) {
818 FillSumPtRegion( maxPtJet1, sumPtRegionPosit/fAreaReg, sumPtRegionNegat/fAreaReg );
819 } else {
820 FillSumPtRegion( maxPtJet1, sumPtRegionNegat/fAreaReg, sumPtRegionPosit/fAreaReg );
821 }
822 if (nTrackRegionPosit > nTrackRegionNegat ) {
823 FillMultRegion( maxPtJet1, nTrackRegionPosit/fAreaReg, nTrackRegionNegat/fAreaReg, sumPtRegionNegat/fAreaReg );
824 } else {
825 FillMultRegion( maxPtJet1, nTrackRegionNegat/fAreaReg, nTrackRegionPosit/fAreaReg, sumPtRegionPosit/fAreaReg );
826 }
827 } else if( fOrdering == 2 ) {
828 if (sumPtRegionPosit > sumPtRegionNegat) {
829 FillSumPtRegion( maxPtJet1, sumPtRegionPosit/fAreaReg, sumPtRegionNegat/fAreaReg );
830 FillMultRegion( maxPtJet1, nTrackRegionPosit/fAreaReg, nTrackRegionNegat/fAreaReg, sumPtRegionNegat/fAreaReg );
831 } else {
832 FillSumPtRegion( maxPtJet1, sumPtRegionNegat/fAreaReg, sumPtRegionPosit/fAreaReg );
833 FillMultRegion( maxPtJet1, nTrackRegionNegat/fAreaReg, nTrackRegionPosit/fAreaReg, sumPtRegionPosit/fAreaReg );
834 }
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 }
843 }
844 fhRegionMaxPartPtMaxVsEt->Fill(maxPtJet1, maxPartPtRegion );
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));
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);
864
865}
866
867//____________________________________________________________________
868void AliAnalysisTaskUE::FillSumPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin )
869{
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);
878 fhMinRegSumPtJetPtBin->Fill(leadingE, ptMin);
879 fhMaxRegSumPtJetPtBin->Fill(leadingE, ptMax);
880}
881
882//____________________________________________________________________
883void AliAnalysisTaskUE::FillAvePartPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin )
884{
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);
893}
894
895//____________________________________________________________________
896void AliAnalysisTaskUE::FillMultRegion( Double_t leadingE, Double_t nTrackPtmax, Double_t nTrackPtmin, Double_t ptMin )
897{
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);
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
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
928
929 } else if( fRegionType == 2 ) {
930 // Cone regions
931 Double_t deltaR = 0.;
932
933 TVector3 positVect,negatVect;
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 }
950 if (TMath::Abs(positVect.DeltaPhi(*partVect)) < fConeRadius ) {
951 region = 1;
952 deltaR = positVect.DrEtaPhi(*partVect);
953 } else if (TMath::Abs(negatVect.DeltaPhi(*partVect)) < fConeRadius) {
954 region = -1;
955 deltaR = negatVect.DrEtaPhi(*partVect);
956 }
957
958 if (deltaR > fConeRadius) region = 0;
959
960 } else
961 AliError("Unknow region type");
962
963 // For debug (to be removed)
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) );
965
966 return region;
967}
968
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}
991
992//____________________________________________________________________
993TObjArray* AliAnalysisTaskUE::FindChargedParticleJets()
994{
995 // Return a TObjArray of "charged particle jets"
996 //
997 // Charged particle jet definition from reference:
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 //
1003 // We defined "jets" as circular regions in eta-phi space with
1004 // radius defined by R = sqrt( (eta-eta0)^2 +(phi-phi0)^2 ).
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.
1015 // We defined the transverse momentum of the jet to be
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;
1035
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();
1050 jet->SetPtEtaPhiE( 1., jet->Eta(), jet->Phi(), pt );
1051 // 3- Go to the next highest pT particle not already included...
1052 AliAODTrack* track1;
1053 Double_t fPt = jet->E();
1054 while ( (track1 = (AliAODTrack*)(itrack.Next())) ) {
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);
1058 Double_t r = TMath::Sqrt( (jet->Eta()-track1->Eta())*(jet->Eta()-track1->Eta()) +
1059 dphi*dphi );
1060 if( r < fConeRadius ) {
1061 fPt = jet->E()+track1->Pt(); // Scalar sum of Pt
1062 // recalculating the centroid
1063 Double_t eta = jet->Eta()*jet->E()/fPt + track1->Eta()*track1->Pt()/fPt;
1064 Double_t phi = jet->Phi()*jet->E()/fPt + tphi*track1->Pt()/fPt;
1065 jet->SetPtEtaPhiE( 1., eta, phi, fPt );
1066 tracks.Remove( track1 );
1067 }
1068 }
1069 tracks.Compress();
1070 nTracks = tracks.GetEntries();
1071 // 4- Continue until all particles are in a jet.
1072 itrack.Reset();
1073 } // end while nTracks
1074
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);
1081 if (jet->E() < fChJetPtMin) continue;
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);
1087
1088 aodjets->AddLast( new AliAODJet(px, py, pz, en) );
1089 }
1090 jets->Delete();
1091 delete jets;
1092
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");
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));
1159}
1160
1161//____________________________________________________________________
1162void AliAnalysisTaskUE::CreateHistos()
1163{
1164 fListOfHistos = new TList();
1165
1166
1167 fhNJets = new TH1F("fhNJets", "n Jet", 20, 0, 20);
1168 fhNJets->SetXTitle("# of jets");
1169 fhNJets->Sumw2();
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)
1199
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);
1211
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);
1214 fhdNdEtaPhiDist->SetXTitle("#Delta#phi");
1215 fhdNdEtaPhiDist->SetYTitle("Leading Jet P_{T}");
1216 fhdNdEtaPhiDist->Sumw2();
1217 fListOfHistos->Add( fhdNdEtaPhiDist ); // At(8)
1218
1219 // Can be use to get part pt distribution for differente Jet Pt bins
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);
1222 fhFullRegPartPtDistVsEt->SetYTitle("Leading Jet P_{T}");
1223 fhFullRegPartPtDistVsEt->SetXTitle("p_{T}");
1224 fhFullRegPartPtDistVsEt->Sumw2();
1225 fListOfHistos->Add( fhFullRegPartPtDistVsEt ); // At(9)
1226
1227 // Can be use to get part pt distribution for differente Jet Pt bins
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);
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)
1240
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)
1250
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)
1255
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)
1260
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)
1265
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)
1275
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)
1280
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)
1285
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
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
1371 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
1372 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
1373 fh1Xsec->Sumw2();
1374 fListOfHistos->Add( fh1Xsec ); //At(35)
1375
1376 fh1Trials = new TH1F("fh1Trials","trials from pyxsec.root",1,0,1);
1377 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
1378 fh1Trials->Sumw2();
1379 fListOfHistos->Add( fh1Trials ); //At(36)
1380
1381 fSettingsTree = new TTree("UEAnalysisSettings","Analysis Settings in UE estimation");
1382 fSettingsTree->Branch("fFilterBit", &fFilterBit,"FilterBit/I");
1383 fSettingsTree->Branch("fTrigger", &fTrigger,"TriggerFlag/I");
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");
1395 fSettingsTree->Branch("fUseChargeHadrons", &fUseChargeHadrons,"UseChHadrons/O");
1396 fSettingsTree->Branch("fUseSingleCharge", &fUseSingleCharge,"UseSingleCh/O");
1397 fSettingsTree->Branch("fUsePositiveCharge", &fUsePositiveCharge,"UsePositiveCh/O");
1398 fSettingsTree->Fill();
1399
1400
1401 fListOfHistos->Add( fSettingsTree ); // At(37)
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 */
1411}
1412
1413
1414
1415//____________________________________________________________________
1416void AliAnalysisTaskUE::Terminate(Option_t */*option*/)
1417{
1418 // Terminate analysis
1419 //
1420
1421 //Save Analysis Settings
1422 // WriteSettings();
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)
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!!!!!
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 }
1440 fhNJets = (TH1F*)fListOfHistos->At(0);
1441 fhEleadingPt = (TH1F*)fListOfHistos->At(1);
1442 fhdNdEtaPhiDist = (TH2F*)fListOfHistos->At(8);
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);
1448 fhRegForwardSumPtVsEt = (TH1F*)fListOfHistos->At(21);
1449 fhRegBackwardSumPtVsEt = (TH1F*)fListOfHistos->At(23);
1450
1451 //fhValidRegion = (TH2F*)fListOfHistos->At(21);
1452
1453 // Canvas
1454 TCanvas* c1 = new TCanvas("c1",Form("sumPt dist (%s)", GetTitle()),60,60,1100,700);
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);
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);
1481 h4r->Divide(fhRegionAveSumPtVsEt,fhEleadingPt,1,1);
1482 //h4r->Scale(2.); // make average
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");
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");
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
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;
1522 if(fDebug > 1)Printf("xSec %f nTrials %f Norm %f \n",xsec,ntrials,normFactor);
1523
1524
1525
1526 TCanvas* c2 = new TCanvas("c2","Jet Pt dist",160,160,1200,800);
1527 TH1 * copy = 0;
1528 c2->Divide(2,2);
1529 c2->cd(1);
1530 fhEleadingPt->SetMarkerStyle(20);
1531 fhEleadingPt->SetMarkerColor(2);
1532 //if( normFactor > 0.) fhEleadingPt->Scale(normFactor);
1533 //fhEleadingPt->Draw("p");
1534 copy = fhEleadingPt->DrawCopy("p");
1535 if( normFactor > 0.) copy->Scale(normFactor);
1536 gPad->SetLogy();
1537
1538 c2->cd(2);
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");
1545 gPad->SetLogy();
1546
1547 c2->cd(3);
1548 fhNJets->DrawCopy();
1549
1550 //c2->cd(4);
1551 //fhValidRegion->DrawCopy("p");
1552
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
1560
1561 // Canvas
1562 TCanvas* c3 = new TCanvas("c3"," pT dist",160,160,1200,800);
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}
1612
1613void AliAnalysisTaskUE::WriteSettings()
1614{
1615 if (fDebug>5){
1616 AliInfo(" All Analysis Settings in Saved Tree");
1617 fSettingsTree->Scan();
1618 }
1619}