PWG4/AliAnalysisTaskPHOSExample.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 //_________________________________________________________________________
  19 // A basic analysis task to analyse photon detected by PHOS
  20 // A basic analysis task to analyse photon detected by PHOS
  21 // A basic analysis task to analyse photon detected by PHOS
  22 // Adapted for AliAnalysisTaskSE and AOD production
  23 // by Gustavo Conesa
  24 //
  25 //*-- Yves Schutz
  26 //////////////////////////////////////////////////////////////////////////////
  27
  28 #include <TCanvas.h>
  29 #include <TH1.h>
  30 #include <TROOT.h>
  31 #include <TLegend.h>
  32 #include <TNtuple.h>
  33 #include <TVector3.h>
  34
  35 #include "AliAnalysisTaskPHOSExample.h"
  36 #include "AliAnalysisManager.h"
  37 #include "AliESDEvent.h"
  38 #include "AliESDCaloCluster.h"
  39 #include "AliAODEvent.h"
  40 #include "AliAODPhoton.h"
  41 #include "AliLog.h"
  42 #include "AliESDVertex.h"
  43
  44 //______________________________________________________________________________
  45 AliAnalysisTaskPHOSExample::AliAnalysisTaskPHOSExample() :
  46   fDebug(0),
  47   fAODPhotons(0x0),
  48   fPhotonsInPhos(0),
  49   fPhotonId(1.0),
  50   fOutputList(0x0),
  51   fhPHOSPos(0),
  52   fhPHOS(0),
  53   fhPHOSEnergy(0),
  54   fhPHOSDigits(0),
  55   fhPHOSRecParticles(0),
  56   fhPHOSPhotons(0),
  57   fhPHOSInvariantMass(0),
  58   fhPHOSDigitsEvent(0)
  59 {
  60   //Default constructor
  61 }
  62 //______________________________________________________________________________
  63 AliAnalysisTaskPHOSExample::AliAnalysisTaskPHOSExample(const char *name) :
  64   AliAnalysisTaskSE(name),
  65   fDebug(0),
  66   fAODPhotons(0x0),
  67   fPhotonsInPhos(0),
  68   fPhotonId(1.0),
  69   fOutputList(0x0),
  70   fhPHOSPos(0),
  71   fhPHOS(0),
  72   fhPHOSEnergy(0),
  73   fhPHOSDigits(0),
  74   fhPHOSRecParticles(0),
  75   fhPHOSPhotons(0),
  76   fhPHOSInvariantMass(0),
  77   fhPHOSDigitsEvent(0)
  78 {
  79   // Constructor.
  80
  81   // Output slots
  82   DefineOutput(1,  TList::Class()) ;
  83 }
  84
  85 //____________________________________________________________________________
  86 AliAnalysisTaskPHOSExample::AliAnalysisTaskPHOSExample(const AliAnalysisTaskPHOSExample& ap) :
  87   AliAnalysisTaskSE(ap.GetName()),
  88   fDebug(ap.fDebug),
  89   fAODPhotons(ap.fAODPhotons),
  90   fPhotonsInPhos(ap.fPhotonsInPhos),
  91   fPhotonId(ap.fPhotonId),
  92   fOutputList(ap.fOutputList),
  93   fhPHOSPos(ap.fhPHOSPos),
  94   fhPHOS(ap.fhPHOS),
  95   fhPHOSEnergy(ap.fhPHOSEnergy),
  96   fhPHOSDigits(ap.fhPHOSDigits),
  97   fhPHOSRecParticles(ap.fhPHOSRecParticles),
  98   fhPHOSPhotons(ap.fhPHOSPhotons),
  99   fhPHOSInvariantMass(ap.fhPHOSInvariantMass),
 100   fhPHOSDigitsEvent(ap.fhPHOSDigitsEvent)
 101 {
 102   // cpy ctor
 103 }
 104
 105 //_____________________________________________________________________________
 106 AliAnalysisTaskPHOSExample& AliAnalysisTaskPHOSExample::operator = (const AliAnalysisTaskPHOSExample& ap)
 107 {
 108 // assignment operator
 109
 110   this->~AliAnalysisTaskPHOSExample();
 111   new(this) AliAnalysisTaskPHOSExample(ap);
 112   return *this;
 113 }
 114
 115 //______________________________________________________________________________
 116 AliAnalysisTaskPHOSExample::~AliAnalysisTaskPHOSExample()
 117 {
 118   // dtor
 119   if(fOutputList) {
 120     fOutputList->Clear() ;
 121     delete fOutputList ;
 122   }
 123 }
 124
 125
 126 //________________________________________________________________________
 127 void AliAnalysisTaskPHOSExample::UserCreateOutputObjects()
 128 {
 129   // Create the outputs containers
 130    //AODs
 131   fAODPhotons = new TClonesArray("AliAODPhoton", 0);
 132   fAODPhotons->SetName("Photons");
 133   AddAODBranch("TClonesArray", &fAODPhotons);
 134
 135   OpenFile(1) ;
 136
 137   fhPHOSPos            = new TNtuple("PHOSPos"         , "Position in PHOS"  , "x:y:z");
 138   fhPHOS               = new TNtuple("PHOS"            , "PHOS"  , "event:digits:clusters:photons");
 139   fhPHOSEnergy         = new TH1D("PHOSEnergy"         , "PHOSEnergy"        , 100, 0., 100. ) ;
 140   fhPHOSDigits         = new TH1I("PHOSDigitsCluster"  , "PHOSDigits"        , 20 , 0 , 20  ) ;
 141   fhPHOSRecParticles   = new TH1D("PHOSRecParticles"   , "PHOSRecParticles" , 20 , 0., 20. ) ;
 142   fhPHOSPhotons        = new TH1I("PHOSPhotons"        , "PHOSPhotons"       , 20 , 0 , 20  ) ;
 143   fhPHOSInvariantMass  = new TH1D("PHOSInvariantMass"  , "PHOSInvariantMass" , 400, 0., 400.) ;
 144   fhPHOSDigitsEvent    = new TH1I("PHOSDigitsEvent"    , "PHOSDigitsEvent"   , 30 , 0 , 30  ) ;
 145
 146   // create output container
 147
 148   fOutputList = new TList() ;
 149   fOutputList->SetName(GetName()) ;
 150
 151   fOutputList->AddAt(fhPHOSPos,             0) ;
 152   fOutputList->AddAt(fhPHOS,                1) ;
 153   fOutputList->AddAt(fhPHOSEnergy,          2) ;
 154   fOutputList->AddAt(fhPHOSDigits,          3) ;
 155   fOutputList->AddAt(fhPHOSRecParticles,    4) ;
 156   fOutputList->AddAt(fhPHOSPhotons,         5) ;
 157   fOutputList->AddAt(fhPHOSInvariantMass,   6) ;
 158   fOutputList->AddAt(fhPHOSDigitsEvent,     7) ;
 159
 160
 161 }
 162
 163 //______________________________________________________________________________
 164 void AliAnalysisTaskPHOSExample::UserExec(Option_t *)
 165 {
 166   // Processing of one event
 167
 168   // if ( !((Entry()-1)%100) )
 169   AliInfo(Form(" Processing event # %lld",  Entry())) ;
 170   AliESDEvent* esd = (AliESDEvent*)InputEvent();
 171
 172   //************************  PHOS *************************************
 173       TRefArray * caloClustersArr  = new TRefArray();
 174   esd->GetPHOSClusters(caloClustersArr);
 175
 176   const Int_t kNumberOfPhosClusters   = caloClustersArr->GetEntries() ;
 177
 178   TVector3 ** phosVector       = new TVector3*[kNumberOfPhosClusters] ;
 179   Float_t  * phosPhotonsEnergy = new Float_t[kNumberOfPhosClusters] ;
 180   Int_t      phosCluster ;
 181   Int_t      numberOfDigitsInPhos   = 0 ;
 182   Double_t v[3] ; //vertex ;
 183   esd->GetVertex()->GetXYZ(v) ;
 184
 185   fPhotonsInPhos  = 0 ;
 186   // loop over the PHOS Cluster
 187   for(phosCluster = 0 ; phosCluster < kNumberOfPhosClusters ; phosCluster++) {
 188     AliESDCaloCluster * caloCluster = (AliESDCaloCluster *) caloClustersArr->At(phosCluster) ;
 189
 190     //AliESDCaloCluster * caloCluster = fESD->GetCaloCluster(phosCluster) ;
 191     if (caloCluster) {
 192       Float_t pos[3] ;
 193       caloCluster->GetPosition( pos ) ;
 194       fhPHOSEnergy->Fill( caloCluster->E() ) ;
 195       fhPHOSPos->Fill( pos[0], pos[1], pos[2] ) ;
 196       fhPHOSDigits->Fill(Entry(), caloCluster->GetNCells() ) ;
 197       numberOfDigitsInPhos += caloCluster->GetNCells() ;
 198       Double_t * pid = caloCluster->GetPid() ;
 199       if(pid[AliPID::kPhoton] > GetPhotonId() ) {
 200         phosVector[fPhotonsInPhos] = new TVector3(pos[0],pos[1],pos[2]) ;
 201         phosPhotonsEnergy[fPhotonsInPhos]=caloCluster->E() ;
 202         TLorentzVector momentum ;
 203         caloCluster->GetMomentum(momentum, v);
 204         new ((*fAODPhotons)[fPhotonsInPhos++]) AliAODPhoton (momentum);
 205       }
 206     }
 207   } //PHOS clusters
 208
 209   fhPHOSRecParticles->Fill(kNumberOfPhosClusters);
 210   fhPHOSPhotons->Fill(fPhotonsInPhos);
 211   fhPHOSDigitsEvent->Fill(numberOfDigitsInPhos);
 212   fhPHOS->Fill(Entry(), numberOfDigitsInPhos, kNumberOfPhosClusters, fPhotonsInPhos) ;
 213
 214   // invariant Mass
 215   if (fPhotonsInPhos > 1 ) {
 216     Int_t phosPhoton1, phosPhoton2 ;
 217     for(phosPhoton1 = 0 ; phosPhoton1 < fPhotonsInPhos ; phosPhoton1++) {
 218       for(phosPhoton2 = phosPhoton1 + 1 ; phosPhoton2 < fPhotonsInPhos ; phosPhoton2++) {
 219         Float_t tempMass = TMath::Sqrt( 2 * phosPhotonsEnergy[phosPhoton1] * phosPhotonsEnergy[phosPhoton2] *
 220                                         ( 1 - TMath::Cos(phosVector[phosPhoton1]->Angle(*phosVector[phosPhoton2])) )
 221                                         );
 222         fhPHOSInvariantMass->Fill(tempMass*1000.);
 223       }
 224     }
 225   }
 226
 227   PostData(1, fOutputList);
 228
 229   delete [] phosVector ;
 230   delete [] phosPhotonsEnergy ;
 231
 232 }
 233
 234
 235 //______________________________________________________________________________
 236 void AliAnalysisTaskPHOSExample::Init()
 237 {
 238   // Intialisation of parameters
 239   AliInfo("Doing initialisation") ;
 240   SetPhotonId(0.9) ;
 241 }
 242
 243 //______________________________________________________________________________
 244 void AliAnalysisTaskPHOSExample::Terminate(Option_t *)
 245 {
 246   // Processing when the event loop is ended
 247
 248   Bool_t problem = kFALSE ;
 249   AliInfo(Form(" *** %s Report:", GetName())) ;
 250   printf("        PHOSEnergy Mean         : %5.3f , RMS : %5.3f \n", fhPHOSEnergy->GetMean(),         fhPHOSEnergy->GetRMS()         ) ;
 251   printf("        PHOSDigits Mean         : %5.3f , RMS : %5.3f \n", fhPHOSDigits->GetMean(),         fhPHOSDigits->GetRMS()         ) ;
 252   printf("        PHOSRecParticles Mean   : %5.3f , RMS : %5.3f \n", fhPHOSRecParticles->GetMean(),   fhPHOSRecParticles->GetRMS()   ) ;
 253   printf("        PHOSPhotons Mean        : %5.3f , RMS : %5.3f \n", fhPHOSPhotons->GetMean(),        fhPHOSPhotons->GetRMS()        ) ;
 254   printf("        PHOSInvariantMass Mean  : %5.3f , RMS : %5.3f \n", fhPHOSInvariantMass->GetMean(),  fhPHOSInvariantMass->GetRMS()  ) ;
 255   printf("        PHOSDigitsEvent Mean    : %5.3f , RMS : %5.3f \n", fhPHOSDigitsEvent->GetMean(),    fhPHOSDigitsEvent->GetRMS()    ) ;
 256
 257   TCanvas  * cPHOS = new TCanvas("cPHOS", "PHOS ESD Test", 400, 10, 600, 700) ;
 258   cPHOS->Divide(3, 2);
 259
 260   cPHOS->cd(1) ;
 261   if ( fhPHOSEnergy->GetMaximum() > 0. )
 262     gPad->SetLogy();
 263   fhPHOSEnergy->SetAxisRange(0, 25.);
 264   fhPHOSEnergy->SetLineColor(2);
 265   fhPHOSEnergy->Draw();
 266
 267   cPHOS->cd(2) ;
 268   fhPHOSDigits->SetAxisRange(0,25.);
 269   fhPHOSDigits->SetLineColor(2);
 270   fhPHOSDigits->Draw();
 271
 272   cPHOS->cd(3) ;
 273   if ( fhPHOSRecParticles->GetMaximum() > 0. )
 274     gPad->SetLogy();
 275   fhPHOSRecParticles->SetAxisRange(0, 25.);
 276   fhPHOSRecParticles->SetLineColor(2);
 277   fhPHOSRecParticles->Draw();
 278
 279   cPHOS->cd(4) ;
 280   if ( fhPHOSPhotons->GetMaximum() > 0. )
 281     gPad->SetLogy();
 282   fhPHOSPhotons->SetAxisRange(0,25.);
 283   fhPHOSPhotons->SetLineColor(2);
 284   fhPHOSPhotons->Draw();
 285
 286   cPHOS->cd(5) ;
 287   fhPHOSInvariantMass->SetLineColor(2);
 288   fhPHOSInvariantMass->Draw();
 289
 290   cPHOS->cd(6) ;
 291   if ( fhPHOSDigitsEvent->GetMaximum() > 0. )
 292     gPad->SetLogy();
 293   fhPHOSDigitsEvent->SetAxisRange(0,40.);
 294   fhPHOSDigitsEvent->SetLineColor(2);
 295   fhPHOSDigitsEvent->Draw();
 296
 297   cPHOS->Print("PHOS.eps");
 298
 299   char line[1024] ;
 300   sprintf(line, ".!tar -zcf %s.tar.gz *.eps", GetName()) ;
 301   gROOT->ProcessLine(line);
 302   sprintf(line, ".!rm -fR *.eps");
 303   gROOT->ProcessLine(line);
 304
 305   AliInfo(Form("!!! All the eps files are in %s.tar.gz !!!", GetName())) ;
 306
 307   char * report ;
 308   if(problem)
 309     report="Problems found, please check!!!";
 310   else
 311     report="OK";
 312
 313   AliInfo(Form("*** %s Summary Report: %s \n",GetName(), report)) ;
 314 }