1)Terminate() method implemented in the frame. Simple examples on what to do with...
[u/mrichter/AliRoot.git] / PWG4 / PartCorrDep / AliAnaCaloTrigger.cxx
index 140ce55cd516de1bfa78d86e917c862a1b2cfa66..a941a8d9ed25e1c48ade9cd2d1ff28816e93b0fd 100755 (executable)
 // Creates an ntuple for 2x2 and NxN triggers
 // Each ntuple connects the maximum trigger amplitudes 
 // and its positions with reconstructed clusters
+// and if MC stack available, with pt of parent.
 //
 //*-- Yves Schutz (CERN) & Gustavo Conesa Balbastre (INFN-LNF)
 //////////////////////////////////////////////////////////////////////////////
-
-#include <TChain.h>
-#include <TFile.h> 
+//Root
 #include <TNtuple.h>
 #include <TVector3.h> 
 
+//Aliroot
 #include "AliAnaCaloTrigger.h" 
-#include "AliAnalysisManager.h"
-#include "AliESDEvent.h" 
+#include "AliStack.h"
 #include "AliLog.h"
 #include "AliESDCaloCluster.h"
+#include "AliMCEvent.h"
+#include "AliESDEvent.h"
 
 //______________________________________________________________________________
 AliAnaCaloTrigger::AliAnaCaloTrigger() :  
-  fChain(0),
-  fESD(0), 
   fOutputContainer(0),
   fCalorimeter("PHOS"),
   fNtTrigger22(0), 
@@ -49,9 +48,7 @@ AliAnaCaloTrigger::AliAnaCaloTrigger() :
 
 //______________________________________________________________________________
 AliAnaCaloTrigger::AliAnaCaloTrigger(const char *name) : 
-  AliAnalysisTask(name,"AnaCaloTrigger"),  
-  fChain(0),
-  fESD(0), 
+  AliAnalysisTaskSE(name),  
   fOutputContainer(0),
   fCalorimeter("PHOS"),
   fNtTrigger22(0), 
@@ -59,22 +56,19 @@ AliAnaCaloTrigger::AliAnaCaloTrigger(const char *name) :
 
 {
   // Constructor.
-  // Input slot #0 works with an Ntuple
-  DefineInput(0, TChain::Class());
-  // Output slot #0 writes into a TH1 container
-  DefineOutput(0,  TObjArray::Class()) ; 
+  // Output slot 
+  DefineOutput(1,  TList::Class()) ; 
 }
+
 //____________________________________________________________________________
 AliAnaCaloTrigger::AliAnaCaloTrigger(const AliAnaCaloTrigger & ct) : 
-  AliAnalysisTask(ct), fChain(ct.fChain), fESD(ct.fESD),
+  AliAnalysisTaskSE(ct.GetName()),
   fOutputContainer(ct.fOutputContainer), fCalorimeter(ct. fCalorimeter),
   fNtTrigger22(ct.fNtTrigger22), fNtTriggerNN(ct.fNtTriggerNN)
 {
 
   // cpy ctor
-  SetName (ct.GetName()) ; 
-  SetTitle(ct.GetTitle()) ;
+   
 }
 
 //_________________________________________________________________________
@@ -82,10 +76,10 @@ AliAnaCaloTrigger & AliAnaCaloTrigger::operator = (const AliAnaCaloTrigger & sou
 {
   // assignment operator
   
-  if(&source == this) return *this;
+  //if(&source == this) return *this;
+  this->~AliAnaCaloTrigger();
+  new(this) AliAnaCaloTrigger(source);
 
-  fChain = source.fChain ; 
-  fESD = source.fESD ;
   fOutputContainer = source.fOutputContainer ;
   fCalorimeter = source. fCalorimeter ;
   fNtTrigger22 = source.fNtTrigger22 ;
@@ -99,46 +93,28 @@ AliAnaCaloTrigger & AliAnaCaloTrigger::operator = (const AliAnaCaloTrigger & sou
 AliAnaCaloTrigger::~AliAnaCaloTrigger()
 {
   // dtor
-  //fOutputContainer->Clear() ; 
-  //delete fOutputContainer ;
-
+       if(fOutputContainer){
+               fOutputContainer->Clear() ; 
+               delete fOutputContainer ;
+       }
 }
 
 
-//______________________________________________________________________________
-void AliAnaCaloTrigger::ConnectInputData(const Option_t*)
-{
-  // Initialisation of branch container and histograms 
-    
-  AliInfo(Form("*** Initialization of %s", GetName())) ; 
-  
-  // Get input data
-  fChain = dynamic_cast<TChain *>(GetInputData(0)) ;
-  if (!fChain) {
-    AliError(Form("Input 0 for %s not found\n", GetName()));
-    return ;
-  }
-  
-  fESD = new AliESDEvent();
-  fESD->ReadFromTree(fChain);
-
-}
-
 //________________________________________________________________________
 
-void AliAnaCaloTrigger::CreateOutputObjects()
+void AliAnaCaloTrigger::UserCreateOutputObjects()
 {  
 
   // Create the outputs containers
-  OpenFile(0) ;
+  OpenFile(1) ;
 
   // create histograms 
-  fNtTrigger22 = new TNtuple(fCalorimeter+"trigger22", "Trigger data 2x2 patch", "a22:a220:enMax:phEnMax:eta22:phi22:etaMax:phiMax:phEtaMax:phPhiMax");
-  fNtTriggerNN = new TNtuple(fCalorimeter+"triggerNN", "Trigger data NxN patch", "aNN:aNN0:enMax:phEnMax:etaNN:phiNN:etaMax:phiMax:phEtaMax:phPhiMax");
+  fNtTrigger22 = new TNtuple(fCalorimeter+"trigger22", "Trigger data 2x2 patch", "a22:a220:ptGen:enMax:phEnMax:eta22:phi22:etaMax:phiMax:phEtaMax:phPhiMax");
+  fNtTriggerNN = new TNtuple(fCalorimeter+"triggerNN", "Trigger data NxN patch", "aNN:aNN0:ptGen:enMax:phEnMax:etaNN:phiNN:etaMax:phiMax:phEtaMax:phPhiMax");
   
   // create output container
   
-  fOutputContainer = new TObjArray(2) ; 
+  fOutputContainer = new TList() ; 
   fOutputContainer->SetName(GetName()) ; 
   
   fOutputContainer->AddAt(fNtTrigger22,             0) ; 
@@ -147,107 +123,118 @@ void AliAnaCaloTrigger::CreateOutputObjects()
 }
 
 //______________________________________________________________________________
-void AliAnaCaloTrigger::Exec(Option_t *) 
+void AliAnaCaloTrigger::UserExec(Option_t *) 
 {
-  // Processing of one event
-  
-  Long64_t entry = fChain->GetReadEntry() ;
-  
-  if (!fESD) {
-    AliError("fESD is not connected to the input!") ; 
-    return ; 
-  }
-  
-  if ( !((entry-1)%100) ) 
-    AliInfo(Form("%s ----> Processing event # %lld",  (dynamic_cast<TChain *>(fChain))->GetFile()->GetName(), entry)) ; 
-  
-  // Get trigger information of fCalorimeter 
-  TArrayF * triggerAmplitudes = 0x0 ;
-  TArrayF * triggerPosition   = 0x0 ;
-  Int_t numberOfCaloClusters  =  fESD->GetNumberOfCaloClusters() ;
-
-  if(fCalorimeter == "PHOS"){
-    triggerAmplitudes      = fESD->GetPHOSTriggerAmplitudes();
-    triggerPosition        = fESD->GetPHOSTriggerPosition();
-  }
-  else if(fCalorimeter == "EMCAL"){
-    triggerAmplitudes    = fESD->GetEMCALTriggerAmplitudes();
-    triggerPosition      = fESD->GetEMCALTriggerPosition();
-  }
-
-  if( triggerAmplitudes && triggerPosition ){
-  // trigger amplitudes
-  const Float_t ka22    = static_cast<Float_t>(triggerAmplitudes->At(0)) ; 
-  const Float_t ka22O   = static_cast<Float_t>(triggerAmplitudes->At(1)) ; 
-  const Float_t kaNN    = static_cast<Float_t>(triggerAmplitudes->At(2)) ; 
-  const Float_t kaNNO   = static_cast<Float_t>(triggerAmplitudes->At(3)) ; 
-
-  // trigger position
-  const Float_t kx22  =  static_cast<Float_t>(triggerPosition->At(0)) ; 
-  const Float_t ky22  =  static_cast<Float_t>(triggerPosition->At(1)) ;
-  const Float_t kz22  =  static_cast<Float_t>(triggerPosition->At(2)) ;
-  const Float_t kxNN  =  static_cast<Float_t>(triggerPosition->At(3)) ; 
-  const Float_t kyNN  =  static_cast<Float_t>(triggerPosition->At(4)) ;
-  const Float_t kzNN  =  static_cast<Float_t>(triggerPosition->At(5)) ; 
-  
-  Float_t enMax       = 0. ;
-  Float_t phEnMax     = 0. ;
-  Float_t etaMax      = 0.5 ;
-  Float_t phiMax      = 0. ; 
-  Float_t phEtaMax    = 0.5 ;
-  Float_t phPhiMax    = 0. ; 
-  
-  TVector3 vpos22(kx22, ky22, kz22) ;
-  TVector3 vposNN(kxNN, kyNN, kzNN) ;
-  Float_t eta22 = vpos22.Eta() ; 
-  Float_t phi22 = vpos22.Phi() * TMath::RadToDeg() + 360. ; 
-  Float_t etaNN = vposNN.Eta() ; 
-  Float_t phiNN = vposNN.Phi() * TMath::RadToDeg() + 360. ; 
-
-  Int_t      icaloCluster ; 
-  
-  // loop over the Calorimeters Clusters
-  
-  for(icaloCluster = 0 ; icaloCluster < numberOfCaloClusters ; icaloCluster++) {
-    
-    AliESDCaloCluster * cluster = fESD->GetCaloCluster(icaloCluster) ;
-    
-    if (cluster && ( (fCalorimeter == "PHOS" && cluster->IsPHOS())  ||  
-                    (fCalorimeter == "EMCAL" && cluster->IsEMCAL()))) {
-         
-      Float_t cluEnergy = cluster->E() ; 
-      Float_t pos[3] ;
-      TVector3 vpos ;
-      
-      cluster->GetPosition( pos ) ;
-      
-      if ( cluEnergy > enMax) { 
-       enMax = cluEnergy ; 
-       vpos.SetXYZ(pos[0], pos[1], pos[2]) ; 
-       etaMax = vpos.Eta() ; 
-       phiMax = vpos.Phi() ; 
-      }
-
-      Double_t * pid = cluster->GetPid() ;
-      
-      if(pid[AliPID::kPhoton] > 0.9) {
-       if ( cluEnergy > phEnMax) { 
-         phEnMax = cluEnergy ; 
-         vpos.SetXYZ(pos[0], pos[1], pos[2]) ; 
-         phEtaMax = vpos.Eta() ; 
-         phPhiMax = vpos.Phi() ; 
+       // Processing of one event
+       
+       if ( !((Entry()-1)%100) ) 
+               printf(" Processing event # %lld\n",  Entry()) ; 
+       AliESDEvent* esd = (AliESDEvent*)InputEvent();
+       
+       //Get MC data, if available
+       AliStack* stack = 0x0; 
+       if(MCEvent())
+               stack = MCEvent()->Stack();
+       
+       // Get trigger information of fCalorimeter 
+       TArrayF * triggerAmplitudes = 0x0 ;
+       TArrayF * triggerPosition   = 0x0 ;
+       Int_t numberOfCaloClusters  =  esd->GetNumberOfCaloClusters() ;
+       
+       if(fCalorimeter == "PHOS"){
+               triggerAmplitudes      = esd->GetPHOSTriggerAmplitudes();
+               triggerPosition        = esd->GetPHOSTriggerPosition();
        }
-      }
-    }//if cluster
-    
-    fNtTrigger22->Fill(ka22, ka22O, enMax, phEnMax, eta22, phi22, etaMax, phiMax * TMath::RadToDeg() + 360., phEtaMax, phPhiMax * TMath::RadToDeg() + 360.);
-    fNtTriggerNN->Fill(kaNN, kaNNO, enMax, phEnMax, etaNN, phiNN, etaMax, phiMax * TMath::RadToDeg() + 360., phEtaMax, phPhiMax * TMath::RadToDeg() + 360.);
-  }//CaloCluster loop
-  
-  }//If trigger arrays filled
-  
-  PostData(0, fOutputContainer);
-  
+       else if(fCalorimeter == "EMCAL"){
+               triggerAmplitudes    = esd->GetEMCALTriggerAmplitudes();
+               triggerPosition      = esd->GetEMCALTriggerPosition();
+       }
+       
+       if( triggerAmplitudes && triggerPosition ){
+               // trigger amplitudes
+               const Float_t ka22    = static_cast<Float_t>(triggerAmplitudes->At(0)) ; 
+               const Float_t ka22O   = static_cast<Float_t>(triggerAmplitudes->At(1)) ; 
+               const Float_t kaNN    = static_cast<Float_t>(triggerAmplitudes->At(2)) ; 
+               const Float_t kaNNO   = static_cast<Float_t>(triggerAmplitudes->At(3)) ; 
+               
+               // trigger position
+               const Float_t kx22  =  static_cast<Float_t>(triggerPosition->At(0)) ; 
+               const Float_t ky22  =  static_cast<Float_t>(triggerPosition->At(1)) ;
+               const Float_t kz22  =  static_cast<Float_t>(triggerPosition->At(2)) ;
+               const Float_t kxNN  =  static_cast<Float_t>(triggerPosition->At(3)) ; 
+               const Float_t kyNN  =  static_cast<Float_t>(triggerPosition->At(4)) ;
+               const Float_t kzNN  =  static_cast<Float_t>(triggerPosition->At(5)) ; 
+               
+               //printf("ka22 %f, ka220 %f, kaNN %f, kaNN0 %f\n",ka22,ka22O,kaNN,kaNNO);
+               //printf("kx22 %f, ky22 %f, kz22 %f, kxNN %f, kyNN %f, kzNN %f \n",kx22,ky22,kz22,kxNN,kyNN,kzNN);
+               
+               Float_t enMax       = 0. ;
+               Float_t phEnMax     = 0. ;
+               Float_t etaMax      = 0.5 ;
+               Float_t phiMax      = 0. ; 
+               Float_t phEtaMax    = 0.5 ;
+               Float_t phPhiMax    = 0. ; 
+               
+               TVector3 vpos22(kx22, ky22, kz22) ;
+               TVector3 vposNN(kxNN, kyNN, kzNN) ;
+               Float_t eta22 = vpos22.Eta() ; 
+               Float_t phi22 = vpos22.Phi() * TMath::RadToDeg() + 360. ; 
+               Float_t etaNN = vposNN.Eta() ; 
+               Float_t phiNN = vposNN.Phi() * TMath::RadToDeg() + 360. ; 
+               
+               
+               Int_t      icaloCluster = 0 ; 
+               Int_t      labelmax     = -1 ;
+               // loop over the Calorimeters Clusters
+               
+               for(icaloCluster = 0 ; icaloCluster < numberOfCaloClusters ; icaloCluster++) {
+                       
+                       AliESDCaloCluster * cluster = esd->GetCaloCluster(icaloCluster) ;
+                       
+                       if (cluster && ( (fCalorimeter == "PHOS" && cluster->IsPHOS())  ||  
+                                                       (fCalorimeter == "EMCAL" && cluster->IsEMCAL()))) {
+                               
+                               Float_t cluEnergy = cluster->E() ; 
+                               Float_t pos[3] ;
+                               TVector3 vpos ;
+                               
+                               cluster->GetPosition( pos ) ;
+                               
+                               if ( cluEnergy > enMax) { 
+                                       enMax    = cluEnergy ; 
+                                       vpos.SetXYZ(pos[0], pos[1], pos[2]) ; 
+                                       etaMax   = vpos.Eta() ; 
+                                       phiMax   = vpos.Phi() ; 
+                                       labelmax = cluster->GetLabel();
+                               }
+                               
+                               Double_t * pid = cluster->GetPid() ;
+                               
+                               if(pid[AliPID::kPhoton] > 0.9) {
+                                       if ( cluEnergy > phEnMax) { 
+                                               phEnMax   = cluEnergy ; 
+                                               vpos.SetXYZ(pos[0], pos[1], pos[2]) ; 
+                                               phEtaMax = vpos.Eta() ; 
+                                               phPhiMax = vpos.Phi() ; 
+                                       }
+                               }
+                       }//if cluster
+                       
+                   Float_t ptGen = -1;
+                       if(stack && labelmax < stack->GetNtrack() && labelmax >= 0 ){
+                               TParticle * particle = stack->Particle(labelmax); 
+                               ptGen = particle->Energy();
+                       }
+                       
+                       fNtTrigger22->Fill(ka22, ka22O, ptGen, enMax, phEnMax, eta22, phi22, etaMax, phiMax * TMath::RadToDeg() + 360., phEtaMax, phPhiMax * TMath::RadToDeg() + 360.);
+                       fNtTriggerNN->Fill(kaNN, kaNNO, ptGen, enMax, phEnMax, etaNN, phiNN, etaMax, phiMax * TMath::RadToDeg() + 360., phEtaMax, phPhiMax * TMath::RadToDeg() + 360.);
+               
+               }//CaloCluster loop
+               
+       }//If trigger arrays filled
+       
+       PostData(1, fOutputContainer);
+       
 }
 
 //______________________________________________________________________________