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