- AliITSPlaneEffSPD.cxx .h

[u/mrichter/AliRoot.git] / PHOS / AliPHOSAnalyze.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$ */

//_________________________________________________________________________
// Algorythm class to analyze PHOS events. In this class we demostrate, 
// how to handle reconstructed objects with AliPHSOIndexToObject.
// As an example we propose sulotions for four most frequently used tasks:
//    DrawRecon(...) - to draw reconstructed objects in the PHOS plane,
//                     very usefull in the debuging
//    InvarianMass(...) - to calculate "REAL" and "MIXED" photon pairs 
//                        invariant mass distributions
//    EnergyResoluition(...) -\ Energy and position resolutions of the 
//    PositionResolution(...)-/ reconstructed photons
//    Contamination(...) - calculates contamination of the photon spectrum and 
//                         pobability of reconstruction of several primaries as 
//                         kGAMMA,kELECTRON etc.
////    User Case:
//    root [0] AliPHOSAnalyze * a = new AliPHOSAnalyze("galice.root")
//                    // set the file you want to analyse
//    root [1] a->DrawRecon(1,3)
//                    // plot RecObjects, made in event 1, PHOS module 3 
//    root [2] a->DrawRecon(1,3,"PHOSRP","another PID")
//                    // plot RecObjets made in the event 1, PHOS module 3,
//                    // produced in the another reconstruction pass,
//                    // which produced PHOS RecParticles ("PHOSRP") with 
//                    // title "another PID".
//    root [3] a->InvariantMass()
//                    // Calculates "REAL" and "MIXED" invariant mass 
//                    // distributions of kGAMMA and (kGAMMA+kNEUTRALEM)
//                    // and APPENDS this to the file "invmass.root"
//    root [4] a->PositionResolution()
//                    // calculates two dimentional histos: energy of the primary
//                    // photon vs distance betwin incedence point and reconstructed 
//                    // poisition. One can analyse the produced file position.root 
//                    // with macro PhotonPosition.C
//    root [5] a->EnergyResolution()
//                    // calculates two dimentional histos: energy of the primary
//                    // photon vs energy of the reconstructed particle. One can 
//                    // analyse the produced file energy.root 
//                    // with macro PhotonEnergy.C
//    root [6] a->Contamination()
//                    // fills spectra of primary photons and several kinds of 
//                    // reconstructed particles, so that analyzing them one can 
//                    // estimate conatmination, efficiency of registration etc.
//*--
//*-- Author: Dmitri Peressounko (SUBATECH & RRC Kurchatov Institute)
//////////////////////////////////////////////////////////////////////////////


// --- ROOT system ---

#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TH2.h"
#include "TParticle.h"
#include "TClonesArray.h"
#include "TMath.h"
#include "TROOT.h"

// --- Standard library ---

// --- AliRoot header files ---

#include "AliLog.h"
#include "AliStack.h"
#include "AliPHOSGeometry.h"
#include "AliPHOSAnalyze.h"
#include "AliPHOSDigit.h"
#include "AliPHOSSDigitizer.h"
#include "AliPHOSEmcRecPoint.h"
#include "AliPHOSCpvRecPoint.h"
#include "AliPHOSTrackSegment.h"
#include "AliPHOSRecParticle.h"
#include "AliPHOSLoader.h"


ClassImp(AliPHOSAnalyze)

//____________________________________________________________________________
AliPHOSAnalyze::AliPHOSAnalyze():
  fCorrection(1.2),  //Value calculated for default parameters of reconstruction
  fEvt(0),
  ffileName(),
  fRunLoader(0)
{
  // default ctor (useless)
}

//____________________________________________________________________________
AliPHOSAnalyze::AliPHOSAnalyze(Text_t * fileName):
  fCorrection(1.05),  //Value calculated for default parameters of reconstruction   
  fEvt(0),
  ffileName(fileName),
  fRunLoader(0)
{
  // ctor: analyze events from root file "name"
  fRunLoader = AliRunLoader::Open(fileName,"AliPHOSAnalyze");
  if (fRunLoader == 0x0)
   {
     AliError(Form("Error Loading session"));
   }
}

//____________________________________________________________________________
AliPHOSAnalyze::AliPHOSAnalyze(const AliPHOSAnalyze & ana): 
  TObject(ana),
  fCorrection(0.),
  fEvt(0),
  ffileName(),
  fRunLoader(0)
{
  // copy ctor
  ( (AliPHOSAnalyze &)ana ).Copy(*this) ;
}

//____________________________________________________________________________
AliPHOSAnalyze::~AliPHOSAnalyze()
{
  // dtor

}
//____________________________________________________________________________
void AliPHOSAnalyze::DrawRecon(Int_t Nevent,Int_t Nmod){
  //Draws pimary particles and reconstructed 
  //digits, RecPoints, RecPartices etc 
  //for event Nevent in the module Nmod.

  //========== Create ObjectLoader
  if (fRunLoader == 0x0)
   {
     AliError(Form("Error Loading session"));
     return;
   }
  
  AliPHOSLoader* gime = dynamic_cast<AliPHOSLoader*>(fRunLoader->GetLoader("PHOSLoader"));
  if ( gime == 0 ) 
   {
     AliError(Form("Could not obtain the Loader object !")); 
     return ;
   } 
  
  
  if(Nevent >= fRunLoader->GetNumberOfEvents() ) {
    AliError(Form("There is no event %d only %d events available", Nevent, fRunLoader->GetNumberOfEvents() )) ;
    return ;
  }
  const AliPHOSGeometry * phosgeom = gime->PHOSGeometry() ; 
  fRunLoader->GetEvent(Nevent);

  Int_t nx = phosgeom->GetNPhi() ;
  Int_t nz = phosgeom->GetNZ() ;
  Float_t * cri= phosgeom->GetEMCAGeometry()->GetCrystalHalfSize() ;
  Float_t x = nx*cri[0] ;
  Float_t z = nz*cri[2] ;
  Int_t nxCPV = (Int_t) (nx*phosgeom->GetPadSizePhi()/(2.*cri[0])) ;
  Int_t nzCPV = (Int_t) (nz*phosgeom->GetPadSizeZ()/(2.*cri[2])) ;
  
  TH2F * emcDigits = (TH2F*) gROOT->FindObject("emcDigits") ;
  if(emcDigits)
    emcDigits->Delete() ;
  emcDigits = new TH2F("emcDigits","EMC digits",  nx,-x,x,nz,-z,z);
  TH2F * emcSdigits =(TH2F*) gROOT->FindObject("emcSdigits") ;  
  if(emcSdigits)
    emcSdigits->Delete() ;
  emcSdigits = new TH2F("emcSdigits","EMC sdigits", nx,-x,x,nz,-z,z);
  TH2F * emcRecPoints = (TH2F*)gROOT->FindObject("emcRecPoints") ; 
  if(emcRecPoints)
    emcRecPoints->Delete() ;
  emcRecPoints = new TH2F("emcRecPoints","EMC RecPoints",nx,-x,x,nz,-z,z);
  TH2F * cpvSdigits =(TH2F*) gROOT->FindObject("cpvSdigits") ;
  if(cpvSdigits)
    cpvSdigits->Delete() ;
  cpvSdigits = new TH2F("cpvSdigits","CPV sdigits", nx,-x,x,nz,-z,z);
  TH2F * cpvDigits = (TH2F*)gROOT->FindObject("cpvDigits") ;
  if(cpvDigits)
    cpvDigits->Delete() ;
  cpvDigits = new TH2F("cpvDigits","CPV digits",   nxCPV,-x,x,nzCPV,-z,z) ;
  TH2F * cpvRecPoints= (TH2F*)gROOT->FindObject("cpvRecPoints") ; 
  if(cpvRecPoints)
    cpvRecPoints->Delete() ;
  cpvRecPoints = new TH2F("cpvRecPoints","CPV RecPoints",    nxCPV,-x,x,nzCPV,-z,z) ;

  TH2F * phot = (TH2F*)gROOT->FindObject("phot") ;
  if(phot)
    phot->Delete() ;
  phot = new TH2F("phot","Primary Photon",  nx,-x,x,nz,-z,z);
  TH2F * recPhot = (TH2F*)gROOT->FindObject("recPhot") ; 
  if(recPhot)
    recPhot->Delete() ;
  recPhot = new TH2F("recPhot","RecParticles with primary Photon",nx,-x,x,nz,-z,z);
  
  //Get Vertex
  Double_t vtx[3]={0.,0.,0.} ;  
//DP: extract vertex either from Generator or from data

  
  //Plot Primary Particles
  
  if (fRunLoader->Stack() == 0x0) fRunLoader->LoadKinematics("READ");
  

  const TParticle * primary ;
  Int_t iPrimary ;
  for ( iPrimary = 0 ; iPrimary < fRunLoader->Stack()->GetNprimary() ; iPrimary++)
    {
      primary = fRunLoader->Stack()->Particle(iPrimary);
      
      Int_t primaryType = primary->GetPdgCode();
//       if( (primaryType == 211)||(primaryType == -211)||(primaryType == 2212)||(primaryType == -2212)
//          ||(primaryType == 11)||(primaryType == -11) ) {
//         Int_t moduleNumber ;
//         Double_t primX, primZ ;
//         phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
//         if(moduleNumber==Nmod)
//           charg->Fill(primZ,primX,primary->Energy()) ;
//       }
      if( primaryType == 22 ) {
        Int_t moduleNumber ;
        Double_t primX, primZ ;
        phosgeom->ImpactOnEmc(vtx,primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
        if(moduleNumber==Nmod) 
          phot->Fill(primZ,primX,primary->Energy()) ;
      }
//       else{
//         if( primaryType == -2112 ) {
//           Int_t moduleNumber ;
//           Double_t primX, primZ ;
//           phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
//           if(moduleNumber==Nmod)
//             nbar->Fill(primZ,primX,primary->Energy()) ;
//         }
//       }
    }  

  
  Int_t iSDigit ;
  AliPHOSDigit * sdigit ;
  const TClonesArray * sdigits = gime->SDigits() ;
  Int_t nsdig[5] = {0,0,0,0,0} ;
  if(sdigits){
    for(iSDigit = 0; iSDigit < sdigits->GetEntriesFast() ; iSDigit++)
      {
       sdigit = (AliPHOSDigit *) sdigits->At(iSDigit) ;
       Int_t relid[4];
       phosgeom->AbsToRelNumbering(sdigit->GetId(), relid) ;
       Float_t x,z ;
       phosgeom->RelPosInModule(relid,x,z);
       AliPHOSSDigitizer* sd = dynamic_cast<AliPHOSSDigitizer*>(gime->SDigitizer());
       Float_t e = sd->Calibrate(sdigit->GetAmp()) ;
       nsdig[relid[0]-1]++ ;
       if(relid[0]==Nmod){
         if(relid[1]==0)  //EMC
           emcSdigits->Fill(x,z,e) ;
         if( relid[1]!=0 )
           cpvSdigits->Fill(x,z,e) ;
       }
      }
  }
  TString message ; 
  message  = "Number of EMC + CPV SDigits per module: \n" ;
  message += "%d %d %d %d %d\n"; 
  AliInfo(Form(message.Data(), nsdig[0], nsdig[1], nsdig[2], nsdig[3], nsdig[4] )) ;

  //Plot digits
  Int_t iDigit ;
  AliPHOSDigit * digit ;
  const TClonesArray * digits = gime->Digits(); 
  if(digits) {
    for(iDigit = 0; iDigit < digits->GetEntriesFast(); iDigit++)
      {
       digit = (AliPHOSDigit *) digits->At(iDigit) ;
       Int_t relid[4];
       phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
       Float_t x,z ;
       phosgeom->RelPosInModule(relid,x,z) ;
       AliPHOSSDigitizer* sd = dynamic_cast<AliPHOSSDigitizer*>(gime->SDigitizer());
       Float_t e = sd->Calibrate(digit->GetAmp()) ;
       if(relid[0]==Nmod){
         if(relid[1]==0)  //EMC
           emcDigits->Fill(x,z,e) ;
         if( relid[1]!=0 )
           cpvDigits->Fill(x,z,e) ;
       }
      }
  }
  
  
  //Plot RecPoints
  Int_t irecp ;
  TVector3 pos ;
  TObjArray * emcrp = gime->EmcRecPoints() ;
  if(emcrp) {
    for(irecp = 0; irecp < emcrp->GetEntriesFast() ; irecp ++){
      AliPHOSEmcRecPoint * emc = (AliPHOSEmcRecPoint *) emcrp->At(irecp) ;
      if(emc->GetPHOSMod()==Nmod){
       emc->GetLocalPosition(pos) ;
       emcRecPoints->Fill(pos.X(),pos.Z(),emc->GetEnergy());
      }
    }
  }
  
  TObjArray * cpvrp = gime->CpvRecPoints() ;
  if(cpvrp) {
    for(irecp = 0; irecp < cpvrp->GetEntriesFast() ; irecp ++){
      AliPHOSRecPoint * cpv = (AliPHOSCpvRecPoint *) cpvrp->At(irecp) ;
      if(cpv->GetPHOSMod()==Nmod){
       cpv->GetLocalPosition(pos) ;
       cpvRecPoints->Fill(pos.X(),pos.Z(),cpv->GetEnergy());
      }
    }
  }
    
  //Plot RecParticles
  AliPHOSRecParticle * recParticle ;
  Int_t iRecParticle ;
  TClonesArray * rp = gime->RecParticles() ;
  TClonesArray * ts = gime->TrackSegments() ;
  if(rp && ts && emcrp) {
    for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast() ; iRecParticle++ )
      {
       recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
       Int_t moduleNumberRec ;
       Double_t recX, recZ ;
       phosgeom->ImpactOnEmc(vtx,recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
       if(moduleNumberRec == Nmod){
         
         Double_t minDistance = 5. ;
         Int_t closestPrimary = -1 ;       

         //extract list of primaries: it is stored at EMC RecPoints
         Int_t emcIndex = ((AliPHOSTrackSegment *) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
         Int_t numberofprimaries ;
         Int_t * listofprimaries  = ((AliPHOSRecPoint*) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries)  ;
         Int_t index ;
         const TParticle * primary ;
         Double_t distance = minDistance ;
         
         for ( index = 0 ; index < numberofprimaries ; index++){
           primary = fRunLoader->Stack()->Particle(listofprimaries[index]) ;
           Int_t moduleNumber ;
           Double_t primX, primZ ;
           phosgeom->ImpactOnEmc(vtx,primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
           if(moduleNumberRec == moduleNumber)
             distance = TMath::Sqrt((recX-primX)*(recX-primX)+(recZ-primZ)*(recZ-primZ) ) ;
           if(minDistance > distance)
             {
              minDistance = distance ;
              closestPrimary = listofprimaries[index] ;
             }
         }
         
         if(closestPrimary >=0 ){
           
           Int_t primaryType = fRunLoader->Stack()->Particle(closestPrimary)->GetPdgCode() ;
           
           if(primaryType==22)
             recPhot->Fill(recZ,recX,recParticle->Energy()) ;
//            else
//              if(primaryType==-2112)
//               recNbar->Fill(recZ,recX,recParticle->Energy()) ; 
         }
       }
      }

  }
  
  //Plot made histograms
  emcSdigits->Draw("box") ;
  emcDigits->SetLineColor(5) ;
  emcDigits->Draw("boxsame") ;
  emcRecPoints->SetLineColor(2) ;
  emcRecPoints->Draw("boxsame") ;
  cpvSdigits->SetLineColor(1) ;
  cpvSdigits->Draw("boxsame") ;
  
}
//____________________________________________________________________________
void AliPHOSAnalyze::Ls(){
  //lists branches and titles of PHOS-related branches of TreeR, TreeD, TreeS
  
  if (fRunLoader == 0x0)
   {
     AliError(Form("Error Loading session"));
     return;
   }
  
  AliPHOSLoader* gime = dynamic_cast<AliPHOSLoader*>(fRunLoader->GetLoader("PHOSLoader"));
  if ( gime == 0 ) 
   {
     AliError(Form("Could not obtain the Loader object !")); 
     return ;
   } 


  Int_t ibranch;
  TObjArray * branches; 
  
  if (gime->TreeS() == 0x0) 
   {
     if (gime->LoadSDigits("READ"))
      {
        AliError(Form("Problems with loading summable digits"));
        return;
      }
   }
  branches = gime->TreeS()->GetListOfBranches() ;
 
  TString message ; 
  message  = "TreeS:\n" ;
  for(ibranch = 0;ibranch <branches->GetEntries();ibranch++){
    TBranch * branch=(TBranch *) branches->At(ibranch) ;
    if(strstr(branch->GetName(),"PHOS") ){
      message += "       " ; 
      message += branch->GetName() ; 
      message += "     " ; 
      message += branch->GetTitle() ;
      message += "\n" ; 
    }
  }
  if (gime->TreeD() == 0x0) 
   {
     if (gime->LoadDigits("READ"))
      {
        AliError(Form("Problems with loading digits"));
        return;
      }
   }

  branches = gime->TreeD()->GetListOfBranches() ;
  
  message += "TreeD:\n" ;
  for(ibranch = 0;ibranch <branches->GetEntries();ibranch++){
    TBranch * branch=(TBranch *) branches->At(ibranch) ;
    if(strstr(branch->GetName(),"PHOS") ) {
      message += "       "; 
      message += branch->GetName() ; 
      message += "     " ; 
      message += branch->GetTitle() ; 
      message +="\n" ;
    }
  }
  
  if (gime->TreeR() == 0x0) 
   {
     if (gime->LoadRecPoints("READ"))
      {
        AliError(Form("Problems with loading rec points"));
        return;
      }
   }

  branches = gime->TreeR()->GetListOfBranches() ;
  
  message += "TreeR: \n" ;
  for(ibranch = 0;ibranch <branches->GetEntries();ibranch++){
    TBranch * branch=(TBranch *) branches->At(ibranch) ;
    if(strstr(branch->GetName(),"PHOS") ) {
      message += "       " ; 
      message += branch->GetName() ; 
      message += "     " ;
      message += branch->GetTitle() ; 
      message += "\n" ;
    }
  }
  AliInfo(Form(message.Data())) ;  
}
//____________________________________________________________________________
 void AliPHOSAnalyze::InvariantMass()
{
  // Calculates Real and Mixed invariant mass distributions
  if (fRunLoader == 0x0)
   {
     AliError(Form("Error Loading session"));
     return;
   }
  
  AliPHOSLoader* gime = dynamic_cast<AliPHOSLoader*>(fRunLoader->GetLoader("PHOSLoader"));
  if ( gime == 0 ) 
   {
     AliError(Form("Could not obtain the Loader object !")); 
     return ;
   } 
  
  gime->LoadRecParticles("READ");
  
  Int_t nMixedEvents = 4 ; //# of events used for calculation of 'mixed' distribution 
  
  
  //opening file
  TFile * mfile = new TFile("invmass.root","update");
  
  //========== Reading /Booking Histograms
  TH2D * hRealEM = 0 ;
  hRealEM = (TH2D*) mfile->Get("hRealEM") ;
  if(hRealEM == 0) 
    hRealEM = new TH2D("hRealEM",   "Real for EM particles",      250,0.,1.,40,0.,4.) ;
  TH2D * hRealPhot = 0 ;

  hRealPhot = (TH2D*)mfile->Get("hRealPhot");
  if(hRealPhot == 0)
    hRealPhot = new TH2D("hRealPhot", "Real for kPhoton particles", 250,0.,1.,40,0.,4.) ;

  TH2D * hMixedEM = 0 ;
  hMixedEM = (TH2D*) mfile->Get("hMixedEM") ;
  if(hMixedEM == 0)
    hMixedEM = new TH2D("hMixedEM",  "Mixed for EM particles",     250,0.,1.,40,0.,4.) ;

  TH2D * hMixedPhot = 0 ;
  hMixedPhot = (TH2D*) mfile->Get("hMixedPhot") ;
  if(hMixedPhot == 0)
    hMixedPhot = new TH2D("hMixedPhot","Mixed for kPhoton particles",250,0.,1.,40,0.,4.) ;
  

  //reading event and copyng it to TConesArray of all photons

  TClonesArray * allRecParticleList  = new TClonesArray("AliPHOSRecParticle", 1000) ;
  Int_t * nRecParticles = new Int_t[nMixedEvents] ; // to mark boundaries of each event in the total list
  for(Int_t index = 0; index < nMixedEvents; index ++)
    nRecParticles[index] = 0 ;
  Int_t iRecPhot = 0 ;                // number of EM particles in total list
  
  //scan over all events
  Int_t event ;
  

  if (fRunLoader->TreeE() == 0x0) fRunLoader->LoadHeader();
  
  
  Int_t maxevent = (Int_t)fRunLoader->TreeE()->GetEntries();
  //  for(event = 0; event < gime->MaxEvent(); event++  ){
  
  
  
  for(event = 0; event < maxevent; event++  ){
    fRunLoader->GetEvent(event);  //will read only TreeR 
    
    //copy EM RecParticles to the "total" list        
    const AliPHOSRecParticle * recParticle ;
    Int_t iRecParticle ;
    TClonesArray * rp = gime->RecParticles() ;
    if(!rp){
      AliError(Form("Can't find RecParticles")) ; 
      return ;
    }

    for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast(); iRecParticle++ )
      {
       recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
       if((recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)||
          (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW))
         new( (*allRecParticleList)[iRecPhot++] ) AliPHOSRecParticle(*recParticle) ;
      }
    
    Int_t mevent = event%nMixedEvents ; //event number in the "mixed" cicle
    nRecParticles[mevent] = iRecPhot-1 ;  
    
    //check, if it is time to calculate invariant mass?
    Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ; 
    if((mevent == 0) && (event +1 == maxevent)){
      
      //   if((mevent == 0) && (event +1 == gime->MaxEvent())){
      
      //calculate invariant mass:
      Int_t irp1,irp2 ;
      Int_t nCurEvent = 0 ;
      
      for(irp1 = 0; irp1 < allRecParticleList->GetEntries()-1; irp1++){
       AliPHOSRecParticle * rp1 = (AliPHOSRecParticle *)allRecParticleList->At(irp1) ;
       
       for(irp2 = irp1+1; irp2 < allRecParticleList->GetEntries(); irp2++){
         AliPHOSRecParticle * rp2 = (AliPHOSRecParticle *)allRecParticleList->At(irp2) ;
         
         Double_t invMass ;
         invMass = (rp1->Energy()+rp2->Energy())*(rp1->Energy()+rp2->Energy())-
           (rp1->Px()+rp2->Px())*(rp1->Px()+rp2->Px())-
           (rp1->Py()+rp2->Py())*(rp1->Py()+rp2->Py())-
           (rp1->Pz()+rp2->Pz())*(rp1->Pz()+rp2->Pz()) ;
         
         if(invMass> 0)
           invMass = TMath::Sqrt(invMass);
         
         Double_t pt ; 
         pt = TMath::Sqrt((rp1->Px()+rp2->Px() )*( rp1->Px()+rp2->Px() ) + 
                        (rp1->Py()+rp2->Py() )*( rp1->Py()+rp2->Py() ) );
         
         if(irp1 > nRecParticles[nCurEvent])
           nCurEvent++;
         
         if(irp2 <= nRecParticles[nCurEvent]){ //'Real' event
           hRealEM->Fill(invMass,pt);
           if((rp1->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)&&
              (rp2->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST) )
             hRealPhot->Fill(invMass,pt);
         }
         else{
           hMixedEM->Fill(invMass,pt);
           if((rp1->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)&&
              (rp2->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST) )
             hMixedPhot->Fill(invMass,pt);
         } //real-mixed
         
      } //loop over second rp
      }//loop over first rp
      
      //Make some cleanings 
      for(Int_t index = 0; index < nMixedEvents; index ++)
       nRecParticles[index] = 0 ;
      iRecPhot = 0 ;              
      allRecParticleList->Clear() ;
      
    }
  }
  delete allRecParticleList ;
  
  //writing output
  mfile->cd();
  
  hRealEM->Write(0,kOverwrite) ;
  hRealPhot->Write(0,kOverwrite) ;
  hMixedEM->Write(0,kOverwrite) ;
  hMixedPhot->Write(0,kOverwrite) ;
  
  mfile->Write();
  mfile->Close();
  delete mfile ;
  delete nRecParticles;

}

//____________________________________________________________________________
 void AliPHOSAnalyze::EnergyResolution()
{
  //fills two dimentional histo: energy of primary vs. energy of reconstructed

  TH2F * hAllEnergy = 0 ;  //all reconstructed with primary photon   
  TH2F * hPhotEnergy= 0 ;  //kGamma  with primary photon   
  TH2F * hEMEnergy  = 0 ;  //electromagnetic with primary photon   

  //opening file and reading histograms if any
  TFile * efile = new TFile("energy.root","update");

  hAllEnergy = (TH2F*)efile->Get("hAllEnergy") ;
  if(hAllEnergy == 0)
    hAllEnergy  = new TH2F("hAllEnergy",  "Energy of any RP with primary photon",100, 0., 5., 100, 0., 5.);

  hPhotEnergy =(TH2F*) efile->Get("hPhotEnergy") ;
  if(hPhotEnergy == 0)
    hPhotEnergy = new TH2F("hPhotEnergy", "Energy of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.);

  hEMEnergy =(TH2F*) efile->Get("hEMEnergy");
  if(hEMEnergy == 0)
    hEMEnergy   = new TH2F("hEMEnergy",   "Energy of EM with primary photon",    100, 0., 5., 100, 0., 5.);


  if (fRunLoader == 0x0)
   {
     AliError(Form("Error Loading session"));
     return;
   }
  
  AliPHOSLoader* gime = dynamic_cast<AliPHOSLoader*>(fRunLoader->GetLoader("PHOSLoader"));
  if ( gime == 0 ) 
   {
     AliError(Form("Could not obtain the Loader object !")); 
     return ;
   } 


  const AliPHOSGeometry * phosgeom = gime->PHOSGeometry();

  Int_t ievent;
  Int_t maxevent = (Int_t)fRunLoader->TreeE()->GetEntries();

  fRunLoader->LoadKinematics("READ");
  gime->LoadTracks("READ");
  
  for ( ievent=0; ievent < maxevent ; ievent++){

    //read the current event
    fRunLoader->GetEvent(ievent) ;

    Double_t vtx[3]={0.,0.,0.} ;  

    const AliPHOSRecParticle * recParticle ;
    Int_t iRecParticle ;
    TClonesArray * rp = gime->RecParticles() ;
    if(!rp) {
      AliError(Form("Event %d,  Can't find RecParticles ", ievent)) ;  
      return ;
    }
    TClonesArray * ts = gime->TrackSegments() ;
    if(!ts) {
      AliError(Form("Event %d,  Can't find TrackSegments", ievent)) ;  
      return ;
    }
    TObjArray * emcrp = gime->EmcRecPoints() ;
    if(!emcrp){
      AliError(Form("Event %d,  Can't find EmcRecPoints")) ; 
      return ;
    }
      
    for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast() ;iRecParticle++ ){
      recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
      
      //find the closest primary
      Int_t moduleNumberRec ;
      Double_t recX, recZ ;
      phosgeom->ImpactOnEmc(vtx,recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
      
      Double_t minDistance  = 100. ;
      Int_t closestPrimary = -1 ;
      
      //extract list of primaries: it is stored at EMC RecPoints
      Int_t emcIndex = ((AliPHOSTrackSegment*) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
      Int_t numberofprimaries ;
      Int_t * listofprimaries  = ((AliPHOSEmcRecPoint*) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries)  ;
      
      Int_t index ;
      const TParticle * primary ;
      Double_t distance = minDistance ;
      Double_t dX, dZ; 
      Double_t dXmin = 0.; 
      Double_t dZmin = 0. ;
      for ( index = 0 ; index < numberofprimaries ; index++){

       primary = fRunLoader->Stack()->Particle(listofprimaries[index]) ;

       Int_t moduleNumber ;
       Double_t primX, primZ ;
       phosgeom->ImpactOnEmc(vtx,primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
       if(moduleNumberRec == moduleNumber) {
         dX = recX - primX;
         dZ = recZ - primZ;
         distance = TMath::Sqrt(dX*dX + dZ*dZ) ;
         if(minDistance > distance) {
           minDistance = distance ;
           dXmin = dX;
           dZmin = dZ;
           closestPrimary = listofprimaries[index] ;
         }
       }
      }

      //if found primary, fill histograms
      if(closestPrimary >=0 ){
       const TParticle * primary = fRunLoader->Stack()->Particle(closestPrimary) ;
       if(primary->GetPdgCode() == 22){
         hAllEnergy->Fill(primary->Energy(), recParticle->Energy()) ;
         if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST){
           hPhotEnergy->Fill(primary->Energy(), recParticle->Energy() ) ; 
           hEMEnergy->Fill(primary->Energy(), recParticle->Energy() ) ; 
         }
         else
           if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW)
             hEMEnergy->Fill(primary->Energy(), recParticle->Energy() ) ; 
       }
      }
    }
  }

  //write filled histograms
  efile->cd() ;
  hAllEnergy->Write(0,kOverwrite) ;
  hPhotEnergy->Write(0,kOverwrite) ;
  hEMEnergy->Write(0,kOverwrite)  ;
  //  efile->Write() ;
  efile->Close() ;
  delete efile ;

}
//____________________________________________________________________________
void AliPHOSAnalyze::PositionResolution()
{
  //fills two dimentional histo: energy vs. primary - reconstructed distance  



  TH2F * hAllPosition  = 0;     // Position of any RP with primary photon
  TH2F * hPhotPosition = 0;    // Position of kGAMMA with primary photon
  TH2F * hEMPosition   = 0;      // Position of EM with primary photon

  TH1F * hAllPositionX = 0;    // X-Position Resolution of photons with photon primary
  TH1F * hAllPositionZ = 0;    // Z-Position Resolution of photons with photon primary


  //opening file and reading histograms if any
  TFile * pfile = new TFile("position.root","update");

  hAllPosition = (TH2F*)pfile->Get("hAllPosition");
  if(hAllPosition == 0)
    hAllPosition  = new TH2F("hAllPosition",  
                          "Position of any RP with primary photon",100, 0., 5., 100, 0., 5.);
  hPhotPosition= (TH2F*)pfile->Get("hPhotPosition");
  if(hPhotPosition == 0)
    hPhotPosition = new TH2F("hPhotPosition", 
                          "Position of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.);
  hEMPosition= (TH2F*)pfile->Get("hEMPosition") ;
  if(hEMPosition == 0)
    hEMPosition   = new TH2F("hEMPosition",   
                          "Position of EM with primary photon",    100, 0., 5., 100, 0., 5.);
  hAllPositionX = (TH1F*)pfile->Get("hAllPositionX") ;                        
  if(hAllPositionX == 0)
    hAllPositionX = new TH1F("hAllPositionX", 
                          "Delta X of any RP with primary photon",100, -2., 2.);
  hAllPositionZ =(TH1F*) pfile->Get("hAllPositionZ") ;
  if(hAllPositionZ == 0)
    hAllPositionZ = new TH1F("hAllPositionZ", 
                          "Delta X of any RP with primary photon",100, -2., 2.);

  if (fRunLoader == 0x0)
   {
     AliError(Form("Error Loading session"));
     return;
   }
  
  AliPHOSLoader* gime = dynamic_cast<AliPHOSLoader*>(fRunLoader->GetLoader("PHOSLoader"));
  if ( gime == 0 ) 
   {
     AliError(Form("Could not obtain the Loader object !")); 
     return ;
   } 
  
  if (fRunLoader->TreeE() == 0x0) fRunLoader->LoadHeader();

  const AliPHOSGeometry * phosgeom = gime->PHOSGeometry() ; 

  Int_t ievent;
  Int_t maxevent = (Int_t)fRunLoader->TreeE()->GetEntries() ; 
  for ( ievent=0; ievent < maxevent ; ievent++){
    
    //read the current event
    fRunLoader->GetEvent(ievent) ;

    //DP:Extract vertex position
    Double_t vtx[3]={0.,0.,0.} ;  

    TClonesArray * rp = gime->RecParticles() ;
    if(!rp) {
      AliError(Form("Event %d,  Can't find RecParticles", ievent)) ;
      return ;
    }
    TClonesArray * ts = gime->TrackSegments() ;
    if(!ts) {
      AliError(Form("Event %d,  Can't find TrackSegments", ievent)) ;
      return ;
    }
    TObjArray * emcrp = gime->EmcRecPoints() ;
    if(!emcrp){
      AliError(Form("Event %d,  Can't find EmcRecPoints", ievent)) ;
      return ;
    }
    
 
    const AliPHOSRecParticle * recParticle ;
    Int_t iRecParticle ;
    for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast(); iRecParticle++ ){
      recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
      
      //find the closest primary
      Int_t moduleNumberRec ;
      Double_t recX, recZ ;
      phosgeom->ImpactOnEmc(vtx,recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
      
      Double_t minDistance  = 100. ;
      Int_t closestPrimary = -1 ;
      
      //extract list of primaries: it is stored at EMC RecPoints
      Int_t emcIndex = ((AliPHOSTrackSegment*) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
      Int_t numberofprimaries ;
      Int_t * listofprimaries  = ((AliPHOSEmcRecPoint *) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries)  ;

      Int_t index ;
      const TParticle * primary ;
      Double_t distance = minDistance ;
      Double_t dX = 1000; // incredible number
      Double_t dZ = 1000; // for the case if no primary will be found
      Double_t dXmin = 0.; 
      Double_t dZmin = 0. ;
      for ( index = 0 ; index < numberofprimaries ; index++){
       primary = fRunLoader->Stack()->Particle(listofprimaries[index]) ;
       Int_t moduleNumber ;
       Double_t primX, primZ ;
       phosgeom->ImpactOnEmc(vtx,primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
       if(moduleNumberRec == moduleNumber) {
         dX = recX - primX;
         dZ = recZ - primZ;
         distance = TMath::Sqrt(dX*dX + dZ*dZ) ;
         if(minDistance > distance) {
           minDistance = distance ;
           dXmin = dX;
           dZmin = dZ;
           closestPrimary = listofprimaries[index] ;
         }
       }
      }
      
      //if found primary, fill histograms
      if(closestPrimary >=0 ){
       const TParticle * primary = fRunLoader->Stack()->Particle(closestPrimary) ;
       if(primary->GetPdgCode() == 22){
         hAllPosition->Fill(primary->Energy(), minDistance) ;
         hAllPositionX->Fill(primary->Energy(), dX) ;
         hAllPositionZ->Fill(primary->Energy(), dZ) ;
         if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST){
           hPhotPosition->Fill(primary->Energy(), minDistance ) ; 
           hEMPosition->Fill(primary->Energy(), minDistance ) ; 
         }
         else
           if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW)
             hEMPosition->Fill(primary->Energy(), minDistance ) ; 
       }
      }
    }
  }
  
  //Write output histgrams
  pfile->cd() ;
  hAllPosition->Write(0,kOverwrite) ;
  hAllPositionX->Write(0,kOverwrite) ;
  hAllPositionZ->Write(0,kOverwrite) ;
  hPhotPosition->Write(0,kOverwrite) ;
  hEMPosition->Write(0,kOverwrite) ;
  pfile->Write() ;
  pfile->Close() ;
  delete pfile ;


}
//____________________________________________________________________________
void AliPHOSAnalyze::Contamination(){
// fills spectra of primary photons and several kinds of 
// reconstructed particles, so that analyzing them one can 
// estimate conatmination, efficiency of registration etc.

  //define several general histograms
  TH1F * hPrimary = 0;   //spectrum (P_t distribution) of primary photons         
  TH1F * hAllRP   = 0;   //spectrum of all RecParticles in PHOS
  TH1F * hPhot    = 0;   //spectrum of kGAMMA RecParticles
  TH1F * hShape   = 0;   //spectrum of all EM RecParticles
  TH1F * hVeto    = 0;   //spectrum of all neutral RecParticles

  //Now separate histograms in accoradance with primary
  //primary - photon
  TH1F * hPhotReg = 0;   //Registeres as photon
  TH1F * hPhotEM  = 0;   //Registered as EM       

  //primary - n
  TH1F * hNReg = 0;   //Registeres as photon          
  TH1F * hNEM  = 0;   //Registered as EM            

  //primary - nBar
  TH1F * hNBarReg = 0;   //Registeres as photon
  TH1F * hNBarEM  = 0;   //Registered as EM          

  //primary - charged hadron (pBar excluded)
  TH1F * hChargedReg = 0;   //Registeres as photon  
  TH1F * hChargedEM  = 0;   //Registered as EM           

  //primary - pBar
  TH1F * hPbarReg = 0;   //Registeres as photon  
  TH1F * hPbarEM  = 0;   //Registered as EM 


  //Reading histograms from the file
  TFile * cfile = new TFile("contamination.root","update") ;

  //read general histograms
  hPrimary = (TH1F*) cfile->Get("hPrimary") ;
  if(hPrimary == 0)
    hPrimary= new TH1F("hPrimary", "Primary photon spectrum", 100, 0., 5.);
  hAllRP = (TH1F*)cfile->Get("hAllRP") ;
  if(hAllRP == 0)
    hAllRP = new TH1F("hAllRP","All Reconstructed particles", 100, 0., 5.);
  hPhot  = (TH1F*)cfile->Get("hPhot") ;
  if(hPhot == 0)
    hPhot = new TH1F("hPhot","All kGAMMA RecParticles",100, 0., 5.);
  hShape = (TH1F*) cfile->Get("hShape") ;
  if(hShape == 0)
    hShape = new TH1F("hShape","All particles with EM shower",100, 0., 5.);
  hVeto= (TH1F*)cfile->Get("hVeto") ;
  if(hVeto == 0) 
    hVeto  = new TH1F("hVeto", "All uncharged particles",      100, 0., 5.);


  //primary - photon
  hPhotReg = (TH1F*)cfile->Get("hPhotReg");
  if(hPhotReg == 0)
    hPhotReg   = new TH1F("hPhotReg","Photon registered as photon",100, 0., 5.);
  hPhotEM  =(TH1F*)cfile->Get("hPhotEM");
  if(hPhotEM== 0)
    hPhotEM   = new TH1F("hPhotEM",  "Photon registered as EM", 100, 0., 5.);

  //primary - n
  hNReg = (TH1F*)cfile->Get("hNReg");
  if(hNReg== 0)
   hNReg      = new TH1F("hNReg",   "N registered as photon",              100, 0., 5.);
  hNEM  = (TH1F*)cfile->Get("hNEM"); 
  if(hNEM== 0)
    hNEM      = new TH1F("hNEM",    "N registered as EM",      100, 0., 5.);

  //primary - nBar
  hNBarReg =(TH1F*)cfile->Get("hNBarReg");
  if(hNBarReg== 0)
   hNBarReg   = new TH1F("hNBarReg", "NBar registered as photon",           100, 0., 5.);
  hNBarEM  =(TH1F*)cfile->Get("hNBarEM"); 
  if(hNBarEM== 0)
    hNBarEM   = new TH1F("hNBarEM",  "NBar registered as EM",   100, 0., 5.);

  //primary - charged hadron (pBar excluded)
  hChargedReg = (TH1F*)cfile->Get("hChargedReg");
  if(hChargedReg== 0)
    hChargedReg= new TH1F("hChargedReg", "Charged hadron registered as photon",100, 0., 5.);
  hChargedEM  = (TH1F*)cfile->Get("hChargedEM"); 
  if(hChargedEM== 0)
    hChargedEM= new TH1F("hChargedEM","Charged registered as EM",100, 0., 5.);
 
  //primary - pBar
  hPbarReg = (TH1F*)cfile->Get("hPbarReg");
  if(hPbarReg== 0)
    hPbarReg= new TH1F("hPbarReg", "pBar registered as photon",100, 0., 5.);
  hPbarEM  = (TH1F*)cfile->Get("hPbarEM");
  if(hPbarEM== 0)
    hPbarEM= new TH1F("hPbarEM","Pbar registered as EM",100, 0., 5.);
  

  //Now make some initializations

  Int_t counter[8][5] ;      //# of registered particles 
  Int_t i1,i2 ;
  for(i1 = 0; i1<8; i1++)
    for(i2 = 0; i2<5; i2++)
      counter[i1][i2] = 0 ;



  if (fRunLoader == 0x0)
   {
     AliError(Form("Error Loading session"));
     return;
   }
  
  AliPHOSLoader* gime = dynamic_cast<AliPHOSLoader*>(fRunLoader->GetLoader("PHOSLoader"));
  if ( gime == 0 ) 
   {
     AliError(Form("Could not obtain the Loader object !")); 
     return ;
   } 
  
  if (fRunLoader->TreeE() == 0x0) fRunLoader->LoadHeader();
  const AliPHOSGeometry * phosgeom = gime->PHOSGeometry() ; 
  
  Int_t ievent;
  Int_t maxevent = (Int_t)fRunLoader->TreeE()->GetEntries() ; 
  for ( ievent=0; ievent < maxevent ; ievent++){
    
    fRunLoader->GetEvent(ievent) ;
    
    //DP:Extract vertex position
    Double_t vtx[3]={0.,0.,0.} ;

    TClonesArray * rp = gime->RecParticles() ;
    if(!rp) {
      AliError(Form("Event %d,  Can't find RecParticles", ievent)) ;
      return ;
    }
    TClonesArray * ts = gime->TrackSegments() ;
    if(!ts) {
      AliError(Form("Event %d,  Can't find TrackSegments", ievent)) ;
      return ;
    }
    TObjArray * emcrp = gime->EmcRecPoints() ;
    if(!emcrp){
      AliError(Form("Event %d,  Can't find EmcRecPoints", ievent)) ;
      return ;
    }
    
    
    //=========== Make spectrum of the primary photons
    const TParticle * primary ;
    Int_t iPrimary ;
    for( iPrimary = 0 ; iPrimary < fRunLoader->Stack()->GetNprimary() ; iPrimary++){
      primary = fRunLoader->Stack()->Particle(iPrimary) ;
      Int_t primaryType = primary->GetPdgCode() ;
      if( primaryType == 22 ) {
       //check, if photons folls onto PHOS
       Int_t moduleNumber ;
       Double_t primX, primZ ;
       phosgeom->ImpactOnEmc(vtx,primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
       if(moduleNumber)
         hPrimary->Fill(primary->Energy()) ;
       
      }
      
    }
    
    //========== Now scan over RecParticles            
    const AliPHOSRecParticle * recParticle ;
    Int_t iRecParticle ;
    for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast(); iRecParticle++ ){
      recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
      //fill histo spectrum of all RecParticles
      hAllRP->Fill(CorrectedEnergy(recParticle->Energy())) ;
      
      //==========find the closest primary       
      Int_t moduleNumberRec ;
      Double_t recX, recZ ;
      phosgeom->ImpactOnEmc(vtx,recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
      
      Double_t minDistance  = 100. ;
      Int_t closestPrimary = -1 ;
      
      //extract list of primaries: it is stored at EMC RecPoints
      Int_t emcIndex = ((AliPHOSTrackSegment *) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
      Int_t numberofprimaries ;
      Int_t * listofprimaries  = ((AliPHOSEmcRecPoint *) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries)  ;
      Int_t index ;
      const TParticle * primary ;
      Double_t distance = minDistance ;
      Double_t dX, dZ; 
      Double_t dXmin = 0.; 
      Double_t dZmin = 0. ;
      for ( index = 0 ; index < numberofprimaries ; index++){
       primary = fRunLoader->Stack()->Particle(listofprimaries[index]) ;
       Int_t moduleNumber ;
       Double_t primX, primZ ;
       phosgeom->ImpactOnEmc(vtx,primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
       if(moduleNumberRec == moduleNumber) {
         dX = recX - primX;
         dZ = recZ - primZ;
         distance = TMath::Sqrt(dX*dX + dZ*dZ) ;
         if(minDistance > distance) {
           minDistance = distance ;
           dXmin = dX;
           dZmin = dZ;
           closestPrimary = listofprimaries[index] ;
         }
       }
      }
      
      //===========define the "type" of closest primary
      if(closestPrimary >=0 ){
       Int_t primaryCode = -1;
       const TParticle * primary = fRunLoader->Stack()->Particle(closestPrimary) ;
       Int_t primaryType = primary->GetPdgCode() ;
       if(primaryType == 22) // photon ?
         primaryCode = 0 ;
       else
         if(primaryType == 2112) // neutron
           primaryCode = 1 ; 
         else
           if(primaryType == -2112) // Anti neutron
             primaryCode = 2 ;
           else
             if(primaryType == -2122) //Anti proton
              primaryCode = 4 ;
             else {
              TParticle tempo(*primary) ; 
              if(tempo.GetPDG()->Charge())
                primaryCode = 3 ;
             }

       //==========Now look at the type of RecParticle
       Float_t energy = CorrectedEnergy(recParticle->Energy()) ;
       if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST){
         hPhot->Fill(energy ) ;        
         switch(primaryCode){
         case 0:
           hPhotReg->Fill(energy ) ; 
           break ;
         case 1:
           hNReg->Fill(energy ) ; 
           break ;
         case 2:
           hNBarReg->Fill(energy ) ; 
           break ;
         case 3:
           hChargedReg->Fill(energy ) ;
           break ;
         case 4:
           hPbarReg->Fill(energy ) ;
           break ;
         default:
           break ;
         }
       }
       if((recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)||
          (recParticle->GetType() == AliPHOSFastRecParticle::kCHARGEDEMFAST)||
          (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW)||
          (recParticle->GetType() == AliPHOSFastRecParticle::kCHARGEDEMSLOW) ){ //with EM shower
         hShape->Fill(energy ) ;
         switch(primaryCode){
         case 0:
           hPhotEM->Fill(energy ) ; 
           break ;
         case 1:
           hNEM->Fill(energy ) ; 
           break ;
         case 2:
           hNBarEM->Fill(energy ) ; 
           break ;
         case 3:
           hChargedEM->Fill(energy ) ; 
           break ;
         case 4:
           hPbarEM->Fill(energy ) ; 
           break ;
         default:
           break ;
         }
       }
       
       if((recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)||
          (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALHAFAST) ||
          (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW) ||
          (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALHASLOW) ) //nuetral
         hVeto->Fill(energy ) ;
       
       //fill number of primaries identified as ...
       if(primaryCode >= 0) // Primary code defined
         counter[recParticle->GetType()][primaryCode]++ ; 
       
      }
      
    } // no closest primary found
  }     
  
  
  //===================  SaveHistograms
  cfile->cd() ;
  hPrimary->Write(0,kOverwrite); 
  hAllRP->Write(0,kOverwrite);  
  hPhot->Write(0,kOverwrite);  
  hShape->Write(0,kOverwrite); 
  hVeto->Write(0,kOverwrite);  
  hPhotReg->Write(0,kOverwrite); 
  hPhotEM->Write(0,kOverwrite);   
  hNReg ->Write(0,kOverwrite);  
  hNEM  ->Write(0,kOverwrite); 
  hNBarReg ->Write(0,kOverwrite); 
  hNBarEM  ->Write(0,kOverwrite); 
  hChargedReg ->Write(0,kOverwrite); 
  hChargedEM  ->Write(0,kOverwrite); 
  hPbarReg ->Write(0,kOverwrite); 
  hPbarEM  ->Write(0,kOverwrite); 
  
  cfile->Write(0,kOverwrite); 
  cfile->Close();
  delete cfile ;
 
  
  //print Final Table
  maxevent = (Int_t)gAlice->TreeE()->GetEntries() ; 

  TString message ; 
  message  = "Resolutions: Analyzed %d event(s)\n" ; 
 
  message += "        Primary:    Photon  Neutron  Antineutron  Charged hadron   AntiProton\n" ; 
  message += "--------------------------------------------------------------------------------\n" ;
  message += "         kGAMMA: " ; 
  message += "%d %d %d %d %d\n" ; 
  message += "       kGAMMAHA: " ;
  message += "%d %d %d %d %d\n" ; 
  message += "     kNEUTRALEM: " ; 
  message += "%d %d %d %d %d\n" ; 
  message += "     kNEUTRALHA: " ; 
  message += "%d %d %d %d %d\n" ;
  message += "      kABSURDEM: ";
  message += "%d %d %d %d %d\n" ;
  message += "      kABSURDHA: " ;
  message += "%d %d %d %d %d\n" ;
  message += "      kELECTRON: " ;
  message += "%d %d %d %d %d\n" ;
  message += "     kCHARGEDHA: " ;  
  message += "%d %d %d %d %d\n" ;
   
  message += "--------------------------------------------------------------------------------" ;

 
  Int_t totalInd = 0 ;
  for(i1 = 0; i1<8; i1++)
    for(i2 = 0; i2<5; i2++)
      totalInd+=counter[i1][i2] ;
  message += "Indentified particles: %d" ; 
  
 AliInfo(Form(message.Data(), maxevent, 
      counter[2][0], counter[2][1], counter[2][2], counter[2][3], counter[2][4], 
      counter[3][0], counter[3][1], counter[3][2], counter[3][3], counter[3][4], 
      counter[0][0], counter[0][1], counter[0][2], counter[0][3], counter[0][4], 
      counter[1][0], counter[1][1], counter[1][2], counter[1][3], counter[1][4], 
      counter[4][0], counter[4][1], counter[4][2], counter[4][3], counter[4][4], 
      counter[5][0], counter[5][1], counter[5][2], counter[5][3], counter[5][4], 
      counter[6][0], counter[6][1], counter[6][2], counter[6][3], counter[6][4], 
      counter[7][0], counter[7][1], counter[7][2], counter[7][3], counter[7][4], 
      totalInd )) ;

}