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