[u/mrichter/AliRoot.git] / PHOS / AnaESD.C

/**************************************************************************
 * 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$ */
//_________________________________________________________________________
// Macros analyzing the ESD file
// Use Case : 
//          root> .L AnaESD.C++
//          root> ana() --> prints the objects stored in ESD
//                                              
// author  : Yves Schutz (CERN/SUBATECH)
// February 2004
//_________________________________________________________________________
#if !defined(__CINT__) || defined(__MAKECINT__)
#include "TFile.h"
#include "TMath.h"
#include "TH1D.h"
#include "TROOT.h"
#include "TBrowser.h"
#include "AliPHOSGetter.h"
#include "AliPHOSGeometry.h"
#include "Riostream.h"
#include "AliESD.h"
#include "AliESDtrack.h"
#include "AliESDCaloTrack.h"
#include "AliEMCALRecParticle.h"
#include "AliPHOSRecParticle.h"
#include "AliKalmanTrack.h"
#include "AliPHOSGridFile.h"
#endif

Bool_t AnaESD(TString filename) ; 
void   Match(AliESDtrack * ct, AliESDtrack * cp, Double_t * dist) ; 


//__________________________________________________________________________
Bool_t Ana(const TString type = "per5", const Int_t run = 1, const Int_t nOfEvt = 1) 
{ 
  // Analyzes data from the AliEn data catalog

  // Data Challenge identification
  const TString kYear("2004") ; 
  const TString kProd("02") ; 
  const TString kVers("V4.01.Rev.00") ; 

  // get the LFN file name in the Grid catalogue ; 
  AliPHOSGridFile lfn ; 
  if (!lfn.IsConnected()) 
    return kFALSE ; 
  lfn.SetPath(kYear, kProd, kVers, type) ;  
  lfn.SetRun(run) ; 

  //loop over the events 
  Int_t nevt, evt = 0 ; 
  for (nevt = 0 ; nevt < nOfEvt ; nevt++) { 
    evt++ ; 
    lfn.SetEvt(evt) ;
    TString fileName = lfn.GetLFN() ; 
    
    if (fileName.IsNull()) {
      nevt-- ; 
      continue ; 
    }

    printf(">>>>>>>>>>>> Processing %s-%s/%s/%s : run # %d event # %d \n", 
	   kYear.Data(), kProd.Data(), kVers.Data(), type.Data(), run, evt) ;
    AnaESD(fileName) ;
  } 
  return kTRUE ; 
}

//__________________________________________________________________________
Bool_t AnaESD(TString fileName) 
{
  // Analyze ESD from file fileName
  // calo track histograms ;
  
  TFile * out = new TFile("AOD.root", "RECREATE") ; 
  TH1D * hCaloEnergyA = 0 ; 
  TH1D * hCaloEnergyE = 0 ; 
  TH1D * hCaloEnergyG = 0 ; 
  TH1D * hCaloEnergyP = 0 ; 
  TH1D * heta = 0 ;
  TH1D * hphi = 0 ; 
  Double_t dist[3] ; 

  AliPHOSGetter * gime = AliPHOSGetter::Instance(fileName) ; 
  
  Int_t nEvent = gime->MaxEvent() ;  
  Int_t event ; 
  AliESD * esd = 0 ;
  for (event = 0 ; event < nEvent; event++) {
    cout << "AnaESD : Processing event # " << event << endl ; 
    esd = gime->ESD(event) ; 
    if (!esd) 
      return kFALSE ; 
    
    //esd->Print();  
    // Calorimeter tracks 
    AliESDtrack * ct ; 
    Int_t caloindex ; 
    for (caloindex = 0 ; caloindex < esd->GetNumberOfTracks() ; caloindex++) {
      // get the tracks and check if it is from PHOS 
      ct = esd->GetTrack(caloindex) ;
      if ( !ct->IsPHOS() ) 
	continue ; 
      Double_t energy = ct->GetPHOSsignal() ; 
      // check if CPV bit is set (see AliPHOSFastRecParticle) 
      Double_t type[AliESDtrack::kSPECIES+4] ;
      ct->GetPHOSpid(type) ; 
      
      if ( (type[AliESDtrack::kElectron] == 1.0) || (type[AliESDtrack::kPhoton] == 1.0) ) {
	if(!hCaloEnergyA) {
	  out->cd() ; 
	  hCaloEnergyA = new TH1D("CaloEnergyA", "Energy in calorimeter Electron/Photon like", 500, 0., 50.) ;
	}
	hCaloEnergyA->Fill(energy) ;  
      }
      
      if ( type[AliESDtrack::kElectron] == 1.0 ) {
	if(!hCaloEnergyE) {
	  out->cd() ; 
	  hCaloEnergyE = new TH1D("CaloEnergyE", "Energy in calorimeter Electron like", 500, 0., 50.) ;
	}
	hCaloEnergyE->Fill(energy) ;  
      } 
  
      if ( type[AliESDtrack::kPhoton] == 1.0 ) {
	if(!hCaloEnergyG) {
	  out->cd() ; 
	  hCaloEnergyG = new TH1D("CaloEnergyG", "Energy in calorimeter Gamma like", 500, 0., 50.) ;
	}
	hCaloEnergyG->Fill(energy) ;
      }  
     
      if ( type[AliESDtrack::kPi0] == 1.0 ) {
	if(!hCaloEnergyP) {
	  out->cd() ; 
	  hCaloEnergyP = new TH1D("CaloEnergyP", "Energy in calorimeter Pi0 like", 500, 0., 100.) ;
	}
	hCaloEnergyP->Fill(energy) ;
      }  
      
      AliESDtrack * cp ; 
      Int_t cpindex ; 
      for (cpindex = 0 ; cpindex < esd->GetNumberOfTracks() ; cpindex++) {
	// get the charged tracks from central tracking
	cp = esd->GetTrack(cpindex) ;
	if ( cp->IsPHOS() ) 
	  continue ; 
	Match(ct, cp, dist) ; 
	if (!heta && !hphi) {
	  heta = new TH1D("Correta", "neutral-charged correlation in eta" , 100, 0., 360.) ; 
	  hphi = new TH1D("Corrphi", "neutral-charged correlation in phi" , 100, 0., 360.) ; 
	}
	heta->Fill(dist[1]) ; 
	hphi->Fill(dist[2]) ;
      }
    }
  }
  TBrowser * bs = new TBrowser("Root Memory Bowser", gROOT->FindObjectAny("AOD.root") ) ;
  bs->Show() ;
  out->Write() ; 
  return kTRUE ; 
}

//__________________________________________________________________________
void Match(AliESDtrack * ct, AliESDtrack * cp, Double_t * dist) 
{
  // Calculates the distance (x,z) between  the particle detected by PHOS and 
  // the charged particle reconstructed by the global tracking 

      
 //  Int_t phN ; 
//   Double_t phZ, phX ; 
  
//   if (ct->IsPHOS()) { // it is a PHOS particle 
//     Double_t cpTheta,  cpPhi ;  
//     Double_t phTheta,  phPhi ; 
//     cpTheta = cpPhi = phTheta = phPhi = 0. ; 
//    //    cout << "PHOS particle # " << " pos (" 
//     // 	 << pp->GetPos().X() << ", " << pp->GetPos().Y() << ", " << pp->GetPos().Z() << ")" << endl ;
    
//     AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
//     gime->PHOSGeometry()->ImpactOnEmc(*pp, phN, phZ, phX) ; 
//     Double_t xyzAtPHOS[3] ; 
//     cp->GetOuterXYZ(xyzAtPHOS) ; 
//     if ( (xyzAtPHOS[0] +  xyzAtPHOS[1] + xyzAtPHOS[2]) != 0.) { //it has reached PHOS
//       //the next check are only if we want high quality tracks 
//       //       ULong_t status = cp->GetStatus() ;  
//       //       if ((status & AliESDtrack::kTRDput)==0) 
//       // 	do not continue;
//       //       if ((status & AliESDtrack::kTRDStop)!=0) 
//       // 	do not continue;  
//       //       cout << "Charged particle # " << " pos (" 
//       // 	   << xyzAtPHOS[0] << ", " << xyzAtPHOS[1] << ", " << xyzAtPHOS[2] << ")" <<  endl ;     
//       TVector3 poscp(xyzAtPHOS[0], xyzAtPHOS[1], xyzAtPHOS[2]) ;
//       Int_t cpN ;
//       Double_t cpZ,cpX ; 
//       gime->PHOSGeometry()->ImpactOnEmc(poscp, cpN, cpZ, cpX) ; 
//       if (cpN) {// we are inside the PHOS acceptance 
// 	// 	cout << "Charged Matching 1: " << cpN << " " << cpZ << " " << cpX << endl ; 
// 	// 	cout << "Charged Matching 2: " << phN << " " << phZ << " " << phX << endl ; 
// 	dist[0] = TMath::Sqrt( (cpZ-phZ)*(cpZ-phZ) + (cpX-phX)*(cpX-phX)) ;  
//       } 
//       phTheta = pp->Theta() ; 
//       phPhi   = pp->Phi() ;
//       TParticle tempo ; 
//       tempo.SetMomentum(xyzAtPHOS[0], xyzAtPHOS[1], xyzAtPHOS[2], 0.) ;  
//       cpTheta = tempo.Theta() ; 
//       cpPhi   = tempo.Phi() ;
//       //cout << phTheta << " " << phPhi << " " << endl 
//       //cout <<	 cpTheta << " " << cpPhi-phPhi << " " << endl ; 
//     }
//     dist[1] = (phTheta - cpTheta)*TMath::RadToDeg() ; 
//     dist[2] = (phPhi - cpPhi)*TMath::RadToDeg() ; 
//   }
  
//   if (ep) {
//     //cout << "EMCAL particle # " << endl ; 
//   }
}
Commit	Line	Data
	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	/* $Id$ */
	16	//_________________________________________________________________________
	17	// Macros analyzing the ESD file
	18	// Use Case :
	19	// root> .L AnaESD.C++
	20	// root> ana() --> prints the objects stored in ESD
	21	//
	22	// author : Yves Schutz (CERN/SUBATECH)
	23	// February 2004
	24	//_________________________________________________________________________
	25	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
	26	#include "TFile.h"
	27	#include "TMath.h"
	28	#include "TH1D.h"
	29	#include "TROOT.h"
	30	#include "TBrowser.h"
	31	#include "AliPHOSGetter.h"
	32	#include "AliPHOSGeometry.h"
	33	#include "Riostream.h"
	34	#include "AliESD.h"
	35	#include "AliESDtrack.h"
	36	#include "AliESDCaloTrack.h"
	37	#include "AliEMCALRecParticle.h"
	38	#include "AliPHOSRecParticle.h"
	39	#include "AliKalmanTrack.h"
	40	#include "AliPHOSGridFile.h"
	41	#endif
	42
	43	Bool_t AnaESD(TString filename) ;
	44	void Match(AliESDtrack * ct, AliESDtrack * cp, Double_t * dist) ;
	45
	46
	47	//__________________________________________________________________________
	48	Bool_t Ana(const TString type = "per5", const Int_t run = 1, const Int_t nOfEvt = 1)
	49	{
	50	// Analyzes data from the AliEn data catalog
	51
	52	// Data Challenge identification
	53	const TString kYear("2004") ;
	54	const TString kProd("02") ;
	55	const TString kVers("V4.01.Rev.00") ;
	56
	57	// get the LFN file name in the Grid catalogue ;
	58	AliPHOSGridFile lfn ;
	59	if (!lfn.IsConnected())
	60	return kFALSE ;
	61	lfn.SetPath(kYear, kProd, kVers, type) ;
	62	lfn.SetRun(run) ;
	63
	64	//loop over the events
	65	Int_t nevt, evt = 0 ;
	66	for (nevt = 0 ; nevt < nOfEvt ; nevt++) {
	67	evt++ ;
	68	lfn.SetEvt(evt) ;
	69	TString fileName = lfn.GetLFN() ;
	70
	71	if (fileName.IsNull()) {
	72	nevt-- ;
	73	continue ;
	74	}
	75
	76	printf(">>>>>>>>>>>> Processing %s-%s/%s/%s : run # %d event # %d \n",
	77	kYear.Data(), kProd.Data(), kVers.Data(), type.Data(), run, evt) ;
	78	AnaESD(fileName) ;
	79	}
	80	return kTRUE ;
	81	}
	82
	83	//__________________________________________________________________________
	84	Bool_t AnaESD(TString fileName)
	85	{
	86	// Analyze ESD from file fileName
	87	// calo track histograms ;
	88
	89	TFile * out = new TFile("AOD.root", "RECREATE") ;
	90	TH1D * hCaloEnergyA = 0 ;
	91	TH1D * hCaloEnergyE = 0 ;
	92	TH1D * hCaloEnergyG = 0 ;
	93	TH1D * hCaloEnergyP = 0 ;
	94	TH1D * heta = 0 ;
	95	TH1D * hphi = 0 ;
	96	Double_t dist[3] ;
	97
	98	AliPHOSGetter * gime = AliPHOSGetter::Instance(fileName) ;
	99
	100	Int_t nEvent = gime->MaxEvent() ;
	101	Int_t event ;
	102	AliESD * esd = 0 ;
	103	for (event = 0 ; event < nEvent; event++) {
	104	cout << "AnaESD : Processing event # " << event << endl ;
	105	esd = gime->ESD(event) ;
	106	if (!esd)
	107	return kFALSE ;
	108
	109	//esd->Print();
	110	// Calorimeter tracks
	111	AliESDtrack * ct ;
	112	Int_t caloindex ;
	113	for (caloindex = 0 ; caloindex < esd->GetNumberOfTracks() ; caloindex++) {
	114	// get the tracks and check if it is from PHOS
	115	ct = esd->GetTrack(caloindex) ;
	116	if ( !ct->IsPHOS() )
	117	continue ;
	118	Double_t energy = ct->GetPHOSsignal() ;
	119	// check if CPV bit is set (see AliPHOSFastRecParticle)
	120	Double_t type[AliESDtrack::kSPECIES+4] ;
	121	ct->GetPHOSpid(type) ;
	122
	123	if ( (type[AliESDtrack::kElectron] == 1.0) \|\| (type[AliESDtrack::kPhoton] == 1.0) ) {
	124	if(!hCaloEnergyA) {
	125	out->cd() ;
	126	hCaloEnergyA = new TH1D("CaloEnergyA", "Energy in calorimeter Electron/Photon like", 500, 0., 50.) ;
	127	}
	128	hCaloEnergyA->Fill(energy) ;
	129	}
	130
	131	if ( type[AliESDtrack::kElectron] == 1.0 ) {
	132	if(!hCaloEnergyE) {
	133	out->cd() ;
	134	hCaloEnergyE = new TH1D("CaloEnergyE", "Energy in calorimeter Electron like", 500, 0., 50.) ;
	135	}
	136	hCaloEnergyE->Fill(energy) ;
	137	}
	138
	139	if ( type[AliESDtrack::kPhoton] == 1.0 ) {
	140	if(!hCaloEnergyG) {
	141	out->cd() ;
	142	hCaloEnergyG = new TH1D("CaloEnergyG", "Energy in calorimeter Gamma like", 500, 0., 50.) ;
	143	}
	144	hCaloEnergyG->Fill(energy) ;
	145	}
	146
	147	if ( type[AliESDtrack::kPi0] == 1.0 ) {
	148	if(!hCaloEnergyP) {
	149	out->cd() ;
	150	hCaloEnergyP = new TH1D("CaloEnergyP", "Energy in calorimeter Pi0 like", 500, 0., 100.) ;
	151	}
	152	hCaloEnergyP->Fill(energy) ;
	153	}
	154
	155	AliESDtrack * cp ;
	156	Int_t cpindex ;
	157	for (cpindex = 0 ; cpindex < esd->GetNumberOfTracks() ; cpindex++) {
	158	// get the charged tracks from central tracking
	159	cp = esd->GetTrack(cpindex) ;
	160	if ( cp->IsPHOS() )
	161	continue ;
	162	Match(ct, cp, dist) ;
	163	if (!heta && !hphi) {
	164	heta = new TH1D("Correta", "neutral-charged correlation in eta" , 100, 0., 360.) ;
	165	hphi = new TH1D("Corrphi", "neutral-charged correlation in phi" , 100, 0., 360.) ;
	166	}
	167	heta->Fill(dist[1]) ;
	168	hphi->Fill(dist[2]) ;
	169	}
	170	}
	171	}
	172	TBrowser * bs = new TBrowser("Root Memory Bowser", gROOT->FindObjectAny("AOD.root") ) ;
	173	bs->Show() ;
	174	out->Write() ;
	175	return kTRUE ;
	176	}
	177
	178	//__________________________________________________________________________
	179	void Match(AliESDtrack * ct, AliESDtrack * cp, Double_t * dist)
	180	{
	181	// Calculates the distance (x,z) between the particle detected by PHOS and
	182	// the charged particle reconstructed by the global tracking
	183
	184
	185	// Int_t phN ;
	186	// Double_t phZ, phX ;
	187
	188	// if (ct->IsPHOS()) { // it is a PHOS particle
	189	// Double_t cpTheta, cpPhi ;
	190	// Double_t phTheta, phPhi ;
	191	// cpTheta = cpPhi = phTheta = phPhi = 0. ;
	192	// // cout << "PHOS particle # " << " pos ("
	193	// // << pp->GetPos().X() << ", " << pp->GetPos().Y() << ", " << pp->GetPos().Z() << ")" << endl ;
	194
	195	// AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
	196	// gime->PHOSGeometry()->ImpactOnEmc(*pp, phN, phZ, phX) ;
	197	// Double_t xyzAtPHOS[3] ;
	198	// cp->GetOuterXYZ(xyzAtPHOS) ;
	199	// if ( (xyzAtPHOS[0] + xyzAtPHOS[1] + xyzAtPHOS[2]) != 0.) { //it has reached PHOS
	200	// //the next check are only if we want high quality tracks
	201	// // ULong_t status = cp->GetStatus() ;
	202	// // if ((status & AliESDtrack::kTRDput)==0)
	203	// // do not continue;
	204	// // if ((status & AliESDtrack::kTRDStop)!=0)
	205	// // do not continue;
	206	// // cout << "Charged particle # " << " pos ("
	207	// // << xyzAtPHOS[0] << ", " << xyzAtPHOS[1] << ", " << xyzAtPHOS[2] << ")" << endl ;
	208	// TVector3 poscp(xyzAtPHOS[0], xyzAtPHOS[1], xyzAtPHOS[2]) ;
	209	// Int_t cpN ;
	210	// Double_t cpZ,cpX ;
	211	// gime->PHOSGeometry()->ImpactOnEmc(poscp, cpN, cpZ, cpX) ;
	212	// if (cpN) {// we are inside the PHOS acceptance
	213	// // cout << "Charged Matching 1: " << cpN << " " << cpZ << " " << cpX << endl ;
	214	// // cout << "Charged Matching 2: " << phN << " " << phZ << " " << phX << endl ;
	215	// dist[0] = TMath::Sqrt( (cpZ-phZ)(cpZ-phZ) + (cpX-phX)(cpX-phX)) ;
	216	// }
	217	// phTheta = pp->Theta() ;
	218	// phPhi = pp->Phi() ;
	219	// TParticle tempo ;
	220	// tempo.SetMomentum(xyzAtPHOS[0], xyzAtPHOS[1], xyzAtPHOS[2], 0.) ;
	221	// cpTheta = tempo.Theta() ;
	222	// cpPhi = tempo.Phi() ;
	223	// //cout << phTheta << " " << phPhi << " " << endl
	224	// //cout << cpTheta << " " << cpPhi-phPhi << " " << endl ;
	225	// }
	226	// dist[1] = (phTheta - cpTheta)*TMath::RadToDeg() ;
	227	// dist[2] = (phPhi - cpPhi)*TMath::RadToDeg() ;
	228	// }
	229
	230	// if (ep) {
	231	// //cout << "EMCAL particle # " << endl ;
	232	// }
	233	}