PWG4/AliAnaCaloTriggerMC.cxx

   1
   2 /**************************************************************************
   3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   4  *                                                                        *
   5  * Author: The ALICE Off-line Project.                                    *
   6  * Contributors are mentioned in the code where appropriate.              *
   7  *                                                                        *
   8  * Permission to use, copy, modify and distribute this software and its   *
   9  * documentation strictly for non-commercial purposes is hereby granted   *
  10  * without fee, provided that the above copyright notice appears in all   *
  11  * copies and that both the copyright notice and this permission notice   *
  12  * appear in the supporting documentation. The authors make no claims     *
  13  * about the suitability of this software for any purpose. It is          *
  14  * provided "as is" without express or implied warranty.                  *
  15  **************************************************************************/
  16
  17 //_________________________________________________________________________
  18 // An analysis task to check the trigger data in ESD with MC data
  19 // Creates an ntuple for 2x2 and NxN triggers
  20 // Each ntuple connects the maximum trigger amplitudes
  21 // and its positions with reconstructed clusters and MC
  22 //
  23 //*-- Yves Schutz (CERN) & Gustavo Conesa Balbastre (INFN-LNF)
  24 //////////////////////////////////////////////////////////////////////////////
  25
  26 #include <TChain.h>
  27 #include <TFile.h>
  28 #include <TNtuple.h>
  29 #include <TVector3.h>
  30
  31 #include "AliAnalysisManager.h"
  32 #include "AliMCEventHandler.h"
  33 #include "AliMCEvent.h"
  34 #include "AliAnaCaloTriggerMC.h"
  35 #include "AliESDEvent.h"
  36 #include "AliESDCaloCluster.h"
  37 #include "AliLog.h"
  38 #include "AliStack.h"
  39 #include "TParticle.h"
  40
  41 //______________________________________________________________________________
  42 AliAnaCaloTriggerMC::AliAnaCaloTriggerMC() :
  43   fChain(0),
  44   fESD(0),
  45   fOutputContainer(0),
  46   fCalorimeter("PHOS"),
  47   fNtTrigger22(0),
  48   fNtTriggerNN(0)
  49
  50 {
  51   // Default constructor.
  52
  53 }
  54 //______________________________________________________________________________
  55 AliAnaCaloTriggerMC::AliAnaCaloTriggerMC(const char *name) :
  56   AliAnalysisTask(name, "AnaCaloTriggerMC"),
  57   fChain(0),
  58   fESD(0),
  59   fOutputContainer(0),
  60   fCalorimeter("PHOS"),
  61   fNtTrigger22(0),
  62   fNtTriggerNN(0)
  63
  64 {
  65   // Constructor.
  66   // Input slot #0 works with an Ntuple
  67   DefineInput(0, TChain::Class());
  68   // Output slot #0 writes into a TH1 container
  69   DefineOutput(0,  TObjArray::Class()) ;
  70 }
  71 //____________________________________________________________________________
  72 AliAnaCaloTriggerMC::AliAnaCaloTriggerMC(const AliAnaCaloTriggerMC & ct) :
  73   AliAnalysisTask(ct),fChain(ct.fChain), fESD(ct.fESD),
  74   fOutputContainer(ct.fOutputContainer), fCalorimeter(ct. fCalorimeter),
  75   fNtTrigger22(ct.fNtTrigger22), fNtTriggerNN(ct.fNtTriggerNN)
  76 {
  77
  78   // cpy ctor
  79   SetName (ct.GetName()) ;
  80   SetTitle(ct.GetTitle()) ;
  81
  82 }
  83
  84 //_________________________________________________________________________
  85 AliAnaCaloTriggerMC & AliAnaCaloTriggerMC::operator = (const AliAnaCaloTriggerMC & source)
  86 {
  87   // assignment operator
  88
  89   if(&source == this) return *this;
  90
  91   fChain = source.fChain ;
  92   fESD = source.fESD ;
  93   fOutputContainer = source.fOutputContainer ;
  94   fCalorimeter = source. fCalorimeter ;
  95   fNtTrigger22 = source.fNtTrigger22 ;
  96   fNtTriggerNN = source.fNtTriggerNN ;
  97
  98   return *this;
  99
 100 }
 101
 102 //______________________________________________________________________________
 103 AliAnaCaloTriggerMC::~AliAnaCaloTriggerMC()
 104 {
 105   // dtor
 106   //  fOutputContainer->Clear() ;
 107   //  delete fOutputContainer ;
 108
 109 }
 110
 111
 112 //______________________________________________________________________________
 113 void AliAnaCaloTriggerMC::ConnectInputData(const Option_t*)
 114 {
 115   // Initialisation of branch container and histograms
 116
 117   AliInfo(Form("*** Initialization of %s", GetName())) ;
 118
 119   // Get input data
 120   fChain = dynamic_cast<TChain *>(GetInputData(0)) ;
 121   if (!fChain) {
 122     AliError(Form("Input 0 for %s not found\n", GetName()));
 123     return ;
 124   }
 125
 126   fESD = new AliESDEvent();
 127   fESD->ReadFromTree(fChain);
 128
 129 }
 130
 131 //________________________________________________________________________
 132
 133 void AliAnaCaloTriggerMC::CreateOutputObjects()
 134 {
 135
 136   // Create the output container
 137   OpenFile(0);
 138
 139   // create histograms
 140   fNtTrigger22 = new TNtuple(fCalorimeter+"trigger22", "Trigger data 2x2 patch", "a22:a220:enMax:phEnMax:eta22:phi22:etaMax:phiMax:phEtaMax:phPhiMax");
 141   fNtTriggerNN = new TNtuple(fCalorimeter+"triggerNN", "Trigger data NxN patch", "aNN:aNN0:enMax:phEnMax:etaNN:phiNN:etaMax:phiMax:phEtaMax:phPhiMax");
 142
 143   // create output container
 144
 145   fOutputContainer = new TObjArray(2) ;
 146   fOutputContainer->SetName(GetName()) ;
 147
 148   fOutputContainer->AddAt(fNtTrigger22,             0) ;
 149   fOutputContainer->AddAt(fNtTriggerNN,             1) ;
 150
 151 }
 152
 153 //______________________________________________________________________________
 154 void AliAnaCaloTriggerMC::Exec(Option_t *)
 155 {
 156
 157   // Processing of one event
 158
 159   Long64_t entry = fChain->GetReadEntry() ;
 160
 161   if (!fESD) {
 162     AliError("fESD is not connected to the input!") ;
 163     return ;
 164   }
 165
 166   if ( !((entry-1)%100) )
 167     AliInfo(Form("%s ----> Processing event # %lld",  (dynamic_cast<TChain *>(fChain))->GetFile()->GetName(), entry)) ;
 168
 169   //Get MC data
 170   AliStack* stack = 0x0;
 171   AliMCEventHandler*    mctruth = (AliMCEventHandler*)
 172     ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
 173
 174   if(mctruth)
 175     stack = mctruth->MCEvent()->Stack();
 176
 177   if (!stack) {
 178     AliError("Stack not found") ;
 179     return ;
 180   }
 181
 182   // Get trigger information of fCalorimeter
 183   TArrayF * triggerAmplitudes = 0x0 ;
 184   TArrayF * triggerPosition   = 0x0 ;
 185   Int_t numberOfCaloClusters  = 0 ;
 186
 187   if(fCalorimeter == "PHOS"){
 188     triggerAmplitudes      = fESD->GetPHOSTriggerAmplitudes();
 189     triggerPosition        = fESD->GetPHOSTriggerPosition();
 190   }
 191   else if(fCalorimeter == "EMCAL"){
 192     triggerAmplitudes    = fESD->GetEMCALTriggerAmplitudes();
 193     triggerPosition      = fESD->GetEMCALTriggerPosition();
 194   }
 195
 196   if( triggerAmplitudes && triggerPosition ){
 197     // trigger amplitudes
 198     const Float_t ka22    = static_cast<Float_t>(triggerAmplitudes->At(0)) ;
 199     const Float_t ka22O   = static_cast<Float_t>(triggerAmplitudes->At(1)) ;
 200     const Float_t kaNN    = static_cast<Float_t>(triggerAmplitudes->At(2)) ;
 201     const Float_t kaNNO   = static_cast<Float_t>(triggerAmplitudes->At(3)) ;
 202
 203     // trigger position
 204     const Float_t kx22  =  static_cast<Float_t>(triggerPosition->At(0)) ;
 205     const Float_t ky22  =  static_cast<Float_t>(triggerPosition->At(1)) ;
 206     const Float_t kz22  =  static_cast<Float_t>(triggerPosition->At(2)) ;
 207     const Float_t kxNN  =  static_cast<Float_t>(triggerPosition->At(3)) ;
 208     const Float_t kyNN  =  static_cast<Float_t>(triggerPosition->At(4)) ;
 209     const Float_t kzNN  =  static_cast<Float_t>(triggerPosition->At(5)) ;
 210
 211     Float_t enMax       = 0. ;
 212     Float_t phEnMax     = 0. ;
 213     Float_t etaMax      = 0.5 ;
 214     Float_t phiMax      = 0. ;
 215     Float_t phEtaMax    = 0.5 ;
 216     Float_t phPhiMax    = 0. ;
 217
 218     TVector3 vpos22(kx22, ky22, kz22) ;
 219     TVector3 vposNN(kxNN, kyNN, kzNN) ;
 220     Float_t eta22 = vpos22.Eta() ;
 221     Float_t phi22 = vpos22.Phi() * TMath::RadToDeg() + 360. ;
 222     Float_t etaNN = vposNN.Eta() ;
 223     Float_t phiNN = vposNN.Phi() * TMath::RadToDeg() + 360. ;
 224
 225     Int_t      icaloCluster ;
 226
 227     // loop over the Calorimeters Clusters
 228     Float_t cluEnergy = 0;
 229     Int_t labelmax = -5;
 230     for(icaloCluster = 0 ; icaloCluster < numberOfCaloClusters ; icaloCluster++) {
 231
 232       AliESDCaloCluster * cluster = fESD->GetCaloCluster(icaloCluster) ;
 233
 234       if (cluster && ( (fCalorimeter == "PHOS" && cluster->IsPHOS())  ||
 235                        (fCalorimeter == "EMCAL" && cluster->IsEMCAL()))) {
 236
 237         cluEnergy = cluster->E() ;
 238         Float_t pos[3] ;
 239         TVector3 vpos ;
 240
 241         cluster->GetPosition( pos ) ;
 242
 243         if ( cluEnergy > enMax) {
 244           enMax = cluEnergy ;
 245           vpos.SetXYZ(pos[0], pos[1], pos[2]) ;
 246           etaMax = vpos.Eta() ;
 247           phiMax = vpos.Phi() ;
 248           labelmax = cluster->GetLabel();
 249         }
 250
 251         Double_t * pid = cluster->GetPid() ;
 252
 253         if(pid[AliPID::kPhoton] > 0.9) {
 254           if ( cluEnergy > phEnMax) {
 255             phEnMax = cluEnergy ;
 256             vpos.SetXYZ(pos[0], pos[1], pos[2]) ;
 257             phEtaMax = vpos.Eta() ;
 258             phPhiMax = vpos.Phi() ;
 259           }
 260         }
 261       }//if cluster
 262     }//CaloCluster loop
 263
 264     if(labelmax < stack->GetNtrack() && labelmax >= 0 ){
 265       TParticle * particle = stack->Particle(labelmax);
 266       Float_t ptgen = particle->Energy();
 267       fNtTrigger22->Fill(ka22, ka22O, ptgen, enMax, phEnMax, eta22, phi22, etaMax, phiMax * TMath::RadToDeg() + 360., phEtaMax, phPhiMax * TMath::RadToDeg() + 360.);
 268       fNtTriggerNN->Fill(kaNN, kaNNO, ptgen, enMax, phEnMax, etaNN, phiNN, etaMax, phiMax * TMath::RadToDeg() + 360., phEtaMax, phPhiMax * TMath::RadToDeg() + 360.);
 269     }
 270     else AliDebug(1, Form("Wrong label %d, ntrack %d, Emax %f ",labelmax, stack->GetNtrack(), phEnMax));
 271   }//If trigger arrays filled
 272
 273   PostData(0, fOutputContainer);
 274
 275 }
 276
 277
 278 //______________________________________________________________________________
 279 void AliAnaCaloTriggerMC::Terminate(Option_t *)
 280 {
 281   // Processing when the event loop is ended
 282
 283 }