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 "AliPHOSPpsdRecPoint.h"
  93 #include "AliPHOSGetter.h"
  94
  95
  96 ClassImp(AliPHOSAnalyze)
  97
  98 //____________________________________________________________________________
  99   AliPHOSAnalyze::AliPHOSAnalyze()
 100 {
 101   // default ctor (useless)
 102   fCorrection = 1.2 ;  //Value calculated for default parameters of reconstruction
 103 }
 104
 105 //____________________________________________________________________________
 106 AliPHOSAnalyze::AliPHOSAnalyze(Text_t * fileName)
 107 {
 108   // ctor: analyze events from root file "name"
 109   ffileName = fileName ;
 110   fCorrection = 1.05 ;  //Value calculated for default parameters of reconstruction
 111 }
 112
 113 //____________________________________________________________________________
 114 AliPHOSAnalyze::AliPHOSAnalyze(const AliPHOSAnalyze & ana)
 115 {
 116   // copy ctor
 117   ( (AliPHOSAnalyze &)ana ).Copy(*this) ;
 118 }
 119
 120 //____________________________________________________________________________
 121 AliPHOSAnalyze::~AliPHOSAnalyze()
 122 {
 123   // dtor
 124
 125 }
 126 //____________________________________________________________________________
 127 void AliPHOSAnalyze::DrawRecon(Int_t Nevent,Int_t Nmod,const char * branchName,const char* branchTitle){
 128   //Draws pimary particles and reconstructed
 129   //digits, RecPoints, RecPartices etc
 130   //for event Nevent in the module Nmod.
 131
 132   TH2F * digitOccupancy  = new TH2F("digitOccupancy","EMC digits", 64,-71.,71.,64,-71.,71.);
 133   TH2F * sdigitOccupancy = new TH2F("sdigitOccupancy","EMC sdigits", 64,-71.,71.,64,-71.,71.);
 134   TH2F * emcOccupancy    = new TH2F("emcOccupancy","EMC RecPoints",64,-71.,71.,64,-71.,71.);
 135   TH2F * ppsdUp          = new TH2F("ppsdUp","PPSD Up digits",     128,-71.,71.,128,-71.,71.) ;
 136   TH2F * ppsdUpCl        = new TH2F("ppsdUpCl","PPSD Up RecPoints",128,-71.,71.,128,-71.,71.) ;
 137   TH2F * ppsdLow         = new TH2F("ppsdLow","PPSD Low digits",     128,-71.,71.,128,-71.,71.) ;
 138   TH2F * ppsdLowCl       = new TH2F("ppsdLowCl","PPSD Low RecPoints",128,-71.,71.,128,-71.,71.) ;
 139   TH2F * nbar            = new TH2F("nbar","Primary nbar",    64,-71.,71.,64,-71.,71.);
 140   TH2F * phot            = new TH2F("phot","Primary Photon",  64,-71.,71.,64,-71.,71.);
 141   TH2F * charg           = new TH2F("charg","Primary charged",64,-71.,71.,64,-71.,71.);
 142   TH2F * recPhot         = new TH2F("recPhot","RecParticles with primary Photon",64,-71.,71.,64,-71.,71.);
 143   TH2F * recNbar         = new TH2F("recNbar","RecParticles with primary Nbar",  64,-71.,71.,64,-71.,71.);
 144
 145   //========== Create ObjectGetter
 146   AliPHOSGetter * gime = AliPHOSGetter::GetInstance(ffileName.Data(),branchTitle) ;
 147   const AliPHOSGeometry * phosgeom = gime->PHOSGeometry() ;
 148   gime->Event(Nevent);
 149
 150   //Plot Primary Particles
 151   const TParticle * primary ;
 152   Int_t iPrimary ;
 153   for ( iPrimary = 0 ; iPrimary < gime->NPrimaries() ; iPrimary++)
 154     {
 155       primary = gime->Primary(iPrimary) ;
 156       Int_t primaryType = primary->GetPdgCode() ;
 157       if( (primaryType == 211)||(primaryType == -211)||(primaryType == 2212)||(primaryType == -2212) ) {
 158         Int_t moduleNumber ;
 159         Double_t primX, primZ ;
 160         phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
 161         if(moduleNumber==Nmod)
 162           charg->Fill(primZ,primX,primary->Energy()) ;
 163       }
 164       if( primaryType == 22 ) {
 165         Int_t moduleNumber ;
 166         Double_t primX, primZ ;
 167         phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
 168         if(moduleNumber==Nmod)
 169           phot->Fill(primZ,primX,primary->Energy()) ;
 170       }
 171       else{
 172         if( primaryType == -2112 ) {
 173           Int_t moduleNumber ;
 174           Double_t primX, primZ ;
 175           phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
 176           if(moduleNumber==Nmod)
 177             nbar->Fill(primZ,primX,primary->Energy()) ;
 178         }
 179       }
 180     }
 181
 182   Int_t iSDigit ;
 183   const AliPHOSDigit * sdigit ;
 184
 185   for(iSDigit = 0; iSDigit < gime->NSDigits(); iSDigit++)
 186       {
 187         sdigit = gime->SDigit(iSDigit) ;
 188         Int_t relid[4];
 189         phosgeom->AbsToRelNumbering(sdigit->GetId(), relid) ;
 190         Float_t x,z ;
 191         phosgeom->RelPosInModule(relid,x,z) ;
 192         Float_t e = gime->SDigitizer()->Calibrate(sdigit->GetAmp()) ;
 193         if(relid[0]==Nmod){
 194           if(relid[1]==0)  //EMC
 195             sdigitOccupancy->Fill(x,z,e) ;
 196           if((relid[1]>0)&&(relid[1]<17))
 197             ppsdUp->Fill(x,z,e) ;
 198           if(relid[1]>16)
 199             ppsdLow->Fill(x,z,e) ;
 200         }
 201       }
 202
 203   //Plot digits
 204   Int_t iDigit ;
 205   const AliPHOSDigit * digit ;
 206   for(iDigit = 0; iDigit < gime->NDigits(); iDigit++)
 207     {
 208       digit = gime->Digit(iDigit) ;
 209       Int_t relid[4];
 210       phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
 211       Float_t x,z ;
 212       phosgeom->RelPosInModule(relid,x,z) ;
 213       Float_t e = gime->SDigitizer()->Calibrate(digit->GetAmp()) ;
 214       if(relid[0]==Nmod){
 215         if(relid[1]==0)  //EMC
 216           digitOccupancy->Fill(x,z,e) ;
 217         if((relid[1]>0)&&(relid[1]<17))
 218           ppsdUp->Fill(x,z,e) ;
 219         if(relid[1]>16)
 220           ppsdLow->Fill(x,z,e) ;
 221       }
 222     }
 223
 224
 225   //Plot RecPoints
 226   Int_t irecp ;
 227   TVector3 pos ;
 228
 229   for(irecp = 0; irecp < gime->NEmcRecPoints() ; irecp ++){
 230     const AliPHOSEmcRecPoint * emc= gime->EmcRecPoint(irecp) ;
 231     if(emc->GetPHOSMod()==Nmod){
 232       emc->GetLocalPosition(pos) ;
 233       emcOccupancy->Fill(pos.X(),pos.Z(),emc->GetEnergy());
 234     }
 235   }
 236
 237
 238   for(irecp = 0; irecp < gime->NCpvRecPoints() ; irecp ++){
 239     const AliPHOSRecPoint * cpv = gime->CpvRecPoint(irecp) ;
 240     if((strcmp(cpv->ClassName(),"AliPHOSPpsdRecPoint" )) == 0){  // PPSD Rec Point
 241       AliPHOSPpsdRecPoint * ppsd = (AliPHOSPpsdRecPoint*) cpv ;
 242       ppsd->GetLocalPosition(pos) ;
 243       if(ppsd->GetPHOSMod()==Nmod){
 244         ppsd->GetLocalPosition(pos) ;
 245         if(ppsd->GetUp())
 246           ppsdUpCl->Fill(pos.X(),pos.Z(),ppsd->GetEnergy());
 247         else
 248           ppsdLowCl->Fill(pos.X(),pos.Z(),ppsd->GetEnergy());
 249       }
 250     }
 251     else{
 252       AliPHOSCpvRecPoint * cpv1 = (AliPHOSCpvRecPoint*) cpv ;
 253       if(cpv1->GetPHOSMod()==Nmod){
 254         cpv1->GetLocalPosition(pos) ;
 255         ppsdUpCl->Fill(pos.X(),pos.Z(),cpv1->GetEnergy());
 256       }
 257     }
 258   }
 259
 260
 261   //Plot RecParticles
 262   const AliPHOSRecParticle * recParticle ;
 263   Int_t iRecParticle ;
 264   for(iRecParticle = 0; iRecParticle < gime->NRecParticles() ;iRecParticle++ )
 265     {
 266       recParticle =  gime->RecParticle(iRecParticle) ;
 267       Int_t moduleNumberRec ;
 268       Double_t recX, recZ ;
 269       phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
 270       if(moduleNumberRec == Nmod){
 271
 272         Double_t minDistance = 5. ;
 273         Int_t closestPrimary = -1 ;
 274
 275
 276         //extract list of primaries: it is stored at EMC RecPoints
 277         Int_t emcIndex = gime->TrackSegment(recParticle->GetPHOSTSIndex())->GetEmcIndex() ;
 278         Int_t numberofprimaries ;
 279         Int_t * listofprimaries  = gime->EmcRecPoint(emcIndex)->GetPrimaries(numberofprimaries)  ;
 280         Int_t index ;
 281         const TParticle * primary ;
 282         Double_t distance = minDistance ;
 283
 284         for ( index = 0 ; index < numberofprimaries ; index++){
 285           primary = gime->Primary(listofprimaries[index]) ;
 286           Int_t moduleNumber ;
 287           Double_t primX, primZ ;
 288           phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
 289           if(moduleNumberRec == moduleNumber)
 290             distance = TMath::Sqrt((recX-primX)*(recX-primX)+(recZ-primZ)*(recZ-primZ) ) ;
 291           if(minDistance > distance)
 292             {
 293               minDistance = distance ;
 294               closestPrimary = listofprimaries[index] ;
 295             }
 296         }
 297
 298         if(closestPrimary >=0 ){
 299
 300           Int_t primaryType = gime->Primary(closestPrimary)->GetPdgCode() ;
 301
 302           if(primaryType==22)
 303             recPhot->Fill(recZ,recX,recParticle->Energy()) ;
 304           else
 305             if(primaryType==-2112)
 306               recNbar->Fill(recZ,recX,recParticle->Energy()) ;
 307         }
 308       }
 309     }
 310
 311
 312   //Plot made histograms
 313   digitOccupancy->Draw("box") ;
 314   sdigitOccupancy->SetLineColor(5) ;
 315   sdigitOccupancy->Draw("box") ;
 316   emcOccupancy->SetLineColor(2) ;
 317   emcOccupancy->Draw("boxsame") ;
 318   ppsdUp->SetLineColor(3) ;
 319   ppsdUp->Draw("boxsame") ;
 320   ppsdLow->SetLineColor(4) ;
 321   ppsdLow->Draw("boxsame") ;
 322   phot->SetLineColor(8) ;
 323   phot->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);
 417
 418     //copy EM RecParticles to the "total" list
 419     const AliPHOSRecParticle * recParticle ;
 420     Int_t iRecParticle ;
 421     for(iRecParticle = 0; iRecParticle < gime->NRecParticles()  ;iRecParticle++ )
 422       {
 423         recParticle = gime->RecParticle(iRecParticle) ;
 424         if((recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)||
 425            (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM))
 426           new( (*allRecParticleList)[iRecPhot++] ) AliPHOSRecParticle(*recParticle) ;
 427       }
 428
 429     Int_t mevent = event%nMixedEvents ; //event number in the "mixed" cicle
 430     nRecParticles[mevent] = iRecPhot-1 ;
 431
 432     //check, if it is time to calculate invariant mass?
 433     Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ;
 434     if((mevent == 0) && (event +1 == maxevent)){
 435
 436     //   if((mevent == 0) && (event +1 == gime->MaxEvent())){
 437
 438       //calculate invariant mass:
 439       Int_t irp1,irp2 ;
 440       Int_t nCurEvent = 0 ;
 441
 442       for(irp1 = 0; irp1 < allRecParticleList->GetEntries()-1; irp1++){
 443         AliPHOSRecParticle * rp1 = (AliPHOSRecParticle *)allRecParticleList->At(irp1) ;
 444
 445         for(irp2 = irp1+1; irp2 < allRecParticleList->GetEntries(); irp2++){
 446           AliPHOSRecParticle * rp2 = (AliPHOSRecParticle *)allRecParticleList->At(irp2) ;
 447
 448           Double_t invMass ;
 449           invMass = (rp1->Energy()+rp2->Energy())*(rp1->Energy()+rp2->Energy())-
 450             (rp1->Px()+rp2->Px())*(rp1->Px()+rp2->Px())-
 451             (rp1->Py()+rp2->Py())*(rp1->Py()+rp2->Py())-
 452             (rp1->Pz()+rp2->Pz())*(rp1->Pz()+rp2->Pz()) ;
 453
 454           if(invMass> 0)
 455             invMass = TMath::Sqrt(invMass);
 456
 457           Double_t pt ;
 458           pt = TMath::Sqrt((rp1->Px()+rp2->Px() )*( rp1->Px()+rp2->Px() ) +
 459                            (rp1->Py()+rp2->Py() )*( rp1->Py()+rp2->Py() ) );
 460
 461           if(irp1 > nRecParticles[nCurEvent])
 462             nCurEvent++;
 463
 464           if(irp2 <= nRecParticles[nCurEvent]){ //'Real' event
 465             hRealEM->Fill(invMass,pt);
 466             if((rp1->GetType() == AliPHOSFastRecParticle::kGAMMA)&&
 467                (rp2->GetType() == AliPHOSFastRecParticle::kGAMMA) )
 468               hRealPhot->Fill(invMass,pt);
 469           }
 470           else{
 471             hMixedEM->Fill(invMass,pt);
 472             if((rp1->GetType() == AliPHOSFastRecParticle::kGAMMA)&&
 473                (rp2->GetType() == AliPHOSFastRecParticle::kGAMMA) )
 474               hMixedPhot->Fill(invMass,pt);
 475           } //real-mixed
 476
 477       } //loop over second rp
 478       }//loop over first rp
 479
 480       //Make some cleanings
 481       for(Int_t index = 0; index < nMixedEvents; index ++)
 482         nRecParticles[index] = 0 ;
 483       iRecPhot = 0 ;
 484       allRecParticleList->Clear() ;
 485
 486     }
 487   }
 488   delete allRecParticleList ;
 489
 490   //writing output
 491   mfile->cd();
 492
 493   hRealEM->Write(0,kOverwrite) ;
 494   hRealPhot->Write(0,kOverwrite) ;
 495   hMixedEM->Write(0,kOverwrite) ;
 496   hMixedPhot->Write(0,kOverwrite) ;
 497
 498   mfile->Write();
 499   mfile->Close();
 500   delete mfile ;
 501   delete nRecParticles;
 502
 503 }
 504
 505 //____________________________________________________________________________
 506  void AliPHOSAnalyze::EnergyResolution(const char * branchTitle)
 507 {
 508   //fills two dimentional histo: energy of primary vs. energy of reconstructed
 509
 510   TH2F * hAllEnergy = 0 ;  //all reconstructed with primary photon
 511   TH2F * hPhotEnergy= 0 ;  //kGamma  with primary photon
 512   TH2F * hEMEnergy  = 0 ;  //electromagnetic with primary photon
 513
 514   //opening file and reading histograms if any
 515   TFile * efile = new TFile("energy.root","update");
 516
 517   hAllEnergy = (TH2F*)efile->Get("hAllEnergy") ;
 518   if(hAllEnergy == 0)
 519     hAllEnergy  = new TH2F("hAllEnergy",  "Energy of any RP with primary photon",100, 0., 5., 100, 0., 5.);
 520
 521   hPhotEnergy =(TH2F*) efile->Get("hPhotEnergy") ;
 522   if(hPhotEnergy == 0)
 523     hPhotEnergy = new TH2F("hPhotEnergy", "Energy of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.);
 524
 525   hEMEnergy =(TH2F*) efile->Get("hEMEnergy");
 526   if(hEMEnergy == 0)
 527     hEMEnergy   = new TH2F("hEMEnergy",   "Energy of EM with primary photon",    100, 0., 5., 100, 0., 5.);
 528
 529
 530   AliPHOSGetter * gime = AliPHOSGetter::GetInstance(ffileName.Data(),branchTitle) ;
 531   const AliPHOSGeometry * phosgeom = gime->PHOSGeometry() ;
 532
 533   Int_t ievent;
 534   Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ;
 535   //  for ( ievent=0; ievent < gime->MaxEvent() ; ievent++){
 536   for ( ievent=0; ievent < maxevent ; ievent++){
 537
 538     //read the current event
 539     gime->Event(ievent) ;
 540
 541     const AliPHOSRecParticle * recParticle ;
 542     Int_t iRecParticle ;
 543     for(iRecParticle = 0; iRecParticle < gime->NRecParticles() ;iRecParticle++ ){
 544       recParticle = gime->RecParticle(iRecParticle) ;
 545
 546       //find the closest primary
 547       Int_t moduleNumberRec ;
 548       Double_t recX, recZ ;
 549       phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
 550
 551       Double_t minDistance  = 100. ;
 552       Int_t closestPrimary = -1 ;
 553
 554       //extract list of primaries: it is stored at EMC RecPoints
 555       Int_t emcIndex = gime->TrackSegment(recParticle->GetPHOSTSIndex())->GetEmcIndex() ;
 556       Int_t numberofprimaries ;
 557       Int_t * listofprimaries  = gime->EmcRecPoint(emcIndex)->GetPrimaries(numberofprimaries)  ;
 558
 559       Int_t index ;
 560       const TParticle * primary ;
 561       Double_t distance = minDistance ;
 562       Double_t dX, dZ;
 563       Double_t dXmin = 0.;
 564       Double_t dZmin = 0. ;
 565       for ( index = 0 ; index < numberofprimaries ; index++){
 566         primary = gime->Primary(listofprimaries[index]) ;
 567         Int_t moduleNumber ;
 568         Double_t primX, primZ ;
 569         phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
 570         if(moduleNumberRec == moduleNumber) {
 571           dX = recX - primX;
 572           dZ = recZ - primZ;
 573           distance = TMath::Sqrt(dX*dX + dZ*dZ) ;
 574           if(minDistance > distance) {
 575             minDistance = distance ;
 576             dXmin = dX;
 577             dZmin = dZ;
 578             closestPrimary = listofprimaries[index] ;
 579           }
 580         }
 581       }
 582
 583       //if found primary, fill histograms
 584       if(closestPrimary >=0 ){
 585         const TParticle * primary = gime->Primary(closestPrimary) ;
 586         if(primary->GetPdgCode() == 22){
 587           hAllEnergy->Fill(primary->Energy(), recParticle->Energy()) ;
 588           if(recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA){
 589             hPhotEnergy->Fill(primary->Energy(), recParticle->Energy() ) ;
 590             hEMEnergy->Fill(primary->Energy(), recParticle->Energy() ) ;
 591           }
 592           else
 593             if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM)
 594               hEMEnergy->Fill(primary->Energy(), recParticle->Energy() ) ;
 595         }
 596       }
 597     }
 598   }
 599
 600   //write filled histograms
 601   efile->cd() ;
 602   hAllEnergy->Write(0,kOverwrite) ;
 603   hPhotEnergy->Write(0,kOverwrite) ;
 604   hEMEnergy->Write(0,kOverwrite)  ;
 605   //  efile->Write() ;
 606   efile->Close() ;
 607   delete efile ;
 608
 609 }
 610 //____________________________________________________________________________
 611 void AliPHOSAnalyze::PositionResolution(const char * branchTitle)
 612 {
 613   //fills two dimentional histo: energy vs. primary - reconstructed distance
 614
 615
 616
 617   TH2F * hAllPosition  = 0;     // Position of any RP with primary photon
 618   TH2F * hPhotPosition = 0;    // Position of kGAMMA with primary photon
 619   TH2F * hEMPosition   = 0;      // Position of EM with primary photon
 620
 621   TH1F * hAllPositionX = 0;    // X-Position Resolution of photons with photon primary
 622   TH1F * hAllPositionZ = 0;    // Z-Position Resolution of photons with photon primary
 623
 624
 625   //opening file and reading histograms if any
 626   TFile * pfile = new TFile("position.root","update");
 627
 628   hAllPosition = (TH2F*)pfile->Get("hAllPosition");
 629   if(hAllPosition == 0)
 630     hAllPosition  = new TH2F("hAllPosition",
 631                              "Position of any RP with primary photon",100, 0., 5., 100, 0., 5.);
 632   hPhotPosition= (TH2F*)pfile->Get("hPhotPosition");
 633   if(hPhotPosition == 0)
 634     hPhotPosition = new TH2F("hPhotPosition",
 635                              "Position of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.);
 636   hEMPosition= (TH2F*)pfile->Get("hEMPosition") ;
 637   if(hEMPosition == 0)
 638     hEMPosition   = new TH2F("hEMPosition",
 639                              "Position of EM with primary photon",    100, 0., 5., 100, 0., 5.);
 640   hAllPositionX = (TH1F*)pfile->Get("hAllPositionX") ;
 641   if(hAllPositionX == 0)
 642     hAllPositionX = new TH1F("hAllPositionX",
 643                              "Delta X of any RP with primary photon",100, -2., 2.);
 644   hAllPositionZ =(TH1F*) pfile->Get("hAllPositionZ") ;
 645   if(hAllPositionZ == 0)
 646     hAllPositionZ = new TH1F("hAllPositionZ",
 647                              "Delta X of any RP with primary photon",100, -2., 2.);
 648
 649
 650   AliPHOSGetter * gime = AliPHOSGetter::GetInstance(ffileName.Data(),branchTitle) ;
 651   const AliPHOSGeometry * phosgeom = gime->PHOSGeometry() ;
 652
 653   Int_t ievent;
 654  Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ;
 655   //  for ( ievent=0; ievent < gime->MaxEvent() ; ievent++){
 656   for ( ievent=0; ievent < maxevent ; ievent++){
 657
 658     //read the current event
 659     gime->Event(ievent) ;
 660
 661     const AliPHOSRecParticle * recParticle ;
 662     Int_t iRecParticle ;
 663     for(iRecParticle = 0; iRecParticle < gime->NRecParticles() ;iRecParticle++ ){
 664       recParticle = gime->RecParticle(iRecParticle) ;
 665
 666       //find the closest primary
 667       Int_t moduleNumberRec ;
 668       Double_t recX, recZ ;
 669       phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
 670
 671       Double_t minDistance  = 100. ;
 672       Int_t closestPrimary = -1 ;
 673
 674       //extract list of primaries: it is stored at EMC RecPoints
 675       Int_t emcIndex = gime->TrackSegment(recParticle->GetPHOSTSIndex())->GetEmcIndex() ;
 676       Int_t numberofprimaries ;
 677       Int_t * listofprimaries  = gime->EmcRecPoint(emcIndex)->GetPrimaries(numberofprimaries)  ;
 678
 679       Int_t index ;
 680       const TParticle * primary ;
 681       Double_t distance = minDistance ;
 682       Double_t dX = 1000; // incredible number
 683       Double_t dZ = 1000; // for the case if no primary will be found
 684       Double_t dXmin = 0.;
 685       Double_t dZmin = 0. ;
 686       for ( index = 0 ; index < numberofprimaries ; index++){
 687         primary = gime->Primary(listofprimaries[index]) ;
 688         Int_t moduleNumber ;
 689         Double_t primX, primZ ;
 690         phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
 691         if(moduleNumberRec == moduleNumber) {
 692           dX = recX - primX;
 693           dZ = recZ - primZ;
 694           distance = TMath::Sqrt(dX*dX + dZ*dZ) ;
 695           if(minDistance > distance) {
 696             minDistance = distance ;
 697             dXmin = dX;
 698             dZmin = dZ;
 699             closestPrimary = listofprimaries[index] ;
 700           }
 701         }
 702       }
 703
 704       //if found primary, fill histograms
 705       if(closestPrimary >=0 ){
 706         const TParticle * primary = gime->Primary(closestPrimary) ;
 707         if(primary->GetPdgCode() == 22){
 708           hAllPosition->Fill(primary->Energy(), minDistance) ;
 709           hAllPositionX->Fill(primary->Energy(), dX) ;
 710           hAllPositionZ->Fill(primary->Energy(), dZ) ;
 711           if(recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA){
 712             hPhotPosition->Fill(primary->Energy(), minDistance ) ;
 713             hEMPosition->Fill(primary->Energy(), minDistance ) ;
 714           }
 715           else
 716             if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM)
 717               hEMPosition->Fill(primary->Energy(), minDistance ) ;
 718         }
 719       }
 720     }
 721   }
 722
 723   //Write output histgrams
 724   pfile->cd() ;
 725   hAllPosition->Write(0,kOverwrite) ;
 726   hAllPositionX->Write(0,kOverwrite) ;
 727   hAllPositionZ->Write(0,kOverwrite) ;
 728   hPhotPosition->Write(0,kOverwrite) ;
 729   hEMPosition->Write(0,kOverwrite) ;
 730   pfile->Write() ;
 731   pfile->Close() ;
 732   delete pfile ;
 733
 734
 735 }
 736 //____________________________________________________________________________
 737 void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){
 738 // fills spectra of primary photons and several kinds of
 739 // reconstructed particles, so that analyzing them one can
 740 // estimate conatmination, efficiency of registration etc.
 741
 742   //define several general histograms
 743   TH1F * hPrimary = 0;   //spectrum (P_t distribution) of primary photons
 744   TH1F * hAllRP   = 0;   //spectrum of all RecParticles in PHOS
 745   TH1F * hPhot    = 0;   //spectrum of kGAMMA RecParticles
 746   TH1F * hPPSD    = 0;   //spectrum of all RecParticles with PPSD signal
 747   TH1F * hShape   = 0;   //spectrum of all EM RecParticles
 748   TH1F * hVeto    = 0;   //spectrum of all neutral RecParticles
 749
 750   //Now separate histograms in accoradance with primary
 751   //primary - photon
 752   TH1F * hPhotReg = 0;   //Registeres as photon
 753   TH1F * hPhotEM  = 0;   //Registered as EM
 754   TH1F * hPhotPPSD= 0;   //Registered as RecParticle with PPSD signal
 755
 756   //primary - n
 757   TH1F * hNReg = 0;   //Registeres as photon
 758   TH1F * hNEM  = 0;   //Registered as EM
 759   TH1F * hNPPSD= 0;   //Registered as RecParticle with PPSD signal
 760
 761   //primary - nBar
 762   TH1F * hNBarReg = 0;   //Registeres as photon
 763   TH1F * hNBarEM  = 0;   //Registered as EM
 764   TH1F * hNBarPPSD= 0;   //Registered as RecParticle with PPSD signal
 765
 766   //primary - charged hadron (pBar excluded)
 767   TH1F * hChargedReg = 0;   //Registeres as photon
 768   TH1F * hChargedEM  = 0;   //Registered as EM
 769   TH1F * hChargedPPSD= 0;   //Registered as RecParticle with PPSD signal
 770
 771   //primary - pBar
 772   TH1F * hPbarReg = 0;   //Registeres as photon
 773   TH1F * hPbarEM  = 0;   //Registered as EM
 774   TH1F * hPbarPPSD= 0;   //Registered as RecParticle with PPSD signal
 775
 776
 777   //Reading histograms from the file
 778   TFile * cfile = new TFile("contamination.root","update") ;
 779
 780   //read general histograms
 781   hPrimary = (TH1F*) cfile->Get("hPrimary") ;
 782   if(hPrimary == 0)
 783     hPrimary= new TH1F("hPrimary", "Primary photon spectrum", 100, 0., 5.);
 784   hAllRP = (TH1F*)cfile->Get("hAllRP") ;
 785   if(hAllRP == 0)
 786     hAllRP = new TH1F("hAllRP","All Reconstructed particles", 100, 0., 5.);
 787   hPhot  = (TH1F*)cfile->Get("hPhot") ;
 788   if(hPhot == 0)
 789     hPhot = new TH1F("hPhot","All kGAMMA RecParticles",100, 0., 5.);
 790   hPPSD = (TH1F*)cfile->Get("hPPSD") ;
 791   if(hPPSD == 0)
 792     hPPSD  = new TH1F("hPPSD", "All PPSD Recparticles",    100, 0., 5.);
 793   hShape = (TH1F*) cfile->Get("hShape") ;
 794   if(hShape == 0)
 795     hShape = new TH1F("hShape","All particles with EM shower",100, 0., 5.);
 796   hVeto= (TH1F*)cfile->Get("hVeto") ;
 797   if(hVeto == 0)
 798     hVeto  = new TH1F("hVeto", "All uncharged particles",      100, 0., 5.);
 799
 800
 801   //primary - photon
 802   hPhotReg = (TH1F*)cfile->Get("hPhotReg");
 803   if(hPhotReg == 0)
 804     hPhotReg   = new TH1F("hPhotReg","Photon registered as photon",100, 0., 5.);
 805   hPhotEM  =(TH1F*)cfile->Get("hPhotEM");
 806   if(hPhotEM== 0)
 807     hPhotEM   = new TH1F("hPhotEM",  "Photon registered as EM", 100, 0., 5.);
 808   hPhotPPSD= (TH1F*)cfile->Get("hPhotPPSD");
 809   if(hPhotPPSD== 0)
 810     hPhotPPSD   = new TH1F("hPhotPPSD","Photon registered as PPSD", 100, 0., 5.);
 811
 812   //primary - n
 813   hNReg = (TH1F*)cfile->Get("hNReg");
 814   if(hNReg== 0)
 815    hNReg      = new TH1F("hNReg",   "N registered as photon",              100, 0., 5.);
 816   hNEM  = (TH1F*)cfile->Get("hNEM");
 817   if(hNEM== 0)
 818     hNEM      = new TH1F("hNEM",    "N registered as EM",      100, 0., 5.);
 819   hNPPSD=(TH1F*)cfile->Get("hNPPSD");
 820   if(hNPPSD== 0)
 821     hNPPSD      = new TH1F("hNPPSD","N registered as PPSD",      100, 0., 5.);
 822
 823   //primary - nBar
 824   hNBarReg =(TH1F*)cfile->Get("hNBarReg");
 825   if(hNBarReg== 0)
 826    hNBarReg   = new TH1F("hNBarReg", "NBar registered as photon",           100, 0., 5.);
 827   hNBarEM  =(TH1F*)cfile->Get("hNBarEM");
 828   if(hNBarEM== 0)
 829     hNBarEM   = new TH1F("hNBarEM",  "NBar registered as EM",   100, 0., 5.);
 830   hNBarPPSD=(TH1F*)cfile->Get("hNBarPPSD");
 831   if(hNBarPPSD== 0)
 832     hNBarPPSD   = new TH1F("hNBarPPSD","NBar registered as PPSD",   100, 0., 5.);
 833
 834   //primary - charged hadron (pBar excluded)
 835   hChargedReg = (TH1F*)cfile->Get("hChargedReg");
 836   if(hChargedReg== 0)
 837     hChargedReg= new TH1F("hChargedReg", "Charged hadron registered as photon",100, 0., 5.);
 838   hChargedEM  = (TH1F*)cfile->Get("hChargedEM");
 839   if(hChargedEM== 0)
 840     hChargedEM= new TH1F("hChargedEM","Charged registered as EM",100, 0., 5.);
 841   hChargedPPSD= (TH1F*)cfile->Get("hChargedPPSD");
 842   if(hChargedPPSD== 0)
 843     hChargedPPSD= new TH1F("hChargedPPSD","Charged registered as PPSD",100, 0., 5.);
 844
 845   //primary - pBar
 846   hPbarReg = (TH1F*)cfile->Get("hPbarReg");
 847   if(hPbarReg== 0)
 848     hPbarReg= new TH1F("hPbarReg", "pBar registered as photon",100, 0., 5.);
 849   hPbarEM  = (TH1F*)cfile->Get("hPbarEM");
 850   if(hPbarEM== 0)
 851     hPbarEM= new TH1F("hPbarEM","Pbar registered as EM",100, 0., 5.);
 852   hPbarPPSD= (TH1F*)cfile->Get("hPbarPPSD");
 853   if(hPbarPPSD== 0)
 854     hPbarPPSD= new TH1F("hPbarPPSD","Pbar as PPSD",100, 0., 5.);
 855
 856
 857   //Now make some initializations
 858
 859   Int_t counter[8][5] ;      //# of registered particles
 860   Int_t i1,i2 ;
 861   for(i1 = 0; i1<8; i1++)
 862     for(i2 = 0; i2<5; i2++)
 863       counter[i1][i2] = 0 ;
 864
 865
 866
 867   AliPHOSGetter * gime = AliPHOSGetter::GetInstance(ffileName.Data(),RecPointsTitle) ;
 868   const AliPHOSGeometry * phosgeom = gime->PHOSGeometry() ;
 869
 870   Int_t ievent;
 871   Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ;
 872   //  for ( ievent=0; ievent < gime->MaxEvent() ; ievent++){
 873   for ( ievent=0; ievent < maxevent ; ievent++){
 874
 875     gime->Event(ievent) ;
 876
 877     //=========== Make spectrum of the primary photons
 878     const TParticle * primary ;
 879     Int_t iPrimary ;
 880     for( iPrimary = 0 ; iPrimary < gime->NPrimaries() ; iPrimary++){
 881       primary = gime->Primary(iPrimary) ;
 882       Int_t primaryType = primary->GetPdgCode() ;
 883       if( primaryType == 22 ) {
 884         //check, if photons folls onto PHOS
 885         Int_t moduleNumber ;
 886         Double_t primX, primZ ;
 887         phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
 888         if(moduleNumber)
 889           hPrimary->Fill(primary->Energy()) ;
 890
 891       }
 892
 893     }
 894     //========== Now scan over RecParticles
 895     const AliPHOSRecParticle * recParticle ;
 896     Int_t iRecParticle ;
 897     for(iRecParticle = 0; iRecParticle < gime->NRecParticles(); iRecParticle++ ){
 898       recParticle = gime->RecParticle(iRecParticle) ;
 899       //fill histo spectrum of all RecParticles
 900       hAllRP->Fill(CorrectedEnergy(recParticle->Energy())) ;
 901
 902       //==========find the closest primary
 903       Int_t moduleNumberRec ;
 904       Double_t recX, recZ ;
 905       phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
 906
 907       Double_t minDistance  = 100. ;
 908       Int_t closestPrimary = -1 ;
 909
 910       //extract list of primaries: it is stored at EMC RecPoints
 911       Int_t emcIndex = gime->TrackSegment(recParticle->GetPHOSTSIndex())->GetEmcIndex() ;
 912       Int_t numberofprimaries ;
 913       Int_t * listofprimaries  = gime->EmcRecPoint(emcIndex)->GetPrimaries(numberofprimaries)  ;
 914       Int_t index ;
 915       const TParticle * primary ;
 916       Double_t distance = minDistance ;
 917       Double_t dX, dZ;
 918       Double_t dXmin = 0.;
 919       Double_t dZmin = 0. ;
 920       for ( index = 0 ; index < numberofprimaries ; index++){
 921         primary = gime->Primary(listofprimaries[index]) ;
 922         Int_t moduleNumber ;
 923         Double_t primX, primZ ;
 924         phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
 925         if(moduleNumberRec == moduleNumber) {
 926           dX = recX - primX;
 927           dZ = recZ - primZ;
 928           distance = TMath::Sqrt(dX*dX + dZ*dZ) ;
 929           if(minDistance > distance) {
 930             minDistance = distance ;
 931             dXmin = dX;
 932             dZmin = dZ;
 933             closestPrimary = listofprimaries[index] ;
 934           }
 935         }
 936       }
 937
 938       //===========define the "type" of closest primary
 939       if(closestPrimary >=0 ){
 940         Int_t primaryCode = -1;
 941         const TParticle * primary = gime->Primary(closestPrimary) ;
 942         Int_t primaryType = primary->GetPdgCode() ;
 943         if(primaryType == 22) // photon ?
 944           primaryCode = 0 ;
 945         else
 946           if(primaryType == 2112) // neutron
 947             primaryCode = 1 ;
 948           else
 949             if(primaryType == -2112) // Anti neutron
 950               primaryCode = 2 ;
 951             else
 952               if(primaryType == -2122) //Anti proton
 953                 primaryCode = 4 ;
 954               else {
 955                 TParticle tempo(*primary) ;
 956                 if(tempo.GetPDG()->Charge())
 957                   primaryCode = 3 ;
 958               }
 959
 960         //==========Now look at the type of RecParticle
 961         Float_t energy = CorrectedEnergy(recParticle->Energy()) ;
 962         if(recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA){
 963           hPhot->Fill(energy ) ;
 964           switch(primaryCode){
 965           case 0:
 966             hPhotReg->Fill(energy ) ;
 967             break ;
 968           case 1:
 969             hNReg->Fill(energy ) ;
 970             break ;
 971           case 2:
 972             hNBarReg->Fill(energy ) ;
 973             break ;
 974           case 3:
 975             hChargedReg->Fill(energy ) ;
 976             break ;
 977           case 4:
 978             hPbarReg->Fill(energy ) ;
 979             break ;
 980           default:
 981             break ;
 982           }
 983         }
 984         if((recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)||
 985            (recParticle->GetType() == AliPHOSFastRecParticle::kGAMMAHA)){ //with PPSD signal
 986           hPPSD->Fill(energy ) ;
 987           switch(primaryCode){
 988           case 0:
 989             hPhotPPSD->Fill(energy ) ;
 990             break ;
 991           case 1:
 992             hNPPSD->Fill(energy ) ;
 993             break ;
 994           case 2:
 995             hNBarPPSD->Fill(energy ) ;
 996             break ;
 997           case 3:
 998             hChargedPPSD->Fill(energy ) ;
 999             break ;
1000           case 4:
1001             hPbarPPSD->Fill(energy ) ;
1002             break ;
1003           default:
1004             break ;
1005           }
1006         }
1007         if((recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)||
1008            (recParticle->GetType() == AliPHOSFastRecParticle::kELECTRON)||
1009            (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM)||
1010            (recParticle->GetType() == AliPHOSFastRecParticle::kABSURDEM) ){ //with EM shower
1011           hShape->Fill(energy ) ;
1012           switch(primaryCode){
1013           case 0:
1014             hPhotEM->Fill(energy ) ;
1015             break ;
1016           case 1:
1017             hNEM->Fill(energy ) ;
1018             break ;
1019           case 2:
1020             hNBarEM->Fill(energy ) ;
1021             break ;
1022           case 3:
1023             hChargedEM->Fill(energy ) ;
1024             break ;
1025           case 4:
1026             hPbarEM->Fill(energy ) ;
1027             break ;
1028           default:
1029             break ;
1030           }
1031         }
1032
1033         if((recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)||
1034            (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALHA) ||
1035            (recParticle->GetType() == AliPHOSFastRecParticle::kGAMMAHA) ||
1036            (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM) ) //nuetral
1037           hVeto->Fill(energy ) ;
1038
1039         //fill number of primaries identified as ...
1040         if(primaryCode >= 0) // Primary code defined
1041           counter[recParticle->GetType()][primaryCode]++ ;
1042
1043       }
1044
1045     } // no closest primary found
1046   }
1047
1048
1049   //===================  SaveHistograms
1050   cfile->cd() ;
1051   hPrimary->Write(0,kOverwrite);
1052   hAllRP->Write(0,kOverwrite);
1053   hPhot->Write(0,kOverwrite);
1054   hPPSD->Write(0,kOverwrite);
1055   hShape->Write(0,kOverwrite);
1056   hVeto->Write(0,kOverwrite);
1057   hPhotReg->Write(0,kOverwrite);
1058   hPhotEM->Write(0,kOverwrite);
1059   hPhotPPSD->Write(0,kOverwrite);
1060   hNReg ->Write(0,kOverwrite);
1061   hNEM  ->Write(0,kOverwrite);
1062   hNPPSD->Write(0,kOverwrite);
1063   hNBarReg ->Write(0,kOverwrite);
1064   hNBarEM  ->Write(0,kOverwrite);
1065   hNBarPPSD->Write(0,kOverwrite);
1066   hChargedReg ->Write(0,kOverwrite);
1067   hChargedEM  ->Write(0,kOverwrite);
1068   hChargedPPSD->Write(0,kOverwrite);
1069   hPbarReg ->Write(0,kOverwrite);
1070   hPbarEM  ->Write(0,kOverwrite);
1071   hPbarPPSD->Write(0,kOverwrite);
1072
1073   cfile->Write(0,kOverwrite);
1074   cfile->Close();
1075   delete cfile ;
1076
1077
1078   //print Final Table
1079   maxevent = (Int_t)gAlice->TreeE()->GetEntries() ;
1080  //  cout << "Resolutions: Analyzed " << gime->MaxEvent() << " event(s)" << endl ;
1081
1082   cout << "Resolutions: Analyzed " << maxevent << " event(s)" << endl ;
1083   cout << endl ;
1084
1085   cout << "        Primary:    Photon  Neutron  Antineutron  Charged hadron   AntiProton" << endl ;
1086   cout << "--------------------------------------------------------------------------------" << endl ;
1087   cout << "         kGAMMA: "
1088        << setw(8) << counter[2][0] << setw(9)  << counter[2][1] << setw(13) << counter[2][2]
1089        << setw(15)<< counter[2][3] << setw(13) << counter[2][4] << endl ;
1090   cout << "       kGAMMAHA: "
1091        << setw(8) << counter[3][0] << setw(9)  << counter[3][1] << setw(13) << counter[3][2]
1092        << setw(15)<< counter[3][3] << setw(13) << counter[3][4] << endl ;
1093   cout << "     kNEUTRALEM: "
1094        << setw(8) << counter[0][0] << setw(9)  << counter[0][1] << setw(13) << counter[0][2]
1095        << setw(15)<< counter[0][3] << setw(13) << counter[0][4] << endl ;
1096   cout << "     kNEUTRALHA: "
1097        << setw(8) << counter[1][0] << setw(9)  << counter[1][1] << setw(13) << counter[1][2]
1098        << setw(15)<< counter[1][3] << setw(13) << counter[1][4] << endl ;
1099   cout << "      kABSURDEM: "
1100        << setw(8) << counter[4][0] << setw(9)  << counter[4][1] << setw(13) << counter[4][2]
1101        << setw(15)<< counter[4][3] << setw(13) << counter[4][4] << endl ;
1102   cout << "      kABSURDHA: "
1103        << setw(8) << counter[5][0] << setw(9)  << counter[5][1] << setw(13) << counter[5][2]
1104        << setw(15)<< counter[5][3] << setw(13) << counter[5][4] << endl ;
1105   cout << "      kELECTRON: "
1106        << setw(8) << counter[6][0] << setw(9)  << counter[6][1] << setw(13) << counter[6][2]
1107        << setw(15)<< counter[6][3] << setw(13) << counter[6][4] << endl ;
1108   cout << "     kCHARGEDHA: "
1109        << setw(8) << counter[7][0] << setw(9)  << counter[7][1] << setw(13) << counter[7][2]
1110        << setw(15)<< counter[7][3] << setw(13) << counter[7][4] << endl ;
1111   cout << "--------------------------------------------------------------------------------" << endl ;
1112
1113   Int_t totalInd = 0 ;
1114   for(i1 = 0; i1<8; i1++)
1115     for(i2 = 0; i2<5; i2++)
1116       totalInd+=counter[i1][i2] ;
1117   cout << "Indentified particles: " << totalInd << endl ;
1118
1119 }
1120
1121
1122