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 <TChain.h>
  20 #include <TFile.h>
  21 #include <TList.h>
  22 #include <TMath.h>
  23 #include <TProfile.h>
  24 #include <TTree.h>
  25 #include <TVector3.h>
  26 #include <TH1F.h>
  27
  28 #include "AliAnalyseUE.h"
  29 #include "AliHistogramsUE.h"
  30 #include "AliAnalysisTaskUE.h"
  31 #include "AliHistogramsUE.h"
  32
  33 #include "AliAnalysisManager.h"
  34 #include "AliMCEventHandler.h"
  35
  36 #include "AliAnalysisHelperJetTasks.h"
  37 #include "AliAODHandler.h"
  38 #include "AliAODInputHandler.h"
  39 #include "AliGenPythiaEventHeader.h"
  40 #include "AliLog.h"
  41 #include "AliInputEventHandler.h"
  42
  43
  44 ////////////////////////////////////////////////////////////////////////
  45 //
  46 // Analysis class for Underlying Event studies
  47 //
  48 // Look for correlations on the tranverse regions to
  49 // the leading charged jet
  50 //
  51 // This class needs as input AOD with track and Jets.
  52 // The output is a list of histograms
  53 //
  54 // AOD can be either connected to the InputEventHandler
  55 // for a chain of AOD files
  56 // or
  57 // to the OutputEventHandler
  58 // for a chain of ESD files, so this case class should be
  59 // in the train after the Jet finder
  60 //
  61 //    Arian.Abrahantes.Quintana@cern.ch
  62 //    Ernesto.Lopez.Torres@cern.ch
  63 //    vallero@physi.uni-heidelberg.de
  64 //
  65 ////////////////////////////////////////////////////////////////////////
  66
  67 ClassImp( AliAnalysisTaskUE)
  68
  69 // Define global pointer
  70 AliAnalysisTaskUE* AliAnalysisTaskUE::fgTaskUE=NULL;
  71
  72 //____________________________________________________________________
  73 AliAnalysisTaskUE:: AliAnalysisTaskUE(const char* name):
  74 AliAnalysisTask(name,""),
  75 fHistosUE(0x0),
  76 fAnaUE(0x0),
  77 fAOD(0x0),
  78 fAODBranch("jets"),
  79 fDebug(0),
  80 fListOfHistos(0x0),
  81 //Configuration
  82 fBinsPtInHist(30),
  83 fIsNorm2Area(kTRUE),
  84 fMaxJetPtInHist(300.),
  85 fMinJetPtInHist(0.),
  86 fConstrainDistance(kTRUE),
  87 fMinDistance(0.2),
  88 fSimulateChJetPt(kFALSE),
  89 fUseAliStack(kTRUE),
  90 fUseMCParticleBranch(kFALSE),
  91 fnTracksVertex(3),  // QA tracks pointing to principal vertex (= 3 default)
  92 fZVertex(5.),
  93 fAnaType(1),
  94 fConePosition(1),
  95 fConeRadius(0.7),
  96 fFilterBit(0xFF),
  97 fJetsOnFly(kFALSE),
  98 fRegionType(1),
  99 fUseChargeHadrons(kFALSE),
 100 fUseChPartJet(kFALSE),
 101 fUsePositiveCharge(kTRUE),
 102 fUseSingleCharge(kFALSE),
 103 fOrdering(1),
 104 fChJetPtMin(5.0),
 105 fJet1EtaCut(0.2),
 106 fJet2DeltaPhiCut(2.616),    // 150 degrees
 107 fJet2RatioPtCut(0.8),
 108 fJet3PtCut(15.),
 109 fTrackEtaCut(0.9),
 110 fTrackPtCut(0.),
 111 //For MC
 112 fAvgTrials(1)
 113 {
 114   // Default constructor
 115   // Define input and output slots here
 116   // Input slot #0 works with a TChain
 117   DefineInput(0, TChain::Class());
 118   // Output slot #0 writes into a TList container
 119   DefineOutput(0, TList::Class());
 120
 121 }
 122
 123 //____________________________________________________________________
 124 AliAnalysisTaskUE:: AliAnalysisTaskUE(const AliAnalysisTaskUE & original):
 125 AliAnalysisTask(),
 126 fHistosUE(original.fHistosUE),
 127 fAnaUE(original.fAnaUE),
 128 fAOD(original.fAOD),
 129 fAODBranch(original.fAODBranch),
 130 fDebug(original.fDebug),
 131 fListOfHistos(original.fListOfHistos),
 132 //Configuration
 133 fBinsPtInHist(original.fBinsPtInHist),
 134 fIsNorm2Area(original.fIsNorm2Area),
 135 fMaxJetPtInHist(original.fMaxJetPtInHist),
 136 fMinJetPtInHist(original.fMinJetPtInHist),
 137 fConstrainDistance(original.fConstrainDistance),
 138 fMinDistance(original.fMinDistance),
 139 fSimulateChJetPt(original.fSimulateChJetPt),
 140 fUseAliStack(original.fUseAliStack),
 141 fUseMCParticleBranch(original.fUseMCParticleBranch),
 142 fnTracksVertex(original.fnTracksVertex),  // QA tracks pointing to principal vertex (= 3 default)
 143 fZVertex(original.fZVertex),
 144 fAnaType(original.fAnaType),
 145 fConePosition(original.fConePosition),
 146 fConeRadius(original.fConeRadius),
 147 fFilterBit(original.fFilterBit),
 148 fJetsOnFly(original.fJetsOnFly),
 149 fRegionType(original.fRegionType),
 150 fUseChargeHadrons(original.fUseChargeHadrons),
 151 fUseChPartJet(original.fUseChPartJet),
 152 fUsePositiveCharge(original.fUsePositiveCharge),
 153 fUseSingleCharge(original.fUseSingleCharge),
 154 fOrdering(original.fOrdering),
 155 fChJetPtMin(original.fChJetPtMin),
 156 fJet1EtaCut(original.fJet1EtaCut),
 157 fJet2DeltaPhiCut(original.fJet2DeltaPhiCut),    // 150 degrees
 158 fJet2RatioPtCut(original.fJet2RatioPtCut),
 159 fJet3PtCut(original.fJet3PtCut),
 160 fTrackEtaCut(original.fTrackEtaCut),
 161 fTrackPtCut(original.fTrackPtCut),
 162 //For MC
 163 fAvgTrials(original.fAvgTrials)
 164 {
 165   // Copy constructor
 166 }
 167
 168
 169 //______________________________________________________________
 170 AliAnalysisTaskUE & AliAnalysisTaskUE::operator = (const AliAnalysisTaskUE & /*source*/)
 171 {
 172   // assignment operator
 173   return *this;
 174 }
 175
 176 //______________________________________________________________
 177 Bool_t AliAnalysisTaskUE::Notify()
 178 {
 179   //
 180   // Implemented Notify() to read the cross sections
 181   // and number of trials from pyxsec.root
 182   // Copy from AliAnalysisTaskJFSystematics
 183   fAvgTrials = 1;
 184   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 185   Float_t xsection = 0;
 186   Float_t trials  = 1;
 187   if(tree){
 188         TFile *curfile = tree->GetCurrentFile();
 189         if (!curfile) {
 190                 Error("Notify","No current file");
 191                 return kFALSE;
 192                 }
 193
 194         AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,trials);
 195
 196         fHistosUE->GetXsec()->Fill("<#sigma>",xsection);
 197
 198         // construct average trials
 199         Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
 200         if(trials>=nEntries && nEntries>0.)fAvgTrials = trials/nEntries;
 201         }
 202
 203   return kTRUE;
 204 }
 205
 206 //____________________________________________________________________
 207 void AliAnalysisTaskUE::ConnectInputData(Option_t* /*option*/)
 208 {
 209   // Connect the input data
 210
 211   // We need AODs with tracks and jets.
 212   // Since AODs can either be connected to the InputEventHandler
 213   // or to the OutputEventHandler ( the AOD is created by a previus task in the train )
 214   // we need to check where it is and get the pointer to AODEvent in the right way
 215
 216   // Delta AODs are also accepted
 217
 218
 219   if (fDebug > 1) AliInfo("ConnectInputData() ");
 220
 221   TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
 222
 223   if( handler && handler->InheritsFrom("AliAODInputHandler") ) { // input AOD
 224         fAOD = ((AliAODInputHandler*)handler)->GetEvent();
 225         if(!fJetsOnFly){
 226                 if (fDebug > 1) AliInfo(" ==== Tracks and Jets from AliAODInputHandler");
 227                 }else{
 228                 if (fDebug > 1) AliInfo(" ==== Tracks from AliAODInputHandler / Jets on-the-fly");
 229                 }
 230         } else {  //output AOD
 231         handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
 232         if( handler && handler->InheritsFrom("AliAODHandler") ) {
 233                 fAOD = ((AliAODHandler*)handler)->GetAOD();
 234                 if (!fJetsOnFly){
 235                         if (fDebug > 1) AliInfo(" ==== Tracks and Jets from AliAODHandler");
 236                         } else {
 237                         if (fDebug > 1) AliInfo(" ==== Tracks from AliAODHandler / Jets on-the-fly");
 238                         }
 239                 }else {
 240                 AliFatal("I can't get any AOD Event Handler");
 241                 return;
 242                 }
 243         }
 244
 245    fAnaUE->Initialize( *this );
 246
 247 }
 248
 249 //____________________________________________________________________
 250 void  AliAnalysisTaskUE::CreateOutputObjects()
 251 {
 252   // Create the output container
 253
 254   if (fDebug > 1) AliInfo("CreateOutPutData()");
 255
 256   //Initialize AliAnalysisUE, a class implementing the main analysis algorithms
 257   fAnaUE = new AliAnalyseUE();
 258   fHistosUE = new AliHistogramsUE();
 259
 260   if (fListOfHistos != NULL){
 261         delete fListOfHistos;
 262         fListOfHistos = NULL;
 263         }
 264   if (!fListOfHistos){
 265         fListOfHistos = new TList();
 266         fListOfHistos->SetOwner(kTRUE);
 267         }
 268
 269   //Initialize output histograms
 270   fHistosUE->CreateHistograms(fListOfHistos,fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist, fTrackEtaCut);
 271   AddSettingsTree();
 272
 273   PostData(0,fListOfHistos);
 274 }
 275
 276 //____________________________________________________________________
 277 void  AliAnalysisTaskUE::Exec(Option_t */*option*/)
 278 {
 279   // Trigger selection ************************************************
 280   AliInputEventHandler* inputHandler = (AliInputEventHandler*)
 281          ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
 282   if (!fAnaUE->TriggerSelection(inputHandler)) return;
 283
 284   // Vertex selection *************************************************
 285   if(!fAnaUE->VertexSelection(fAOD,fnTracksVertex,fZVertex)) return;
 286   //if(!fAnaUE->VertexSelectionOld(fAOD)) return; // temporary to compare with old task and to have same cuts for MC !!!
 287
 288   // Execute analysis for current event ******************************
 289
 290   if ( fDebug > 3 ) AliInfo( " Processing event..." );
 291
 292   // fetch the pythia header info and get the trials
 293   AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 294   Float_t nTrials = 1;
 295   if (mcHandler) {
 296         AliMCEvent* mcEvent = mcHandler->MCEvent();
 297         if (mcEvent) {
 298                 AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
 299                 if(pythiaGenHeader){
 300                         nTrials = pythiaGenHeader->Trials();
 301                         }
 302                 }
 303         }
 304   fHistosUE->GetTrials()->Fill("#sum{ntrials}",fAvgTrials);
 305
 306   //analyse the event
 307   AnalyseUE();
 308
 309  PostData(0,fListOfHistos);
 310 }
 311
 312 //____________________________________________________________________
 313 void  AliAnalysisTaskUE::AddSettingsTree()
 314 {
 315   //Write settings to output list
 316   TTree *settingsTree   = new TTree("UEAnalysisSettings","Analysis Settings in UE estimation");
 317   settingsTree->Branch("fFilterBit", &fFilterBit,"FilterBit/I");
 318   settingsTree->Branch("fConeRadius", &fConeRadius,"Rad/D");
 319   settingsTree->Branch("fJet1EtaCut", &fJet1EtaCut, "LeadJetEtaCut/D");
 320   settingsTree->Branch("fJet2DeltaPhiCut", &fJet2DeltaPhiCut, "DeltaPhi/D");
 321   settingsTree->Branch("fJet2RatioPtCut", &fJet2RatioPtCut, "Jet2Ratio/D");
 322   settingsTree->Branch("fJet3PtCut", &fJet3PtCut, "Jet3PtCut/D");
 323   settingsTree->Branch("fTrackPtCut", &fTrackPtCut, "TrackPtCut/D");
 324   settingsTree->Branch("fTrackEtaCut", &fTrackEtaCut, "TrackEtaCut/D");
 325   settingsTree->Branch("fAnaType", &fAnaType, "Ana/I");
 326   settingsTree->Branch("fRegionType", &fRegionType,"Reg/I");
 327   settingsTree->Branch("fOrdering", &fOrdering,"OrderMeth/I");
 328   settingsTree->Branch("fUseChPartJet", &fUseChPartJet,"UseChPart/O");
 329   settingsTree->Branch("fUseChargeHadrons", &fUseChargeHadrons,"UseChHadrons/O");
 330   settingsTree->Branch("fUseSingleCharge", &fUseSingleCharge,"UseSingleCh/O");
 331   settingsTree->Branch("fUsePositiveCharge", &fUsePositiveCharge,"UsePositiveCh/O");
 332   settingsTree->Fill();
 333   fListOfHistos->Add(settingsTree);
 334 }
 335
 336 //____________________________________________________________________
 337 void  AliAnalysisTaskUE::AnalyseUE()
 338 {
 339
 340
 341    // Get jets and order by pT
 342    TVector3 jetVect[3];
 343    *jetVect = fAnaUE->GetOrderedClusters(fAODBranch, fUseChPartJet, fChJetPtMin );
 344
 345    if( jetVect[0].Pt() < 0. ) {
 346         if( fDebug > 1 ) AliInfo("\n   Skipping Event, not jet found...");
 347         return;
 348         } else {
 349                 if (fDebug >1 ) AliInfo(Form("\n   Pt Leading Jet = %6.1f eta=%5.3f ",  jetVect[0].Pt(), jetVect[0].Eta() ));
 350                 }
 351
 352    // Select events according to analysis type ***********************************
 353    if ( ! (fAnaUE->AnaTypeSelection(jetVect))) return;
 354
 355   // Find max and min regions with real tracks
 356   if (!fUseMCParticleBranch){
 357         // Primary vertex distribution
 358         AliAODVertex* vertex = (AliAODVertex*)fAOD->GetPrimaryVertex();
 359         fHistosUE->FillHistogram("hVertexMult",vertex->GetNContributors());
 360
 361         // Fill leading "jet" histogram
 362         fHistosUE->FillHistogram("hEleadingPt",jetVect[0].Pt());
 363
 364         fAnaUE->FindMaxMinRegions( jetVect,  fConePosition, 0, 0 );
 365
 366   }else {
 367
 368     // this is the part we only use when we have MC information
 369     // More than a test for values of it also resumes the reconstruction efficiency of jets
 370     // As commented bellow if available for the data, we try to pair reconstructed jets with simulated ones
 371     // afterwards we kept angular variables of MC jet to perform UE analysis over MC particles
 372     // TODO: Handle Multiple jet environment. 06/2009 just suited for inclusive jet condition ( fAnaType = 1 )
 373
 374       AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 375       if (!mcHandler) {
 376         Printf("ERROR: Could not retrieve MC event handler");
 377         return;
 378       }
 379
 380       AliMCEvent* mcEvent = mcHandler->MCEvent();
 381       if (!mcEvent) {
 382         Printf("ERROR: Could not retrieve MC event");
 383         return;
 384       }
 385     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
 386     if(!pythiaGenHeader){
 387       return;
 388     }
 389
 390     fAnaUE->AnalyseMC(jetVect,mcEvent,pythiaGenHeader, fConePosition, fUseAliStack, fConstrainDistance, fMinDistance);
 391
 392   }
 393
 394   fAnaUE->FillRegions(fIsNorm2Area, jetVect);
 395
 396 }
 397
 398 //____________________________________________________________________
 399 AliAnalysisTaskUE* AliAnalysisTaskUE::Instance()
 400 {
 401
 402   //Create instance
 403   if (fgTaskUE) {
 404         return fgTaskUE;
 405   } else {
 406         fgTaskUE = new AliAnalysisTaskUE();
 407         return fgTaskUE;
 408         }
 409 }
 410
 411
 412 //____________________________________________________________________
 413 void  AliAnalysisTaskUE::Terminate(Option_t */*option*/)
 414 {
 415
 416   // Terminate analysis
 417
 418   if (!gROOT->IsBatch()){
 419   fListOfHistos = dynamic_cast<TList*> (GetOutputData(0));
 420   if (!fListOfHistos){
 421         AliError("Histogram List is not available");
 422         return;
 423         }
 424   //call class AliHistogramsUE
 425   AliHistogramsUE *histos=new AliHistogramsUE(fListOfHistos);
 426   histos->DrawUE(0);
 427   } else {
 428         AliInfo(" Batch mode, not histograms will be shown...");
 429   }
 430
 431   if( fDebug > 1 )
 432     AliInfo("End analysis");
 433
 434 }
 435
 436 void  AliAnalysisTaskUE::WriteSettings()
 437 {
 438
 439   //Print analysis settings on screen
 440   if (fDebug > 5){
 441     AliInfo(" All Analysis Settings in Saved Tree");
 442     fListOfHistos = dynamic_cast<TList*> (GetOutputData(0));
 443     if (!fListOfHistos){
 444         AliError("Histogram List is not available");
 445         return;
 446         }
 447      TTree *tree = (TTree*)(fListOfHistos->FindObject("UEAnalysisSettings"));
 448      tree->Scan();
 449      }
 450 }