PHOS/AliPHOSAnalyze.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 // Algorythm class to analyze PHOS events. In this class we demostrate,
  20 // how to handle reconstructed objects with AliPHSOIndexToObject.
  21 // As an example we propose sulotions for four most frequently used tasks:
  22 //    DrawRecon(...) - to draw reconstructed objects in the PHOS plane,
  23 //                     very usefull in the debuging
  24 //    InvarianMass(...) - to calculate "REAL" and "MIXED" photon pairs
  25 //                        invariant mass distributions
  26 //    EnergyResoluition(...) -\ Energy and position resolutions of the
  27 //    PositionResolution(...)-/ reconstructed photons
  28 //    Contamination(...) - calculates contamination of the photon spectrum and
  29 //                         pobability of reconstruction of several primaries as
  30 //                         kGAMMA,kELECTRON etc.
  31 //
  32 //    User Case:
  33 //    root [0] AliPHOSAnalyze * a = new AliPHOSAnalyze("galice.root")
  34 //                    // set the file you want to analyse
  35 //    root [1] a->DrawRecon(1,3)
  36 //                    // plot RecObjects, made in event 1, PHOS module 3
  37 //    root [2] a->DrawRecon(1,3,"PHOSRP","another PID")
  38 //                    // plot RecObjets made in the event 1, PHOS module 3,
  39 //                    // produced in the another reconstruction pass,
  40 //                    // which produced PHOS RecParticles ("PHOSRP") with
  41 //                    // title "another PID".
  42 //    root [3] a->InvariantMass()
  43 //                    // Calculates "REAL" and "MIXED" invariant mass
  44 //                    // distributions of kGAMMA and (kGAMMA+kNEUTRALEM)
  45 //                    // and APPENDS this to the file "invmass.root"
  46 //    root [4] a->PositionResolution()
  47 //                    // calculates two dimentional histos: energy of the primary
  48 //                    // photon vs distance betwin incedence point and reconstructed
  49 //                    // poisition. One can analyse the produced file position.root
  50 //                    // with macro PhotonPosition.C
  51 //    root [5] a->EnergyResolution()
  52 //                    // calculates two dimentional histos: energy of the primary
  53 //                    // photon vs energy of the reconstructed particle. One can
  54 //                    // analyse the produced file energy.root
  55 //                    // with macro PhotonEnergy.C
  56 //    root [6] a->Contamination()
  57 //                    // fills spectra of primary photons and several kinds of
  58 //                    // reconstructed particles, so that analyzing them one can
  59 //                    // estimate conatmination, efficiency of registration etc.
  60 //*--
  61 //*-- Author: Dmitri Peressounko (SUBATECH & RRC Kurchatov Institute)
  62 //////////////////////////////////////////////////////////////////////////////
  63
  64
  65 // --- ROOT system ---
  66
  67 #include "TFile.h"
  68 #include "TH1.h"
  69 #include "TH2.h"
  70 #include "TH2.h"
  71 #include "TParticle.h"
  72 #include "TClonesArray.h"
  73 #include "TTree.h"
  74 #include "TMath.h"
  75 #include "TROOT.h"
  76 #include "TFolder.h"
  77
  78 // --- Standard library ---
  79
  80 #include <iomanip.h>
  81
  82 // --- AliRoot header files ---
  83
  84 #include "AliRun.h"
  85 #include "AliPHOSv1.h"
  86 #include "AliPHOSAnalyze.h"
  87 #include "AliPHOSDigit.h"
  88 #include "AliPHOSSDigitizer.h"
  89 #include "AliPHOSTrackSegment.h"
  90 #include "AliPHOSRecParticle.h"
  91 #include "AliPHOSCpvRecPoint.h"
  92 #include "AliPHOSGetter.h"
  93
  94
  95 ClassImp(AliPHOSAnalyze)
  96
  97 //____________________________________________________________________________
  98   AliPHOSAnalyze::AliPHOSAnalyze()
  99 {
 100   // default ctor (useless)
 101   fCorrection = 1.2 ;  //Value calculated for default parameters of reconstruction
 102 }
 103
 104 //____________________________________________________________________________
 105 AliPHOSAnalyze::AliPHOSAnalyze(Text_t * fileName)
 106 {
 107   // ctor: analyze events from root file "name"
 108   ffileName = fileName ;
 109   fCorrection = 1.05 ;  //Value calculated for default parameters of reconstruction
 110 }
 111
 112 //____________________________________________________________________________
 113 AliPHOSAnalyze::AliPHOSAnalyze(const AliPHOSAnalyze & ana)
 114 {
 115   // copy ctor
 116   ( (AliPHOSAnalyze &)ana ).Copy(*this) ;
 117 }
 118
 119 //____________________________________________________________________________
 120 AliPHOSAnalyze::~AliPHOSAnalyze()
 121 {
 122   // dtor
 123
 124 }
 125 //____________________________________________________________________________
 126 void AliPHOSAnalyze::DrawRecon(Int_t Nevent,Int_t Nmod,const char * branchName,const char* branchTitle){
 127   //Draws pimary particles and reconstructed
 128   //digits, RecPoints, RecPartices etc
 129   //for event Nevent in the module Nmod.
 130
 131   //========== Create ObjectGetter
 132   AliPHOSGetter * gime = AliPHOSGetter::GetInstance(ffileName.Data(),branchTitle) ;
 133   const AliPHOSGeometry * phosgeom = gime->PHOSGeometry() ;
 134   gime->Event(Nevent);
 135
 136   Int_t nx = phosgeom->GetNPhi() ;
 137   Int_t nz = phosgeom->GetNZ() ;
 138   Float_t * cri= phosgeom->GetEMCAGeometry()->GetCrystalHalfSize() ;
 139   Float_t x = nx*cri[0] ;
 140   Float_t z = nz*cri[2] ;
 141   Int_t nxCPV = (Int_t) (nx*phosgeom->GetPadSizePhi()/(2.*cri[0])) ;
 142   Int_t nzCPV = (Int_t) (nz*phosgeom->GetPadSizeZ()/(2.*cri[2])) ;
 143
 144   TH2F * emcDigits    = new TH2F("emcDigits","EMC digits",  nx,-x,x,nz,-z,z);
 145   TH2F * emcSdigits   = new TH2F("emcSdigits","EMC sdigits", nx,-x,x,nz,-z,z);
 146   TH2F * emcRecPoints = new TH2F("emcRecPoints","EMC RecPoints",nx,-x,x,nz,-z,z);
 147   TH2F * cpvSdigits   = new TH2F("cpvSdigits","CPV sdigits", nx,-x,x,nz,-z,z);
 148   TH2F * cpvDigits    = new TH2F("cpvDigits","CPV digits",   nxCPV,-x,x,nzCPV,-z,z) ;
 149   TH2F * cpvRecPoints = new TH2F("cpvCl","CPV RecPoints",    nxCPV,-x,x,nzCPV,-z,z) ;
 150   TH2F * nbar         = new TH2F("nbar","Primary nbar",    nx,-x,x,nz,-z,z);
 151   TH2F * phot         = new TH2F("phot","Primary Photon",  nx,-x,x,nz,-z,z);
 152   TH2F * charg        = new TH2F("charg","Primary charged",nx,-x,x,nz,-z,z);
 153   TH2F * recPhot      = new TH2F("recPhot","RecParticles with primary Photon",nx,-x,x,nz,-z,z);
 154   TH2F * recNbar      = new TH2F("recNbar","RecParticles with primary Nbar",  nx,-x,x,nz,-z,z);
 155
 156
 157   //Plot Primary Particles
 158   const TParticle * primary ;
 159   Int_t iPrimary ;
 160   for ( iPrimary = 0 ; iPrimary < gime->NPrimaries() ; iPrimary++)
 161     {
 162       primary = gime->Primary(iPrimary) ;
 163       Int_t primaryType = primary->GetPdgCode() ;
 164       if( (primaryType == 211)||(primaryType == -211)||(primaryType == 2212)||(primaryType == -2212)
 165           ||(primaryType == 11)||(primaryType == -11) ) {
 166         Int_t moduleNumber ;
 167         Double_t primX, primZ ;
 168         phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
 169         if(moduleNumber==Nmod)
 170           charg->Fill(primZ,primX,primary->Energy()) ;
 171       }
 172       if( primaryType == 22 ) {
 173         Int_t moduleNumber ;
 174         Double_t primX, primZ ;
 175         phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
 176         if(moduleNumber==Nmod)
 177           phot->Fill(primZ,primX,primary->Energy()) ;
 178       }
 179       else{
 180         if( primaryType == -2112 ) {
 181           Int_t moduleNumber ;
 182           Double_t primX, primZ ;
 183           phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
 184           if(moduleNumber==Nmod)
 185             nbar->Fill(primZ,primX,primary->Energy()) ;
 186         }
 187       }
 188     }
 189
 190
 191   Int_t iSDigit ;
 192   AliPHOSDigit * sdigit ;
 193   TClonesArray * sdigits = gime->SDigits() ;
 194   if(sdigits){
 195     for(iSDigit = 0; iSDigit < sdigits->GetEntriesFast() ; iSDigit++)
 196       {
 197         sdigit = (AliPHOSDigit *) sdigits->At(iSDigit) ;
 198         Int_t relid[4];
 199         phosgeom->AbsToRelNumbering(sdigit->GetId(), relid) ;
 200         Float_t x,z ;
 201         phosgeom->RelPosInModule(relid,x,z) ;
 202         Float_t e = gime->SDigitizer()->Calibrate(sdigit->GetAmp()) ;
 203         if(relid[0]==Nmod){
 204           if(relid[1]==0)  //EMC
 205             emcSdigits->Fill(x,z,e) ;
 206           if( relid[1]!=0 )
 207             cpvSdigits->Fill(x,z,e) ;
 208         }
 209       }
 210   }
 211
 212   //Plot digits
 213   Int_t iDigit ;
 214   AliPHOSDigit * digit ;
 215   TClonesArray * digits = gime->Digits();
 216   if(digits) {
 217     for(iDigit = 0; iDigit < digits->GetEntriesFast(); iDigit++)
 218       {
 219         digit = (AliPHOSDigit *) digits->At(iDigit) ;
 220         Int_t relid[4];
 221         phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
 222         Float_t x,z ;
 223         phosgeom->RelPosInModule(relid,x,z) ;
 224         Float_t e = gime->SDigitizer()->Calibrate(digit->GetAmp()) ;
 225         if(relid[0]==Nmod){
 226           if(relid[1]==0)  //EMC
 227             emcDigits->Fill(x,z,e) ;
 228           if( relid[1]!=0 )
 229             cpvDigits->Fill(x,z,e) ;
 230         }
 231       }
 232   }
 233
 234
 235   //Plot RecPoints
 236   Int_t irecp ;
 237   TVector3 pos ;
 238   TObjArray * emcrp = gime->EmcRecPoints() ;
 239   if(emcrp) {
 240     for(irecp = 0; irecp < emcrp->GetEntriesFast() ; irecp ++){
 241       AliPHOSEmcRecPoint * emc = (AliPHOSEmcRecPoint *) emcrp->At(irecp) ;
 242       if(emc->GetPHOSMod()==Nmod){
 243         emc->GetLocalPosition(pos) ;
 244         emcRecPoints->Fill(pos.X(),pos.Z(),emc->GetEnergy());
 245       }
 246     }
 247   }
 248
 249   TObjArray * cpvrp = gime->CpvRecPoints() ;
 250   if(cpvrp) {
 251     for(irecp = 0; irecp < cpvrp->GetEntriesFast() ; irecp ++){
 252       AliPHOSRecPoint * cpv = (AliPHOSCpvRecPoint *) cpvrp->At(irecp) ;
 253       if(cpv->GetPHOSMod()==Nmod){
 254         cpv->GetLocalPosition(pos) ;
 255         cpvRecPoints->Fill(pos.X(),pos.Z(),cpv->GetEnergy());
 256       }
 257     }
 258   }
 259
 260   //Plot RecParticles
 261   AliPHOSRecParticle * recParticle ;
 262   Int_t iRecParticle ;
 263   TClonesArray * rp = gime->RecParticles() ;
 264   TClonesArray * ts = gime->TrackSegments() ;
 265   if(rp && ts && emcrp) {
 266     for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast() ; iRecParticle++ )
 267       {
 268         recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
 269         Int_t moduleNumberRec ;
 270         Double_t recX, recZ ;
 271         phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
 272         if(moduleNumberRec == Nmod){
 273
 274           Double_t minDistance = 5. ;
 275           Int_t closestPrimary = -1 ;
 276
 277           //extract list of primaries: it is stored at EMC RecPoints
 278           Int_t emcIndex = ((AliPHOSTrackSegment *) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
 279           Int_t numberofprimaries ;
 280           Int_t * listofprimaries  = ((AliPHOSRecPoint*) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries)  ;
 281           Int_t index ;
 282           const TParticle * primary ;
 283           Double_t distance = minDistance ;
 284
 285           for ( index = 0 ; index < numberofprimaries ; index++){
 286             primary = gime->Primary(listofprimaries[index]) ;
 287             Int_t moduleNumber ;
 288             Double_t primX, primZ ;
 289             phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
 290             if(moduleNumberRec == moduleNumber)
 291               distance = TMath::Sqrt((recX-primX)*(recX-primX)+(recZ-primZ)*(recZ-primZ) ) ;
 292             if(minDistance > distance)
 293               {
 294                 minDistance = distance ;
 295                 closestPrimary = listofprimaries[index] ;
 296               }
 297           }
 298
 299           if(closestPrimary >=0 ){
 300
 301             Int_t primaryType = gime->Primary(closestPrimary)->GetPdgCode() ;
 302
 303             if(primaryType==22)
 304               recPhot->Fill(recZ,recX,recParticle->Energy()) ;
 305             else
 306               if(primaryType==-2112)
 307                 recNbar->Fill(recZ,recX,recParticle->Energy()) ;
 308           }
 309         }
 310       }
 311
 312   }
 313
 314   //Plot made histograms
 315   emcSdigits->Draw("box") ;
 316   emcDigits->SetLineColor(5) ;
 317   emcDigits->Draw("box") ;
 318   emcRecPoints->SetLineColor(2) ;
 319   emcRecPoints->Draw("boxsame") ;
 320   cpvSdigits->SetLineColor(1) ;
 321   cpvSdigits->Draw("box") ;
 322   charg->SetLineColor(2) ;
 323   charg->Draw("boxsame") ;
 324   nbar->SetLineColor(6) ;
 325   nbar->Draw("boxsame") ;
 326
 327 }
 328 //____________________________________________________________________________
 329 void AliPHOSAnalyze::Ls(){
 330   //lists branches and titles of PHOS-related branches of TreeR, TreeD, TreeS
 331
 332   AliPHOSGetter::GetInstance(ffileName.Data()) ;
 333
 334   Int_t ibranch;
 335   TObjArray * branches;
 336
 337   branches = gAlice->TreeS()->GetListOfBranches() ;
 338
 339   cout << "TreeS: " << endl ;
 340   for(ibranch = 0;ibranch <branches->GetEntries();ibranch++){
 341     TBranch * branch=(TBranch *) branches->At(ibranch) ;
 342     if(strstr(branch->GetName(),"PHOS") )
 343       cout << "       " << branch->GetName() << "     " << branch->GetTitle() << endl ;
 344   }
 345   cout << endl ;
 346
 347   branches = gAlice->TreeD()->GetListOfBranches() ;
 348
 349   cout << "TreeD: " << endl ;
 350   for(ibranch = 0;ibranch <branches->GetEntries();ibranch++){
 351     TBranch * branch=(TBranch *) branches->At(ibranch) ;
 352     if(strstr(branch->GetName(),"PHOS") )
 353       cout << "       " << branch->GetName() << "     " << branch->GetTitle() << endl ;
 354   }
 355   cout << endl ;
 356
 357
 358   branches = gAlice->TreeR()->GetListOfBranches() ;
 359
 360   cout << "TreeR: " << endl ;
 361   for(ibranch = 0;ibranch <branches->GetEntries();ibranch++){
 362     TBranch * branch=(TBranch *) branches->At(ibranch) ;
 363     if(strstr(branch->GetName(),"PHOS") )
 364       cout << "       " << branch->GetName() << "     " << branch->GetTitle() << endl ;
 365   }
 366   cout << endl ;
 367
 368
 369 }
 370 //____________________________________________________________________________
 371  void AliPHOSAnalyze::InvariantMass(const char* branchTitle)
 372 {
 373   // Calculates Real and Mixed invariant mass distributions
 374   AliPHOSGetter * gime  = AliPHOSGetter::GetInstance(ffileName.Data(),branchTitle) ;
 375
 376   Int_t nMixedEvents = 4 ; //# of events used for calculation of 'mixed' distribution
 377
 378   //opening file
 379   TFile * mfile = new TFile("invmass.root","update");
 380
 381   //========== Reading /Booking Histograms
 382   TH2D * hRealEM = 0 ;
 383   hRealEM = (TH2D*) mfile->Get("hRealEM") ;
 384   if(hRealEM == 0)
 385     hRealEM = new TH2D("hRealEM",   "Real for EM particles",      250,0.,1.,40,0.,4.) ;
 386   TH2D * hRealPhot = 0 ;
 387
 388   hRealPhot = (TH2D*)mfile->Get("hRealPhot");
 389   if(hRealPhot == 0)
 390     hRealPhot = new TH2D("hRealPhot", "Real for kPhoton particles", 250,0.,1.,40,0.,4.) ;
 391
 392   TH2D * hMixedEM = 0 ;
 393   hMixedEM = (TH2D*) mfile->Get("hMixedEM") ;
 394   if(hMixedEM == 0)
 395     hMixedEM = new TH2D("hMixedEM",  "Mixed for EM particles",     250,0.,1.,40,0.,4.) ;
 396
 397   TH2D * hMixedPhot = 0 ;
 398   hMixedPhot = (TH2D*) mfile->Get("hMixedPhot") ;
 399   if(hMixedPhot == 0)
 400     hMixedPhot = new TH2D("hMixedPhot","Mixed for kPhoton particles",250,0.,1.,40,0.,4.) ;
 401
 402
 403   //reading event and copyng it to TConesArray of all photons
 404
 405   TClonesArray * allRecParticleList  = new TClonesArray("AliPHOSRecParticle", 1000) ;
 406   Int_t * nRecParticles = new Int_t[nMixedEvents] ; // to mark boundaries of each event in the total list
 407   for(Int_t index = 0; index < nMixedEvents; index ++)
 408     nRecParticles[index] = 0 ;
 409   Int_t iRecPhot = 0 ;                // number of EM particles in total list
 410
 411   //scan over all events
 412   Int_t event ;
 413   Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ;
 414   //  for(event = 0; event < gime->MaxEvent(); event++  ){
 415   for(event = 0; event < maxevent; event++  ){
 416     gime->Event(event,"R");  //will read only TreeR
 417
 418     //copy EM RecParticles to the "total" list
 419     const AliPHOSRecParticle * recParticle ;
 420     Int_t iRecParticle ;
 421     TClonesArray * rp = gime->RecParticles() ;
 422     if(!rp){
 423       cout << "AliPHOSAnalyze::InvariantMass --> Can't find RecParticles " << endl ;
 424       return ;
 425     }
 426
 427     for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast(); iRecParticle++ )
 428       {
 429         recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
 430         if((recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)||
 431            (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW))
 432           new( (*allRecParticleList)[iRecPhot++] ) AliPHOSRecParticle(*recParticle) ;
 433       }
 434
 435     Int_t mevent = event%nMixedEvents ; //event number in the "mixed" cicle
 436     nRecParticles[mevent] = iRecPhot-1 ;
 437
 438     //check, if it is time to calculate invariant mass?
 439     Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ;
 440     if((mevent == 0) && (event +1 == maxevent)){
 441
 442       //   if((mevent == 0) && (event +1 == gime->MaxEvent())){
 443
 444       //calculate invariant mass:
 445       Int_t irp1,irp2 ;
 446       Int_t nCurEvent = 0 ;
 447
 448       for(irp1 = 0; irp1 < allRecParticleList->GetEntries()-1; irp1++){
 449         AliPHOSRecParticle * rp1 = (AliPHOSRecParticle *)allRecParticleList->At(irp1) ;
 450
 451         for(irp2 = irp1+1; irp2 < allRecParticleList->GetEntries(); irp2++){
 452           AliPHOSRecParticle * rp2 = (AliPHOSRecParticle *)allRecParticleList->At(irp2) ;
 453
 454           Double_t invMass ;
 455           invMass = (rp1->Energy()+rp2->Energy())*(rp1->Energy()+rp2->Energy())-
 456             (rp1->Px()+rp2->Px())*(rp1->Px()+rp2->Px())-
 457             (rp1->Py()+rp2->Py())*(rp1->Py()+rp2->Py())-
 458             (rp1->Pz()+rp2->Pz())*(rp1->Pz()+rp2->Pz()) ;
 459
 460           if(invMass> 0)
 461             invMass = TMath::Sqrt(invMass);
 462
 463           Double_t pt ;
 464           pt = TMath::Sqrt((rp1->Px()+rp2->Px() )*( rp1->Px()+rp2->Px() ) +
 465                            (rp1->Py()+rp2->Py() )*( rp1->Py()+rp2->Py() ) );
 466
 467           if(irp1 > nRecParticles[nCurEvent])
 468             nCurEvent++;
 469
 470           if(irp2 <= nRecParticles[nCurEvent]){ //'Real' event
 471             hRealEM->Fill(invMass,pt);
 472             if((rp1->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)&&
 473                (rp2->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST) )
 474               hRealPhot->Fill(invMass,pt);
 475           }
 476           else{
 477             hMixedEM->Fill(invMass,pt);
 478             if((rp1->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)&&
 479                (rp2->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST) )
 480               hMixedPhot->Fill(invMass,pt);
 481           } //real-mixed
 482
 483       } //loop over second rp
 484       }//loop over first rp
 485
 486       //Make some cleanings
 487       for(Int_t index = 0; index < nMixedEvents; index ++)
 488         nRecParticles[index] = 0 ;
 489       iRecPhot = 0 ;
 490       allRecParticleList->Clear() ;
 491
 492     }
 493   }
 494   delete allRecParticleList ;
 495
 496   //writing output
 497   mfile->cd();
 498
 499   hRealEM->Write(0,kOverwrite) ;
 500   hRealPhot->Write(0,kOverwrite) ;
 501   hMixedEM->Write(0,kOverwrite) ;
 502   hMixedPhot->Write(0,kOverwrite) ;
 503
 504   mfile->Write();
 505   mfile->Close();
 506   delete mfile ;
 507   delete nRecParticles;
 508
 509 }
 510
 511 //____________________________________________________________________________
 512  void AliPHOSAnalyze::EnergyResolution(const char * branchTitle)
 513 {
 514   //fills two dimentional histo: energy of primary vs. energy of reconstructed
 515
 516   TH2F * hAllEnergy = 0 ;  //all reconstructed with primary photon
 517   TH2F * hPhotEnergy= 0 ;  //kGamma  with primary photon
 518   TH2F * hEMEnergy  = 0 ;  //electromagnetic with primary photon
 519
 520   //opening file and reading histograms if any
 521   TFile * efile = new TFile("energy.root","update");
 522
 523   hAllEnergy = (TH2F*)efile->Get("hAllEnergy") ;
 524   if(hAllEnergy == 0)
 525     hAllEnergy  = new TH2F("hAllEnergy",  "Energy of any RP with primary photon",100, 0., 5., 100, 0., 5.);
 526
 527   hPhotEnergy =(TH2F*) efile->Get("hPhotEnergy") ;
 528   if(hPhotEnergy == 0)
 529     hPhotEnergy = new TH2F("hPhotEnergy", "Energy of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.);
 530
 531   hEMEnergy =(TH2F*) efile->Get("hEMEnergy");
 532   if(hEMEnergy == 0)
 533     hEMEnergy   = new TH2F("hEMEnergy",   "Energy of EM with primary photon",    100, 0., 5., 100, 0., 5.);
 534
 535
 536   AliPHOSGetter * gime = AliPHOSGetter::GetInstance(ffileName.Data(),branchTitle) ;
 537   const AliPHOSGeometry * phosgeom = gime->PHOSGeometry() ;
 538
 539   Int_t ievent;
 540   Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ;
 541   for ( ievent=0; ievent < maxevent ; ievent++){
 542
 543     //read the current event
 544     gime->Event(ievent) ;
 545
 546     const AliPHOSRecParticle * recParticle ;
 547     Int_t iRecParticle ;
 548     TClonesArray * rp = gime->RecParticles() ;
 549     if(!rp) {
 550       cout << "AliPHOSAnalyze::EnergyResolution --> Event " <<ievent
 551            << ",  Can't find RecParticles " << endl ;
 552       return ;
 553     }
 554     TClonesArray * ts = gime->TrackSegments() ;
 555     if(!ts) {
 556       cout << "AliPHOSAnalyze::EnergyResolution --> Event " <<ievent
 557            << ",  Can't find TrackSegments " << endl ;
 558       return ;
 559     }
 560     TObjArray * emcrp = gime->EmcRecPoints() ;
 561     if(!emcrp){
 562       cout << "AliPHOSAnalyze::EnergyResolution --> Event " <<ievent
 563            << ",  Can't find EmcRecPoints " << endl ;
 564       return ;
 565     }
 566
 567     for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast() ;iRecParticle++ ){
 568       recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
 569
 570       //find the closest primary
 571       Int_t moduleNumberRec ;
 572       Double_t recX, recZ ;
 573       phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
 574
 575       Double_t minDistance  = 100. ;
 576       Int_t closestPrimary = -1 ;
 577
 578       //extract list of primaries: it is stored at EMC RecPoints
 579       Int_t emcIndex = ((AliPHOSTrackSegment*) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
 580       Int_t numberofprimaries ;
 581       Int_t * listofprimaries  = ((AliPHOSEmcRecPoint*) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries)  ;
 582
 583       Int_t index ;
 584       const TParticle * primary ;
 585       Double_t distance = minDistance ;
 586       Double_t dX, dZ;
 587       Double_t dXmin = 0.;
 588       Double_t dZmin = 0. ;
 589       for ( index = 0 ; index < numberofprimaries ; index++){
 590         primary = gime->Primary(listofprimaries[index]) ;
 591         Int_t moduleNumber ;
 592         Double_t primX, primZ ;
 593         phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
 594         if(moduleNumberRec == moduleNumber) {
 595           dX = recX - primX;
 596           dZ = recZ - primZ;
 597           distance = TMath::Sqrt(dX*dX + dZ*dZ) ;
 598           if(minDistance > distance) {
 599             minDistance = distance ;
 600             dXmin = dX;
 601             dZmin = dZ;
 602             closestPrimary = listofprimaries[index] ;
 603           }
 604         }
 605       }
 606
 607       //if found primary, fill histograms
 608       if(closestPrimary >=0 ){
 609         const TParticle * primary = gime->Primary(closestPrimary) ;
 610         if(primary->GetPdgCode() == 22){
 611           hAllEnergy->Fill(primary->Energy(), recParticle->Energy()) ;
 612           if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST){
 613             hPhotEnergy->Fill(primary->Energy(), recParticle->Energy() ) ;
 614             hEMEnergy->Fill(primary->Energy(), recParticle->Energy() ) ;
 615           }
 616           else
 617             if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW)
 618               hEMEnergy->Fill(primary->Energy(), recParticle->Energy() ) ;
 619         }
 620       }
 621     }
 622   }
 623
 624   //write filled histograms
 625   efile->cd() ;
 626   hAllEnergy->Write(0,kOverwrite) ;
 627   hPhotEnergy->Write(0,kOverwrite) ;
 628   hEMEnergy->Write(0,kOverwrite)  ;
 629   //  efile->Write() ;
 630   efile->Close() ;
 631   delete efile ;
 632
 633 }
 634 //____________________________________________________________________________
 635 void AliPHOSAnalyze::PositionResolution(const char * branchTitle)
 636 {
 637   //fills two dimentional histo: energy vs. primary - reconstructed distance
 638
 639
 640
 641   TH2F * hAllPosition  = 0;     // Position of any RP with primary photon
 642   TH2F * hPhotPosition = 0;    // Position of kGAMMA with primary photon
 643   TH2F * hEMPosition   = 0;      // Position of EM with primary photon
 644
 645   TH1F * hAllPositionX = 0;    // X-Position Resolution of photons with photon primary
 646   TH1F * hAllPositionZ = 0;    // Z-Position Resolution of photons with photon primary
 647
 648
 649   //opening file and reading histograms if any
 650   TFile * pfile = new TFile("position.root","update");
 651
 652   hAllPosition = (TH2F*)pfile->Get("hAllPosition");
 653   if(hAllPosition == 0)
 654     hAllPosition  = new TH2F("hAllPosition",
 655                              "Position of any RP with primary photon",100, 0., 5., 100, 0., 5.);
 656   hPhotPosition= (TH2F*)pfile->Get("hPhotPosition");
 657   if(hPhotPosition == 0)
 658     hPhotPosition = new TH2F("hPhotPosition",
 659                              "Position of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.);
 660   hEMPosition= (TH2F*)pfile->Get("hEMPosition") ;
 661   if(hEMPosition == 0)
 662     hEMPosition   = new TH2F("hEMPosition",
 663                              "Position of EM with primary photon",    100, 0., 5., 100, 0., 5.);
 664   hAllPositionX = (TH1F*)pfile->Get("hAllPositionX") ;
 665   if(hAllPositionX == 0)
 666     hAllPositionX = new TH1F("hAllPositionX",
 667                              "Delta X of any RP with primary photon",100, -2., 2.);
 668   hAllPositionZ =(TH1F*) pfile->Get("hAllPositionZ") ;
 669   if(hAllPositionZ == 0)
 670     hAllPositionZ = new TH1F("hAllPositionZ",
 671                              "Delta X of any RP with primary photon",100, -2., 2.);
 672
 673
 674   AliPHOSGetter * gime = AliPHOSGetter::GetInstance(ffileName.Data(),branchTitle) ;
 675   const AliPHOSGeometry * phosgeom = gime->PHOSGeometry() ;
 676
 677   Int_t ievent;
 678   Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ;
 679   for ( ievent=0; ievent < maxevent ; ievent++){
 680
 681     //read the current event
 682     gime->Event(ievent) ;
 683     TClonesArray * rp = gime->RecParticles() ;
 684     if(!rp) {
 685       cout << "AliPHOSAnalyze::PositionResolution --> Event " <<ievent
 686            << ",  Can't find RecParticles " << endl ;
 687       return ;
 688     }
 689     TClonesArray * ts = gime->TrackSegments() ;
 690     if(!ts) {
 691       cout << "AliPHOSAnalyze::PositionResolution --> Event " <<ievent
 692            << ",  Can't find TrackSegments " << endl ;
 693       return ;
 694     }
 695     TObjArray * emcrp = gime->EmcRecPoints() ;
 696     if(!emcrp){
 697       cout << "AliPHOSAnalyze::PositionResolution --> Event " <<ievent
 698            << ",  Can't find EmcRecPoints " << endl ;
 699       return ;
 700     }
 701
 702
 703     const AliPHOSRecParticle * recParticle ;
 704     Int_t iRecParticle ;
 705     for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast(); iRecParticle++ ){
 706       recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
 707
 708       //find the closest primary
 709       Int_t moduleNumberRec ;
 710       Double_t recX, recZ ;
 711       phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
 712
 713       Double_t minDistance  = 100. ;
 714       Int_t closestPrimary = -1 ;
 715
 716       //extract list of primaries: it is stored at EMC RecPoints
 717       Int_t emcIndex = ((AliPHOSTrackSegment*) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
 718       Int_t numberofprimaries ;
 719       Int_t * listofprimaries  = ((AliPHOSEmcRecPoint *) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries)  ;
 720
 721       Int_t index ;
 722       const TParticle * primary ;
 723       Double_t distance = minDistance ;
 724       Double_t dX = 1000; // incredible number
 725       Double_t dZ = 1000; // for the case if no primary will be found
 726       Double_t dXmin = 0.;
 727       Double_t dZmin = 0. ;
 728       for ( index = 0 ; index < numberofprimaries ; index++){
 729         primary = gime->Primary(listofprimaries[index]) ;
 730         Int_t moduleNumber ;
 731         Double_t primX, primZ ;
 732         phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
 733         if(moduleNumberRec == moduleNumber) {
 734           dX = recX - primX;
 735           dZ = recZ - primZ;
 736           distance = TMath::Sqrt(dX*dX + dZ*dZ) ;
 737           if(minDistance > distance) {
 738             minDistance = distance ;
 739             dXmin = dX;
 740             dZmin = dZ;
 741             closestPrimary = listofprimaries[index] ;
 742           }
 743         }
 744       }
 745
 746       //if found primary, fill histograms
 747       if(closestPrimary >=0 ){
 748         const TParticle * primary = gime->Primary(closestPrimary) ;
 749         if(primary->GetPdgCode() == 22){
 750           hAllPosition->Fill(primary->Energy(), minDistance) ;
 751           hAllPositionX->Fill(primary->Energy(), dX) ;
 752           hAllPositionZ->Fill(primary->Energy(), dZ) ;
 753           if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST){
 754             hPhotPosition->Fill(primary->Energy(), minDistance ) ;
 755             hEMPosition->Fill(primary->Energy(), minDistance ) ;
 756           }
 757           else
 758             if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW)
 759               hEMPosition->Fill(primary->Energy(), minDistance ) ;
 760         }
 761       }
 762     }
 763   }
 764
 765   //Write output histgrams
 766   pfile->cd() ;
 767   hAllPosition->Write(0,kOverwrite) ;
 768   hAllPositionX->Write(0,kOverwrite) ;
 769   hAllPositionZ->Write(0,kOverwrite) ;
 770   hPhotPosition->Write(0,kOverwrite) ;
 771   hEMPosition->Write(0,kOverwrite) ;
 772   pfile->Write() ;
 773   pfile->Close() ;
 774   delete pfile ;
 775
 776
 777 }
 778 //____________________________________________________________________________
 779 void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){
 780 // fills spectra of primary photons and several kinds of
 781 // reconstructed particles, so that analyzing them one can
 782 // estimate conatmination, efficiency of registration etc.
 783
 784   //define several general histograms
 785   TH1F * hPrimary = 0;   //spectrum (P_t distribution) of primary photons
 786   TH1F * hAllRP   = 0;   //spectrum of all RecParticles in PHOS
 787   TH1F * hPhot    = 0;   //spectrum of kGAMMA RecParticles
 788   TH1F * hShape   = 0;   //spectrum of all EM RecParticles
 789   TH1F * hVeto    = 0;   //spectrum of all neutral RecParticles
 790
 791   //Now separate histograms in accoradance with primary
 792   //primary - photon
 793   TH1F * hPhotReg = 0;   //Registeres as photon
 794   TH1F * hPhotEM  = 0;   //Registered as EM
 795
 796   //primary - n
 797   TH1F * hNReg = 0;   //Registeres as photon
 798   TH1F * hNEM  = 0;   //Registered as EM
 799
 800   //primary - nBar
 801   TH1F * hNBarReg = 0;   //Registeres as photon
 802   TH1F * hNBarEM  = 0;   //Registered as EM
 803
 804   //primary - charged hadron (pBar excluded)
 805   TH1F * hChargedReg = 0;   //Registeres as photon
 806   TH1F * hChargedEM  = 0;   //Registered as EM
 807
 808   //primary - pBar
 809   TH1F * hPbarReg = 0;   //Registeres as photon
 810   TH1F * hPbarEM  = 0;   //Registered as EM
 811
 812
 813   //Reading histograms from the file
 814   TFile * cfile = new TFile("contamination.root","update") ;
 815
 816   //read general histograms
 817   hPrimary = (TH1F*) cfile->Get("hPrimary") ;
 818   if(hPrimary == 0)
 819     hPrimary= new TH1F("hPrimary", "Primary photon spectrum", 100, 0., 5.);
 820   hAllRP = (TH1F*)cfile->Get("hAllRP") ;
 821   if(hAllRP == 0)
 822     hAllRP = new TH1F("hAllRP","All Reconstructed particles", 100, 0., 5.);
 823   hPhot  = (TH1F*)cfile->Get("hPhot") ;
 824   if(hPhot == 0)
 825     hPhot = new TH1F("hPhot","All kGAMMA RecParticles",100, 0., 5.);
 826   hShape = (TH1F*) cfile->Get("hShape") ;
 827   if(hShape == 0)
 828     hShape = new TH1F("hShape","All particles with EM shower",100, 0., 5.);
 829   hVeto= (TH1F*)cfile->Get("hVeto") ;
 830   if(hVeto == 0)
 831     hVeto  = new TH1F("hVeto", "All uncharged particles",      100, 0., 5.);
 832
 833
 834   //primary - photon
 835   hPhotReg = (TH1F*)cfile->Get("hPhotReg");
 836   if(hPhotReg == 0)
 837     hPhotReg   = new TH1F("hPhotReg","Photon registered as photon",100, 0., 5.);
 838   hPhotEM  =(TH1F*)cfile->Get("hPhotEM");
 839   if(hPhotEM== 0)
 840     hPhotEM   = new TH1F("hPhotEM",  "Photon registered as EM", 100, 0., 5.);
 841
 842   //primary - n
 843   hNReg = (TH1F*)cfile->Get("hNReg");
 844   if(hNReg== 0)
 845    hNReg      = new TH1F("hNReg",   "N registered as photon",              100, 0., 5.);
 846   hNEM  = (TH1F*)cfile->Get("hNEM");
 847   if(hNEM== 0)
 848     hNEM      = new TH1F("hNEM",    "N registered as EM",      100, 0., 5.);
 849
 850   //primary - nBar
 851   hNBarReg =(TH1F*)cfile->Get("hNBarReg");
 852   if(hNBarReg== 0)
 853    hNBarReg   = new TH1F("hNBarReg", "NBar registered as photon",           100, 0., 5.);
 854   hNBarEM  =(TH1F*)cfile->Get("hNBarEM");
 855   if(hNBarEM== 0)
 856     hNBarEM   = new TH1F("hNBarEM",  "NBar registered as EM",   100, 0., 5.);
 857
 858   //primary - charged hadron (pBar excluded)
 859   hChargedReg = (TH1F*)cfile->Get("hChargedReg");
 860   if(hChargedReg== 0)
 861     hChargedReg= new TH1F("hChargedReg", "Charged hadron registered as photon",100, 0., 5.);
 862   hChargedEM  = (TH1F*)cfile->Get("hChargedEM");
 863   if(hChargedEM== 0)
 864     hChargedEM= new TH1F("hChargedEM","Charged registered as EM",100, 0., 5.);
 865
 866   //primary - pBar
 867   hPbarReg = (TH1F*)cfile->Get("hPbarReg");
 868   if(hPbarReg== 0)
 869     hPbarReg= new TH1F("hPbarReg", "pBar registered as photon",100, 0., 5.);
 870   hPbarEM  = (TH1F*)cfile->Get("hPbarEM");
 871   if(hPbarEM== 0)
 872     hPbarEM= new TH1F("hPbarEM","Pbar registered as EM",100, 0., 5.);
 873
 874
 875   //Now make some initializations
 876
 877   Int_t counter[8][5] ;      //# of registered particles
 878   Int_t i1,i2 ;
 879   for(i1 = 0; i1<8; i1++)
 880     for(i2 = 0; i2<5; i2++)
 881       counter[i1][i2] = 0 ;
 882
 883
 884
 885   AliPHOSGetter * gime = AliPHOSGetter::GetInstance(ffileName.Data(),RecPointsTitle) ;
 886   const AliPHOSGeometry * phosgeom = gime->PHOSGeometry() ;
 887
 888   Int_t ievent;
 889   Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ;
 890   for ( ievent=0; ievent < maxevent ; ievent++){
 891
 892     gime->Event(ievent) ;
 893
 894     TClonesArray * rp = gime->RecParticles() ;
 895     if(!rp) {
 896       cout << "AliPHOSAnalyze::Contamination --> Event " <<ievent
 897            << ",  Can't find RecParticles " << endl ;
 898       return ;
 899     }
 900     TClonesArray * ts = gime->TrackSegments() ;
 901     if(!ts) {
 902       cout << "AliPHOSAnalyze::Contamination --> Event " <<ievent
 903            << ",  Can't find TrackSegments " << endl ;
 904       return ;
 905     }
 906     TObjArray * emcrp = gime->EmcRecPoints() ;
 907     if(!emcrp){
 908       cout << "AliPHOSAnalyze::Contamination --> Event " <<ievent
 909            << ",  Can't find EmcRecPoints " << endl ;
 910       return ;
 911     }
 912
 913
 914     //=========== Make spectrum of the primary photons
 915     const TParticle * primary ;
 916     Int_t iPrimary ;
 917     for( iPrimary = 0 ; iPrimary < gime->NPrimaries() ; iPrimary++){
 918       primary = gime->Primary(iPrimary) ;
 919       Int_t primaryType = primary->GetPdgCode() ;
 920       if( primaryType == 22 ) {
 921         //check, if photons folls onto PHOS
 922         Int_t moduleNumber ;
 923         Double_t primX, primZ ;
 924         phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
 925         if(moduleNumber)
 926           hPrimary->Fill(primary->Energy()) ;
 927
 928       }
 929
 930     }
 931
 932     //========== Now scan over RecParticles
 933     const AliPHOSRecParticle * recParticle ;
 934     Int_t iRecParticle ;
 935     for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast(); iRecParticle++ ){
 936       recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
 937       //fill histo spectrum of all RecParticles
 938       hAllRP->Fill(CorrectedEnergy(recParticle->Energy())) ;
 939
 940       //==========find the closest primary
 941       Int_t moduleNumberRec ;
 942       Double_t recX, recZ ;
 943       phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
 944
 945       Double_t minDistance  = 100. ;
 946       Int_t closestPrimary = -1 ;
 947
 948       //extract list of primaries: it is stored at EMC RecPoints
 949       Int_t emcIndex = ((AliPHOSTrackSegment *) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
 950       Int_t numberofprimaries ;
 951       Int_t * listofprimaries  = ((AliPHOSEmcRecPoint *) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries)  ;
 952       Int_t index ;
 953       const TParticle * primary ;
 954       Double_t distance = minDistance ;
 955       Double_t dX, dZ;
 956       Double_t dXmin = 0.;
 957       Double_t dZmin = 0. ;
 958       for ( index = 0 ; index < numberofprimaries ; index++){
 959         primary = gime->Primary(listofprimaries[index]) ;
 960         Int_t moduleNumber ;
 961         Double_t primX, primZ ;
 962         phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
 963         if(moduleNumberRec == moduleNumber) {
 964           dX = recX - primX;
 965           dZ = recZ - primZ;
 966           distance = TMath::Sqrt(dX*dX + dZ*dZ) ;
 967           if(minDistance > distance) {
 968             minDistance = distance ;
 969             dXmin = dX;
 970             dZmin = dZ;
 971             closestPrimary = listofprimaries[index] ;
 972           }
 973         }
 974       }
 975
 976       //===========define the "type" of closest primary
 977       if(closestPrimary >=0 ){
 978         Int_t primaryCode = -1;
 979         const TParticle * primary = gime->Primary(closestPrimary) ;
 980         Int_t primaryType = primary->GetPdgCode() ;
 981         if(primaryType == 22) // photon ?
 982           primaryCode = 0 ;
 983         else
 984           if(primaryType == 2112) // neutron
 985             primaryCode = 1 ;
 986           else
 987             if(primaryType == -2112) // Anti neutron
 988               primaryCode = 2 ;
 989             else
 990               if(primaryType == -2122) //Anti proton
 991                 primaryCode = 4 ;
 992               else {
 993                 TParticle tempo(*primary) ;
 994                 if(tempo.GetPDG()->Charge())
 995                   primaryCode = 3 ;
 996               }
 997
 998         //==========Now look at the type of RecParticle
 999         Float_t energy = CorrectedEnergy(recParticle->Energy()) ;
1000         if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST){
1001           hPhot->Fill(energy ) ;
1002           switch(primaryCode){
1003           case 0:
1004             hPhotReg->Fill(energy ) ;
1005             break ;
1006           case 1:
1007             hNReg->Fill(energy ) ;
1008             break ;
1009           case 2:
1010             hNBarReg->Fill(energy ) ;
1011             break ;
1012           case 3:
1013             hChargedReg->Fill(energy ) ;
1014             break ;
1015           case 4:
1016             hPbarReg->Fill(energy ) ;
1017             break ;
1018           default:
1019             break ;
1020           }
1021         }
1022         if((recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)||
1023            (recParticle->GetType() == AliPHOSFastRecParticle::kCHARGEDEMFAST)||
1024            (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW)||
1025            (recParticle->GetType() == AliPHOSFastRecParticle::kCHARGEDEMSLOW) ){ //with EM shower
1026           hShape->Fill(energy ) ;
1027           switch(primaryCode){
1028           case 0:
1029             hPhotEM->Fill(energy ) ;
1030             break ;
1031           case 1:
1032             hNEM->Fill(energy ) ;
1033             break ;
1034           case 2:
1035             hNBarEM->Fill(energy ) ;
1036             break ;
1037           case 3:
1038             hChargedEM->Fill(energy ) ;
1039             break ;
1040           case 4:
1041             hPbarEM->Fill(energy ) ;
1042             break ;
1043           default:
1044             break ;
1045           }
1046         }
1047
1048         if((recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)||
1049            (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALHAFAST) ||
1050            (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW) ||
1051            (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALHASLOW) ) //nuetral
1052           hVeto->Fill(energy ) ;
1053
1054         //fill number of primaries identified as ...
1055         if(primaryCode >= 0) // Primary code defined
1056           counter[recParticle->GetType()][primaryCode]++ ;
1057
1058       }
1059
1060     } // no closest primary found
1061   }
1062
1063
1064   //===================  SaveHistograms
1065   cfile->cd() ;
1066   hPrimary->Write(0,kOverwrite);
1067   hAllRP->Write(0,kOverwrite);
1068   hPhot->Write(0,kOverwrite);
1069   hShape->Write(0,kOverwrite);
1070   hVeto->Write(0,kOverwrite);
1071   hPhotReg->Write(0,kOverwrite);
1072   hPhotEM->Write(0,kOverwrite);
1073   hNReg ->Write(0,kOverwrite);
1074   hNEM  ->Write(0,kOverwrite);
1075   hNBarReg ->Write(0,kOverwrite);
1076   hNBarEM  ->Write(0,kOverwrite);
1077   hChargedReg ->Write(0,kOverwrite);
1078   hChargedEM  ->Write(0,kOverwrite);
1079   hPbarReg ->Write(0,kOverwrite);
1080   hPbarEM  ->Write(0,kOverwrite);
1081
1082   cfile->Write(0,kOverwrite);
1083   cfile->Close();
1084   delete cfile ;
1085
1086
1087   //print Final Table
1088   maxevent = (Int_t)gAlice->TreeE()->GetEntries() ;
1089  //  cout << "Resolutions: Analyzed " << gime->MaxEvent() << " event(s)" << endl ;
1090
1091   cout << "Resolutions: Analyzed " << maxevent << " event(s)" << endl ;
1092   cout << endl ;
1093
1094   cout << "        Primary:    Photon  Neutron  Antineutron  Charged hadron   AntiProton" << endl ;
1095   cout << "--------------------------------------------------------------------------------" << endl ;
1096   cout << "         kGAMMA: "
1097        << setw(8) << counter[2][0] << setw(9)  << counter[2][1] << setw(13) << counter[2][2]
1098        << setw(15)<< counter[2][3] << setw(13) << counter[2][4] << endl ;
1099   cout << "       kGAMMAHA: "
1100        << setw(8) << counter[3][0] << setw(9)  << counter[3][1] << setw(13) << counter[3][2]
1101        << setw(15)<< counter[3][3] << setw(13) << counter[3][4] << endl ;
1102   cout << "     kNEUTRALEM: "
1103        << setw(8) << counter[0][0] << setw(9)  << counter[0][1] << setw(13) << counter[0][2]
1104        << setw(15)<< counter[0][3] << setw(13) << counter[0][4] << endl ;
1105   cout << "     kNEUTRALHA: "
1106        << setw(8) << counter[1][0] << setw(9)  << counter[1][1] << setw(13) << counter[1][2]
1107        << setw(15)<< counter[1][3] << setw(13) << counter[1][4] << endl ;
1108   cout << "      kABSURDEM: "
1109        << setw(8) << counter[4][0] << setw(9)  << counter[4][1] << setw(13) << counter[4][2]
1110        << setw(15)<< counter[4][3] << setw(13) << counter[4][4] << endl ;
1111   cout << "      kABSURDHA: "
1112        << setw(8) << counter[5][0] << setw(9)  << counter[5][1] << setw(13) << counter[5][2]
1113        << setw(15)<< counter[5][3] << setw(13) << counter[5][4] << endl ;
1114   cout << "      kELECTRON: "
1115        << setw(8) << counter[6][0] << setw(9)  << counter[6][1] << setw(13) << counter[6][2]
1116        << setw(15)<< counter[6][3] << setw(13) << counter[6][4] << endl ;
1117   cout << "     kCHARGEDHA: "
1118        << setw(8) << counter[7][0] << setw(9)  << counter[7][1] << setw(13) << counter[7][2]
1119        << setw(15)<< counter[7][3] << setw(13) << counter[7][4] << endl ;
1120   cout << "--------------------------------------------------------------------------------" << endl ;
1121
1122   Int_t totalInd = 0 ;
1123   for(i1 = 0; i1<8; i1++)
1124     for(i2 = 0; i2<5; i2++)
1125       totalInd+=counter[i1][i2] ;
1126   cout << "Indentified particles: " << totalInd << endl ;
1127
1128 }
1129
1130
1131