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