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