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