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