PWG4/JetTasks/AliAnalysisTaskUE.cxx

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