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