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   TRefArray * caloClustersArr  = new TRefArray();
 213   fESD->GetPHOSClusters(caloClustersArr);
 214
 215   const Int_t kNumberOfPhosClusters   = caloClustersArr->GetEntries() ;
 216
 217   TVector3 ** phosVector       = new TVector3*[kNumberOfPhosClusters] ;
 218   Float_t  * phosPhotonsEnergy = new Float_t[kNumberOfPhosClusters] ;
 219   Int_t      phosCluster ;
 220   Int_t      numberOfDigitsInPhos   = 0 ;
 221
 222   fPhotonsInPhos  = 0 ;
 223   // loop over the PHOS Cluster
 224   for(phosCluster = 0 ; phosCluster < kNumberOfPhosClusters ; phosCluster++) {
 225     AliESDCaloCluster * caloCluster = (AliESDCaloCluster *) caloClustersArr->At(phosCluster) ;
 226     //AliESDCaloCluster * caloCluster = fESD->GetCaloCluster(phosCluster) ;
 227     if (caloCluster) {
 228       Float_t pos[3] ;
 229       caloCluster->GetPosition( pos ) ;
 230       fhPHOSEnergy->Fill( caloCluster->E() ) ;
 231       fhPHOSPos->Fill( pos[0], pos[1], pos[2] ) ;
 232       fhPHOSDigits->Fill(entry, caloCluster->GetNCells() ) ;
 233       numberOfDigitsInPhos += caloCluster->GetNCells() ;
 234       Double_t * pid = caloCluster->GetPid() ;
 235       if(pid[AliPID::kPhoton] > GetPhotonId() ) {
 236         phosVector[fPhotonsInPhos] = new TVector3(pos[0],pos[1],pos[2]) ;
 237         phosPhotonsEnergy[fPhotonsInPhos]=caloCluster->E() ;
 238         //new ((*fAODPhotons)[fPhotonsInPhos++;]) AliAODPhoton ( );
 239       }
 240     }
 241   } //PHOS clusters
 242
 243   fhPHOSRecParticles->Fill(kNumberOfPhosClusters);
 244   fhPHOSPhotons->Fill(fPhotonsInPhos);
 245   fhPHOSDigitsEvent->Fill(numberOfDigitsInPhos);
 246   fhPHOS->Fill(entry, numberOfDigitsInPhos, kNumberOfPhosClusters, fPhotonsInPhos) ;
 247
 248   // invariant Mass
 249   if (fPhotonsInPhos > 1 ) {
 250     Int_t phosPhoton1, phosPhoton2 ;
 251     for(phosPhoton1 = 0 ; phosPhoton1 < fPhotonsInPhos ; phosPhoton1++) {
 252       for(phosPhoton2 = phosPhoton1 + 1 ; phosPhoton2 < fPhotonsInPhos ; phosPhoton2++) {
 253         Float_t tempMass = TMath::Sqrt( 2 * phosPhotonsEnergy[phosPhoton1] * phosPhotonsEnergy[phosPhoton2] *
 254                                         ( 1 - TMath::Cos(phosVector[phosPhoton1]->Angle(*phosVector[phosPhoton2])) )
 255                                         );
 256         fhPHOSInvariantMass->Fill(tempMass*1000.);
 257       }
 258     }
 259   }
 260
 261   PostData(0, fTreeA) ;
 262   PostData(1, fOutputList);
 263
 264   delete [] phosVector ;
 265   delete [] phosPhotonsEnergy ;
 266
 267 }
 268
 269
 270 //______________________________________________________________________________
 271 void AliAnaGammaPhos::Init()
 272 {
 273   // Intialisation of parameters
 274   AliInfo("Doing initialisation") ;
 275   SetPhotonId(0.9) ;
 276 }
 277
 278 //______________________________________________________________________________
 279 void AliAnaGammaPhos::Terminate(Option_t *)
 280 {
 281   // Processing when the event loop is ended
 282
 283 //   fOutputList = (TList*)GetOutputData(0);
 284 //   fhPHOSPos            = (TNtuple*)fOutputList->At(0);
 285 //   fhPHOS               = (TNtuple*)fOutputList->At(1);
 286 //   fhPHOSEnergy         = (TH1D*)fOutputList->At(2);
 287 //   fhPHOSDigits         = (TH1I*)fOutputList->At(3);
 288 //   fhPHOSRecParticles   = (TH1D*)fOutputList->At(4);
 289 //   fhPHOSPhotons        = (TH1I*)fOutputList->At(5);
 290 //   fhPHOSInvariantMass  = (TH1D*)fOutputList->At(6);
 291 //   fhPHOSDigitsEvent    = (TH1I*)fOutputList->At(7);
 292
 293
 294   Bool_t problem = kFALSE ;
 295   AliInfo(Form(" *** %s Report:", GetName())) ;
 296   printf("        PHOSEnergy Mean         : %5.3f , RMS : %5.3f \n", fhPHOSEnergy->GetMean(),         fhPHOSEnergy->GetRMS()         ) ;
 297   printf("        PHOSDigits Mean         : %5.3f , RMS : %5.3f \n", fhPHOSDigits->GetMean(),         fhPHOSDigits->GetRMS()         ) ;
 298   printf("        PHOSRecParticles Mean   : %5.3f , RMS : %5.3f \n", fhPHOSRecParticles->GetMean(),   fhPHOSRecParticles->GetRMS()   ) ;
 299   printf("        PHOSPhotons Mean        : %5.3f , RMS : %5.3f \n", fhPHOSPhotons->GetMean(),        fhPHOSPhotons->GetRMS()        ) ;
 300   printf("        PHOSInvariantMass Mean  : %5.3f , RMS : %5.3f \n", fhPHOSInvariantMass->GetMean(),  fhPHOSInvariantMass->GetRMS()  ) ;
 301   printf("        PHOSDigitsEvent Mean    : %5.3f , RMS : %5.3f \n", fhPHOSDigitsEvent->GetMean(),    fhPHOSDigitsEvent->GetRMS()    ) ;
 302
 303   TCanvas  * cPHOS = new TCanvas("cPHOS", "PHOS ESD Test", 400, 10, 600, 700) ;
 304   cPHOS->Divide(3, 2);
 305
 306   cPHOS->cd(1) ;
 307   if ( fhPHOSEnergy->GetMaximum() > 0. )
 308     gPad->SetLogy();
 309   fhPHOSEnergy->SetAxisRange(0, 25.);
 310   fhPHOSEnergy->SetLineColor(2);
 311   fhPHOSEnergy->Draw();
 312
 313   cPHOS->cd(2) ;
 314   fhPHOSDigits->SetAxisRange(0,25.);
 315   fhPHOSDigits->SetLineColor(2);
 316   fhPHOSDigits->Draw();
 317
 318   cPHOS->cd(3) ;
 319   if ( fhPHOSRecParticles->GetMaximum() > 0. )
 320     gPad->SetLogy();
 321   fhPHOSRecParticles->SetAxisRange(0, 25.);
 322   fhPHOSRecParticles->SetLineColor(2);
 323   fhPHOSRecParticles->Draw();
 324
 325   cPHOS->cd(4) ;
 326   if ( fhPHOSPhotons->GetMaximum() > 0. )
 327     gPad->SetLogy();
 328   fhPHOSPhotons->SetAxisRange(0,25.);
 329   fhPHOSPhotons->SetLineColor(2);
 330   fhPHOSPhotons->Draw();
 331
 332   cPHOS->cd(5) ;
 333   fhPHOSInvariantMass->SetLineColor(2);
 334   fhPHOSInvariantMass->Draw();
 335
 336   cPHOS->cd(6) ;
 337   if ( fhPHOSDigitsEvent->GetMaximum() > 0. )
 338     gPad->SetLogy();
 339   fhPHOSDigitsEvent->SetAxisRange(0,40.);
 340   fhPHOSDigitsEvent->SetLineColor(2);
 341   fhPHOSDigitsEvent->Draw();
 342
 343   cPHOS->Print("PHOS.eps");
 344
 345   char line[1024] ;
 346   sprintf(line, ".!tar -zcf %s.tar.gz *.eps", GetName()) ;
 347   gROOT->ProcessLine(line);
 348   sprintf(line, ".!rm -fR *.eps");
 349   gROOT->ProcessLine(line);
 350
 351   AliInfo(Form("!!! All the eps files are in %s.tar.gz !!!", GetName())) ;
 352
 353   char * report ;
 354   if(problem)
 355     report="Problems found, please check!!!";
 356   else
 357     report="OK";
 358
 359   AliInfo(Form("*** %s Summary Report: %s \n",GetName(), report)) ;
 360 }