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[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
 511 }
 512
 513 //____________________________________________________________________________
 514  void AliPHOSAnalyze::EnergyResolution(const char * branchTitle)
 515 {
 516   //fills two dimentional histo: energy of primary vs. energy of reconstructed
 517
 518   TH2F * hAllEnergy = 0 ;  //all reconstructed with primary photon
 519   TH2F * hPhotEnergy= 0 ;  //kGamma  with primary photon
 520   TH2F * hEMEnergy  = 0 ;  //electromagnetic with primary photon
 521
 522   //opening file and reading histograms if any
 523   TFile * efile = new TFile("energy.root","update");
 524
 525   hAllEnergy = (TH2F*)efile->Get("hAllEnergy") ;
 526   if(hAllEnergy == 0)
 527     hAllEnergy  = new TH2F("hAllEnergy",  "Energy of any RP with primary photon",100, 0., 5., 100, 0., 5.);
 528
 529   hPhotEnergy =(TH2F*) efile->Get("hPhotEnergy") ;
 530   if(hPhotEnergy == 0)
 531     hPhotEnergy = new TH2F("hPhotEnergy", "Energy of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.);
 532
 533   hEMEnergy =(TH2F*) efile->Get("hEMEnergy");
 534   if(hEMEnergy == 0)
 535     hEMEnergy   = new TH2F("hEMEnergy",   "Energy of EM with primary photon",    100, 0., 5., 100, 0., 5.);
 536
 537
 538   fObjGetter = AliPHOSIndexToObject::GetInstance(ffileName.Data(),"PHOSRP",branchTitle) ;
 539   fGeom = AliPHOSGeometry::GetInstance(((AliPHOS*)gAlice->GetModule("PHOS"))->GetGeometry()->GetName(),
 540                                        ((AliPHOS*)gAlice->GetModule("PHOS"))->GetGeometry()->GetTitle() );
 541
 542   Int_t ievent;
 543   for ( ievent=0; ievent < fObjGetter->GetMaxEvent() ; ievent++){
 544
 545     //read the current event
 546     fObjGetter->GetEvent(ievent) ;
 547
 548     AliPHOSRecParticle * recParticle ;
 549     Int_t iRecParticle ;
 550     for(iRecParticle = 0; iRecParticle < fObjGetter->GimeNRecParticles() ;iRecParticle++ ){
 551       recParticle = fObjGetter->GimeRecParticle(iRecParticle) ;
 552
 553       //find the closest primary
 554       Int_t moduleNumberRec ;
 555       Double_t recX, recZ ;
 556       fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
 557
 558       Double_t minDistance  = 100. ;
 559       Int_t closestPrimary = -1 ;
 560
 561       //extract list of primaries: it is stored at EMC RecPoints
 562       Int_t emcIndex = fObjGetter->GimeTrackSegment(recParticle->GetPHOSTSIndex())->GetEmcIndex() ;
 563       Int_t numberofprimaries ;
 564       Int_t * listofprimaries  = fObjGetter->GimeEmcRecPoint(emcIndex)->GetPrimaries(numberofprimaries)  ;
 565
 566       Int_t index ;
 567       TParticle * primary ;
 568       Double_t distance = minDistance ;
 569       Double_t dX, dZ;
 570       Double_t dXmin = 0.;
 571       Double_t dZmin = 0. ;
 572       for ( index = 0 ; index < numberofprimaries ; index++){
 573         primary = fObjGetter->GimePrimary(listofprimaries[index]) ;
 574         Int_t moduleNumber ;
 575         Double_t primX, primZ ;
 576         fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
 577         if(moduleNumberRec == moduleNumber) {
 578           dX = recX - primX;
 579           dZ = recZ - primZ;
 580           distance = TMath::Sqrt(dX*dX + dZ*dZ) ;
 581           if(minDistance > distance) {
 582             minDistance = distance ;
 583             dXmin = dX;
 584             dZmin = dZ;
 585             closestPrimary = listofprimaries[index] ;
 586           }
 587         }
 588       }
 589
 590       //if found primary, fill histograms
 591       if(closestPrimary >=0 ){
 592         TParticle * primary = fObjGetter->GimePrimary(closestPrimary) ;
 593         if(primary->GetPdgCode() == 22){
 594           hAllEnergy->Fill(primary->Energy(), recParticle->Energy()) ;
 595           if(recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA){
 596             hPhotEnergy->Fill(primary->Energy(), recParticle->Energy() ) ;
 597             hEMEnergy->Fill(primary->Energy(), recParticle->Energy() ) ;
 598           }
 599           else
 600             if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM)
 601               hEMEnergy->Fill(primary->Energy(), recParticle->Energy() ) ;
 602         }
 603       }
 604     }
 605   }
 606
 607   //write filled histograms
 608   efile->cd() ;
 609   hAllEnergy->Write(0,kOverwrite) ;
 610   hPhotEnergy->Write(0,kOverwrite) ;
 611   hEMEnergy->Write(0,kOverwrite)  ;
 612   //  efile->Write() ;
 613   efile->Close() ;
 614   delete efile ;
 615
 616 }
 617 //____________________________________________________________________________
 618 void AliPHOSAnalyze::PositionResolution(const char * branchTitle)
 619 {
 620   //fills two dimentional histo: energy vs. primary - reconstructed distance
 621
 622
 623
 624   TH2F * hAllPosition  = 0;     // Position of any RP with primary photon
 625   TH2F * hPhotPosition = 0;    // Position of kGAMMA with primary photon
 626   TH2F * hEMPosition   = 0;      // Position of EM with primary photon
 627
 628   TH1F * hAllPositionX = 0;    // X-Position Resolution of photons with photon primary
 629   TH1F * hAllPositionZ = 0;    // Z-Position Resolution of photons with photon primary
 630
 631
 632   //opening file and reading histograms if any
 633   TFile * pfile = new TFile("position.root","update");
 634
 635   hAllPosition = (TH2F*)pfile->Get("hAllPosition");
 636   if(hAllPosition == 0)
 637     hAllPosition  = new TH2F("hAllPosition",
 638                              "Position of any RP with primary photon",100, 0., 5., 100, 0., 5.);
 639   hPhotPosition= (TH2F*)pfile->Get("hPhotPosition");
 640   if(hPhotPosition == 0)
 641     hPhotPosition = new TH2F("hPhotPosition",
 642                              "Position of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.);
 643   hEMPosition= (TH2F*)pfile->Get("hEMPosition") ;
 644   if(hEMPosition == 0)
 645     hEMPosition   = new TH2F("hEMPosition",
 646                              "Position of EM with primary photon",    100, 0., 5., 100, 0., 5.);
 647   hAllPositionX = (TH1F*)pfile->Get("hAllPositionX") ;
 648   if(hAllPositionX == 0)
 649     hAllPositionX = new TH1F("hAllPositionX",
 650                              "Delta X of any RP with primary photon",100, -2., 2.);
 651   hAllPositionZ =(TH1F*) pfile->Get("hAllPositionZ") ;
 652   if(hAllPositionZ == 0)
 653     hAllPositionZ = new TH1F("hAllPositionZ",
 654                              "Delta X of any RP with primary photon",100, -2., 2.);
 655
 656
 657   fObjGetter = AliPHOSIndexToObject::GetInstance(ffileName.Data(),"PHOSRP",branchTitle) ;
 658   fGeom = AliPHOSGeometry::GetInstance(((AliPHOS*)gAlice->GetModule("PHOS"))->GetGeometry()->GetName(),
 659                                        ((AliPHOS*)gAlice->GetModule("PHOS"))->GetGeometry()->GetTitle() );
 660
 661   Int_t ievent;
 662   for ( ievent=0; ievent < fObjGetter->GetMaxEvent() ; ievent++){
 663
 664     //read the current event
 665     fObjGetter->GetEvent(ievent) ;
 666
 667     AliPHOSRecParticle * recParticle ;
 668     Int_t iRecParticle ;
 669     for(iRecParticle = 0; iRecParticle < fObjGetter->GimeNRecParticles() ;iRecParticle++ ){
 670       recParticle = fObjGetter->GimeRecParticle(iRecParticle) ;
 671
 672       //find the closest primary
 673       Int_t moduleNumberRec ;
 674       Double_t recX, recZ ;
 675       fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
 676
 677       Double_t minDistance  = 100. ;
 678       Int_t closestPrimary = -1 ;
 679
 680       //extract list of primaries: it is stored at EMC RecPoints
 681       Int_t emcIndex = fObjGetter->GimeTrackSegment(recParticle->GetPHOSTSIndex())->GetEmcIndex() ;
 682       Int_t numberofprimaries ;
 683       Int_t * listofprimaries  = fObjGetter->GimeEmcRecPoint(emcIndex)->GetPrimaries(numberofprimaries)  ;
 684
 685       Int_t index ;
 686       TParticle * primary ;
 687       Double_t distance = minDistance ;
 688       Double_t dX = 1000; // incredible number
 689       Double_t dZ = 1000; // for the case if no primary will be found
 690       Double_t dXmin = 0.;
 691       Double_t dZmin = 0. ;
 692       for ( index = 0 ; index < numberofprimaries ; index++){
 693         primary = fObjGetter->GimePrimary(listofprimaries[index]) ;
 694         Int_t moduleNumber ;
 695         Double_t primX, primZ ;
 696         fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
 697         if(moduleNumberRec == moduleNumber) {
 698           dX = recX - primX;
 699           dZ = recZ - primZ;
 700           distance = TMath::Sqrt(dX*dX + dZ*dZ) ;
 701           if(minDistance > distance) {
 702             minDistance = distance ;
 703             dXmin = dX;
 704             dZmin = dZ;
 705             closestPrimary = listofprimaries[index] ;
 706           }
 707         }
 708       }
 709
 710       //if found primary, fill histograms
 711       if(closestPrimary >=0 ){
 712         TParticle * primary = fObjGetter->GimePrimary(closestPrimary) ;
 713         if(primary->GetPdgCode() == 22){
 714           hAllPosition->Fill(primary->Energy(), minDistance) ;
 715           hAllPositionX->Fill(primary->Energy(), dX) ;
 716           hAllPositionZ->Fill(primary->Energy(), dZ) ;
 717           if(recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA){
 718             hPhotPosition->Fill(primary->Energy(), minDistance ) ;
 719             hEMPosition->Fill(primary->Energy(), minDistance ) ;
 720           }
 721           else
 722             if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM)
 723               hEMPosition->Fill(primary->Energy(), minDistance ) ;
 724         }
 725       }
 726     }
 727   }
 728
 729   //Write output histgrams
 730   pfile->cd() ;
 731   hAllPosition->Write(0,kOverwrite) ;
 732   hAllPositionX->Write(0,kOverwrite) ;
 733   hAllPositionZ->Write(0,kOverwrite) ;
 734   hPhotPosition->Write(0,kOverwrite) ;
 735   hEMPosition->Write(0,kOverwrite) ;
 736   pfile->Write() ;
 737   pfile->Close() ;
 738   delete pfile ;
 739
 740
 741 }
 742 //____________________________________________________________________________
 743 void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){
 744 // fills spectra of primary photons and several kinds of
 745 // reconstructed particles, so that analyzing them one can
 746 // estimate conatmination, efficiency of registration etc.
 747
 748   //define several general histograms
 749   TH1F * hPrimary = 0;   //spectrum (P_t distribution) of primary photons
 750   TH1F * hAllRP   = 0;   //spectrum of all RecParticles in PHOS
 751   TH1F * hPhot    = 0;   //spectrum of kGAMMA RecParticles
 752   TH1F * hPPSD    = 0;   //spectrum of all RecParticles with PPSD signal
 753   TH1F * hShape   = 0;   //spectrum of all EM RecParticles
 754   TH1F * hVeto    = 0;   //spectrum of all neutral RecParticles
 755
 756   //Now separate histograms in accoradance with primary
 757   //primary - photon
 758   TH1F * hPhotReg = 0;   //Registeres as photon
 759   TH1F * hPhotEM  = 0;   //Registered as EM
 760   TH1F * hPhotPPSD= 0;   //Registered as RecParticle with PPSD signal
 761
 762   //primary - n
 763   TH1F * hNReg = 0;   //Registeres as photon
 764   TH1F * hNEM  = 0;   //Registered as EM
 765   TH1F * hNPPSD= 0;   //Registered as RecParticle with PPSD signal
 766
 767   //primary - nBar
 768   TH1F * hNBarReg = 0;   //Registeres as photon
 769   TH1F * hNBarEM  = 0;   //Registered as EM
 770   TH1F * hNBarPPSD= 0;   //Registered as RecParticle with PPSD signal
 771
 772   //primary - charged hadron (pBar excluded)
 773   TH1F * hChargedReg = 0;   //Registeres as photon
 774   TH1F * hChargedEM  = 0;   //Registered as EM
 775   TH1F * hChargedPPSD= 0;   //Registered as RecParticle with PPSD signal
 776
 777   //primary - pBar
 778   TH1F * hPbarReg = 0;   //Registeres as photon
 779   TH1F * hPbarEM  = 0;   //Registered as EM
 780   TH1F * hPbarPPSD= 0;   //Registered as RecParticle with PPSD signal
 781
 782
 783   //Reading histograms from the file
 784   TFile * cfile = new TFile("contamination.root","update") ;
 785
 786   //read general histograms
 787   hPrimary = (TH1F*) cfile->Get("hPrimary") ;
 788   if(hPrimary == 0)
 789     hPrimary= new TH1F("hPrimary", "Primary photon spectrum", 100, 0., 5.);
 790   hAllRP = (TH1F*)cfile->Get("hAllRP") ;
 791   if(hAllRP == 0)
 792     hAllRP = new TH1F("hAllRP","All Reconstructed particles", 100, 0., 5.);
 793   hPhot  = (TH1F*)cfile->Get("hPhot") ;
 794   if(hPhot == 0)
 795     hPhot = new TH1F("hPhot","All kGAMMA RecParticles",100, 0., 5.);
 796   hPPSD = (TH1F*)cfile->Get("hPPSD") ;
 797   if(hPPSD == 0)
 798     hPPSD  = new TH1F("hPPSD", "All PPSD Recparticles",    100, 0., 5.);
 799   hShape = (TH1F*) cfile->Get("hShape") ;
 800   if(hShape == 0)
 801     hShape = new TH1F("hShape","All particles with EM shower",100, 0., 5.);
 802   hVeto= (TH1F*)cfile->Get("hVeto") ;
 803   if(hVeto == 0)
 804     hVeto  = new TH1F("hVeto", "All uncharged particles",      100, 0., 5.);
 805
 806
 807   //primary - photon
 808   hPhotReg = (TH1F*)cfile->Get("hPhotReg");
 809   if(hPhotReg == 0)
 810     hPhotReg   = new TH1F("hPhotReg","Photon registered as photon",100, 0., 5.);
 811   hPhotEM  =(TH1F*)cfile->Get("hPhotEM");
 812   if(hPhotEM== 0)
 813     hPhotEM   = new TH1F("hPhotEM",  "Photon registered as EM", 100, 0., 5.);
 814   hPhotPPSD= (TH1F*)cfile->Get("hPhotPPSD");
 815   if(hPhotPPSD== 0)
 816     hPhotPPSD   = new TH1F("hPhotPPSD","Photon registered as PPSD", 100, 0., 5.);
 817
 818   //primary - n
 819   hNReg = (TH1F*)cfile->Get("hNReg");
 820   if(hNReg== 0)
 821    hNReg      = new TH1F("hNReg",   "N registered as photon",              100, 0., 5.);
 822   hNEM  = (TH1F*)cfile->Get("hNEM");
 823   if(hNEM== 0)
 824     hNEM      = new TH1F("hNEM",    "N registered as EM",      100, 0., 5.);
 825   hNPPSD=(TH1F*)cfile->Get("hNPPSD");
 826   if(hNPPSD== 0)
 827     hNPPSD      = new TH1F("hNPPSD","N registered as PPSD",      100, 0., 5.);
 828
 829   //primary - nBar
 830   hNBarReg =(TH1F*)cfile->Get("hNBarReg");
 831   if(hNBarReg== 0)
 832    hNBarReg   = new TH1F("hNBarReg", "NBar registered as photon",           100, 0., 5.);
 833   hNBarEM  =(TH1F*)cfile->Get("hNBarEM");
 834   if(hNBarEM== 0)
 835     hNBarEM   = new TH1F("hNBarEM",  "NBar registered as EM",   100, 0., 5.);
 836   hNBarPPSD=(TH1F*)cfile->Get("hNBarPPSD");
 837   if(hNBarPPSD== 0)
 838     hNBarPPSD   = new TH1F("hNBarPPSD","NBar registered as PPSD",   100, 0., 5.);
 839
 840   //primary - charged hadron (pBar excluded)
 841   hChargedReg = (TH1F*)cfile->Get("hChargedReg");
 842   if(hChargedReg== 0)
 843     hChargedReg= new TH1F("hChargedReg", "Charged hadron registered as photon",100, 0., 5.);
 844   hChargedEM  = (TH1F*)cfile->Get("hChargedEM");
 845   if(hChargedEM== 0)
 846     hChargedEM= new TH1F("hChargedEM","Charged registered as EM",100, 0., 5.);
 847   hChargedPPSD= (TH1F*)cfile->Get("hChargedPPSD");
 848   if(hChargedPPSD== 0)
 849     hChargedPPSD= new TH1F("hChargedPPSD","Charged registered as PPSD",100, 0., 5.);
 850
 851   //primary - pBar
 852   hPbarReg = (TH1F*)cfile->Get("hPbarReg");
 853   if(hPbarReg== 0)
 854     hPbarReg= new TH1F("hPbarReg", "pBar registered as photon",100, 0., 5.);
 855   hPbarEM  = (TH1F*)cfile->Get("hPbarEM");
 856   if(hPbarEM== 0)
 857     hPbarEM= new TH1F("hPbarEM","Pbar registered as EM",100, 0., 5.);
 858   hPbarPPSD= (TH1F*)cfile->Get("hPbarPPSD");
 859   if(hPbarPPSD== 0)
 860     hPbarPPSD= new TH1F("hPbarPPSD","Pbar as PPSD",100, 0., 5.);
 861
 862
 863   //Now make some initializations
 864
 865   Int_t counter[8][5] ;      //# of registered particles
 866   Int_t i1,i2 ;
 867   for(i1 = 0; i1<8; i1++)
 868     for(i2 = 0; i2<5; i2++)
 869       counter[i1][i2] = 0 ;
 870
 871
 872
 873   fObjGetter = AliPHOSIndexToObject::GetInstance(ffileName.Data(),"PHOSRP",RecPointsTitle) ;
 874   fGeom = AliPHOSGeometry::GetInstance(((AliPHOS*)gAlice->GetModule("PHOS"))->GetGeometry()->GetName(),
 875                                        ((AliPHOS*)gAlice->GetModule("PHOS"))->GetGeometry()->GetTitle() );
 876
 877   Int_t ievent;
 878   for ( ievent=0; ievent < fObjGetter->GetMaxEvent() ; ievent++){
 879
 880     fObjGetter->GetEvent(ievent) ;
 881
 882     //=========== Make spectrum of the primary photons
 883     TParticle * primary ;
 884     Int_t iPrimary ;
 885     for( iPrimary = 0 ; iPrimary < fObjGetter->GimeNPrimaries() ; iPrimary++){
 886       primary = fObjGetter->GimePrimary(iPrimary) ;
 887       Int_t primaryType = primary->GetPdgCode() ;
 888       if( primaryType == 22 ) {
 889         //check, if photons folls onto PHOS
 890         Int_t moduleNumber ;
 891         Double_t primX, primZ ;
 892         fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
 893         if(moduleNumber)
 894           hPrimary->Fill(primary->Energy()) ;
 895
 896       }
 897
 898     }
 899     //========== Now scan over RecParticles
 900     AliPHOSRecParticle * recParticle ;
 901     Int_t iRecParticle ;
 902     for(iRecParticle = 0; iRecParticle < fObjGetter->GimeNRecParticles(); iRecParticle++ ){
 903       recParticle = fObjGetter->GimeRecParticle(iRecParticle) ;
 904       //fill histo spectrum of all RecParticles
 905       hAllRP->Fill(CorrectedEnergy(recParticle->Energy())) ;
 906
 907       //==========find the closest primary
 908       Int_t moduleNumberRec ;
 909       Double_t recX, recZ ;
 910       fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
 911
 912       Double_t minDistance  = 100. ;
 913       Int_t closestPrimary = -1 ;
 914
 915       //extract list of primaries: it is stored at EMC RecPoints
 916       Int_t emcIndex = fObjGetter->GimeTrackSegment(recParticle->GetPHOSTSIndex())->GetEmcIndex() ;
 917       Int_t numberofprimaries ;
 918       Int_t * listofprimaries  = fObjGetter->GimeEmcRecPoint(emcIndex)->GetPrimaries(numberofprimaries)  ;
 919       Int_t index ;
 920       TParticle * primary ;
 921       Double_t distance = minDistance ;
 922       Double_t dX, dZ;
 923       Double_t dXmin = 0.;
 924       Double_t dZmin = 0. ;
 925       for ( index = 0 ; index < numberofprimaries ; index++){
 926         primary = fObjGetter->GimePrimary(listofprimaries[index]) ;
 927         Int_t moduleNumber ;
 928         Double_t primX, primZ ;
 929         fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
 930         if(moduleNumberRec == moduleNumber) {
 931           dX = recX - primX;
 932           dZ = recZ - primZ;
 933           distance = TMath::Sqrt(dX*dX + dZ*dZ) ;
 934           if(minDistance > distance) {
 935             minDistance = distance ;
 936             dXmin = dX;
 937             dZmin = dZ;
 938             closestPrimary = listofprimaries[index] ;
 939           }
 940         }
 941       }
 942
 943       //===========define the "type" of closest primary
 944       if(closestPrimary >=0 ){
 945         Int_t primaryCode = -1;
 946         TParticle * primary = fObjGetter->GimePrimary(closestPrimary) ;
 947         Int_t primaryType = primary->GetPdgCode() ;
 948         if(primaryType == 22) // photon ?
 949           primaryCode = 0 ;
 950         else
 951           if(primaryType == 2112) // neutron
 952             primaryCode = 1 ;
 953           else
 954             if(primaryType == -2112) // Anti neutron
 955               primaryCode = 2 ;
 956             else
 957               if(primaryType == -2122) //Anti proton
 958                 primaryCode = 4 ;
 959               else
 960                 if(fObjGetter->GimePrimary(closestPrimary)->GetPDG()->Charge())
 961                   primaryCode = 3 ;
 962
 963         //==========Now look at the type of RecParticle
 964         Float_t energy = CorrectedEnergy(recParticle->Energy()) ;
 965         if(recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA){
 966           hPhot->Fill(energy ) ;
 967           switch(primaryCode){
 968           case 0:
 969             hPhotReg->Fill(energy ) ;
 970             break ;
 971           case 1:
 972             hNReg->Fill(energy ) ;
 973             break ;
 974           case 2:
 975             hNBarReg->Fill(energy ) ;
 976             break ;
 977           case 3:
 978             hChargedReg->Fill(energy ) ;
 979             break ;
 980           case 4:
 981             hPbarReg->Fill(energy ) ;
 982             break ;
 983           default:
 984             break ;
 985           }
 986         }
 987         if((recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)||
 988            (recParticle->GetType() == AliPHOSFastRecParticle::kGAMMAHA)){ //with PPSD signal
 989           hPPSD->Fill(energy ) ;
 990           switch(primaryCode){
 991           case 0:
 992             hPhotPPSD->Fill(energy ) ;
 993             break ;
 994           case 1:
 995             hNPPSD->Fill(energy ) ;
 996             break ;
 997           case 2:
 998             hNBarPPSD->Fill(energy ) ;
 999             break ;
1000           case 3:
1001             hChargedPPSD->Fill(energy ) ;
1002             break ;
1003           case 4:
1004             hPbarPPSD->Fill(energy ) ;
1005             break ;
1006           default:
1007             break ;
1008           }
1009         }
1010         if((recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)||
1011            (recParticle->GetType() == AliPHOSFastRecParticle::kELECTRON)||
1012            (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM)||
1013            (recParticle->GetType() == AliPHOSFastRecParticle::kABSURDEM) ){ //with EM shower
1014           hShape->Fill(energy ) ;
1015           switch(primaryCode){
1016           case 0:
1017             hPhotEM->Fill(energy ) ;
1018             break ;
1019           case 1:
1020             hNEM->Fill(energy ) ;
1021             break ;
1022           case 2:
1023             hNBarEM->Fill(energy ) ;
1024             break ;
1025           case 3:
1026             hChargedEM->Fill(energy ) ;
1027             break ;
1028           case 4:
1029             hPbarEM->Fill(energy ) ;
1030             break ;
1031           default:
1032             break ;
1033           }
1034         }
1035
1036         if((recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)||
1037            (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALHA) ||
1038            (recParticle->GetType() == AliPHOSFastRecParticle::kGAMMAHA) ||
1039            (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM) ) //nuetral
1040           hVeto->Fill(energy ) ;
1041
1042         //fill number of primaries identified as ...
1043         if(primaryCode >= 0) // Primary code defined
1044           counter[recParticle->GetType()][primaryCode]++ ;
1045
1046       }
1047
1048     } // no closest primary found
1049   }
1050
1051
1052   //===================  SaveHistograms
1053   cfile->cd() ;
1054   hPrimary->Write(0,kOverwrite);
1055   hAllRP->Write(0,kOverwrite);
1056   hPhot->Write(0,kOverwrite);
1057   hPPSD->Write(0,kOverwrite);
1058   hShape->Write(0,kOverwrite);
1059   hVeto->Write(0,kOverwrite);
1060   hPhotReg->Write(0,kOverwrite);
1061   hPhotEM->Write(0,kOverwrite);
1062   hPhotPPSD->Write(0,kOverwrite);
1063   hNReg ->Write(0,kOverwrite);
1064   hNEM  ->Write(0,kOverwrite);
1065   hNPPSD->Write(0,kOverwrite);
1066   hNBarReg ->Write(0,kOverwrite);
1067   hNBarEM  ->Write(0,kOverwrite);
1068   hNBarPPSD->Write(0,kOverwrite);
1069   hChargedReg ->Write(0,kOverwrite);
1070   hChargedEM  ->Write(0,kOverwrite);
1071   hChargedPPSD->Write(0,kOverwrite);
1072   hPbarReg ->Write(0,kOverwrite);
1073   hPbarEM  ->Write(0,kOverwrite);
1074   hPbarPPSD->Write(0,kOverwrite);
1075
1076   cfile->Write(0,kOverwrite);
1077   cfile->Close();
1078   delete cfile ;
1079
1080
1081   //print Final Table
1082
1083   cout << "Resolutions: Analyzed " << fObjGetter->GetMaxEvent() << " event(s)" << endl ;
1084   cout << endl ;
1085
1086   cout << "        Primary:    Photon  Neutron  Antineutron  Charged hadron   AntiProton" << endl ;
1087   cout << "--------------------------------------------------------------------------------" << endl ;
1088   cout << "         kGAMMA: "
1089        << setw(8) << counter[2][0] << setw(9)  << counter[2][1] << setw(13) << counter[2][2]
1090        << setw(15)<< counter[2][3] << setw(13) << counter[2][4] << endl ;
1091   cout << "       kGAMMAHA: "
1092        << setw(8) << counter[3][0] << setw(9)  << counter[3][1] << setw(13) << counter[3][2]
1093        << setw(15)<< counter[3][3] << setw(13) << counter[3][4] << endl ;
1094   cout << "     kNEUTRALEM: "
1095        << setw(8) << counter[0][0] << setw(9)  << counter[0][1] << setw(13) << counter[0][2]
1096        << setw(15)<< counter[0][3] << setw(13) << counter[0][4] << endl ;
1097   cout << "     kNEUTRALHA: "
1098        << setw(8) << counter[1][0] << setw(9)  << counter[1][1] << setw(13) << counter[1][2]
1099        << setw(15)<< counter[1][3] << setw(13) << counter[1][4] << endl ;
1100   cout << "      kABSURDEM: "
1101        << setw(8) << counter[4][0] << setw(9)  << counter[4][1] << setw(13) << counter[4][2]
1102        << setw(15)<< counter[4][3] << setw(13) << counter[4][4] << endl ;
1103   cout << "      kABSURDHA: "
1104        << setw(8) << counter[5][0] << setw(9)  << counter[5][1] << setw(13) << counter[5][2]
1105        << setw(15)<< counter[5][3] << setw(13) << counter[5][4] << endl ;
1106   cout << "      kELECTRON: "
1107        << setw(8) << counter[6][0] << setw(9)  << counter[6][1] << setw(13) << counter[6][2]
1108        << setw(15)<< counter[6][3] << setw(13) << counter[6][4] << endl ;
1109   cout << "     kCHARGEDHA: "
1110        << setw(8) << counter[7][0] << setw(9)  << counter[7][1] << setw(13) << counter[7][2]
1111        << setw(15)<< counter[7][3] << setw(13) << counter[7][4] << endl ;
1112   cout << "--------------------------------------------------------------------------------" << endl ;
1113
1114   Int_t totalInd = 0 ;
1115   for(i1 = 0; i1<8; i1++)
1116     for(i2 = 0; i2<5; i2++)
1117       totalInd+=counter[i1][i2] ;
1118   cout << "Indentified particles: " << totalInd << endl ;
1119
1120 }
1121
1122
1123