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