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