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