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