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 PHOSv1 events:
  20 // Construct histograms and displays them.
  21 // Use the macro EditorBar.C for best access to the functionnalities
  22 //
  23 //*-- Author: Y. Schutz (SUBATECH) & Gines Martinez (SUBATECH)
  24 //////////////////////////////////////////////////////////////////////////////
  25
  26 // --- ROOT system ---
  27
  28 #include "TFile.h"
  29 #include "TH1.h"
  30 #include "TPad.h"
  31 #include "TH2.h"
  32 #include "TParticle.h"
  33 #include "TClonesArray.h"
  34 #include "TTree.h"
  35 #include "TMath.h"
  36 #include "TCanvas.h"
  37
  38 // --- Standard library ---
  39
  40 #include <iostream.h>
  41 #include <stdio.h>
  42
  43 // --- AliRoot header files ---
  44
  45 #include "AliRun.h"
  46 #include "AliPHOSAnalyze.h"
  47 #include "AliPHOSClusterizerv1.h"
  48 #include "AliPHOSTrackSegmentMakerv1.h"
  49 #include "AliPHOSPIDv1.h"
  50 #include "AliPHOSReconstructioner.h"
  51 #include "AliPHOSDigit.h"
  52 #include "AliPHOSTrackSegment.h"
  53 #include "AliPHOSRecParticle.h"
  54 #include "AliPHOSIndexToObject.h"
  55
  56 ClassImp(AliPHOSAnalyze)
  57
  58
  59 //____________________________________________________________________________
  60   AliPHOSAnalyze::AliPHOSAnalyze()
  61 {
  62   // default ctor (useless)
  63
  64   fRootFile = 0 ;
  65 }
  66
  67 //____________________________________________________________________________
  68 AliPHOSAnalyze::AliPHOSAnalyze(Text_t * name)
  69 {
  70   // ctor: analyze events from root file "name"
  71
  72   Bool_t ok = OpenRootFile(name)  ;
  73   if ( !ok ) {
  74     cout << " AliPHOSAnalyze > Error opening " << name << endl ;
  75   }
  76   else {
  77       //========== Get AliRun object from file
  78       gAlice = (AliRun*) fRootFile->Get("gAlice") ;
  79
  80       //=========== Get the PHOS object and associated geometry from the file
  81       fPHOS  = (AliPHOSv1 *)gAlice->GetDetector("PHOS") ;
  82       fGeom  = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() );
  83
  84       //========== Initializes the Index to Object converter
  85       fObjGetter = AliPHOSIndexToObject::GetInstance(fPHOS) ;
  86       //========== Current event number
  87       fEvt = -999 ;
  88
  89   }
  90   fClu = 0 ;
  91   fPID = 0 ;
  92   fTrs = 0 ;
  93   fRec = 0 ;
  94   ResetHistograms() ;
  95 }
  96
  97 //____________________________________________________________________________
  98 AliPHOSAnalyze::AliPHOSAnalyze(const AliPHOSAnalyze & ana)
  99 {
 100   // copy ctor
 101   ( (AliPHOSAnalyze &)ana ).Copy(*this) ;
 102 }
 103
 104 //____________________________________________________________________________
 105 void AliPHOSAnalyze::Copy(TObject & obj)
 106 {
 107   // copy an analysis into an other one
 108   TObject::Copy(obj) ;
 109   // I do nothing more because the copy is silly but the Code checkers requires one
 110 }
 111
 112 //____________________________________________________________________________
 113 AliPHOSAnalyze::~AliPHOSAnalyze()
 114 {
 115   // dtor
 116
 117   if (fRootFile->IsOpen() )
 118     fRootFile->Close() ;
 119   if(fRootFile)
 120     delete fRootFile ;
 121
 122   if(fPHOS)
 123     delete fPHOS ;
 124
 125   if(fClu)
 126     delete fClu ;
 127
 128   if(fPID)
 129     delete fPID ;
 130
 131   if(fRec)
 132     delete fRec ;
 133
 134   if(fTrs)
 135     delete fTrs ;
 136
 137 }
 138
 139 //____________________________________________________________________________
 140 void AliPHOSAnalyze::ActivePPSD(Int_t Nevents=1){
 141
 142   fhEnergyCorrelations  = new TH2F("hEnergyCorrelations","hEnergyCorrelations",40,  0., 0.15, 30, 0., 3.e-5);
 143   //========== Create the Clusterizer
 144   fClu = new AliPHOSClusterizerv1() ;
 145   fClu->SetEmcEnergyThreshold(0.01) ;
 146   fClu->SetEmcClusteringThreshold(0.20) ;
 147   fClu->SetPpsdEnergyThreshold    (0.0000002) ;
 148   fClu->SetPpsdClusteringThreshold(0.0000001) ;
 149   fClu->SetLocalMaxCut(0.02) ;
 150   fClu->SetCalibrationParameters(0., 0.00000001) ;
 151
 152   Int_t ievent;
 153
 154   for ( ievent=0; ievent<Nevents; ievent++)
 155     {
 156
 157       //========== Event Number>
 158       if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. )
 159         cout <<  "AnalyzeResolutions > " << "Event is " << ievent << endl ;
 160
 161       //=========== Connects the various Tree's for evt
 162       gAlice->GetEvent(ievent);
 163
 164       //=========== Gets the Kine TTree
 165       gAlice->TreeK()->GetEvent(0) ;
 166
 167       //=========== Get the Digit Tree
 168       gAlice->TreeD()->GetEvent(0) ;
 169
 170       //========== Creating branches ===================================
 171       AliPHOSRecPoint::RecPointsList ** EmcRecPoints =  fPHOS->EmcRecPoints() ;
 172       gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP", EmcRecPoints ) ;
 173
 174       AliPHOSRecPoint::RecPointsList ** PpsdRecPoints = fPHOS->PpsdRecPoints() ;
 175       gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", PpsdRecPoints ) ;
 176
 177       AliPHOSTrackSegment::TrackSegmentsList **  TrackSegmentsList = fPHOS->TrackSegments() ;
 178       if( (*TrackSegmentsList) )
 179         (*TrackSegmentsList)->Clear() ;
 180       gAlice->TreeR()->SetBranchAddress( "PHOSTS", TrackSegmentsList ) ;
 181
 182       AliPHOSRecParticle::RecParticlesList ** RecParticleList  = fPHOS->RecParticles() ;
 183       if( (*RecParticleList) )
 184         (*RecParticleList)->Clear() ;
 185       gAlice->TreeR()->SetBranchAddress( "PHOSRP", RecParticleList ) ;
 186
 187
 188       //=========== Gets the Reconstraction TTree
 189       gAlice->TreeR()->GetEvent(0) ;
 190
 191       AliPHOSPpsdRecPoint * RecPoint ;
 192       Int_t relid[4] ;
 193       TIter nextRP(*fPHOS->PpsdRecPoints() ) ;
 194       while( ( RecPoint = (AliPHOSPpsdRecPoint *)nextRP() ) )
 195         {
 196           if(!(RecPoint->GetUp()) ) {
 197             AliPHOSDigit *digit ;
 198             Int_t iDigit ;
 199             for(iDigit = 0; iDigit < fPHOS->Digits()->GetEntries(); iDigit++)
 200               {
 201                 digit = (AliPHOSDigit *) fPHOS->Digits()->At(iDigit) ;
 202                 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
 203                 if((relid[2]==1)&&(relid[3]==1)&&(relid[0]==RecPoint->GetPHOSMod())){
 204                   Float_t ConvertorEnergy = fClu->Calibrate(digit->GetAmp()) ;
 205                   fhEnergyCorrelations->Fill(ConvertorEnergy,RecPoint->GetTotalEnergy() );
 206                   break ;
 207                 }
 208               }
 209             break ;
 210           }
 211         }
 212     }
 213   SaveHistograms() ;
 214   fhEnergyCorrelations->Draw("BOX") ;
 215 }
 216
 217
 218 //____________________________________________________________________________
 219 void AliPHOSAnalyze::AnalyzeManyEvents(Int_t Nevents, Int_t module)
 220 {
 221   // analyzes Nevents events in a single PHOS module
 222   // Events should be reconstructed by Reconstruct()
 223
 224   if ( fRootFile == 0 )
 225     cout << "AnalyzeManyEvents > " << "Root File not openned" << endl ;
 226   else
 227     {
 228       //========== Booking Histograms
 229       cout << "AnalyzeManyEvents > " << "Booking Histograms" << endl ;
 230       BookingHistograms();
 231
 232       Int_t ievent;
 233       Int_t relid[4] ;
 234       AliPHOSDigit * digit ;
 235       AliPHOSEmcRecPoint * emc ;
 236       AliPHOSPpsdRecPoint * ppsd ;
 237       //      AliPHOSTrackSegment * tracksegment ;
 238       AliPHOSRecParticle * recparticle;
 239
 240       for ( ievent=0; ievent<Nevents; ievent++)
 241         {
 242           //========== Event Number>
 243           if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. )
 244             cout <<  "AnalyzeManyEvents > " << "Event is " << ievent << endl ;
 245
 246           //=========== Connects the various Tree's for evt
 247           gAlice->GetEvent(ievent);
 248
 249           //=========== Gets the Digit TTree
 250           gAlice->TreeD()->GetEvent(0) ;
 251
 252           //=========== Gets the number of entries in the Digits array
 253           TIter nextdigit(fPHOS->Digits()) ;
 254           while( ( digit = (AliPHOSDigit *)nextdigit() ) )
 255             {
 256               fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
 257               if (fClu->IsInEmc(digit)) fhEmcDigit->Fill(fClu->Calibrate(digit->GetAmp())) ;
 258               else
 259                 {
 260                   if (relid[1]<17) fhVetoDigit->Fill(fClu->Calibrate(digit->GetAmp()));
 261                   if (relid[1]>16) fhConvertorDigit->Fill(fClu->Calibrate(digit->GetAmp()));
 262                 }
 263             }
 264
 265
 266           //=========== Cluster in module
 267           TIter nextEmc(*fPHOS->EmcRecPoints()  ) ;
 268           while((emc = (AliPHOSEmcRecPoint *)nextEmc()))
 269             {
 270               if ( emc->GetPHOSMod() == module )
 271                 {
 272                   fhEmcCluster->Fill(  emc->GetTotalEnergy()  );
 273                   TIter nextPpsd( *fPHOS->PpsdRecPoints()) ;
 274                   while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd()))
 275                     {
 276                       if ( ppsd->GetPHOSMod() == module )
 277                         {
 278                           if (!ppsd->GetUp()) fhConvertorEmc->Fill(emc->GetTotalEnergy(),ppsd->GetTotalEnergy()) ;
 279                         }
 280                     }
 281                 }
 282             }
 283
 284           //=========== Cluster in module PPSD Down
 285           TIter nextPpsd(*fPHOS->PpsdRecPoints() ) ;
 286           while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd()))
 287             {
 288               if ( ppsd->GetPHOSMod() == module )
 289                 {
 290                   if (!ppsd->GetUp()) fhConvertorCluster->Fill(ppsd->GetTotalEnergy()) ;
 291                   if (ppsd->GetUp())  fhVetoCluster     ->Fill(ppsd->GetTotalEnergy()) ;
 292                 }
 293             }
 294
 295           //========== TRackSegments in the event
 296           TIter nextRecParticle(*fPHOS->RecParticles() ) ;
 297           while((recparticle = (AliPHOSRecParticle *)nextRecParticle()))
 298             {
 299               if ( recparticle->GetPHOSTrackSegment()->GetPHOSMod() == module )
 300                 {
 301                   cout << "Particle type is " << recparticle->GetType() << endl ;
 302                   Int_t numberofprimaries = 0 ;
 303                   Int_t * listofprimaries = recparticle->GetPrimaries(numberofprimaries) ;
 304                   cout << "Number of primaries = " << numberofprimaries << endl ;
 305                   Int_t index ;
 306                   for ( index = 0 ; index < numberofprimaries ; index++)
 307                     cout << "    primary # " << index << " =  " << listofprimaries[index] << endl ;
 308                 }
 309             }
 310         }   // endfor
 311       SaveHistograms();
 312     }       // endif
 313 }           // endfunction
 314
 315 //____________________________________________________________________________
 316  void AliPHOSAnalyze::Reconstruct(Int_t Nevents,Int_t FirstEvent )
 317 {
 318   Int_t ievent ;
 319   for ( ievent=FirstEvent; ievent<Nevents; ievent++)
 320     {
 321       if (ievent==FirstEvent)
 322         {
 323           cout << "Analyze > Starting Reconstructing " << endl ;
 324           //========== Create the Clusterizer
 325           fClu = new AliPHOSClusterizerv1() ;
 326           fClu->SetEmcEnergyThreshold(0.03) ;
 327           fClu->SetEmcClusteringThreshold(0.20) ;
 328           fClu->SetPpsdEnergyThreshold    (0.0000001) ;
 329           fClu->SetPpsdClusteringThreshold(0.0000001) ;
 330           fClu->SetLocalMaxCut(0.02) ;
 331           fClu->SetCalibrationParameters(0., 0.00000001) ;
 332
 333           //========== Creates the track segment maker
 334           fTrs = new AliPHOSTrackSegmentMakerv1()  ;
 335           //      fTrs->UnsetUnfoldFlag() ;
 336
 337           //========== Creates the particle identifier
 338           fPID = new AliPHOSPIDv1() ;
 339           fPID->SetShowerProfileCuts(0.3, 1.8, 0.3, 1.8 ) ;
 340           fPID->SetDispersionCutOff(2.0) ;
 341           fPID->SetRelativeDistanceCut(3.) ;
 342
 343           //========== Creates the Reconstructioner
 344           fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;
 345           //              fRec -> SetDebugReconstruction(kTRUE);
 346
 347         }
 348
 349       //========== Event Number>
 350       //      if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. )
 351         cout <<  "Reconstruct > Event is " << ievent << endl ;
 352
 353       //=========== Connects the various Tree's for evt
 354       gAlice->GetEvent(ievent);
 355
 356       //=========== Gets the Digit TTree
 357       gAlice->TreeD()->GetEvent(0) ;
 358
 359       //=========== Do the reconstruction
 360       fPHOS->Reconstruction(fRec);
 361     }
 362
 363     fClu->Delete();
 364     fClu=0 ;
 365     fTrs->Delete();
 366     fTrs = 0 ;
 367     fPID->Delete();
 368     fPID = 0 ;
 369     fRec->Delete();
 370     fRec = 0 ;
 371
 372 }
 373 //____________________________________________________________________________
 374  void AliPHOSAnalyze::InvariantMass(Int_t Nevents )
 375 {
 376   // Calculates Real and Mixed invariant mass distributions
 377   Int_t NMixedEvents = 4 ; //# of events used for calculation of 'mixed' distribution
 378   Int_t MixedLoops = (Int_t )TMath::Ceil(Nevents/NMixedEvents) ;
 379
 380   //========== Booking Histograms
 381   TH2D * hRealEM   = new TH2D("hRealEM",   "Real for EM particles",      250,0.,1.,40,0.,4.) ;
 382   TH2D * hRealPhot = new TH2D("hRealPhot", "Real for kPhoton particles", 250,0.,1.,40,0.,4.) ;
 383   TH2D * hMixedEM  = new TH2D("hMixedEM",  "Mixed for EM particles",     250,0.,1.,40,0.,4.) ;
 384   TH2D * hMixedPhot= new TH2D("hMixedPhot","Mixed for kPhoton particles",250,0.,1.,40,0.,4.) ;
 385
 386   Int_t ievent;
 387   Int_t EventInMixedLoop ;
 388
 389   Int_t NRecParticles[NMixedEvents] ;
 390
 391   AliPHOSRecParticle::RecParticlesList * AllRecParticleList  = new TClonesArray("AliPHOSRecParticle", NMixedEvents*1000) ;
 392
 393   for(EventInMixedLoop = 0; EventInMixedLoop < MixedLoops; EventInMixedLoop++  ){
 394     Int_t iRecPhot = 0 ;
 395
 396     for ( ievent=0; ievent < NMixedEvents; ievent++){
 397
 398       Int_t AbsEventNumber = EventInMixedLoop*NMixedEvents + ievent ;
 399
 400       //=========== Connects the various Tree's for evt
 401       gAlice->GetEvent(AbsEventNumber);
 402
 403       //=========== Get the Digit Tree
 404       gAlice->TreeD()->GetEvent(0) ;
 405
 406       //========== Creating branches ===================================
 407
 408       AliPHOSRecParticle::RecParticlesList ** RecParticleList  = fPHOS->RecParticles() ;
 409       if( (*RecParticleList) )
 410         (*RecParticleList)->Clear() ;
 411       gAlice->TreeR()->SetBranchAddress( "PHOSRP", RecParticleList ) ;
 412
 413       //=========== Gets the Reconstraction TTree
 414       gAlice->TreeR()->GetEvent(0) ;
 415
 416       AliPHOSRecParticle * RecParticle ;
 417       Int_t iRecParticle ;
 418       for(iRecParticle = 0; iRecParticle < (*RecParticleList)->GetEntries() ;iRecParticle++ )
 419         {
 420           RecParticle = (AliPHOSRecParticle *) (*RecParticleList)->At(iRecParticle) ;
 421           if((RecParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)||
 422              (RecParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM)){
 423             new( (*AllRecParticleList)[iRecPhot] ) AliPHOSRecParticle(*RecParticle) ;
 424             iRecPhot++;
 425           }
 426         }
 427
 428         NRecParticles[ievent] = iRecPhot-1 ;
 429     }
 430
 431
 432     //Now calculate invariant mass:
 433     Int_t irp1,irp2 ;
 434     Int_t NCurEvent = 0 ;
 435
 436     for(irp1 = 0; irp1 < AllRecParticleList->GetEntries()-1; irp1++){
 437       AliPHOSRecParticle * rp1 = (AliPHOSRecParticle *)AllRecParticleList->At(irp1) ;
 438
 439       for(irp2 = irp1+1; irp2 < AllRecParticleList->GetEntries(); irp2++){
 440         AliPHOSRecParticle * rp2 = (AliPHOSRecParticle *)AllRecParticleList->At(irp2) ;
 441
 442         Double_t InvMass ;
 443         InvMass = (rp1->Energy()+rp2->Energy())*(rp1->Energy()+rp2->Energy())-
 444           (rp1->Px()+rp2->Px())*(rp1->Px()+rp2->Px())-
 445           (rp1->Py()+rp2->Py())*(rp1->Py()+rp2->Py())-
 446           (rp1->Pz()+rp2->Pz())*(rp1->Pz()+rp2->Pz()) ;
 447
 448         if(InvMass> 0)
 449           InvMass = TMath::Sqrt(InvMass);
 450
 451         Double_t Pt ;
 452         Pt = TMath::Sqrt((rp1->Px()+rp2->Px() )*( rp1->Px()+rp2->Px() ) +(rp1->Py()+rp2->Py())*(rp1->Py()+rp2->Py()));
 453
 454         if(irp1 > NRecParticles[NCurEvent])
 455           NCurEvent++;
 456
 457         if(irp2 <= NRecParticles[NCurEvent]){ //'Real' event
 458           hRealEM->Fill(InvMass,Pt);
 459           if((rp1->GetType() == AliPHOSFastRecParticle::kGAMMA)&&(rp2->GetType() == AliPHOSFastRecParticle::kGAMMA))
 460             hRealPhot->Fill(InvMass,Pt);
 461         }
 462         else{
 463           hMixedEM->Fill(InvMass,Pt);
 464           if((rp1->GetType() == AliPHOSFastRecParticle::kGAMMA)&&(rp2->GetType() == AliPHOSFastRecParticle::kGAMMA))
 465             hMixedPhot->Fill(InvMass,Pt);
 466         } //real-mixed
 467
 468       } //loop over second rp
 469     }//loop over first rp
 470
 471
 472     AllRecParticleList->Delete() ;
 473   } //Loop over events
 474
 475   delete AllRecParticleList ;
 476
 477   //writing output
 478   TFile output("invmass.root","RECREATE");
 479   output.cd();
 480
 481   hRealEM->Write() ;
 482   hRealPhot->Write() ;
 483   hMixedEM->Write() ;
 484   hMixedPhot->Write() ;
 485
 486   output.Write();
 487   output.Close();
 488
 489 }
 490
 491 //____________________________________________________________________________
 492  void AliPHOSAnalyze::AnalyzeResolutions(Int_t Nevents )
 493 {
 494   // analyzes Nevents events and calculate Energy and Position resolution as well as
 495   // probaility of correct indentifiing of the incident particle
 496
 497   //========== Booking Histograms
 498   cout << "AnalyzeResolutions > " << "Booking Histograms" << endl ;
 499   BookResolutionHistograms();
 500
 501   Int_t Counter[9][5] ;
 502   Int_t i1,i2,TotalInd = 0 ;
 503   for(i1 = 0; i1<9; i1++)
 504     for(i2 = 0; i2<5; i2++)
 505       Counter[i1][i2] = 0 ;
 506
 507   Int_t TotalPrimary = 0 ;
 508   Int_t TotalRecPart = 0 ;
 509   Int_t TotalRPwithPrim = 0 ;
 510   Int_t ievent;
 511
 512   cout << "Start Analysing"<< endl ;
 513   for ( ievent=0; ievent<Nevents; ievent++)
 514     {
 515
 516       //========== Event Number>
 517       //      if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. )
 518         cout <<  "AnalyzeResolutions > " << "Event is " << ievent << endl ;
 519
 520       //=========== Connects the various Tree's for evt
 521       gAlice->GetEvent(ievent);
 522
 523       //=========== Gets the Kine TTree
 524       gAlice->TreeK()->GetEvent(0) ;
 525
 526       //=========== Gets the list of Primari Particles
 527       TClonesArray * PrimaryList  = gAlice->Particles();
 528
 529       TParticle * Primary ;
 530       Int_t iPrimary ;
 531       for ( iPrimary = 0 ; iPrimary < PrimaryList->GetEntries() ; iPrimary++)
 532         {
 533           Primary = (TParticle*)PrimaryList->UncheckedAt(iPrimary) ;
 534           Int_t PrimaryType = Primary->GetPdgCode() ;
 535           if( PrimaryType == 22 ) {
 536             Int_t ModuleNumber ;
 537             Double_t PrimX, PrimZ ;
 538             fGeom->ImpactOnEmc(Primary->Theta(), Primary->Phi(), ModuleNumber, PrimX, PrimZ) ;
 539             if(ModuleNumber){
 540               fhPrimary->Fill(Primary->Energy()) ;
 541               if(Primary->Energy() > 0.3)
 542                 TotalPrimary++ ;
 543             }
 544           }
 545         }
 546
 547       //=========== Get the Digit Tree
 548       gAlice->TreeD()->GetEvent(0) ;
 549
 550       //========== Creating branches ===================================
 551       AliPHOSRecPoint::RecPointsList ** EmcRecPoints =  fPHOS->EmcRecPoints() ;
 552       gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP", EmcRecPoints ) ;
 553
 554       AliPHOSRecPoint::RecPointsList ** PpsdRecPoints = fPHOS->PpsdRecPoints() ;
 555       gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", PpsdRecPoints ) ;
 556
 557       AliPHOSTrackSegment::TrackSegmentsList **  TrackSegmentsList = fPHOS->TrackSegments() ;
 558       if( (*TrackSegmentsList) )
 559         (*TrackSegmentsList)->Clear() ;
 560       gAlice->TreeR()->SetBranchAddress( "PHOSTS", TrackSegmentsList ) ;
 561
 562       AliPHOSRecParticle::RecParticlesList ** RecParticleList  = fPHOS->RecParticles() ;
 563       if( (*RecParticleList) )
 564         (*RecParticleList)->Clear() ;
 565       gAlice->TreeR()->SetBranchAddress( "PHOSRP", RecParticleList ) ;
 566
 567       //=========== Gets the Reconstraction TTree
 568       gAlice->TreeR()->GetEvent(0) ;
 569
 570       AliPHOSRecParticle * RecParticle ;
 571       Int_t iRecParticle ;
 572       for(iRecParticle = 0; iRecParticle < (*RecParticleList)->GetEntries() ;iRecParticle++ )
 573         {
 574           RecParticle = (AliPHOSRecParticle *) (*RecParticleList)->At(iRecParticle) ;
 575           fhAllRP->Fill(CorrectEnergy(RecParticle->Energy())) ;
 576
 577           Int_t ModuleNumberRec ;
 578           Double_t RecX, RecZ ;
 579           fGeom->ImpactOnEmc(RecParticle->Theta(), RecParticle->Phi(), ModuleNumberRec, RecX, RecZ) ;
 580
 581           Double_t MinDistance = 5. ;
 582           Int_t ClosestPrimary = -1 ;
 583
 584           Int_t numberofprimaries ;
 585           Int_t * listofprimaries  = RecParticle->GetPrimaries(numberofprimaries)  ;
 586           Int_t index ;
 587           TParticle * Primary ;
 588           Double_t Distance = MinDistance ;
 589           for ( index = 0 ; index < numberofprimaries ; index++){
 590             Primary = (TParticle*)PrimaryList->UncheckedAt(listofprimaries[index]) ;
 591             Int_t ModuleNumber ;
 592             Double_t PrimX, PrimZ ;
 593             fGeom->ImpactOnEmc(Primary->Theta(), Primary->Phi(), ModuleNumber, PrimX, PrimZ) ;
 594             if(ModuleNumberRec == ModuleNumber)
 595               Distance = TMath::Sqrt((RecX-PrimX)*(RecX-PrimX)+(RecZ-PrimZ)*(RecZ-PrimZ) ) ;
 596             if(MinDistance > Distance)
 597               {
 598                 MinDistance = Distance ;
 599                 ClosestPrimary = listofprimaries[index] ;
 600               }
 601           }
 602           TotalRecPart++ ;
 603
 604           if(ClosestPrimary >=0 ){
 605             TotalRPwithPrim++;
 606
 607             Int_t PrimaryType = ((TParticle *)PrimaryList->At(ClosestPrimary))->GetPdgCode() ;
 608             TParticlePDG* PDGparticle = ((TParticle *)PrimaryList->At(ClosestPrimary))->GetPDG();
 609             Double_t charge =  PDGparticle->Charge() ;
 610             Int_t PrimaryCode ;
 611             switch(PrimaryType)
 612               {
 613               case 22:
 614                 PrimaryCode = 0;  //Photon
 615                 fhAllEnergy->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(), RecParticle->Energy()) ;
 616                 fhAllPosition->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(),MinDistance) ;
 617                 break;
 618               case 11 :
 619                 PrimaryCode = 1;  //Electron
 620                 break;
 621               case -11 :
 622                 PrimaryCode = 1;  //positron
 623                 break;
 624               case 321 :
 625                 PrimaryCode = 4;  //K+
 626                 break;
 627               case -321 :
 628                 PrimaryCode = 4;  //K-
 629                 break;
 630               case 310 :
 631                 PrimaryCode = 4;  //K0s
 632                 break;
 633               case 130 :
 634                 PrimaryCode = 4;  //K0l
 635                 break;
 636               default:
 637                 if(charge)
 638                   PrimaryCode = 2; //Charged hadron
 639                 else
 640                   PrimaryCode = 3; //Neutral hadron
 641                 break;
 642               }
 643
 644             switch(RecParticle->GetType())
 645               {
 646               case AliPHOSFastRecParticle::kGAMMA:
 647                 if(PrimaryType == 22){
 648                   fhPhotEnergy->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(), RecParticle->Energy() ) ;
 649                   fhEMEnergy->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(), RecParticle->Energy() ) ;
 650                   fhPPSDEnergy->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(), RecParticle->Energy() ) ;
 651
 652                   fhPhotPosition->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(),MinDistance) ;
 653                   fhEMPosition->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(),MinDistance) ;
 654                   fhPPSDPosition->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(),MinDistance) ;
 655
 656                   fhPhotReg->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 657                   fhPhotEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 658                   fhPhotPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 659
 660                   fhPhotPhot->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 661                 }
 662                 if(PrimaryType == 2112){ //neutron
 663                   fhNReg->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 664                   fhNEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 665                   fhNPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 666                 }
 667
 668                 if(PrimaryType == -2112){ //neutron ~
 669                   fhNBarReg->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 670                   fhNBarEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 671                   fhNBarPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 672
 673                 }
 674                 if(PrimaryCode == 2){
 675                   fhChargedReg->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 676                   fhChargedEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 677                   fhChargedPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 678                 }
 679
 680                 fhAllReg->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 681                 fhAllEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 682                 fhAllPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 683                 fhShape->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 684                 fhVeto->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 685                 fhPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 686                 Counter[0][PrimaryCode]++;
 687                 break;
 688               case  AliPHOSFastRecParticle::kELECTRON:
 689                 if(PrimaryType == 22){
 690                   fhPhotElec->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 691                   fhEMEnergy->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(), RecParticle->Energy() ) ;
 692                   fhEMPosition->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(),MinDistance) ;
 693                   fhPhotEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 694                   fhPhotPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 695                 }
 696                 if(PrimaryType == 2112){ //neutron
 697                   fhNEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 698                   fhNPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 699                 }
 700
 701                 if(PrimaryType == -2112){ //neutron ~
 702                   fhNBarEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 703                   fhNBarPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 704
 705                 }
 706                 if(PrimaryCode == 2){
 707                   fhChargedEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 708                   fhChargedPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 709                 }
 710
 711                 fhAllEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 712                 fhAllPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 713                 fhShape->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 714                 fhPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 715                 Counter[1][PrimaryCode]++;
 716                 break;
 717               case  AliPHOSFastRecParticle::kNEUTRALHA:
 718                 if(PrimaryType == 22)
 719                   fhPhotNeuH->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 720
 721                 fhVeto->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 722                 Counter[2][PrimaryCode]++;
 723                 break ;
 724               case  AliPHOSFastRecParticle::kNEUTRALEM:
 725                 if(PrimaryType == 22){
 726                   fhEMEnergy->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(),RecParticle->Energy() ) ;
 727                   fhEMPosition->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(),MinDistance ) ;
 728
 729                   fhPhotNuEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 730                   fhPhotEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 731                 }
 732                 if(PrimaryType == 2112) //neutron
 733                   fhNEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 734
 735                 if(PrimaryType == -2112) //neutron ~
 736                   fhNBarEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 737
 738                 if(PrimaryCode == 2)
 739                   fhChargedEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 740
 741                 fhAllEM->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 742                 fhShape->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 743                 fhVeto->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 744
 745                 Counter[3][PrimaryCode]++;
 746                 break ;
 747               case  AliPHOSFastRecParticle::kCHARGEDHA:
 748                 if(PrimaryType == 22) //photon
 749                   fhPhotChHa->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 750
 751                 Counter[4][PrimaryCode]++ ;
 752                 break ;
 753               case  AliPHOSFastRecParticle::kGAMMAHA:
 754                   if(PrimaryType == 22){ //photon
 755                     fhPhotGaHa->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 756                     fhPPSDEnergy->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(), RecParticle->Energy() ) ;
 757                     fhPPSDPosition->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(),MinDistance) ;
 758                     fhPhotPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 759                   }
 760                   if(PrimaryType == 2112){ //neutron
 761                     fhNPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 762                   }
 763
 764                   if(PrimaryType == -2112){ //neutron ~
 765                     fhNBarPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 766                   }
 767                   if(PrimaryCode == 2){
 768                     fhChargedPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 769                   }
 770
 771                   fhAllPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 772                   fhVeto->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 773                   fhPPSD->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 774
 775                   Counter[5][PrimaryCode]++ ;
 776                   break ;
 777               case  AliPHOSFastRecParticle::kABSURDEM:
 778                 Counter[6][PrimaryCode]++ ;
 779                 fhShape->Fill(CorrectEnergy(RecParticle->Energy()) ) ;
 780                 break;
 781               case  AliPHOSFastRecParticle::kABSURDHA:
 782                 Counter[7][PrimaryCode]++ ;
 783                 break;
 784               default:
 785                 Counter[8][PrimaryCode]++ ;
 786                 break;
 787               }
 788           }
 789         }
 790     }   // endfor
 791   SaveHistograms();
 792   cout << "Resolutions: Analyzed " << Nevents << " event(s)" << endl ;
 793   cout << "Resolutions: Total primary       " << TotalPrimary << endl ;
 794   cout << "Resoluitons: Total reconstracted " << TotalRecPart << endl ;
 795   cout << "TotalReconstructed with Primarie " << TotalRPwithPrim << endl ;
 796   cout << "                        Primary:   Photon   Electron   Ch. Hadr.  Neutr. Hadr  Kaons" << endl ;
 797   cout << "             Detected as photon       " << Counter[0][0] << "          " << Counter[0][1] << "          " << Counter[0][2] << "          " <<Counter[0][3] << "          " << Counter[0][4] << endl ;
 798   cout << "           Detected as electron       " << Counter[1][0] << "          " << Counter[1][1] << "          " << Counter[1][2] << "          " <<Counter[1][3] << "          " << Counter[1][4] << endl ;
 799   cout << "     Detected as neutral hadron       " << Counter[2][0] << "          " << Counter[2][1] << "          " << Counter[2][2] << "          " <<Counter[2][3] << "          " << Counter[2][4] << endl ;
 800   cout << "         Detected as neutral EM       " << Counter[3][0] << "          " << Counter[3][1] << "          " << Counter[3][2] << "          " <<Counter[3][3] << "          " << Counter[3][4] << endl ;
 801   cout << "     Detected as charged hadron       " << Counter[4][0] << "          " << Counter[4][1] << "          " << Counter[4][2] << "          " <<Counter[4][3] << "          " << Counter[4][4] << endl ;
 802   cout << "       Detected as gamma-hadron       " << Counter[5][0] << "          " << Counter[5][1] << "          " << Counter[5][2] << "          " <<Counter[5][3] << "          " << Counter[5][4] << endl ;
 803   cout << "          Detected as Absurd EM       " << Counter[6][0] << "          " << Counter[6][1] << "          " << Counter[6][2] << "          " <<Counter[6][3] << "          " << Counter[6][4] << endl ;
 804   cout << "      Detected as absurd hadron       " << Counter[7][0] << "          " << Counter[7][1] << "          " << Counter[7][2] << "          " <<Counter[7][3] << "          " << Counter[7][4] << endl ;
 805   cout << "          Detected as undefined       " << Counter[8][0] << "          " << Counter[8][1] << "          " << Counter[8][2] << "          " <<Counter[8][3] << "          " << Counter[8][4] << endl ;
 806
 807       for(i1 = 0; i1<9; i1++)
 808         for(i2 = 0; i2<5; i2++)
 809           TotalInd+=Counter[i1][i2] ;
 810       cout << "Indentified particles            " << TotalInd << endl ;
 811
 812 }           // endfunction
 813
 814
 815 //____________________________________________________________________________
 816 void  AliPHOSAnalyze::BookingHistograms()
 817 {
 818   // Books the histograms where the results of the analysis are stored (to be changed)
 819
 820   delete fhEmcDigit  ;
 821   delete fhVetoDigit  ;
 822   delete fhConvertorDigit   ;
 823   delete  fhEmcCluster   ;
 824   delete fhVetoCluster   ;
 825   delete fhConvertorCluster  ;
 826   delete fhConvertorEmc  ;
 827
 828   fhEmcDigit                = new TH1F("hEmcDigit",      "hEmcDigit",         1000,  0. ,  25.);
 829   fhVetoDigit               = new TH1F("hVetoDigit",     "hVetoDigit",         500,  0. ,  3.e-5);
 830   fhConvertorDigit          = new TH1F("hConvertorDigit","hConvertorDigit",    500,  0. ,  3.e-5);
 831   fhEmcCluster              = new TH1F("hEmcCluster",    "hEmcCluster",       1000,  0. ,  30.);
 832   fhVetoCluster             = new TH1F("hVetoCluster",   "hVetoCluster",       500,  0. ,  3.e-5);
 833   fhConvertorCluster        = new TH1F("hConvertorCluster","hConvertorCluster",500,  0. ,  3.e-5);
 834   fhConvertorEmc            = new TH2F("hConvertorEmc",  "hConvertorEmc",      200,  1. ,  3., 200, 0., 3.e-5);
 835
 836 }
 837 //____________________________________________________________________________
 838 void  AliPHOSAnalyze::BookResolutionHistograms()
 839 {
 840   // Books the histograms where the results of the Resolution analysis are stored
 841
 842 //   if(fhAllEnergy)
 843 //     delete fhAllEnergy ;
 844 //   if(fhPhotEnergy)
 845 //     delete fhPhotEnergy ;
 846 //   if(fhEMEnergy)
 847 //     delete fhEMEnergy ;
 848 //   if(fhPPSDEnergy)
 849 //     delete fhPPSDEnergy ;
 850
 851
 852   fhAllEnergy  = new TH2F("hAllEnergy",  "Energy of any RP with primary photon",100, 0., 5., 100, 0., 5.);
 853   fhPhotEnergy = new TH2F("hPhotEnergy", "Energy of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.);
 854   fhEMEnergy   = new TH2F("hEMEnergy",   "Energy of EM with primary photon",    100, 0., 5., 100, 0., 5.);
 855   fhPPSDEnergy = new TH2F("hPPSDEnergy", "Energy of PPSD with primary photon",  100, 0., 5., 100, 0., 5.);
 856
 857 //   if(fhAllPosition)
 858 //     delete fhAllPosition ;
 859 //   if(fhPhotPosition)
 860 //     delete fhPhotPosition ;
 861 //   if(fhEMPosition)
 862 //     delete fhEMPosition ;
 863 //   if(fhPPSDPosition)
 864 //     delete fhPPSDPosition ;
 865
 866
 867   fhAllPosition  = new TH2F("hAllPosition",  "Position of any RP with primary photon",100, 0., 5., 100, 0., 5.);
 868   fhPhotPosition = new TH2F("hPhotPosition", "Position of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.);
 869   fhEMPosition   = new TH2F("hEMPosition",   "Position of EM with primary photon",    100, 0., 5., 100, 0., 5.);
 870   fhPPSDPosition = new TH2F("hPPSDPosition", "Position of PPSD with primary photon",  100, 0., 5., 100, 0., 5.);
 871
 872 //   if(fhAllReg)
 873 //     delete fhAllReg ;
 874 //   if(fhPhotReg)
 875 //     delete fhPhotReg ;
 876 //   if(fhNReg)
 877 //     delete fhNReg ;
 878 //   if(fhNBarReg)
 879 //     delete fhNBarReg ;
 880 //   if(fhChargedReg)
 881 //     delete fhChargedReg ;
 882
 883   fhAllReg    = new TH1F("hAllReg",    "All primaries registered as photon",  100, 0., 5.);
 884   fhPhotReg   = new TH1F("hPhotReg",   "Photon registered as photon",         100, 0., 5.);
 885   fhNReg      = new TH1F("hNReg",      "N registered as photon",              100, 0., 5.);
 886   fhNBarReg   = new TH1F("hNBarReg",   "NBar registered as photon",           100, 0., 5.);
 887   fhChargedReg= new TH1F("hChargedReg", "Charged hadron registered as photon",100, 0., 5.);
 888
 889 //   if(fhAllEM)
 890 //     delete fhAllEM ;
 891 //   if(fhPhotEM)
 892 //     delete fhPhotEM ;
 893 //   if(fhNEM)
 894 //     delete fhNEM ;
 895 //   if(fhNBarEM)
 896 //     delete fhNBarEM ;
 897 //   if(fhChargedEM)
 898 //     delete fhChargedEM ;
 899
 900   fhAllEM    = new TH1F("hAllEM",    "All primary registered as EM",100, 0., 5.);
 901   fhPhotEM   = new TH1F("hPhotEM",   "Photon registered as EM", 100, 0., 5.);
 902   fhNEM      = new TH1F("hNEM",      "N registered as EM",      100, 0., 5.);
 903   fhNBarEM   = new TH1F("hNBarEM",   "NBar registered as EM",   100, 0., 5.);
 904   fhChargedEM= new TH1F("hChargedEM","Charged registered as EM",100, 0., 5.);
 905
 906 //   if(fhAllPPSD)
 907 //     delete fhAllPPSD ;
 908 //   if(fhPhotPPSD)
 909 //     delete fhPhotPPSD ;
 910 //   if(fhNPPSD)
 911 //     delete fhNPPSD ;
 912 //   if(fhNBarPPSD)
 913 //     delete fhNBarPPSD ;
 914 //   if(fhChargedPPSD)
 915 //     delete fhChargedPPSD ;
 916
 917   fhAllPPSD    = new TH1F("hAllPPSD",    "All primary registered as PPSD",100, 0., 5.);
 918   fhPhotPPSD   = new TH1F("hPhotPPSD",   "Photon registered as PPSD", 100, 0., 5.);
 919   fhNPPSD      = new TH1F("hNPPSD",      "N registered as PPSD",      100, 0., 5.);
 920   fhNBarPPSD   = new TH1F("hNBarPPSD",   "NBar registered as PPSD",   100, 0., 5.);
 921   fhChargedPPSD= new TH1F("hChargedPPSD","Charged registered as PPSD",100, 0., 5.);
 922
 923 //   if(fhPrimary)
 924 //     delete fhPrimary ;
 925   fhPrimary= new TH1F("hPrimary", "hPrimary",  100, 0., 5.);
 926
 927 //   if(fhAllRP)
 928 //     delete fhAllRP ;
 929 //   if(fhVeto)
 930 //     delete fhVeto ;
 931 //   if(fhShape)
 932 //     delete fhShape ;
 933 //   if(fhPPSD)
 934 //     delete fhPPSD ;
 935
 936   fhAllRP = new TH1F("hAllRP","All Reconstructed particles",  100, 0., 5.);
 937   fhVeto  = new TH1F("hVeto", "All uncharged particles",      100, 0., 5.);
 938   fhShape = new TH1F("hShape","All particles with EM shaower",100, 0., 5.);
 939   fhPPSD  = new TH1F("hPPSD", "All PPSD photon particles",    100, 0., 5.);
 940
 941
 942 //   if(fhPhotPhot)
 943 //     delete fhPhotPhot ;
 944 //   if(fhPhotElec)
 945 //     delete fhPhotElec ;
 946 //   if(fhPhotNeuH)
 947 //     delete fhPhotNeuH ;
 948 //   if(fhPhotNuEM)
 949 //     delete fhPhotNuEM ;
 950 //   if(fhPhotChHa)
 951 //     delete fhPhotChHa ;
 952 //   if(fhPhotGaHa)
 953 //     delete fhPhotGaHa ;
 954
 955   fhPhotPhot = new TH1F("hPhotPhot","hPhotPhot", 100, 0., 5.);   //Photon registered as photon
 956   fhPhotElec = new TH1F("hPhotElec","hPhotElec", 100, 0., 5.);   //Photon registered as Electron
 957   fhPhotNeuH = new TH1F("hPhotNeuH","hPhotNeuH", 100, 0., 5.);   //Photon registered as Neutral Hadron
 958   fhPhotNuEM = new TH1F("hPhotNuEM","hPhotNuEM", 100, 0., 5.);   //Photon registered as Neutral EM
 959   fhPhotChHa = new TH1F("hPhotChHa","hPhotChHa", 100, 0., 5.);   //Photon registered as Charged Hadron
 960   fhPhotGaHa = new TH1F("hPhotGaHa","hPhotGaHa", 100, 0., 5.);   //Photon registered as Gamma-Hadron
 961
 962
 963 }
 964 //____________________________________________________________________________
 965 Bool_t AliPHOSAnalyze::Init(Int_t evt)
 966 {
 967   // Do a few initializations: open the root file
 968   //                           get the AliRun object
 969   //                           defines the clusterizer, tracksegment maker and particle identifier
 970   //                           sets the associated parameters
 971
 972   Bool_t ok = kTRUE ;
 973
 974    //========== Open galice root file
 975
 976   if ( fRootFile == 0 ) {
 977     Text_t * name  = new Text_t[80] ;
 978     cout << "AnalyzeOneEvent > Enter file root file name : " ;
 979     cin >> name ;
 980     Bool_t ok = OpenRootFile(name) ;
 981     if ( !ok )
 982       cout << " AliPHOSAnalyze > Error opening " << name << endl ;
 983     else {
 984       //========== Get AliRun object from file
 985
 986       gAlice = (AliRun*) fRootFile->Get("gAlice") ;
 987
 988       //=========== Get the PHOS object and associated geometry from the file
 989
 990       fPHOS  = (AliPHOSv1 *)gAlice->GetDetector("PHOS") ;
 991       fGeom = fPHOS->GetGeometry();
 992       //      fGeom  = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() );
 993
 994     } // else !ok
 995   } // if fRootFile
 996
 997   if ( ok ) {
 998
 999     //========== Create the Clusterizer
1000
1001     fClu =  new AliPHOSClusterizerv1() ;
1002     fClu->SetEmcEnergyThreshold(0.030) ;
1003     fClu->SetEmcClusteringThreshold(0.20) ;
1004     fClu->SetPpsdEnergyThreshold    (0.0000002) ;
1005     fClu->SetPpsdClusteringThreshold(0.0000001) ;
1006     fClu->SetLocalMaxCut(0.03) ;
1007     fClu->SetCalibrationParameters(0., 0.00000001) ;
1008     cout <<  "AnalyzeOneEvent > using clusterizer " << fClu->GetName() << endl ;
1009     fClu->PrintParameters() ;
1010
1011     //========== Creates the track segment maker
1012
1013     fTrs = new AliPHOSTrackSegmentMakerv1() ;
1014     cout <<  "AnalyzeOneEvent > using tack segment maker " << fTrs->GetName() << endl ;
1015     //   fTrs->UnsetUnfoldFlag() ;
1016
1017     //========== Creates the particle identifier
1018
1019     fPID = new AliPHOSPIDv1() ;
1020     cout <<  "AnalyzeOneEvent > using particle identifier " << fPID->GetName() << endl ;
1021     //fPID->SetShowerProfileCuts(Float_t l1m, Float_t l1M, Float_t l2m, Float_t l2M) ;
1022     fPID->SetShowerProfileCuts(0.7, 2.0 , 0.6 , 1.5) ;
1023
1024     //========== Creates the Reconstructioner
1025
1026     fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;
1027     fRec -> SetDebugReconstruction(kFALSE);
1028
1029     //=========== Connect the various Tree's for evt
1030
1031     if ( evt == -999 ) {
1032       cout <<  "AnalyzeOneEvent > Enter event number : " ;
1033       cin >> evt ;
1034       cout << evt << endl ;
1035     }
1036     fEvt = evt ;
1037     gAlice->GetEvent(evt);
1038
1039     //=========== Get the Digit TTree
1040
1041     gAlice->TreeD()->GetEvent(0) ;
1042
1043   } // ok
1044
1045   return ok ;
1046 }
1047
1048
1049 //____________________________________________________________________________
1050 void AliPHOSAnalyze::DisplayKineEvent(Int_t evt)
1051 {
1052   // Display particles from the Kine Tree in global Alice (theta, phi) coordinates.
1053   // One PHOS module at the time.
1054   // The particle type can be selected.
1055
1056   if (evt == -999)
1057     evt = fEvt ;
1058
1059   Int_t module ;
1060   cout <<  "DisplayKineEvent > which module (1-5,  -1: all) ? " ;
1061   cin >> module ; cout << module << endl ;
1062
1063   Int_t testparticle ;
1064   cout << " 22      : PHOTON " << endl
1065        << " (-)11   : (POSITRON)ELECTRON " << endl
1066        << " (-)2112 : (ANTI)NEUTRON " << endl
1067        << " -999    : Everything else " << endl ;
1068   cout  <<  "DisplayKineEvent > enter PDG particle code to display " ;
1069   cin >> testparticle ; cout << testparticle << endl ;
1070
1071   Text_t histoname[80] ;
1072   sprintf(histoname,"Event %d: Incident particles in module %d", evt, module) ;
1073
1074   Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
1075   fGeom->EmcModuleCoverage(module, tm, tM, pm, pM, AliPHOSGeometry::Degre() ) ;
1076
1077   Double_t theta, phi ;
1078   fGeom->EmcXtalCoverage(theta, phi, AliPHOSGeometry::Degre() ) ;
1079
1080   Int_t tdim = (Int_t)( (tM - tm) / theta ) ;
1081   Int_t pdim = (Int_t)( (pM - pm) / phi ) ;
1082
1083   tm -= theta ;
1084   tM += theta ;
1085   pm -= phi ;
1086   pM += phi ;
1087
1088   TH2F * histoparticle = new TH2F("histoparticle",  histoname,
1089                                           pdim, pm, pM, tdim, tm, tM) ;
1090   histoparticle->SetStats(kFALSE) ;
1091
1092   // Get pointers to Alice Particle TClonesArray
1093
1094   TParticle * particle;
1095   TClonesArray * particlearray  = gAlice->Particles();
1096
1097   Text_t canvasname[80];
1098   sprintf(canvasname,"Particles incident in PHOS/EMC module # %d",module) ;
1099   TCanvas * kinecanvas = new TCanvas("kinecanvas", canvasname, 650, 500) ;
1100
1101   // get the KINE Tree
1102
1103   TTree * kine =  gAlice->TreeK() ;
1104   Stat_t nParticles =  kine->GetEntries() ;
1105   cout << "DisplayKineEvent > events in kine " << nParticles << endl ;
1106
1107   // loop over particles
1108
1109   Double_t kRADDEG = 180. / TMath::Pi() ;
1110   Int_t index1 ;
1111   Int_t nparticlein = 0 ;
1112   for (index1 = 0 ; index1 < nParticles ; index1++){
1113     Int_t nparticle = particlearray->GetEntriesFast() ;
1114     Int_t index2 ;
1115     for( index2 = 0 ; index2 < nparticle ; index2++) {
1116       particle            = (TParticle*)particlearray->UncheckedAt(index2) ;
1117       Int_t  particletype = particle->GetPdgCode() ;
1118       if (testparticle == -999 || testparticle == particletype) {
1119         Double_t phi        = particle->Phi() ;
1120         Double_t theta      = particle->Theta() ;
1121         Int_t mod ;
1122         Double_t x, z ;
1123         fGeom->ImpactOnEmc(theta, phi, mod, z, x) ;
1124         if ( mod == module ) {
1125           nparticlein++ ;
1126           if (particle->Energy() >  fClu->GetEmcClusteringThreshold()  )
1127             histoparticle->Fill(phi*kRADDEG, theta*kRADDEG, particle->Energy() ) ;
1128         }
1129       }
1130     }
1131   }
1132   kinecanvas->Draw() ;
1133   histoparticle->Draw("color") ;
1134   TPaveText *  pavetext = new TPaveText(294, 100, 300, 101);
1135   Text_t text[40] ;
1136   sprintf(text, "Particles: %d ", nparticlein) ;
1137   pavetext->AddText(text) ;
1138   pavetext->Draw() ;
1139   kinecanvas->Update();
1140
1141 }
1142 //____________________________________________________________________________
1143 void AliPHOSAnalyze::DisplayRecParticles()
1144 {
1145   // Display reconstructed particles in global Alice(theta, phi) coordinates.
1146   // One PHOS module at the time.
1147   // Click on symbols indicate the reconstructed particle type.
1148
1149   if (fEvt == -999) {
1150     cout << "DisplayRecParticles > Analyze an event first ... (y/n) " ;
1151     Text_t answer[1] ;
1152     cin >> answer ; cout << answer ;
1153 //     if ( answer == "y" )
1154 //       AnalyzeOneEvent() ;
1155   }
1156     if (fEvt != -999) {
1157
1158       Int_t module ;
1159       cout <<  "DisplayRecParticles > which module (1-5,  -1: all) ? " ;
1160       cin >> module ; cout << module << endl ;
1161       Text_t histoname[80] ;
1162       sprintf(histoname,"Event %d: Reconstructed particles in module %d", fEvt, module) ;
1163       Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
1164       fGeom->EmcModuleCoverage(module, tm, tM, pm, pM, AliPHOSGeometry::Degre() ) ;
1165       Double_t theta, phi ;
1166       fGeom->EmcXtalCoverage(theta, phi, AliPHOSGeometry::Degre() ) ;
1167       Int_t tdim = (Int_t)( (tM - tm) / theta ) ;
1168       Int_t pdim = (Int_t)( (pM - pm) / phi ) ;
1169       tm -= theta ;
1170       tM += theta ;
1171       pm -= phi ;
1172       TH2F * histoRparticle = new TH2F("histoRparticle",  histoname,
1173                                        pdim, pm, pM, tdim, tm, tM) ;
1174       histoRparticle->SetStats(kFALSE) ;
1175       Text_t canvasname[80] ;
1176       sprintf(canvasname, "Reconstructed particles in PHOSmodule # %d", module) ;
1177       TCanvas * rparticlecanvas = new TCanvas("RparticleCanvas", canvasname, 650, 500) ;
1178       AliPHOSRecParticle::RecParticlesList * rpl = *fPHOS->RecParticles() ;
1179       Int_t nRecParticles = rpl->GetEntries() ;
1180       Int_t nRecParticlesInModule = 0 ;
1181       TIter nextRecPart(rpl) ;
1182       AliPHOSRecParticle * rp ;
1183       cout << "DisplayRecParticles > " << nRecParticles << " reconstructed particles " << endl ;
1184       Double_t kRADDEG = 180. / TMath::Pi() ;
1185       while ( (rp = (AliPHOSRecParticle *)nextRecPart() ) ) {
1186         AliPHOSTrackSegment * ts = rp->GetPHOSTrackSegment() ;
1187         if ( ts->GetPHOSMod() == module ) {
1188           Int_t numberofprimaries = 0 ;
1189           Int_t * listofprimaries = 0;
1190           rp->GetPrimaries(numberofprimaries) ;
1191           cout << "Number of primaries = " << numberofprimaries << endl ;
1192           Int_t index ;
1193           for ( index = 0 ; index < numberofprimaries ; index++)
1194             cout << "    primary # " << index << " =  " << listofprimaries[index] << endl ;
1195
1196           nRecParticlesInModule++ ;
1197           Double_t theta = rp->Theta() * kRADDEG ;
1198           Double_t phi   = rp->Phi() * kRADDEG ;
1199           Double_t energy = rp->Energy() ;
1200           histoRparticle->Fill(phi, theta, energy) ;
1201         }
1202       }
1203       histoRparticle->Draw("color") ;
1204
1205       nextRecPart.Reset() ;
1206       while ( (rp = (AliPHOSRecParticle *)nextRecPart() ) ) {
1207         AliPHOSTrackSegment * ts = rp->GetPHOSTrackSegment() ;
1208         if ( ts->GetPHOSMod() == module )
1209           rp->Draw("P") ;
1210       }
1211
1212       Text_t text[80] ;
1213       sprintf(text, "reconstructed particles: %d", nRecParticlesInModule) ;
1214       TPaveText *  pavetext = new TPaveText(292, 100, 300, 101);
1215       pavetext->AddText(text) ;
1216       pavetext->Draw() ;
1217       rparticlecanvas->Update() ;
1218     }
1219 }
1220
1221 //____________________________________________________________________________
1222 void AliPHOSAnalyze::DisplayRecPoints()
1223 {
1224   // Display reconstructed points in local PHOS-module (x, z) coordinates.
1225   // One PHOS module at the time.
1226   // Click on symbols displays the EMC cluster, or PPSD information.
1227
1228   if (fEvt == -999) {
1229     cout << "DisplayRecPoints > Analyze an event first ... (y/n) " ;
1230     Text_t answer[1] ;
1231     cin >> answer ; cout << answer ;
1232 //     if ( answer == "y" )
1233 //       AnalyzeOneEvent() ;
1234   }
1235     if (fEvt != -999) {
1236
1237       Int_t module ;
1238       cout <<  "DisplayRecPoints > which module (1-5,  -1: all) ? " ;
1239       cin >> module ; cout << module << endl ;
1240
1241       Text_t canvasname[80];
1242       sprintf(canvasname,"Digits in PHOS/EMC module # %d",module) ;
1243       TCanvas * modulecanvas = new TCanvas("module", canvasname, 650, 500) ;
1244       modulecanvas->Draw() ;
1245
1246       //=========== Creating 2d-histogram of the PHOS module
1247       // a little bit junkie but is used to test Geom functinalities
1248
1249       Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
1250
1251       fGeom->EmcModuleCoverage(module, tm, tM, pm, pM);
1252       // convert angles into coordinates local to the EMC module of interest
1253
1254       Int_t emcModuleNumber ;
1255       Double_t emcModulexm, emcModulezm ; // minimum local coordinate in a given EMCA module
1256       Double_t emcModulexM, emcModulezM ; // maximum local coordinate in a given EMCA module
1257       fGeom->ImpactOnEmc(tm, pm, emcModuleNumber, emcModulezm, emcModulexm) ;
1258       fGeom->ImpactOnEmc(tM, pM, emcModuleNumber, emcModulezM, emcModulexM) ;
1259       Int_t xdim = (Int_t)( ( emcModulexM - emcModulexm ) / fGeom->GetCrystalSize(0) ) ;
1260       Int_t zdim = (Int_t)( ( emcModulezM - emcModulezm ) / fGeom->GetCrystalSize(2) ) ;
1261       Float_t xmin = emcModulexm - fGeom->GetCrystalSize(0) ;
1262       Float_t xMax = emcModulexM + fGeom->GetCrystalSize(0) ;
1263       Float_t zmin = emcModulezm - fGeom->GetCrystalSize(2) ;
1264       Float_t zMax = emcModulezM + fGeom->GetCrystalSize(2) ;
1265       Text_t histoname[80];
1266       sprintf(histoname,"Event %d: Digits and RecPoints in module %d", fEvt, module) ;
1267       TH2F * hModule = new TH2F("HistoReconstructed", histoname,
1268                                 xdim, xmin, xMax, zdim, zmin, zMax) ;
1269       hModule->SetMaximum(2.0);
1270       hModule->SetMinimum(0.0);
1271       hModule->SetStats(kFALSE);
1272
1273       TIter next(fPHOS->Digits()) ;
1274       Float_t energy, y, z;
1275       Float_t etot=0.;
1276       Int_t relid[4]; Int_t nDigits = 0 ;
1277       AliPHOSDigit * digit ;
1278
1279       // Making 2D histogram of the EMC module
1280       while((digit = (AliPHOSDigit *)next()))
1281         {
1282           fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
1283           if (relid[0] == module && relid[1] == 0)
1284             {
1285               energy = fClu->Calibrate(digit->GetAmp()) ;
1286               cout << "Energy is " << energy << " and threshold is " << fClu->GetEmcEnergyThreshold() << endl;
1287               if (energy >  fClu->GetEmcEnergyThreshold()  ){
1288                 nDigits++ ;
1289                 etot += energy ;
1290                 fGeom->RelPosInModule(relid,y,z) ;
1291                 hModule->Fill(y, z, energy) ;
1292               }
1293             }
1294         }
1295       cout <<"DrawRecPoints >  Found in module "
1296            << module << " " << nDigits << "  digits with total energy " << etot << endl ;
1297       hModule->Draw("col2") ;
1298
1299       //=========== Cluster in module
1300
1301       //      TClonesArray * emcRP = fPHOS->EmcClusters() ;
1302       TObjArray * emcRP = *(fPHOS->EmcRecPoints()) ;
1303
1304       etot = 0.;
1305       Int_t totalnClusters = 0 ;
1306       Int_t nClusters = 0 ;
1307       TIter nextemc(emcRP) ;
1308       AliPHOSEmcRecPoint * emc ;
1309       while((emc = (AliPHOSEmcRecPoint *)nextemc()))
1310         {
1311           //      Int_t numberofprimaries ;
1312           //      Int_t * primariesarray = new Int_t[10] ;
1313           //      emc->GetPrimaries(numberofprimaries, primariesarray) ;
1314           totalnClusters++ ;
1315           if ( emc->GetPHOSMod() == module )
1316             {
1317               nClusters++ ;
1318               energy = emc->GetTotalEnergy() ;
1319               etot+= energy ;
1320               emc->Draw("M") ;
1321             }
1322         }
1323       cout << "DrawRecPoints > Found " << totalnClusters << " EMC Clusters in PHOS" << endl ;
1324       cout << "DrawRecPoints > Found in module " << module << "  " << nClusters << " EMC Clusters " << endl ;
1325       cout << "DrawRecPoints > total energy  " << etot << endl ;
1326
1327       TPaveText *  pavetext = new TPaveText(22, 80, 83, 90);
1328       Text_t text[40] ;
1329       sprintf(text, "digits: %d;  clusters: %d", nDigits, nClusters) ;
1330       pavetext->AddText(text) ;
1331       pavetext->Draw() ;
1332       modulecanvas->Update();
1333
1334       //=========== Cluster in module PPSD Down
1335
1336       //      TClonesArray * ppsdRP = fPHOS->PpsdClusters() ;
1337       TObjArray * ppsdRP = *(fPHOS->PpsdRecPoints() );
1338
1339       etot = 0.;
1340       TIter nextPpsd(ppsdRP) ;
1341       AliPHOSPpsdRecPoint * ppsd ;
1342       while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd()))
1343         {
1344           totalnClusters++ ;
1345           if ( ppsd->GetPHOSMod() == module )
1346             {
1347               nClusters++ ;
1348               energy = ppsd->GetEnergy() ;
1349               etot+=energy ;
1350               if (!ppsd->GetUp()) ppsd->Draw("P") ;
1351             }
1352         }
1353       cout << "DrawRecPoints > Found " << totalnClusters << " Ppsd Down Clusters in PHOS" << endl ;
1354       cout << "DrawRecPoints > Found in module " << module << "  " << nClusters << " Ppsd Down Clusters " << endl ;
1355       cout << "DrawRecPoints > total energy  " << etot << endl ;
1356
1357       //=========== Cluster in module PPSD Up
1358
1359       ppsdRP = *(fPHOS->PpsdRecPoints()) ;
1360
1361       etot = 0.;
1362       TIter nextPpsdUp(ppsdRP) ;
1363       while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsdUp()))
1364         {
1365           totalnClusters++ ;
1366           if ( ppsd->GetPHOSMod() == module )
1367             {
1368               nClusters++ ;
1369               energy = ppsd->GetEnergy() ;
1370               etot+=energy ;
1371               if (ppsd->GetUp()) ppsd->Draw("P") ;
1372             }
1373         }
1374   cout << "DrawRecPoints > Found " << totalnClusters << " Ppsd Up Clusters in PHOS" << endl ;
1375   cout << "DrawRecPoints > Found in module " << module << "  " << nClusters << " Ppsd Up Clusters " << endl ;
1376   cout << "DrawRecPoints > total energy  " << etot << endl ;
1377
1378     } // if !-999
1379 }
1380
1381 //____________________________________________________________________________
1382 void AliPHOSAnalyze::DisplayTrackSegments()
1383 {
1384   // Display track segments in local PHOS-module (x, z) coordinates.
1385   // One PHOS module at the time.
1386   // One symbol per PHOS subsystem: EMC, upper PPSD, lower PPSD.
1387
1388   if (fEvt == -999) {
1389     cout << "DisplayTrackSegments > Analyze an event first ... (y/n) " ;
1390     Text_t answer[1] ;
1391     cin >> answer ; cout << answer ;
1392 //     if ( answer == "y" )
1393 //       AnalyzeOneEvent() ;
1394   }
1395     if (fEvt != -999) {
1396
1397       Int_t module ;
1398       cout <<  "DisplayTrackSegments > which module (1-5,  -1: all) ? " ;
1399       cin >> module ; cout << module << endl ;
1400       //=========== Creating 2d-histogram of the PHOS module
1401       // a little bit junkie but is used to test Geom functinalities
1402
1403       Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
1404
1405       fGeom->EmcModuleCoverage(module, tm, tM, pm, pM);
1406       // convert angles into coordinates local to the EMC module of interest
1407
1408       Int_t emcModuleNumber ;
1409       Double_t emcModulexm, emcModulezm ; // minimum local coordinate in a given EMCA module
1410       Double_t emcModulexM, emcModulezM ; // maximum local coordinate in a given EMCA module
1411       fGeom->ImpactOnEmc(tm, pm, emcModuleNumber, emcModulezm, emcModulexm) ;
1412       fGeom->ImpactOnEmc(tM, pM, emcModuleNumber, emcModulezM, emcModulexM) ;
1413       Int_t xdim = (Int_t)( ( emcModulexM - emcModulexm ) / fGeom->GetCrystalSize(0) ) ;
1414       Int_t zdim = (Int_t)( ( emcModulezM - emcModulezm ) / fGeom->GetCrystalSize(2) ) ;
1415       Float_t xmin = emcModulexm - fGeom->GetCrystalSize(0) ;
1416       Float_t xMax = emcModulexM + fGeom->GetCrystalSize(0) ;
1417       Float_t zmin = emcModulezm - fGeom->GetCrystalSize(2) ;
1418       Float_t zMax = emcModulezM + fGeom->GetCrystalSize(2) ;
1419       Text_t histoname[80];
1420       sprintf(histoname,"Event %d: Track Segments in module %d", fEvt, module) ;
1421       TH2F * histotrack = new TH2F("histotrack",  histoname,
1422                                    xdim, xmin, xMax, zdim, zmin, zMax) ;
1423       histotrack->SetStats(kFALSE);
1424       Text_t canvasname[80];
1425       sprintf(canvasname,"Track segments in PHOS/EMC-PPSD module # %d", module) ;
1426       TCanvas * trackcanvas = new TCanvas("TrackSegmentCanvas", canvasname, 650, 500) ;
1427       histotrack->Draw() ;
1428
1429       AliPHOSTrackSegment::TrackSegmentsList * trsegl = *(fPHOS->TrackSegments()) ;
1430       AliPHOSTrackSegment * trseg ;
1431
1432       Int_t nTrackSegments = trsegl->GetEntries() ;
1433       Int_t index ;
1434       Float_t etot = 0 ;
1435       Int_t nTrackSegmentsInModule = 0 ;
1436       for(index = 0; index < nTrackSegments ; index++){
1437         trseg = (AliPHOSTrackSegment * )trsegl->At(index) ;
1438         etot+= trseg->GetEnergy() ;
1439         if ( trseg->GetPHOSMod() == module ) {
1440           nTrackSegmentsInModule++ ;
1441           trseg->Draw("P");
1442         }
1443       }
1444       Text_t text[80] ;
1445       sprintf(text, "track segments: %d", nTrackSegmentsInModule) ;
1446       TPaveText *  pavetext = new TPaveText(22, 80, 83, 90);
1447       pavetext->AddText(text) ;
1448       pavetext->Draw() ;
1449       trackcanvas->Update() ;
1450       cout << "DisplayTrackSegments > Found " << trsegl->GetEntries() << " Track segments with total energy "<< etot << endl ;
1451
1452    }
1453 }
1454 //____________________________________________________________________________
1455 Bool_t AliPHOSAnalyze::OpenRootFile(Text_t * name)
1456 {
1457   // Open the root file named "name"
1458
1459   fRootFile   = new TFile(name, "update") ;
1460   return  fRootFile->IsOpen() ;
1461 }
1462 //____________________________________________________________________________
1463 void AliPHOSAnalyze::SaveHistograms()
1464 {
1465   // Saves the histograms in a root file named "name.analyzed"
1466
1467   Text_t outputname[80] ;
1468   sprintf(outputname,"%s.analyzed",fRootFile->GetName());
1469   TFile output(outputname,"RECREATE");
1470   output.cd();
1471
1472   if (fhAllEnergy)
1473     fhAllEnergy->Write() ;
1474   if (fhPhotEnergy)
1475     fhPhotEnergy->Write() ;
1476   if(fhEMEnergy)
1477     fhEMEnergy->Write()  ;
1478   if(fhPPSDEnergy)
1479     fhPPSDEnergy->Write() ;
1480   if(fhAllPosition)
1481     fhAllPosition->Write() ;
1482   if(fhPhotPosition)
1483     fhPhotPosition->Write() ;
1484   if(fhEMPosition)
1485     fhEMPosition->Write() ;
1486   if(fhPPSDPosition)
1487     fhPPSDPosition->Write() ;
1488   if (fhAllReg)
1489     fhAllReg->Write() ;
1490   if (fhPhotReg)
1491     fhPhotReg->Write() ;
1492   if(fhNReg)
1493     fhNReg->Write() ;
1494   if(fhNBarReg)
1495     fhNBarReg->Write() ;
1496   if(fhChargedReg)
1497     fhChargedReg->Write() ;
1498   if (fhAllEM)
1499     fhAllEM->Write() ;
1500   if (fhPhotEM)
1501     fhPhotEM->Write() ;
1502   if(fhNEM)
1503     fhNEM->Write() ;
1504   if(fhNBarEM)
1505     fhNBarEM->Write() ;
1506   if(fhChargedEM)
1507     fhChargedEM->Write() ;
1508   if (fhAllPPSD)
1509     fhAllPPSD->Write() ;
1510   if (fhPhotPPSD)
1511     fhPhotPPSD->Write() ;
1512   if(fhNPPSD)
1513     fhNPPSD->Write() ;
1514   if(fhNBarPPSD)
1515     fhNBarPPSD->Write() ;
1516   if(fhChargedPPSD)
1517     fhChargedPPSD->Write() ;
1518   if(fhPrimary)
1519     fhPrimary->Write() ;
1520   if(fhAllRP)
1521     fhAllRP->Write()  ;
1522   if(fhVeto)
1523     fhVeto->Write()  ;
1524   if(fhShape)
1525     fhShape->Write()  ;
1526   if(fhPPSD)
1527     fhPPSD->Write()  ;
1528   if(fhPhotPhot)
1529     fhPhotPhot->Write() ;
1530   if(fhPhotElec)
1531     fhPhotElec->Write() ;
1532   if(fhPhotNeuH)
1533     fhPhotNeuH->Write() ;
1534   if(fhPhotNuEM)
1535     fhPhotNuEM->Write() ;
1536   if(fhPhotNuEM)
1537     fhPhotNuEM->Write() ;
1538   if(fhPhotChHa)
1539     fhPhotChHa->Write() ;
1540   if(fhPhotGaHa)
1541     fhPhotGaHa->Write() ;
1542   if(fhEnergyCorrelations)
1543     fhEnergyCorrelations->Write() ;
1544
1545   output.Write();
1546   output.Close();
1547 }
1548 //____________________________________________________________________________
1549 Float_t AliPHOSAnalyze::CorrectEnergy(Float_t ERecPart)
1550 {
1551   return ERecPart/0.8783 ;
1552 }
1553
1554 //____________________________________________________________________________
1555 void AliPHOSAnalyze::ResetHistograms()
1556 {
1557    fhEnergyCorrelations = 0 ;     //Energy correlations between Eloss in Convertor and PPSD(2)
1558
1559    fhEmcDigit = 0 ;               // Histo of digit energies in the Emc
1560    fhVetoDigit = 0 ;              // Histo of digit energies in the Veto
1561    fhConvertorDigit = 0 ;         // Histo of digit energies in the Convertor
1562    fhEmcCluster = 0 ;             // Histo of Cluster energies in Emc
1563    fhVetoCluster = 0 ;            // Histo of Cluster energies in Veto
1564    fhConvertorCluster = 0 ;       // Histo of Cluster energies in Convertor
1565    fhConvertorEmc = 0 ;           // 2d Convertor versus Emc energies
1566
1567    fhAllEnergy = 0 ;
1568    fhPhotEnergy = 0 ;        // Total spectrum of detected photons
1569    fhEMEnergy = 0 ;         // Spectrum of detected electrons with electron primary
1570    fhPPSDEnergy = 0 ;
1571    fhAllPosition = 0 ;
1572    fhPhotPosition = 0 ;
1573    fhEMPosition = 0 ;
1574    fhPPSDPosition = 0 ;
1575
1576    fhPhotReg = 0 ;
1577    fhAllReg = 0 ;
1578    fhNReg = 0 ;
1579    fhNBarReg = 0 ;
1580    fhChargedReg = 0 ;
1581    fhPhotEM = 0 ;
1582    fhAllEM = 0 ;
1583    fhNEM = 0 ;
1584    fhNBarEM = 0 ;
1585    fhChargedEM = 0 ;
1586    fhPhotPPSD = 0 ;
1587    fhAllPPSD = 0 ;
1588    fhNPPSD = 0 ;
1589    fhNBarPPSD = 0 ;
1590    fhChargedPPSD = 0 ;
1591
1592    fhPrimary = 0 ;
1593
1594    fhPhotPhot = 0 ;
1595    fhPhotElec = 0 ;
1596    fhPhotNeuH = 0 ;
1597    fhPhotNuEM = 0 ;
1598    fhPhotChHa = 0 ;
1599    fhPhotGaHa = 0 ;
1600
1601
1602 }
1603
1604