[u/mrichter/AliRoot.git] / PWG4 / AliAnaGammaPhos.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id: */

//_________________________________________________________________________
// A basic analysis task to analyse photon detected by PHOS
// A basic analysis task to analyse photon detected by PHOS
// A basic analysis task to analyse photon detected by PHOS
// A basic analysis task to analyse photon detected by PHOS
// A basic analysis task to analyse photon detected by PHOS
//
//*-- Yves Schutz 
//////////////////////////////////////////////////////////////////////////////

#include <TCanvas.h>
#include <TChain.h>
#include <TFile.h> 
#include <TH1.h>
#include <TLegend.h> 
#include <TNtuple.h>
#include <TROOT.h> 
#include <TVector3.h> 

#include "AliAnaGammaPhos.h" 
#include "AliAnalysisManager.h"
#include "AliESDEvent.h" 
#include "AliESDCaloCluster.h" 
#include "AliAODEvent.h"
#include "AliAODHandler.h"
#include "AliLog.h"

//______________________________________________________________________________
AliAnaGammaPhos::AliAnaGammaPhos() : 
  fChain(0x0),
  fDebug(0),
  fESD(0x0),
  fAOD(0x0),
  fAODPhotons(0x0), 
  fPhotonsInPhos(0),
  fTreeA(0x0),
  fPhotonId(1.0),
  fOutputList(0x0), 
  fhPHOSPos(0),
  fhPHOS(0),
  fhPHOSEnergy(0),
  fhPHOSDigits(0),
  fhPHOSRecParticles(0),
  fhPHOSPhotons(0),
  fhPHOSInvariantMass(0),
  fhPHOSDigitsEvent(0)
{
  //Default constructor
}
//______________________________________________________________________________
AliAnaGammaPhos::AliAnaGammaPhos(const char *name) : 
  AliAnalysisTask(name,""),  
  fChain(0x0),
  fDebug(0),
  fESD(0x0),
  fAOD(0x0),
  fAODPhotons(0x0), 
  fPhotonsInPhos(0),
  fTreeA(0x0),
  fPhotonId(1.0),
  fOutputList(0x0), 
  fhPHOSPos(0),
  fhPHOS(0),
  fhPHOSEnergy(0),
  fhPHOSDigits(0),
  fhPHOSRecParticles(0),
  fhPHOSPhotons(0),
  fhPHOSInvariantMass(0),
  fhPHOSDigitsEvent(0)
{
  // Constructor.
  // Input slot #0 
  DefineInput(0, TChain::Class());
  // Output slots 
  DefineOutput(0,  TTree::Class()) ; 
  DefineOutput(1,  TList::Class()) ; 
}

//____________________________________________________________________________
AliAnaGammaPhos::AliAnaGammaPhos(const AliAnaGammaPhos& ap) :
  AliAnalysisTask(ap.GetName(),""),  
  fChain(ap.fChain),
  fDebug(ap.fDebug),
  fESD(ap.fESD),
  fAOD(ap.fAOD),
  fAODPhotons(ap.fAODPhotons), 
  fPhotonsInPhos(ap.fPhotonsInPhos),
  fTreeA(ap.fTreeA),
  fPhotonId(ap.fPhotonId),
  fOutputList(ap.fOutputList), 
  fhPHOSPos(ap.fhPHOSPos),
  fhPHOS(ap.fhPHOS),
  fhPHOSEnergy(ap.fhPHOSEnergy),
  fhPHOSDigits(ap.fhPHOSDigits),
  fhPHOSRecParticles(ap.fhPHOSRecParticles),
  fhPHOSPhotons(ap.fhPHOSPhotons),
  fhPHOSInvariantMass(ap.fhPHOSInvariantMass),
  fhPHOSDigitsEvent(ap.fhPHOSDigitsEvent)
{ 
  // cpy ctor
}

//_____________________________________________________________________________
AliAnaGammaPhos& AliAnaGammaPhos::operator = (const AliAnaGammaPhos& ap)
{
// assignment operator

  this->~AliAnaGammaPhos();
  new(this) AliAnaGammaPhos(ap);
  return *this;
}

//______________________________________________________________________________
AliAnaGammaPhos::~AliAnaGammaPhos()
{
  // dtor
  //  fOutputList->Clear() ; 
  //delete fOutputList ;
}


//______________________________________________________________________________
void AliAnaGammaPhos::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 AliAnaGammaPhos::CreateOutputObjects()
{  
  // Create the outputs containers
 
  OpenFile(0) ;
  AliAODHandler* handler = (AliAODHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());  
  fTreeA = handler->GetTree() ; 
  fAOD   = handler->GetAOD();
  fAODPhotons = fAOD->GetCaloClusters() ; 
  

  OpenFile(1) ; 

  fhPHOSPos            = new TNtuple("PHOSPos"         , "Position in PHOS"  , "x:y:z");
  fhPHOS               = new TNtuple("PHOS"            , "PHOS"  , "event:digits:clusters:photons");
  fhPHOSEnergy         = new TH1D("PHOSEnergy"         , "PHOSEnergy"        , 100, 0., 100. ) ;
  fhPHOSDigits         = new TH1I("PHOSDigitsCluster"  , "PHOSDigits"        , 20 , 0 , 20  ) ;
  fhPHOSRecParticles   = new TH1D("PHOSRecParticles"   , "PHOSRecParticles" , 20 , 0., 20. ) ;
  fhPHOSPhotons        = new TH1I("PHOSPhotons"        , "PHOSPhotons"       , 20 , 0 , 20  ) ;
  fhPHOSInvariantMass  = new TH1D("PHOSInvariantMass"  , "PHOSInvariantMass" , 400, 0., 400.) ;
  fhPHOSDigitsEvent    = new TH1I("PHOSDigitsEvent"    , "PHOSDigitsEvent"   , 30 , 0 , 30  ) ;
  
  // create output container
  
  fOutputList = new TList() ; 
  fOutputList->SetName(GetName()) ; 

  fOutputList->AddAt(fhPHOSPos,             0) ; 
  fOutputList->AddAt(fhPHOS,                1) ; 
  fOutputList->AddAt(fhPHOSEnergy,          2) ; 
  fOutputList->AddAt(fhPHOSDigits,          3) ; 
  fOutputList->AddAt(fhPHOSRecParticles,    4) ; 
  fOutputList->AddAt(fhPHOSPhotons,         5) ; 
  fOutputList->AddAt(fhPHOSInvariantMass,   6) ; 
  fOutputList->AddAt(fhPHOSDigitsEvent,     7) ; 
}

//______________________________________________________________________________
void AliAnaGammaPhos::Exec(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)) ; 
  
  //************************  PHOS *************************************
  TRefArray * caloClustersArr  = new TRefArray();  
  fESD->GetPHOSClusters(caloClustersArr);

  const Int_t kNumberOfPhosClusters   = caloClustersArr->GetEntries() ;  
  
  TVector3 ** phosVector       = new TVector3*[kNumberOfPhosClusters] ;
  Float_t  * phosPhotonsEnergy = new Float_t[kNumberOfPhosClusters] ;
  Int_t      phosCluster ; 
  Int_t      numberOfDigitsInPhos   = 0 ;
  
  fPhotonsInPhos  = 0 ;
  // loop over the PHOS Cluster
  for(phosCluster = 0 ; phosCluster < kNumberOfPhosClusters ; phosCluster++) {
    AliESDCaloCluster * caloCluster = (AliESDCaloCluster *) caloClustersArr->At(phosCluster) ;
    //AliESDCaloCluster * caloCluster = fESD->GetCaloCluster(phosCluster) ;
    if (caloCluster) {
      Float_t pos[3] ;
      caloCluster->GetPosition( pos ) ;
      fhPHOSEnergy->Fill( caloCluster->E() ) ;
      fhPHOSPos->Fill( pos[0], pos[1], pos[2] ) ;
      fhPHOSDigits->Fill(entry, caloCluster->GetNCells() ) ;
      numberOfDigitsInPhos += caloCluster->GetNCells() ;
      Double_t * pid = caloCluster->GetPid() ;
      if(pid[AliPID::kPhoton] > GetPhotonId() ) {
	phosVector[fPhotonsInPhos] = new TVector3(pos[0],pos[1],pos[2]) ;
	phosPhotonsEnergy[fPhotonsInPhos]=caloCluster->E() ;
        //new ((*fAODPhotons)[fPhotonsInPhos++;]) AliAODPhoton ( );
      }
    }
  } //PHOS clusters
    
  fhPHOSRecParticles->Fill(kNumberOfPhosClusters);
  fhPHOSPhotons->Fill(fPhotonsInPhos);
  fhPHOSDigitsEvent->Fill(numberOfDigitsInPhos);
  fhPHOS->Fill(entry, numberOfDigitsInPhos, kNumberOfPhosClusters, fPhotonsInPhos) ; 

  // invariant Mass
  if (fPhotonsInPhos > 1 ) {
    Int_t phosPhoton1, phosPhoton2 ; 
    for(phosPhoton1 = 0 ; phosPhoton1 < fPhotonsInPhos ; phosPhoton1++) {
      for(phosPhoton2 = phosPhoton1 + 1 ; phosPhoton2 < fPhotonsInPhos ; phosPhoton2++) {      
	Float_t tempMass = TMath::Sqrt( 2 * phosPhotonsEnergy[phosPhoton1] * phosPhotonsEnergy[phosPhoton2] *
					( 1 - TMath::Cos(phosVector[phosPhoton1]->Angle(*phosVector[phosPhoton2])) ) 
					);
	fhPHOSInvariantMass->Fill(tempMass*1000.);
      }
    }    
  }
    
  PostData(0, fTreeA) ; 
  PostData(1, fOutputList);

  delete [] phosVector ; 
  delete [] phosPhotonsEnergy ; 
  
}


//______________________________________________________________________________
void AliAnaGammaPhos::Init()
{
  // Intialisation of parameters
  AliInfo("Doing initialisation") ; 
  SetPhotonId(0.9) ; 
}

//______________________________________________________________________________
void AliAnaGammaPhos::Terminate(Option_t *)
{
  // Processing when the event loop is ended

//   fOutputList = (TList*)GetOutputData(0);  
//   fhPHOSPos            = (TNtuple*)fOutputList->At(0);
//   fhPHOS               = (TNtuple*)fOutputList->At(1);
//   fhPHOSEnergy         = (TH1D*)fOutputList->At(2);
//   fhPHOSDigits         = (TH1I*)fOutputList->At(3);
//   fhPHOSRecParticles   = (TH1D*)fOutputList->At(4);
//   fhPHOSPhotons        = (TH1I*)fOutputList->At(5);
//   fhPHOSInvariantMass  = (TH1D*)fOutputList->At(6);
//   fhPHOSDigitsEvent    = (TH1I*)fOutputList->At(7);

  
  Bool_t problem = kFALSE ; 
  AliInfo(Form(" *** %s Report:", GetName())) ; 
  printf("        PHOSEnergy Mean         : %5.3f , RMS : %5.3f \n", fhPHOSEnergy->GetMean(),         fhPHOSEnergy->GetRMS()         ) ;
  printf("        PHOSDigits Mean         : %5.3f , RMS : %5.3f \n", fhPHOSDigits->GetMean(),         fhPHOSDigits->GetRMS()         ) ;
  printf("        PHOSRecParticles Mean   : %5.3f , RMS : %5.3f \n", fhPHOSRecParticles->GetMean(),   fhPHOSRecParticles->GetRMS()   ) ;
  printf("        PHOSPhotons Mean        : %5.3f , RMS : %5.3f \n", fhPHOSPhotons->GetMean(),        fhPHOSPhotons->GetRMS()        ) ;
  printf("        PHOSInvariantMass Mean  : %5.3f , RMS : %5.3f \n", fhPHOSInvariantMass->GetMean(),  fhPHOSInvariantMass->GetRMS()  ) ;
  printf("        PHOSDigitsEvent Mean    : %5.3f , RMS : %5.3f \n", fhPHOSDigitsEvent->GetMean(),    fhPHOSDigitsEvent->GetRMS()    ) ;

  TCanvas  * cPHOS = new TCanvas("cPHOS", "PHOS ESD Test", 400, 10, 600, 700) ;
  cPHOS->Divide(3, 2);

  cPHOS->cd(1) ; 
  if ( fhPHOSEnergy->GetMaximum() > 0. ) 
    gPad->SetLogy();
  fhPHOSEnergy->SetAxisRange(0, 25.);
  fhPHOSEnergy->SetLineColor(2);
  fhPHOSEnergy->Draw();

  cPHOS->cd(2) ; 
  fhPHOSDigits->SetAxisRange(0,25.);
  fhPHOSDigits->SetLineColor(2);
  fhPHOSDigits->Draw();

  cPHOS->cd(3) ; 
  if ( fhPHOSRecParticles->GetMaximum() > 0. ) 
    gPad->SetLogy();
  fhPHOSRecParticles->SetAxisRange(0, 25.);
  fhPHOSRecParticles->SetLineColor(2);
  fhPHOSRecParticles->Draw();

  cPHOS->cd(4) ; 
  if ( fhPHOSPhotons->GetMaximum() > 0. ) 
    gPad->SetLogy();
  fhPHOSPhotons->SetAxisRange(0,25.);
  fhPHOSPhotons->SetLineColor(2);
  fhPHOSPhotons->Draw();

  cPHOS->cd(5) ; 
  fhPHOSInvariantMass->SetLineColor(2);
  fhPHOSInvariantMass->Draw();
 
  cPHOS->cd(6) ; 
  if ( fhPHOSDigitsEvent->GetMaximum() > 0. ) 
    gPad->SetLogy();
  fhPHOSDigitsEvent->SetAxisRange(0,40.);
  fhPHOSDigitsEvent->SetLineColor(2);
  fhPHOSDigitsEvent->Draw();
 
  cPHOS->Print("PHOS.eps");
 
  char line[1024] ; 
  sprintf(line, ".!tar -zcf %s.tar.gz *.eps", GetName()) ; 
  gROOT->ProcessLine(line);
  sprintf(line, ".!rm -fR *.eps"); 
  gROOT->ProcessLine(line);
 
  AliInfo(Form("!!! All the eps files are in %s.tar.gz !!!", GetName())) ;

  char * report ; 
  if(problem)
    report="Problems found, please check!!!";  
  else 
    report="OK";

  AliInfo(Form("*** %s Summary Report: %s \n",GetName(), report)) ; 
}
Commit	Line	Data
ce9d0223	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
d40a7452	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
ce9d0223	24	//
	25	//*-- Yves Schutz
	26	//////////////////////////////////////////////////////////////////////////////
	27
	28	#include <TCanvas.h>
	29	#include <TChain.h>
	30	#include <TFile.h>
	31	#include <TH1.h>
ce9d0223	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"
9ea45617	39	#include "AliESDEvent.h"
9ea45617	40	#include "AliESDCaloCluster.h"
ce9d0223	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),
3bb2c538	56	fhPHOSPos(0),
ce9d0223	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),
3bb2c538	79	fhPHOSPos(0),
ce9d0223	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
d40a7452	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
ce9d0223	130	//______________________________________________________________________________
	131	AliAnaGammaPhos::~AliAnaGammaPhos()
	132	{
	133	// dtor
c0acc236	134	// fOutputList->Clear() ;
c0acc236	135	//delete fOutputList ;
ce9d0223	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
9ea45617	153	fESD = new AliESDEvent();
	154	fESD->ReadFromTree(fChain);
	155
ce9d0223	156	}
	157
	158	//________________________________________________________________________
	159	void AliAnaGammaPhos::CreateOutputObjects()
	160	{
	161	// Create the outputs containers
	162
	163	OpenFile(0) ;
691685d6	164	AliAODHandler* handler = (AliAODHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
c0acc236	165	fTreeA = handler->GetTree() ;
ce9d0223	166	fAOD = handler->GetAOD();
4b707925	167	fAODPhotons = fAOD->GetCaloClusters() ;
c0acc236	168
ce9d0223	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");
bdcfac30	174	fhPHOSEnergy = new TH1D("PHOSEnergy" , "PHOSEnergy" , 100, 0., 100. ) ;
ce9d0223	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
c427aac7	200
ce9d0223	201	Long64_t entry = fChain->GetReadEntry() ;
	202
	203	if (!fESD) {
	204	AliError("fESD is not connected to the input!") ;
	205	return ;
	206	}
	207
c427aac7	208	// if ( !((entry-1)%100) )
c427aac7	209	AliInfo(Form("%s ----> Processing event # %lld", (dynamic_cast<TChain *>(fChain))->GetFile()->GetName(), entry)) ;
ce9d0223	210
ce9d0223	211	//********************** PHOS ***********************************
c427aac7	212	TRefArray * caloClustersArr = new TRefArray();
	213	fESD->GetPHOSClusters(caloClustersArr);
	214
	215	const Int_t kNumberOfPhosClusters = caloClustersArr->GetEntries() ;
ce9d0223	216
d40a7452	217	TVector3 ** phosVector = new TVector3*[kNumberOfPhosClusters] ;
d40a7452	218	Float_t * phosPhotonsEnergy = new Float_t[kNumberOfPhosClusters] ;
ce9d0223	219	Int_t phosCluster ;
	220	Int_t numberOfDigitsInPhos = 0 ;
	221
	222	fPhotonsInPhos = 0 ;
	223	// loop over the PHOS Cluster
c427aac7	224	for(phosCluster = 0 ; phosCluster < kNumberOfPhosClusters ; phosCluster++) {
	225	AliESDCaloCluster * caloCluster = (AliESDCaloCluster *) caloClustersArr->At(phosCluster) ;
	226	//AliESDCaloCluster * caloCluster = fESD->GetCaloCluster(phosCluster) ;
ce9d0223	227	if (caloCluster) {
ce9d0223	228	Float_t pos[3] ;
c0acc236	229	caloCluster->GetPosition( pos ) ;
c0acc236	230	fhPHOSEnergy->Fill( caloCluster->E() ) ;
ce9d0223	231	fhPHOSPos->Fill( pos[0], pos[1], pos[2] ) ;
c427aac7	232	fhPHOSDigits->Fill(entry, caloCluster->GetNCells() ) ;
c427aac7	233	numberOfDigitsInPhos += caloCluster->GetNCells() ;
4b707925	234	Double_t * pid = caloCluster->GetPid() ;
ce9d0223	235	if(pid[AliPID::kPhoton] > GetPhotonId() ) {
ce9d0223	236	phosVector[fPhotonsInPhos] = new TVector3(pos[0],pos[1],pos[2]) ;
c0acc236	237	phosPhotonsEnergy[fPhotonsInPhos]=caloCluster->E() ;
ce9d0223	238	//new ((*fAODPhotons)[fPhotonsInPhos++;]) AliAODPhoton ( );
	239	}
	240	}
	241	} //PHOS clusters
	242
d40a7452	243	fhPHOSRecParticles->Fill(kNumberOfPhosClusters);
ce9d0223	244	fhPHOSPhotons->Fill(fPhotonsInPhos);
ce9d0223	245	fhPHOSDigitsEvent->Fill(numberOfDigitsInPhos);
d40a7452	246	fhPHOS->Fill(entry, numberOfDigitsInPhos, kNumberOfPhosClusters, fPhotonsInPhos) ;
ce9d0223	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	}
32392f88	260
ce9d0223	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
75480d65	278	//______________________________________________________________________________
ce9d0223	279	void AliAnaGammaPhos::Terminate(Option_t *)
	280	{
	281	// Processing when the event loop is ended
	282
c0acc236	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);
ce9d0223	292
32392f88	293
ce9d0223	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	}