Added PID task, some fixes for coding conventions
[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>
26#include <TCanvas.h>
27#include <TVector3.h>
28#include <TLorentzVector.h>
29#include <TMath.h>
6ef5bfa4 30#include <TTree.h>
31#include <TBranch.h>
f3050824 32
33#include "AliAnalysisTaskUE.h"
34#include "AliAnalysisManager.h"
35#include "AliMCEventHandler.h"
36#include "AliAODEvent.h"
37#include "AliAODInputHandler.h"
38#include "AliAODHandler.h"
39#include "AliStack.h"
40#include "AliAODJet.h"
41#include "AliAODTrack.h"
42
43#include "AliLog.h"
44
45//
46// Analysis class for Underlying Event studies
47//
48// Look for correlations on the tranverse regions to
49// the leading charged jet
50//
51// This class needs as input AOD with track and Jets
52// the output is a list of histograms
53//
54// AOD can be either connected to the InputEventHandler
55// for a chain of AOD files
56// or
57// to the OutputEventHandler
58// for a chain of ESD files, so this case class should be
59// in the train after the Jet finder
60//
61// Arian.Abrahantes.Quintana@cern.ch
62// Ernesto.Lopez.Torres@cern.ch
63//
64
65
66ClassImp( AliAnalysisTaskUE)
67
68////////////////////////////////////////////////////////////////////////
69
70
71//____________________________________________________________________
72AliAnalysisTaskUE:: AliAnalysisTaskUE(const char* name):
6ef5bfa4 73AliAnalysisTask(name, ""),
74fDebug(kFALSE),
75fAOD(0x0),
76fAODjets(0x0),
77fListOfHistos(0x0),
78fBinsPtInHist(30),
79fMinJetPtInHist(0.),
80fMaxJetPtInHist(300.),
81fAnaType(1),
82fRegionType(1),
83fConeRadius(0.7),
84fAreaReg(1.5393), // Pi*0.7*0.7
85fUseChPartJet(kFALSE),
86fUseSingleCharge(kFALSE),
87fUsePositiveCharge(kTRUE),
88fOrdering(1),
89fFilterBit(0xFF),
90fJetsOnFly(kFALSE),
91fJet1EtaCut(0.2),
92fJet2DeltaPhiCut(2.616), // 150 degrees
93fJet2RatioPtCut(0.8),
94fJet3PtCut(15.),
95fTrackPtCut(0.),
96fTrackEtaCut(0.9),
97fhNJets(0x0),
98fhEleadingPt(0x0),
99fhMinRegPtDist(0x0),
100fhRegionMultMin(0x0),
101fhMinRegAvePt(0x0),
102fhMinRegSumPt(0x0),
103fhMinRegMaxPtPart(0x0),
104fhMinRegSumPtvsMult(0x0),
de6f8090 105fhdNdEtaPhiDist(0x0),
6ef5bfa4 106fhFullRegPartPtDistVsEt(0x0),
107fhTransRegPartPtDistVsEt(0x0),
108fhRegionSumPtMaxVsEt(0x0),
109fhRegionMultMax(0x0),
110fhRegionMultMaxVsEt(0x0),
111fhRegionSumPtMinVsEt(0x0), //fhRegionMultMin(0x0),
112fhRegionMultMinVsEt(0x0),
113fhRegionAveSumPtVsEt(0x0),
114fhRegionDiffSumPtVsEt(0x0),
115fhRegionAvePartPtMaxVsEt(0x0),
116fhRegionAvePartPtMinVsEt(0x0),
117fhRegionMaxPartPtMaxVsEt(0x0),
118fSettingsTree(0x0)//, fhValidRegion(0x0)
f3050824 119{
120 // Default constructor
121 // Define input and output slots here
122 // Input slot #0 works with a TChain
123 DefineInput(0, TChain::Class());
124 // Output slot #0 writes into a TList container
125 DefineOutput(0, TList::Class());
126}
127
f3050824 128//____________________________________________________________________
129void AliAnalysisTaskUE::ConnectInputData(Option_t* /*option*/)
130{
6ef5bfa4 131 // Connect the input data
f3050824 132
6ef5bfa4 133 // We need AOD with tracks and jets.
134 // Since AOD can be either connected to the InputEventHandler (input chain fron AOD files)
135 // or to the OutputEventHandler ( AOD is created by a previus task in the train )
136 // we need to check where it is and get the pointer to AODEvent in the right way
f3050824 137
6ef5bfa4 138 if (fDebug > 1) AliInfo("ConnectInputData() \n");
f3050824 139
140 TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
141
142 if( handler && handler->InheritsFrom("AliAODInputHandler") ) {
6ef5bfa4 143 fAOD = ((AliAODInputHandler*)handler)->GetEvent();
144 if (fDebug > 1) AliInfo(" ==== Tracks from AliAODInputHandler");
145 // Case when jets are reconstructed on the fly from AOD tracks
146 // (the Jet Finder is using the AliJetAODReader) of InputEventHandler
147 // and put in the OutputEventHandler AOD. Useful whe you want to reconstruct jets with
148 // different parameters to default ones stored in the AOD or to use a different algorithm
149 if( fJetsOnFly ) {
f3050824 150 handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
151 if( handler && handler->InheritsFrom("AliAODHandler") ) {
6ef5bfa4 152 fAODjets = ((AliAODHandler*)handler)->GetAOD();
153 if (fDebug > 1) AliInfo(" ==== Jets from AliAODHandler");
f3050824 154 }
6ef5bfa4 155 } else {
156 fAODjets = fAOD;
157 if (fDebug > 1) AliInfo(" ==== Jets from AliAODInputHandler");
158 }
159 } else {
160 handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
161 if( handler && handler->InheritsFrom("AliAODHandler") ) {
162 fAOD = ((AliAODHandler*)handler)->GetAOD();
163 fAODjets = fAOD;
164 if (fDebug > 1) AliInfo(" ==== Tracks and Jets from AliAODHandler");
165 } else {
166 AliFatal("I can't get any AOD Event Handler");
167 return;
168 }
169 }
f3050824 170}
171
172//____________________________________________________________________
173void AliAnalysisTaskUE::CreateOutputObjects()
174{
6ef5bfa4 175 // Create the output container
176 //
f3050824 177 if (fDebug > 1) AliInfo("CreateOutPutData()");
178 //
179 // Histograms
180 CreateHistos();
181 fListOfHistos->SetOwner(kTRUE);
6ef5bfa4 182 // OpenFile(0);
f3050824 183}
184
185
186//____________________________________________________________________
187void AliAnalysisTaskUE::Exec(Option_t */*option*/)
188{
6ef5bfa4 189 // Execute analysis for current event
190 //
191 if ( fDebug > 3 ) AliInfo( " Processing event..." );
192 AnalyseUE();
193
f3050824 194 // Post the data
195 PostData(0, fListOfHistos);
196}
f3050824 197
198//____________________________________________________________________
199void AliAnalysisTaskUE::AnalyseUE()
200{
201 //
202 // Look for correlations on the tranverse regions to
203 // the leading charged jet
204 //
6ef5bfa4 205
f3050824 206 // ------------------------------------------------
207 // Find Leading Jets 1,2,3
208 // (could be skipped if Jets are sort by Pt...)
209 Double_t maxPtJet1 = 0.;
210 Int_t index1 = -1;
211 Double_t maxPtJet2 = 0.; // jet 2 need for back to back inclusive
212 Int_t index2 = -1;
213 Double_t maxPtJet3 = 0.; // jet 3 need for back to back exclusive
214 Int_t index3 = -1;
6ef5bfa4 215 TVector3 jetVect[3];
216 Int_t nJets = 0;
217
218 if( !fUseChPartJet ) {
219
220 // Use AOD Jets
221 nJets = fAODjets->GetNJets();
222 // printf("AOD %d jets \n", nJets);
223 for( Int_t i=0; i<nJets; ++i ) {
224 AliAODJet* jet = fAODjets->GetJet(i);
225 Double_t jetPt = jet->Pt();//*1.666; // FIXME Jet Pt Correction ?????!!!
226 if( jetPt > maxPtJet1 ) {
227 maxPtJet3 = maxPtJet2; index3 = index2;
228 maxPtJet2 = maxPtJet1; index2 = index1;
229 maxPtJet1 = jetPt; index1 = i;
230 } else if( jetPt > maxPtJet2 ) {
231 maxPtJet3 = maxPtJet2; index3 = index2;
232 maxPtJet2 = jetPt; index2 = i;
233 } else if( jetPt > maxPtJet3 ) {
234 maxPtJet3 = jetPt; index3 = i;
235 }
236 }
237 if( index1 != -1 ) {
238 AliAODJet* jet = fAODjets->GetJet(index1);
239 jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
240 }
241 if( index2 != -1 ) {
242 AliAODJet* jet = fAODjets->GetJet(index2);
243 jetVect[1].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
244 }
245 if( index3 != -1 ) {
246 AliAODJet* jet = fAODjets->GetJet(index3);
247 jetVect[2].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
248 }
249
250 } else {
251
252 // Use "Charged Particle Jets"
253 TObjArray* jets = FindChargedParticleJets();
254 if( jets ) {
255 nJets = jets->GetEntriesFast();
256 if( nJets > 0 ) {
257 index1 = 0;
258 AliAODJet* jet = (AliAODJet*)jets->At(0);
259 maxPtJet1 = jet->Pt();
260 jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
261 }
262 if( nJets > 1 ) {
263 index2 = 1;
264 AliAODJet* jet = (AliAODJet*)jets->At(1);
265 maxPtJet1 = jet->Pt();
266 jetVect[1].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
267 }
268 if( nJets > 2 ) {
269 index3 = 2;
270 AliAODJet* jet = (AliAODJet*)jets->At(2);
271 maxPtJet1 = jet->Pt();
272 jetVect[2].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
273 }
274
275 jets->Delete();
276 delete jets;
f3050824 277 }
278 }
6ef5bfa4 279
f3050824 280
281 fhNJets->Fill(nJets);
282
6ef5bfa4 283 if( fDebug > 1 ) {
284 if( index1 < 0 ) {
285 AliInfo("\n Skipping Event, not jet found...");
286 return;
287 } else
288 AliInfo(Form("\n Pt Leading Jet = %6.1f eta=%5.3f ", maxPtJet1, jetVect[0].Eta() ));
289 }
290
f3050824 291
292 // ----------------------------------------------
6ef5bfa4 293 // Cut events by jets topology
f3050824 294 // fAnaType:
6ef5bfa4 295 // 1 = inclusive,
f3050824 296 // - Jet1 |eta| < fJet1EtaCut
297 // 2 = back to back inclusive
298 // - fulfill case 1
299 // - |Jet1.Phi - Jet2.Phi| > fJet2DeltaPhiCut
300 // - Jet2.Pt/Jet1Pt > fJet2RatioPtCut
6ef5bfa4 301 // 3 = back to back exclusive
f3050824 302 // - fulfill case 2
303 // - Jet3.Pt < fJet3PtCut
304
6ef5bfa4 305 if( index1 < 0 || TMath::Abs(jetVect[0].Eta()) > fJet1EtaCut) {
306 if( fDebug > 1 ) AliInfo("\n Skipping Event...Jet1 |eta| > fJet1EtaCut");
f3050824 307 return;
308 }
f3050824 309 // back to back inclusive
6ef5bfa4 310 if( fAnaType > 1 && index2 == -1 ) {
311 if( fDebug > 1 ) AliInfo("\n Skipping Event... no second Jet found");
312 return;
313 }
314 if( fAnaType > 1 && index2 > -1 ) {
f3050824 315 if( TMath::Abs(jetVect[0].DeltaPhi(jetVect[1])) < fJet2DeltaPhiCut ||
6ef5bfa4 316 maxPtJet2/maxPtJet1 < fJet2RatioPtCut ) {
317 if( fDebug > 1 ) AliInfo("\n Skipping Event... |Jet1.Phi - Jet2.Phi| < fJet2DeltaPhiCut");
318 return;
f3050824 319 }
f3050824 320 }
f3050824 321 // back to back exclusive
6ef5bfa4 322 if( fAnaType > 2 && index3 > -1 ) {
323 if( maxPtJet3 > fJet3PtCut ) {
324 if( fDebug > 1 ) AliInfo("\n Skipping Event... Jet3.Pt > fJet3PtCut ");
f3050824 325 return;
326 }
f3050824 327 }
f3050824 328
f3050824 329 fhEleadingPt->Fill( maxPtJet1 );
6ef5bfa4 330 //Area for Normalization Purpose at Display histos
331 SetRegionArea(jetVect);
332
f3050824 333 // ----------------------------------------------
334 // Find max and min regions
335 Double_t sumPtRegionPosit = 0.;
336 Double_t sumPtRegionNegat = 0.;
337 Double_t maxPartPtRegion = 0.;
338 Int_t nTrackRegionPosit = 0;
339 Int_t nTrackRegionNegat = 0;
340 static const Double_t k270rad = 270.*TMath::Pi()/180.;
341
342 Int_t nTracks = fAOD->GetNTracks();
343 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
344 AliAODTrack* part = fAOD->GetTrack( ipart );
6ef5bfa4 345 if ( !part->TestFilterBit(fFilterBit) ) continue; // track cut selection
346 // PID Selection: Reject everything but hadrons
347 Bool_t isHadron = part->GetMostProbablePID()==AliAODTrack::kPion ||
348 part->GetMostProbablePID()==AliAODTrack::kKaon ||
349 part->GetMostProbablePID()==AliAODTrack::kProton;
350 if ( !isHadron ) continue;
351
352 if ( !part->Charge() ) continue; //Only charged
353 if ( fUseSingleCharge ) { // Charge selection
354 if ( fUsePositiveCharge && part->Charge() < 0.) continue; // keep Positives
355 if ( !fUsePositiveCharge && part->Charge() > 0.) continue; // keep Negatives
356 }
357
f3050824 358 if ( part->Pt() < fTrackPtCut ) continue;
6ef5bfa4 359
f3050824 360 TVector3 partVect(part->Px(), part->Py(), part->Pz());
361
362 Double_t deltaPhi = jetVect[0].DeltaPhi(partVect)+k270rad;
363 if( deltaPhi > 2.*TMath::Pi() ) deltaPhi-= 2.*TMath::Pi();
de6f8090 364 fhdNdEtaPhiDist->Fill( deltaPhi );
f3050824 365 fhFullRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 );
6ef5bfa4 366
f3050824 367 Int_t region = IsTrackInsideRegion( jetVect, &partVect );
368
369 if (region > 0) {
6ef5bfa4 370 if( maxPartPtRegion < part->Pt() ) maxPartPtRegion = part->Pt();
371 sumPtRegionPosit += part->Pt();
372 nTrackRegionPosit++;
373 fhTransRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 );
f3050824 374 }
375 if (region < 0) {
6ef5bfa4 376 if( maxPartPtRegion < part->Pt() ) maxPartPtRegion = part->Pt();
377 sumPtRegionNegat += part->Pt();
378 nTrackRegionNegat++;
379 fhTransRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 );
f3050824 380 }
381 }
6ef5bfa4 382
f3050824 383 //How quantities will be sorted before Fill Min and Max Histogram
384 // 1=Plots will be CDF-like
385 // 2=Plots will be Marchesini-like
6ef5bfa4 386 if( fOrdering == 1 ) {
387 if( sumPtRegionPosit > sumPtRegionNegat ) {
388 FillSumPtRegion( maxPtJet1, sumPtRegionPosit/fAreaReg, sumPtRegionNegat/fAreaReg );
f3050824 389 } else {
6ef5bfa4 390 FillSumPtRegion( maxPtJet1, sumPtRegionNegat/fAreaReg, sumPtRegionPosit/fAreaReg );
f3050824 391 }
392 if (nTrackRegionPosit > nTrackRegionNegat ) {
6ef5bfa4 393 FillMultRegion( maxPtJet1, nTrackRegionPosit/fAreaReg, nTrackRegionNegat/fAreaReg, sumPtRegionNegat/fAreaReg );
f3050824 394 } else {
6ef5bfa4 395 FillMultRegion( maxPtJet1, nTrackRegionNegat/fAreaReg, nTrackRegionPosit/fAreaReg, sumPtRegionPosit/fAreaReg );
f3050824 396 }
6ef5bfa4 397 } else if( fOrdering == 2 ) {
f3050824 398 if (sumPtRegionPosit > sumPtRegionNegat) {
6ef5bfa4 399 FillSumPtRegion( maxPtJet1, sumPtRegionPosit/fAreaReg, sumPtRegionNegat/fAreaReg );
400 FillMultRegion( maxPtJet1, nTrackRegionPosit/fAreaReg, nTrackRegionNegat/fAreaReg, sumPtRegionNegat/fAreaReg );
f3050824 401 } else {
6ef5bfa4 402 FillSumPtRegion( maxPtJet1, sumPtRegionNegat/fAreaReg, sumPtRegionPosit/fAreaReg );
403 FillMultRegion( maxPtJet1, nTrackRegionNegat/fAreaReg, nTrackRegionPosit/fAreaReg, sumPtRegionPosit/fAreaReg );
f3050824 404 }
405 }
406
407 Double_t avePosRegion = (nTrackRegionPosit) ? sumPtRegionPosit/nTrackRegionPosit : 0.;
408 Double_t aveNegRegion = (nTrackRegionNegat) ? sumPtRegionNegat/nTrackRegionNegat : 0.;
6ef5bfa4 409 if( avePosRegion > aveNegRegion ) {
410 FillAvePartPtRegion( maxPtJet1, avePosRegion/fAreaReg, aveNegRegion/fAreaReg );
f3050824 411 } else {
6ef5bfa4 412 FillAvePartPtRegion( maxPtJet1, aveNegRegion/fAreaReg, avePosRegion/fAreaReg );
f3050824 413 }
6ef5bfa4 414
f3050824 415 fhRegionMaxPartPtMaxVsEt->Fill(maxPtJet1, maxPartPtRegion );
6ef5bfa4 416
417 // Compute pedestal like magnitudes
418 fhRegionDiffSumPtVsEt->Fill(maxPtJet1, TMath::Abs(sumPtRegionPosit-sumPtRegionNegat)/(2.0*fAreaReg));
419 fhRegionAveSumPtVsEt->Fill(maxPtJet1, (sumPtRegionPosit+sumPtRegionNegat)/(2.0*fAreaReg));
420
f3050824 421}
422
f3050824 423//____________________________________________________________________
424void AliAnalysisTaskUE::FillSumPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin )
425{
6ef5bfa4 426 // Fill sumPt of control regions
427
428 // Max cone
429 fhRegionSumPtMaxVsEt->Fill( leadingE, ptMax );
430 // Min cone
431 fhRegionSumPtMinVsEt->Fill( leadingE, ptMin );
432 // MAke distributions for UE comparison with MB data
433 fhMinRegSumPt->Fill(ptMin);
f3050824 434}
435
436//____________________________________________________________________
437void AliAnalysisTaskUE::FillAvePartPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin )
438{
6ef5bfa4 439 // Fill average particle Pt of control regions
440
441 // Max cone
442 fhRegionAvePartPtMaxVsEt->Fill( leadingE, ptMax );
443 // Min cone
444 fhRegionAvePartPtMinVsEt->Fill( leadingE, ptMin );
445 // MAke distributions for UE comparison with MB data
446 fhMinRegAvePt->Fill(ptMin);
f3050824 447}
448
449//____________________________________________________________________
450void AliAnalysisTaskUE::FillMultRegion( Double_t leadingE, Double_t nTrackPtmax, Double_t nTrackPtmin, Double_t ptMin )
451{
6ef5bfa4 452 // Fill Nch multiplicity of control regions
453
454 // Max cone
455 fhRegionMultMaxVsEt->Fill( leadingE, nTrackPtmax );
456 fhRegionMultMax->Fill( nTrackPtmax );
457 // Min cone
458 fhRegionMultMinVsEt->Fill( leadingE, nTrackPtmin );
459 fhRegionMultMin->Fill( nTrackPtmin );
460 // MAke distributions for UE comparison with MB data
461 fhMinRegSumPtvsMult->Fill(nTrackPtmin,ptMin);
f3050824 462}
463
464//____________________________________________________________________
465Int_t AliAnalysisTaskUE::IsTrackInsideRegion(TVector3 *jetVect, TVector3 *partVect)
466{
467 // return de region in delta phi
468 // -1 negative delta phi
469 // 1 positive delta phi
470 // 0 outside region
471 static const Double_t k60rad = 60.*TMath::Pi()/180.;
472 static const Double_t k120rad = 120.*TMath::Pi()/180.;
473
474 Int_t region = 0;
475 if( fRegionType == 1 ) {
476 if( TMath::Abs(partVect->Eta()) > fTrackEtaCut ) return 0;
477 // transverse regions
478 if (jetVect[0].DeltaPhi(*partVect) < -k60rad && jetVect[0].DeltaPhi(*partVect) > -k120rad ) region = -1;
479 if (jetVect[0].DeltaPhi(*partVect) > k60rad && jetVect[0].DeltaPhi(*partVect) < k120rad ) region = 1;
6ef5bfa4 480
f3050824 481 } else if( fRegionType == 2 ) {
482 // Cone regions
483 Double_t deltaR = 0.;
484
485 TVector3 positVect,negatVect;
486 positVect.SetMagThetaPhi(1, 2.*atan(exp(-jetVect[0].Eta())), jetVect[0].Phi()+TMath::PiOver2());
487 negatVect.SetMagThetaPhi(1, 2.*atan(exp(-jetVect[0].Eta())), jetVect[0].Phi()-TMath::PiOver2());
488
489 if (TMath::Abs(positVect.DeltaPhi(*partVect)) < fConeRadius ) {
6ef5bfa4 490 region = 1;
491 deltaR = positVect.DrEtaPhi(*partVect);
f3050824 492 } else if (TMath::Abs(negatVect.DeltaPhi(*partVect)) < fConeRadius) {
6ef5bfa4 493 region = -1;
494 deltaR = negatVect.DrEtaPhi(*partVect);
f3050824 495 }
496
497 if (deltaR > fConeRadius) region = 0;
6ef5bfa4 498
f3050824 499 } else
500 AliError("Unknow region type");
6ef5bfa4 501
f3050824 502 // For debug (to be removed)
503 //if( region != 0 ) fhValidRegion->Fill( partVect->Eta()-jetVect[0].Eta(), jetVect[0].DeltaPhi(*partVect) );
504
505 return region;
506}
507
6ef5bfa4 508
509//____________________________________________________________________
510TObjArray* AliAnalysisTaskUE::FindChargedParticleJets()
511{
512 // Return a TObjArray of "charged particle jets"
513 //
514 // Charged particle jet definition from reference:
515 //
516 // "Charged jet evolution and the underlying event
517 // in proton-antiproton collisions at 1.8 TeV"
518 // PHYSICAL REVIEW D 65 092002, CDF Collaboration
519 //
520 // We define "jets" as circular regions in eta-phi space with
521 // radius defined by R = sqrt( (eta-eta0)^2 +(phi-phi0)^2 ).
522 // Our jet algorithm is as follows:
523 // 1- Order all charged particles according to their pT .
524 // 2- Start with the highest pT particle and include in the jet all
525 // particles within the radius R=0.7 considering each particle
526 // in the order of decreasing pT and recalculating the centroid
527 // of the jet after each new particle is added to the jet .
528 // 3- Go to the next highest pT particle not already included in
529 // a jet and add to the jet all particles not already included in
530 // a jet within R=0.7.
531 // 4- Continue until all particles are in a jet.
532 // We define the transverse momentum of the jet to be
533 // the scalar pT sum of all the particles within the jet, where pT
534 // is measured with respect to the beam axis
535
536 // 1 - Order all charged particles according to their pT .
537 Int_t nTracks = fAOD->GetNTracks();
538 if( !nTracks ) return 0;
539 TObjArray tracks(nTracks);
540
541 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
542 AliAODTrack* part = fAOD->GetTrack( ipart );
543 if( !part->TestFilterBit(fFilterBit) ) continue; // track cut selection
544 if( !part->Charge() ) continue;
545 if( part->Pt() < fTrackPtCut ) continue;
546 tracks.AddLast(part);
547 }
548 QSortTracks( tracks, 0, tracks.GetEntriesFast() );
549
550 nTracks = tracks.GetEntriesFast();
551 if( !nTracks ) return 0;
552 TObjArray *jets = new TObjArray(nTracks);
553 TIter itrack(&tracks);
554 while( nTracks ) {
555 // 2- Start with the highest pT particle ...
556 Float_t px,py,pz,pt;
557 AliAODTrack* track = (AliAODTrack*)itrack.Next();
558 if( !track ) continue;
559 px = track->Px();
560 py = track->Py();
561 pz = track->Pz();
562 pt = track->Pt(); // Use the energy member to store Pt
563 jets->AddLast( new TLorentzVector(px, py, pz, pt) );
564 tracks.Remove( track );
565 TLorentzVector* jet = (TLorentzVector*)jets->Last();
566 // 3- Go to the next highest pT particle not already included...
567 AliAODTrack* track1;
568 while ( (track1 = (AliAODTrack*)(itrack.Next())) ) {
569 Double_t dphi = TVector2::Phi_mpi_pi(jet->Phi()-track1->Phi());
570 Double_t r = TMath::Sqrt( (jet->Eta()-track1->Eta())*(jet->Eta()-track1->Eta()) +
571 dphi*dphi );
572 if( r < fConeRadius ) {
de6f8090 573 Double_t fPt = jet->E()+track1->Pt(); // Scalar sum of Pt
6ef5bfa4 574 // recalculating the centroid
de6f8090 575 Double_t eta = jet->Eta()*jet->E()/fPt + track1->Eta()/fPt;
576 Double_t phi = jet->Phi()*jet->E()/fPt + track1->Phi()/fPt;
577 jet->SetPtEtaPhiE( 1., eta, phi, fPt );
6ef5bfa4 578 tracks.Remove( track1 );
579 }
580 }
581
582 tracks.Compress();
583 nTracks = tracks.GetEntries();
584 // 4- Continue until all particles are in a jet.
585 itrack.Reset();
586 } // end while nTracks
587
588 // Convert to AODjets....
589 Int_t njets = jets->GetEntriesFast();
590 TObjArray* aodjets = new TObjArray(njets);
591 aodjets->SetOwner(kTRUE);
592 for(Int_t ijet=0; ijet<njets; ++ijet) {
593 TLorentzVector* jet = (TLorentzVector*)jets->At(ijet);
594 Float_t px, py,pz,en; // convert to 4-vector
595 px = jet->E() * TMath::Cos(jet->Phi()); // Pt * cos(phi)
596 py = jet->E() * TMath::Sin(jet->Phi()); // Pt * sin(phi)
597 pz = jet->E() / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-jet->Eta())));
598 en = TMath::Sqrt(px * px + py * py + pz * pz);
599 aodjets->AddLast( new AliAODJet(px, py, pz, en) );
600 }
601 // Order jets according to their pT .
602 QSortTracks( *aodjets, 0, aodjets->GetEntriesFast() );
603
604 // debug
605 if (fDebug>3) AliInfo(Form(" %d Charged jets found\n",njets));
606
607 return aodjets;
608}
609
610//____________________________________________________________________
611void AliAnalysisTaskUE::QSortTracks(TObjArray &a, Int_t first, Int_t last)
612{
613 // Sort array of TObjArray of tracks by Pt using a quicksort algorithm.
614
615 static TObject *tmp;
616 static int i; // "static" to save stack space
617 int j;
618
619 while (last - first > 1) {
620 i = first;
621 j = last;
622 for (;;) {
623 while (++i < last && ((AliVParticle*)a[i])->Pt() > ((AliVParticle*)a[first])->Pt() )
624 ;
625 while (--j > first && ((AliVParticle*)a[j])->Pt() < ((AliVParticle*)a[first])->Pt() )
626 ;
627 if (i >= j)
628 break;
629
630 tmp = a[i];
631 a[i] = a[j];
632 a[j] = tmp;
633 }
634 if (j == first) {
635 ++first;
636 continue;
637 }
638 tmp = a[first];
639 a[first] = a[j];
640 a[j] = tmp;
641 if (j - first < last - (j + 1)) {
642 QSortTracks(a, first, j);
643 first = j + 1; // QSortTracks(j + 1, last);
644 } else {
645 QSortTracks(a, j + 1, last);
646 last = j; // QSortTracks(first, j);
647 }
648 }
649}
650
651//____________________________________________________________________
652void AliAnalysisTaskUE::SetRegionArea(TVector3 *jetVect)
653{
654 Double_t fAreaCorrFactor=0.;
655 Double_t deltaEta = 0.;
656 if (fRegionType==1) fAreaReg = 2.*fTrackEtaCut*60.*TMath::Pi()/180.;
657 else if (fRegionType==2){
658 deltaEta = 0.9-TMath::Abs(jetVect[0].Eta());
659 if (deltaEta>fConeRadius) fAreaReg = TMath::Pi()*fConeRadius*fConeRadius;
660 else{
661 fAreaCorrFactor = fConeRadius*fConeRadius*TMath::ACos(deltaEta/fConeRadius) -
662 (deltaEta*fConeRadius)*TMath::Sqrt( 1. - deltaEta*deltaEta/(fConeRadius*fConeRadius));
663 fAreaReg=TMath::Pi()*fConeRadius*fConeRadius-fAreaCorrFactor;
664 }
665 }else AliWarning("Unknown Rgion Type");
666 if (fDebug>5) 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));
667}
668
f3050824 669//____________________________________________________________________
670void AliAnalysisTaskUE::CreateHistos()
671{
672 fListOfHistos = new TList();
673
674
675 fhNJets = new TH1F("fhNJets", "n Jet", 10, 0, 10);
676 fhNJets->SetXTitle("# of jets");
677 fhNJets->Sumw2();
f3050824 678 fListOfHistos->Add( fhNJets ); // At(0)
679
680 fhEleadingPt = new TH1F("hEleadingPt", "Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
681 fhEleadingPt->SetXTitle("P_{T} (GeV/c)");
682 fhEleadingPt->SetYTitle("dN/dP_{T}");
683 fhEleadingPt->Sumw2();
684 fListOfHistos->Add( fhEleadingPt ); // At(1)
685
686 fhMinRegPtDist = new TH1F("hMinRegPtDist", "P_{T} distribution in Min zone", 50,0.,20.);
687 fhMinRegPtDist->SetXTitle("P_{T} (GeV/c)");
688 fhMinRegPtDist->SetYTitle("dN/dP_{T}");
689 fhMinRegPtDist->Sumw2();
690 fListOfHistos->Add( fhMinRegPtDist ); // At(2)
691
692 fhRegionMultMin = new TH1F("hRegionMultMin", "N_{ch}^{90, min}", 21, -0.5, 20.5);
693 fhRegionMultMin->SetXTitle("N_{ch tracks}");
694 fhRegionMultMin->Sumw2();
695 fListOfHistos->Add( fhRegionMultMin ); // At(3)
696
697 fhMinRegAvePt = new TH1F("hMinRegAvePt", "#LTp_{T}#GT", 50, 0., 20.);
698 fhMinRegAvePt->SetXTitle("P_{T} (GeV/c)");
699 fhMinRegAvePt->Sumw2();
700 fListOfHistos->Add( fhMinRegAvePt ); // At(4)
701
702 fhMinRegSumPt = new TH1F("hMinRegSumPt", "#Sigma p_{T} ", 50, 0., 20.);
703 fhMinRegSumPt->SetYTitle("Ed^{3}N_{tracks}/dp^{3} (c^{3}/GeV^{2})");
704 fhMinRegSumPt->SetXTitle("#Sigma p_{T} (GeV/c)");
705 fhMinRegSumPt->Sumw2();
706 fListOfHistos->Add( fhMinRegSumPt ); // At(5)
6ef5bfa4 707
f3050824 708 fhMinRegMaxPtPart = new TH1F("hMinRegMaxPtPart", "max(p_{T})|_{event} ", 50, 0., 20.);
709 fhMinRegMaxPtPart->SetYTitle("Ed^{3}N_{tracks}/dp^{3} (c^{3}/GeV^{2})");
710 fhMinRegMaxPtPart->SetXTitle("p_{T} (GeV/c)");
711 fhMinRegMaxPtPart->Sumw2();
712 fListOfHistos->Add( fhMinRegMaxPtPart ); // At(6)
713
714 fhMinRegSumPtvsMult = new TH1F("hMinRegSumPtvsMult", "#Sigma p_{T} vs. Multiplicity ", 21, -0.5, 20.5);
715 fhMinRegSumPtvsMult->SetYTitle("#Sigma p_{T} (GeV/c)");
716 fhMinRegSumPtvsMult->SetXTitle("N_{charge}");
717 fhMinRegSumPtvsMult->Sumw2();
718 fListOfHistos->Add( fhMinRegSumPtvsMult ); // At(7);
6ef5bfa4 719
de6f8090 720 fhdNdEtaPhiDist = new TH1F("hdNdEtaPhiDist", "Charge particle density |#eta|< 1 vs #Delta#phi", 120, 0., 2.*TMath::Pi());
721 fhdNdEtaPhiDist->SetXTitle("#Delta#phi");
722 fhdNdEtaPhiDist->SetYTitle("dN_{ch}/d#etad#phi");
723 fhdNdEtaPhiDist->Sumw2();
724 fListOfHistos->Add( fhdNdEtaPhiDist ); // At(8)
6ef5bfa4 725
f3050824 726 // Can be use to get part pt distribution for differente Jet Pt bins
727 fhFullRegPartPtDistVsEt = new TH2F("hFullRegPartPtDistVsEt", "dN/dP_{T} |#eta|<1 vs Leading Jet P_{T}",
6ef5bfa4 728 50,0.,50., fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
f3050824 729 fhFullRegPartPtDistVsEt->SetYTitle("Leading Jet P_{T}");
730 fhFullRegPartPtDistVsEt->SetXTitle("p_{T}");
731 fhFullRegPartPtDistVsEt->Sumw2();
732 fListOfHistos->Add( fhFullRegPartPtDistVsEt ); // At(9)
6ef5bfa4 733
734 // Can be use to get part pt distribution for differente Jet Pt bins
f3050824 735 fhTransRegPartPtDistVsEt = new TH2F("hTransRegPartPtDistVsEt", "dN/dP_{T} in tranvese regions |#eta|<1 vs Leading Jet P_{T}",
6ef5bfa4 736 50,0.,50., fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
f3050824 737 fhTransRegPartPtDistVsEt->SetYTitle("Leading Jet P_{T}");
738 fhTransRegPartPtDistVsEt->SetXTitle("p_{T}");
739 fhTransRegPartPtDistVsEt->Sumw2();
740 fListOfHistos->Add( fhTransRegPartPtDistVsEt ); // At(10)
741
742
743 fhRegionSumPtMaxVsEt = new TH1F("hRegionSumPtMaxVsEt", "P_{T}^{90, max} vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
744 fhRegionSumPtMaxVsEt->SetXTitle("P_{T} (GeV/c)");
745 fhRegionSumPtMaxVsEt->Sumw2();
746 fListOfHistos->Add( fhRegionSumPtMaxVsEt ); // At(11)
6ef5bfa4 747
f3050824 748 fhRegionSumPtMinVsEt = new TH1F("hRegionSumPtMinVsEt", "P_{T}^{90, min} vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
749 fhRegionSumPtMinVsEt->SetXTitle("P_{T} (GeV/c)");
750 fhRegionSumPtMinVsEt->Sumw2();
751 fListOfHistos->Add( fhRegionSumPtMinVsEt ); // At(12)
752
753 fhRegionMultMax = new TH1I("hRegionMultMax", "N_{ch}^{90, max}", 21, -0.5, 20.5);
754 fhRegionMultMax->SetXTitle("N_{ch tracks}");
755 fhRegionMultMax->Sumw2();
756 fListOfHistos->Add( fhRegionMultMax ); // At(13)
6ef5bfa4 757
f3050824 758 fhRegionMultMaxVsEt = new TH1F("hRegionMultMaxVsEt", "N_{ch}^{90, max} vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
759 fhRegionMultMaxVsEt->SetXTitle("E (GeV hRegionAveSumPtVsEt/c)");
760 fhRegionMultMaxVsEt->Sumw2();
761 fListOfHistos->Add( fhRegionMultMaxVsEt ); // At(14)
6ef5bfa4 762
f3050824 763 fhRegionMultMinVsEt = new TH1F("hRegionMultMinVsEt", "N_{ch}^{90, min} vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
764 fhRegionMultMinVsEt->SetXTitle("E (GeV/c)");
765 fhRegionMultMinVsEt->Sumw2();
766 fListOfHistos->Add( fhRegionMultMinVsEt ); // At(15)
6ef5bfa4 767
f3050824 768 fhRegionAveSumPtVsEt = new TH1F("hRegionAveSumPtVsEt", "(P_{T}^{90, max} + P_{T}^{90, min})/2 vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
769 fhRegionAveSumPtVsEt->SetXTitle("P_{T} (GeV/c)");
770 fhRegionAveSumPtVsEt->Sumw2();
771 fListOfHistos->Add( fhRegionAveSumPtVsEt ); // At(16)
6ef5bfa4 772
f3050824 773 fhRegionDiffSumPtVsEt= new TH1F("hRegionPtDiffVsEt", "(P_{T}^{90, max} - P_{T}^{90, min}) vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
774 fhRegionDiffSumPtVsEt->SetXTitle("P_{T} (GeV/c)");
775 fhRegionDiffSumPtVsEt->Sumw2();
776 fListOfHistos->Add( fhRegionDiffSumPtVsEt ); // At(17)
777
778 fhRegionAvePartPtMaxVsEt = new TH1F("hRegionAvePartPtMaxVsEt", "#LTp_{T}#GT^{90, max} vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
779 fhRegionAvePartPtMaxVsEt->SetXTitle("P_{T} (GeV/c)");
780 fhRegionAvePartPtMaxVsEt->Sumw2();
781 fListOfHistos->Add( fhRegionAvePartPtMaxVsEt ); // At(18)
6ef5bfa4 782
f3050824 783 fhRegionAvePartPtMinVsEt = new TH1F("hRegionAvePartPtMinVsEt", "#LTp_{T}#GT^{90, min} vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
784 fhRegionAvePartPtMinVsEt->SetXTitle("P_{T} (GeV/c)");
785 fhRegionAvePartPtMinVsEt->Sumw2();
786 fListOfHistos->Add( fhRegionAvePartPtMinVsEt ); // At(19)
6ef5bfa4 787
f3050824 788 fhRegionMaxPartPtMaxVsEt = new TH1F("hRegionMaxPartPtMaxVsEt", "max(p_{T})^{90} vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
789 fhRegionMaxPartPtMaxVsEt->SetXTitle("P_{T} (GeV/c)");
790 fhRegionMaxPartPtMaxVsEt->Sumw2();
791 fListOfHistos->Add( fhRegionMaxPartPtMaxVsEt ); // At(20)
6ef5bfa4 792
793 fSettingsTree = new TTree("UEAnalysisSettings","Analysis Settings in UE estimation");
794 fListOfHistos->Add( fSettingsTree ); // At(21)
795
796 /*
797 // For debug region selection
798 fhValidRegion = new TH2F("hValidRegion", "dN/d#eta/d#phi",
799 20, -2.,2., 62, -TMath::Pi(), TMath::Pi());
800 fhValidRegion->SetYTitle("#Delta#phi");
801 fhValidRegion->Sumw2();
802 fListOfHistos->Add( fhValidRegion ); // At(15)
803 */
f3050824 804}
f3050824 805
6ef5bfa4 806//____________________________________________________________________
807void AliAnalysisTaskUE::Terminate(Option_t */*option*/)
808{
809 // Terminate analysis
810 //
811
812 //Save Analysis Settings
813 WriteSettings();
814
815 // Normalize histos to region area TODO:
816 // Normalization done at Analysis, taking into account
817 // area variations on per-event basis (cone case)
818
819 // Draw some Test plot to the screen
820 if (!gROOT->IsBatch()) {
821
822 // Update pointers reading them from the output slot
823 fListOfHistos = dynamic_cast<TList*> (GetOutputData(0));
824 if( !fListOfHistos ) {
825 AliError("Histogram List is not available");
826 return;
827 }
828 fhEleadingPt = (TH1F*)fListOfHistos->At(1);
de6f8090 829 fhdNdEtaPhiDist = (TH1F*)fListOfHistos->At(8);
6ef5bfa4 830 fhRegionSumPtMaxVsEt = (TH1F*)fListOfHistos->At(11);
831 fhRegionSumPtMinVsEt = (TH1F*)fListOfHistos->At(12);
832 fhRegionMultMaxVsEt = (TH1F*)fListOfHistos->At(14);
833 fhRegionMultMinVsEt = (TH1F*)fListOfHistos->At(15);
834 fhRegionAveSumPtVsEt = (TH1F*)fListOfHistos->At(16);
835
836 fhNJets = (TH1F*)fListOfHistos->At(0);
837
838 // Canvas
839 TCanvas* c1 = new TCanvas("c1",Form("sumPt dist (%s)", GetTitle()),60,60,1200,800);
840 c1->Divide(2,2);
841 c1->cd(1);
842 TH1F *h1r = new TH1F("hRegionEtvsSumPtMax" , "", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
843 h1r->Divide(fhRegionSumPtMaxVsEt,fhEleadingPt,1,1);
844 //h1r->Scale( areafactor );
845 h1r->SetMarkerStyle(20);
846 h1r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
847 h1r->SetYTitle("P_{T}^{90, max}");
848 h1r->DrawCopy("p");
849
850 c1->cd(2);
851
852 TH1F *h2r = new TH1F("hRegionEtvsSumPtMin" , "", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
853 h2r->Divide(fhRegionSumPtMinVsEt,fhEleadingPt,1,1);
854 //h2r->Scale( areafactor );
855 h2r->SetMarkerStyle(20);
856 h2r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
857 h2r->SetYTitle("P_{T}^{90, min}");
858 h2r->DrawCopy("p");
859
860 c1->cd(3);
861 TH1F *h4r = new TH1F("hRegionEtvsDiffPt" , "", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
862 h4r->Divide(fhRegionAveSumPtVsEt,fhEleadingPt,1,1);
863 h4r->Scale(2.); // make average
864 //h4r->Scale( areafactor );
865 h4r->SetYTitle("#DeltaP_{T}^{90}");
866 h4r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
867 h4r->SetMarkerStyle(20);
868 h4r->DrawCopy("p");
869
870 c1->cd(4);
871 TH1F *h5r = new TH1F("hRegionMultMaxVsEtleading", "", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
872 TH1F *h6r = new TH1F("hRegionMultMinVsEtleading", "", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
873
874 h5r->Divide(fhRegionMultMaxVsEt,fhEleadingPt,1,1);
875 h6r->Divide(fhRegionMultMinVsEt,fhEleadingPt,1,1);
876 //h5r->Scale( areafactor );
877 //h6r->Scale( areafactor );
878 h5r->SetYTitle("N_{Tracks}^{90}");
879 h5r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
880 h5r->SetMarkerStyle(20);
881 h5r->DrawCopy("p");
882 h6r->SetMarkerStyle(21);
883 h6r->SetMarkerColor(2);
884 h6r->SetYTitle("N_{Tracks}^{90}");
885 h6r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
886 h6r->DrawCopy("p same");
887 c1->Update();
888
889 TCanvas* c2 = new TCanvas("c2","Jet Pt dist",160,160,1200,800);
890 c2->Divide(2,2);
891 c2->cd(1);
892 fhEleadingPt->SetMarkerStyle(20);
893 fhEleadingPt->SetMarkerColor(2);
894 fhEleadingPt->DrawCopy("p");
895 gPad->SetLogy();
896
897 c2->cd(2);
de6f8090 898 fhdNdEtaPhiDist->SetMarkerStyle(20);
899 fhdNdEtaPhiDist->SetMarkerColor(2);
900 fhdNdEtaPhiDist->DrawCopy("p");
6ef5bfa4 901 gPad->SetLogy();
902
903 c2->cd(3);
904 fhNJets->DrawCopy();
905
906 //fhTransRegPartPtDist = (TH1F*)fListOfHistos->At(2);
907 fhRegionMultMin = (TH1F*)fListOfHistos->At(3);
908 fhMinRegAvePt = (TH1F*)fListOfHistos->At(4);
909 fhMinRegSumPt = (TH1F*)fListOfHistos->At(5);
910 //fhMinRegMaxPtPart = (TH1F*)fListOfHistos->At(6);
911 fhMinRegSumPtvsMult = (TH1F*)fListOfHistos->At(7);
912
913 // Canvas
914 TCanvas* c3 = new TCanvas("c3"," p_{T} dist",160,160,1200,800);
915 c3->Divide(2,2);
916 c3->cd(1);
917 /*fhTransRegPartPtDist->SetMarkerStyle(20);
918 fhTransRegPartPtDist->SetMarkerColor(2);
919 fhTransRegPartPtDist->Scale(areafactor/fhTransRegPartPtDist->GetEntries());
920 fhTransRegPartPtDist->DrawCopy("p");
921 gPad->SetLogy();
922 */
923 c3->cd(2);
924 fhMinRegSumPt->SetMarkerStyle(20);
925 fhMinRegSumPt->SetMarkerColor(2);
926 //fhMinRegSumPt->Scale(areafactor);
927 fhMinRegSumPt->DrawCopy("p");
928 gPad->SetLogy();
929
930 c3->cd(3);
931 fhMinRegAvePt->SetMarkerStyle(20);
932 fhMinRegAvePt->SetMarkerColor(2);
933 //fhMinRegAvePt->Scale(areafactor);
934 fhMinRegAvePt->DrawCopy("p");
935 gPad->SetLogy();
936
937 c3->cd(4);
938
939 TH1F *h7r = new TH1F("hRegionMultMinVsMult", "", 21, -0.5, 20.5);
940
941 h7r->Divide(fhMinRegSumPtvsMult,fhRegionMultMin,1,1);
942
943 h7r->SetMarkerStyle(20);
944 h7r->SetMarkerColor(2);
945 h7r->DrawCopy("p");
946
947 c3->Update();
948
949
950 /* c2->cd(3);
951 fhValidRegion->SetMarkerStyle(20);
952 fhValidRegion->SetMarkerColor(2);
953 fhValidRegion->DrawCopy("p");
954 */
955 c2->Update();
956 } else {
957 AliInfo(" Batch mode, not histograms will be shown...");
958 }
959
960 if( fDebug > 1 )
961 AliInfo("End analysis");
962
963}
f3050824 964
6ef5bfa4 965void AliAnalysisTaskUE::WriteSettings()
966{
967
968 fSettingsTree->Branch("fAnaType", &fAnaType, "Ana/I");
969 fSettingsTree->Branch("fRegionType", &fRegionType,"Reg/I");
970 fSettingsTree->Branch("fConeRadius", &fConeRadius,"Rad/D");
971 fSettingsTree->Branch("fUseChPartJet", &fUseChPartJet,"UseChPart/O");
972 fSettingsTree->Branch("fUseSingleCharge", &fUseSingleCharge,"UseSingleCh/O");
973 fSettingsTree->Branch("fUsePositiveCharge", &fUsePositiveCharge,"UsePositiveCh/O");
974 fSettingsTree->Branch("fOrdering", &fOrdering,"OrderMeth/I");
975 fSettingsTree->Branch("fJet1EtaCut", &fJet1EtaCut, "LeadJetEtaCut/D");
976 fSettingsTree->Branch("fJet2DeltaPhiCut", &fJet2DeltaPhiCut, "DeltaPhi/D");
977 fSettingsTree->Branch("fJet2RatioPtCut", &fJet2RatioPtCut, "Jet2Ratio/D");
978 fSettingsTree->Branch("fJet3PtCut", &fJet3PtCut, "Jet3PtCut/D");
979 fSettingsTree->Branch("fTrackPtCut", &fTrackPtCut, "TrackPtCut/D");
980 fSettingsTree->Branch("fTrackEtaCut", &fTrackEtaCut, "TrackEtaCut/D");
981 fSettingsTree->Fill();
982
983 if (fDebug>5){
984 AliInfo(" All Analysis Settings in Saved Tree");
985 fSettingsTree->Scan();
986 }
987}