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