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