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