PWG4/JetTasks/AliAnalysisTaskUE.cxx

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