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