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 #include "TStyle.h"
  38
  39 // --- Standard library ---
  40
  41 #include <iostream.h>
  42 #include <stdio.h>
  43
  44 // --- AliRoot header files ---
  45
  46 #include "AliRun.h"
  47 #include "AliPHOSAnalyze.h"
  48 #include "AliPHOSClusterizerv1.h"
  49 #include "AliPHOSTrackSegmentMakerv1.h"
  50 #include "AliPHOSPIDv1.h"
  51 #include "AliPHOSReconstructioner.h"
  52 #include "AliPHOSDigit.h"
  53 #include "AliPHOSTrackSegment.h"
  54 #include "AliPHOSRecParticle.h"
  55 #include "AliPHOSIndexToObject.h"
  56 #include "AliPHOSHit.h"
  57 #include "AliPHOSCPVHit.h"
  58 #include "AliPHOSCpvRecPoint.h"
  59
  60 ClassImp(AliPHOSAnalyze)
  61
  62 //____________________________________________________________________________
  63   AliPHOSAnalyze::AliPHOSAnalyze()
  64 {
  65   // default ctor (useless)
  66
  67   fRootFile = 0 ;
  68 }
  69
  70 //____________________________________________________________________________
  71 AliPHOSAnalyze::AliPHOSAnalyze(Text_t * name)
  72 {
  73   // ctor: analyze events from root file "name"
  74
  75   Bool_t ok = OpenRootFile(name)  ;
  76   if ( !ok ) {
  77     cout << " AliPHOSAnalyze > Error opening " << name << endl ;
  78   }
  79   else {
  80       //========== Get AliRun object from file
  81       gAlice = (AliRun*) fRootFile->Get("gAlice") ;
  82
  83       //=========== Get the PHOS object and associated geometry from the file
  84       fPHOS  = (AliPHOSv1 *)gAlice->GetDetector("PHOS") ;
  85       fGeom  = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() );
  86
  87       //========== Initializes the Index to Object converter
  88       fObjGetter = AliPHOSIndexToObject::GetInstance(fPHOS) ;
  89       //========== Current event number
  90       fEvt = -999 ;
  91
  92   }
  93   fDebugLevel = 0;
  94   fClu = 0 ;
  95   fPID = 0 ;
  96   fTrs = 0 ;
  97   fRec = 0 ;
  98   ResetHistograms() ;
  99 }
 100
 101 //____________________________________________________________________________
 102 AliPHOSAnalyze::AliPHOSAnalyze(const AliPHOSAnalyze & ana)
 103 {
 104   // copy ctor
 105   ( (AliPHOSAnalyze &)ana ).Copy(*this) ;
 106 }
 107
 108 //____________________________________________________________________________
 109 void AliPHOSAnalyze::Copy(TObject & obj)
 110 {
 111   // copy an analysis into an other one
 112   TObject::Copy(obj) ;
 113   // I do nothing more because the copy is silly but the Code checkers requires one
 114 }
 115
 116 //____________________________________________________________________________
 117 AliPHOSAnalyze::~AliPHOSAnalyze()
 118 {
 119   // dtor
 120
 121   if(fRootFile->IsOpen()) fRootFile->Close() ;
 122   if(fRootFile)   {delete fRootFile ; fRootFile=0 ;}
 123   if(fPHOS)       {delete fPHOS     ; fPHOS    =0 ;}
 124   if(fClu)        {delete fClu      ; fClu     =0 ;}
 125   if(fPID)        {delete fPID      ; fPID     =0 ;}
 126   if(fRec)        {delete fRec      ; fRec     =0 ;}
 127   if(fTrs)        {delete fTrs      ; fTrs     =0 ;}
 128
 129 }
 130 //____________________________________________________________________________
 131 void AliPHOSAnalyze::DrawRecon(Int_t Nevent,Int_t Nmod){
 132   //Draws pimary particles and reconstructed
 133   //digits, RecPoints, RecPartices etc
 134   //for event Nevent in the module Nmod.
 135
 136   TH2F * digitOccupancy  = new TH2F("digitOccupancy","EMC digits", 64,-71.,71.,64,-71.,71.);
 137   TH2F * emcOccupancy    = new TH2F("emcOccupancy","EMC RecPoints",64,-71.,71.,64,-71.,71.);
 138   TH2F * ppsdUp          = new TH2F("ppsdUp","PPSD Up digits",     128,-71.,71.,128,-71.,71.) ;
 139   TH2F * ppsdUpCl        = new TH2F("ppsdUpCl","PPSD Up RecPoints",128,-71.,71.,128,-71.,71.) ;
 140   TH2F * ppsdLow         = new TH2F("ppsdLow","PPSD Low digits",     128,-71.,71.,128,-71.,71.) ;
 141   TH2F * ppsdLowCl       = new TH2F("ppsdLowCl","PPSD Low RecPoints",128,-71.,71.,128,-71.,71.) ;
 142   TH2F * nbar            = new TH2F("nbar","Primary nbar",    64,-71.,71.,64,-71.,71.);
 143   TH2F * phot            = new TH2F("phot","Primary Photon",  64,-71.,71.,64,-71.,71.);
 144   TH2F * charg           = new TH2F("charg","Primary charged",64,-71.,71.,64,-71.,71.);
 145   TH2F * recPhot         = new TH2F("recPhot","RecParticles with primary Photon",64,-71.,71.,64,-71.,71.);
 146   TH2F * recNbar         = new TH2F("recNbar","RecParticles with primary Nbar",  64,-71.,71.,64,-71.,71.);
 147
 148   //========== Create the Clusterizer
 149   fClu = new AliPHOSClusterizerv1() ;
 150
 151   gAlice->GetEvent(Nevent);
 152
 153   TParticle * primary ;
 154   Int_t iPrimary ;
 155   for ( iPrimary = 0 ; iPrimary < gAlice->GetNtrack() ; iPrimary++)
 156     {
 157       primary = gAlice->Particle(iPrimary) ;
 158       Int_t primaryType = primary->GetPdgCode() ;
 159       if( (primaryType == 211)||(primaryType == -211)||(primaryType == 2212)||(primaryType == -2212) ) {
 160         Int_t moduleNumber ;
 161         Double_t primX, primZ ;
 162         fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
 163         if(moduleNumber==Nmod)
 164           charg->Fill(primZ,primX,primary->Energy()) ;
 165       }
 166       if( primaryType == 22 ) {
 167         Int_t moduleNumber ;
 168         Double_t primX, primZ ;
 169         fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
 170         if(moduleNumber==Nmod)
 171           phot->Fill(primZ,primX,primary->Energy()) ;
 172       }
 173       else{
 174         if( primaryType == -2112 ) {
 175           Int_t moduleNumber ;
 176           Double_t primX, primZ ;
 177           fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
 178           if(moduleNumber==Nmod)
 179             nbar->Fill(primZ,primX,primary->Energy()) ;
 180         }
 181       }
 182     }
 183
 184   fPHOS->SetTreeAddress() ;
 185
 186   gAlice->TreeD()->GetEvent(0) ;
 187   gAlice->TreeR()->GetEvent(0) ;
 188
 189   TObjArray ** emcRecPoints =  fPHOS->EmcRecPoints() ;
 190   TObjArray ** ppsdRecPoints = fPHOS->PpsdRecPoints() ;
 191   TClonesArray ** recParticleList  = fPHOS->RecParticles() ;
 192
 193   Int_t iDigit ;
 194   AliPHOSDigit * digit ;
 195
 196   for(iDigit = 0; iDigit < fPHOS->Digits()->GetEntries(); iDigit++)
 197     {
 198       digit = (AliPHOSDigit *) fPHOS->Digits()->At(iDigit) ;
 199       Int_t relid[4];
 200       fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
 201       Float_t x,z ;
 202       fGeom->RelPosInModule(relid,x,z) ;
 203       Float_t e = fClu->Calibrate(digit->GetAmp()) ;
 204       if(relid[0]==Nmod){
 205         if(relid[1]==0)  //EMC
 206           digitOccupancy->Fill(x,z,e) ;
 207         if((relid[1]>0)&&(relid[1]<17))
 208           ppsdUp->Fill(x,z,e) ;
 209         if(relid[1]>16)
 210           ppsdLow->Fill(x,z,e) ;
 211       }
 212     }
 213
 214   Int_t irecp ;
 215   TVector3 pos ;
 216
 217   for(irecp = 0; irecp < (*emcRecPoints)->GetEntries() ; irecp ++){
 218     AliPHOSEmcRecPoint * emc= (AliPHOSEmcRecPoint*)(*emcRecPoints)->At(irecp) ;
 219     if(emc->GetPHOSMod()==Nmod){
 220       emc->GetLocalPosition(pos) ;
 221       emcOccupancy->Fill(pos.X(),pos.Z(),emc->GetEnergy());
 222     }
 223   }
 224
 225   for(irecp = 0; irecp < (*ppsdRecPoints)->GetEntries() ; irecp ++){
 226     AliPHOSPpsdRecPoint * ppsd= (AliPHOSPpsdRecPoint *)(*ppsdRecPoints)->At(irecp) ;
 227     if(ppsd->GetPHOSMod()==Nmod){
 228       ppsd->GetLocalPosition(pos) ;
 229       if(ppsd->GetUp())
 230         ppsdUpCl->Fill(pos.X(),pos.Z(),ppsd->GetEnergy());
 231       else
 232         ppsdLowCl->Fill(pos.X(),pos.Z(),ppsd->GetEnergy());
 233     }
 234   }
 235
 236   AliPHOSRecParticle * recParticle ;
 237   Int_t iRecParticle ;
 238   for(iRecParticle = 0; iRecParticle < (*recParticleList)->GetEntries() ;iRecParticle++ )
 239     {
 240       recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(iRecParticle) ;
 241
 242       Int_t moduleNumberRec ;
 243       Double_t recX, recZ ;
 244       fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
 245       if(moduleNumberRec == Nmod){
 246
 247         Double_t minDistance = 5. ;
 248         Int_t closestPrimary = -1 ;
 249
 250         Int_t numberofprimaries ;
 251         Int_t * listofprimaries  = recParticle->GetPrimaries(numberofprimaries)  ;
 252         Int_t index ;
 253         TParticle * primary ;
 254         Double_t distance = minDistance ;
 255
 256         for ( index = 0 ; index < numberofprimaries ; index++){
 257           primary = gAlice->Particle(listofprimaries[index]) ;
 258           Int_t moduleNumber ;
 259           Double_t primX, primZ ;
 260           fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
 261           if(moduleNumberRec == moduleNumber)
 262             distance = TMath::Sqrt((recX-primX)*(recX-primX)+(recZ-primZ)*(recZ-primZ) ) ;
 263           if(minDistance > distance)
 264             {
 265               minDistance = distance ;
 266               closestPrimary = listofprimaries[index] ;
 267             }
 268         }
 269
 270         if(closestPrimary >=0 ){
 271
 272           Int_t primaryType = gAlice->Particle(closestPrimary)->GetPdgCode() ;
 273
 274           if(primaryType==22)
 275             recPhot->Fill(recZ,recX,recParticle->Energy()) ;
 276           else
 277             if(primaryType==-2112)
 278               recNbar->Fill(recZ,recX,recParticle->Energy()) ;
 279         }
 280       }
 281     }
 282
 283
 284   digitOccupancy->Draw("box") ;
 285   emcOccupancy->SetLineColor(2) ;
 286   emcOccupancy->Draw("boxsame") ;
 287   ppsdUp->SetLineColor(3) ;
 288   ppsdUp->Draw("boxsame") ;
 289   ppsdLow->SetLineColor(4) ;
 290   ppsdLow->Draw("boxsame") ;
 291   phot->SetLineColor(8) ;
 292   phot->Draw("boxsame") ;
 293   nbar->SetLineColor(6) ;
 294   nbar->Draw("boxsame") ;
 295
 296 }
 297 //____________________________________________________________________________
 298  void AliPHOSAnalyze::Reconstruct(Int_t nevents,Int_t firstEvent )
 299 {
 300
 301   // Performs reconstruction of EMC and CPV (GPS2, IHEP or MIXT)
 302   // for events from FirstEvent to Nevents
 303
 304   Int_t ievent ;
 305   for ( ievent=firstEvent; ievent<nevents; ievent++) {
 306     if (ievent==firstEvent) {
 307       cout << "Analyze > Starting Reconstructing " << endl ;
 308       //========== Create the Clusterizer
 309       fClu = new AliPHOSClusterizerv1() ;
 310
 311       //========== Creates the track segment maker
 312       fTrs = new AliPHOSTrackSegmentMakerv1()  ;
 313       //          fTrs->UnsetUnfoldFlag() ;
 314
 315       //========== Creates the particle identifier
 316       fPID = new AliPHOSPIDv1() ;
 317       fPID->SetShowerProfileCuts(0.3, 1.8, 0.3, 1.8 ) ;
 318
 319       //========== Creates the Reconstructioner
 320       fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;
 321       if (fDebugLevel != 0) fRec -> SetDebugReconstruction(kTRUE);
 322     }
 323
 324     if (fDebugLevel != 0 ||
 325         (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
 326       cout <<  "======= Analyze ======> Event " << ievent+1 << endl ;
 327
 328     fPHOS->Enable() ;
 329
 330     gAlice->Hits2Digits() ;
 331
 332     //=========== Do the reconstruction
 333     fPHOS->Reconstruction(fRec);
 334
 335   }
 336
 337   if(fClu)      {delete fClu      ; fClu     =0 ;}
 338   if(fPID)      {delete fPID      ; fPID     =0 ;}
 339   if(fRec)      {delete fRec      ; fRec     =0 ;}
 340   if(fTrs)      {delete fTrs      ; fTrs     =0 ;}
 341
 342 }
 343
 344 //-------------------------------------------------------------------------------------
 345 void AliPHOSAnalyze::ReadAndPrintCPV(Int_t EvFirst, Int_t EvLast)
 346 {
 347 //   //
 348 //   // Read and print generated and reconstructed hits in CPV
 349 //   // for events from EvFirst to Nevent.
 350 //   // If only EvFirst is defined, print only this one event.
 351 //   // Author: Yuri Kharlov
 352 //   // 12 October 2000
 353 //   //
 354
 355 //   if (EvFirst!=0 && EvLast==0) EvLast=EvFirst;
 356 //   for ( Int_t ievent=EvFirst; ievent<=EvLast; ievent++) {
 357
 358 //     //========== Event Number>
 359 //     cout << endl <<  "==== ReadAndPrintCPV ====> Event is " << ievent+1 << endl ;
 360
 361 //     //=========== Connects the various Tree's for evt
 362 //     Int_t ntracks = gAlice->GetEvent(ievent);
 363
 364 //     //========== Creating branches ===================================
 365 //     AliPHOSRecPoint::RecPointsList ** emcRecPoints = fPHOS->EmcRecPoints() ;
 366 //     gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP" , emcRecPoints  ) ;
 367
 368 //     AliPHOSRecPoint::RecPointsList ** cpvRecPoints = fPHOS->PpsdRecPoints() ;
 369 //     gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", cpvRecPoints ) ;
 370
 371 //     // Read and print CPV hits
 372
 373 //     AliPHOSCPVModule cpvModule;
 374 //     TClonesArray    *cpvHits;
 375 //     Int_t           nCPVhits;
 376 //     AliPHOSCPVHit   *cpvHit;
 377 //     TLorentzVector   p;
 378 //     Float_t          xgen, zgen;
 379 //     Int_t            ipart;
 380 //     Int_t            nGenHits = 0;
 381 //     for (Int_t itrack=0; itrack<ntracks; itrack++) {
 382 //       //=========== Get the Hits Tree for the Primary track itrack
 383 //       gAlice->ResetHits();
 384 //       gAlice->TreeH()->GetEvent(itrack);
 385 //       Int_t iModule = 0 ;
 386 //       for (iModule=0; iModule < fGeom->GetNCPVModules(); iModule++) {
 387 //      cpvModule = fPHOS->GetCPVModule(iModule);
 388 //      cpvHits   = cpvModule.Hits();
 389 //      nCPVhits  = cpvHits->GetEntriesFast();
 390 //      for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
 391 //        nGenHits++;
 392 //        cpvHit = (AliPHOSCPVHit*)cpvHits->UncheckedAt(ihit);
 393 //        p      = cpvHit->GetMomentum();
 394 //        xgen   = cpvHit->X();
 395 //        zgen   = cpvHit->Y();
 396 //        ipart  = cpvHit->GetIpart();
 397 //        printf("CPV hit in module %d: ",iModule+1);
 398 //        printf(" p = (%f, %f, %f, %f) GeV,\n",
 399 //               p.Px(),p.Py(),p.Pz(),p.Energy());
 400 //        printf("                  (X,Z) = (%8.4f, %8.4f) cm, ipart = %d\n",
 401 //               xgen,zgen,ipart);
 402 //      }
 403 //       }
 404 //     }
 405
 406 //     // Read and print CPV reconstructed points
 407
 408 //     //=========== Gets the Reconstruction TTree
 409 //     gAlice->TreeR()->GetEvent(0) ;
 410 //     printf("Recpoints: %d\n",(*fPHOS->CpvRecPoints())->GetEntries());
 411 //     TIter nextRP(*fPHOS->CpvRecPoints() ) ;
 412 //     AliPHOSCpvRecPoint *cpvRecPoint ;
 413 //     Int_t nRecPoints = 0;
 414 //     while( ( cpvRecPoint = (AliPHOSCpvRecPoint *)nextRP() ) ) {
 415 //       nRecPoints++;
 416 //       TVector3  locpos;
 417 //       cpvRecPoint->GetLocalPosition(locpos);
 418 //       Int_t phosModule = cpvRecPoint->GetPHOSMod();
 419 //       printf("CPV recpoint in module %d: (X,Z) = (%f,%f) cm\n",
 420 //           phosModule,locpos.X(),locpos.Z());
 421 //     }
 422 //     printf("This event has %d generated hits and %d reconstructed points\n",
 423 //         nGenHits,nRecPoints);
 424 //   }
 425 }
 426
 427 //____________________________________________________________________________
 428 void AliPHOSAnalyze::AnalyzeCPV(Int_t Nevents)
 429 {
 430 //   //
 431 //   // Analyzes CPV characteristics
 432 //   // Author: Yuri Kharlov
 433 //   // 9 October 2000
 434 //   //
 435
 436 //   // Book histograms
 437
 438 //   TH1F *hDx   = new TH1F("hDx"  ,"CPV x-resolution@reconstruction",100,-5. , 5.);
 439 //   TH1F *hDz   = new TH1F("hDz"  ,"CPV z-resolution@reconstruction",100,-5. , 5.);
 440 //   TH1F *hDr   = new TH1F("hDr"  ,"CPV r-resolution@reconstruction",100, 0. , 5.);
 441 //   TH1S *hNrp  = new TH1S("hNrp" ,"CPV rec.point multiplicity",      21,-0.5,20.5);
 442 //   TH1S *hNrpX = new TH1S("hNrpX","CPV rec.point Phi-length"  ,      21,-0.5,20.5);
 443 //   TH1S *hNrpZ = new TH1S("hNrpZ","CPV rec.point Z-length"    ,      21,-0.5,20.5);
 444
 445 //   cout << "Start CPV Analysis"<< endl ;
 446 //   for ( Int_t ievent=0; ievent<Nevents; ievent++) {
 447
 448 //     //========== Event Number>
 449 // //      if ( (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
 450 //       cout << endl <<  "==== AnalyzeCPV ====> Event is " << ievent+1 << endl ;
 451
 452 //     //=========== Connects the various Tree's for evt
 453 //     Int_t ntracks = gAlice->GetEvent(ievent);
 454
 455 //     //========== Creating branches ===================================
 456 //     AliPHOSRecPoint::RecPointsList ** emcRecPoints = fPHOS->EmcRecPoints() ;
 457 //     gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP" , emcRecPoints  ) ;
 458
 459 //     AliPHOSRecPoint::RecPointsList ** cpvRecPoints = fPHOS->PpsdRecPoints() ;
 460 //     gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", cpvRecPoints ) ;
 461
 462 //     // Create and fill arrays of hits for each CPV module
 463
 464 //     Int_t nOfModules = fGeom->GetNModules();
 465 //     TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
 466 //     Int_t iModule = 0;
 467 //     for (iModule=0; iModule < nOfModules; iModule++)
 468 //       hitsPerModule[iModule] = new TClonesArray("AliPHOSCPVHit",100);
 469
 470 //     AliPHOSCPVModule cpvModule;
 471 //     TClonesArray    *cpvHits;
 472 //     Int_t           nCPVhits;
 473 //     AliPHOSCPVHit   *cpvHit;
 474 //     TLorentzVector   p;
 475 //     Float_t          xzgen[2];
 476 //     Int_t            ipart;
 477
 478 //     // First go through all primary tracks and fill the arrays
 479 //     // of hits per each CPV module
 480
 481 //     for (Int_t itrack=0; itrack<ntracks; itrack++) {
 482 //       // Get the Hits Tree for the Primary track itrack
 483 //       gAlice->ResetHits();
 484 //       gAlice->TreeH()->GetEvent(itrack);
 485 //       for (Int_t iModule=0; iModule < nOfModules; iModule++) {
 486 //      cpvModule = fPHOS->GetCPVModule(iModule);
 487 //      cpvHits   = cpvModule.Hits();
 488 //      nCPVhits  = cpvHits->GetEntriesFast();
 489 //      for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
 490 //        cpvHit   = (AliPHOSCPVHit*)cpvHits->UncheckedAt(ihit);
 491 //        p        = cpvHit->GetMomentum();
 492 //        xzgen[0] = cpvHit->X();
 493 //        xzgen[1] = cpvHit->Y();
 494 //        ipart    = cpvHit->GetIpart();
 495 //        TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
 496 //        new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSCPVHit(*cpvHit);
 497 //      }
 498 //      cpvModule.Clear();
 499 //       }
 500 //     }
 501 //     for (iModule=0; iModule < nOfModules; iModule++) {
 502 //       Int_t nsum = hitsPerModule[iModule]->GetEntriesFast();
 503 //       printf("Module %d has %d hits\n",iModule,nsum);
 504 //     }
 505
 506 //     // Then go through reconstructed points and for each find
 507 //     // the closeset hit
 508 //     // The distance from the rec.point to the closest hit
 509 //     // gives the coordinate resolution of the CPV
 510
 511 //     // Get the Reconstruction Tree
 512 //     gAlice->TreeR()->GetEvent(0) ;
 513 //     TIter nextRP(*fPHOS->PpsdRecPoints() ) ;
 514 //     AliPHOSCpvRecPoint *cpvRecPoint ;
 515 //     Float_t xgen, zgen;
 516 //     while( ( cpvRecPoint = (AliPHOSCpvRecPoint *)nextRP() ) ) {
 517 //       TVector3  locpos;
 518 //       cpvRecPoint->GetLocalPosition(locpos);
 519 //       Int_t phosModule = cpvRecPoint->GetPHOSMod();
 520 //       Int_t rpMult     = cpvRecPoint->GetDigitsMultiplicity();
 521 //       Int_t rpMultX, rpMultZ;
 522 //       cpvRecPoint->GetClusterLengths(rpMultX,rpMultZ);
 523 //       Float_t xrec  = locpos.X();
 524 //       Float_t zrec  = locpos.Z();
 525 //       Float_t dxmin = 1.e+10;
 526 //       Float_t dzmin = 1.e+10;
 527 //       Float_t r2min = 1.e+10;
 528 //       Float_t r2;
 529
 530 //       cpvHits = hitsPerModule[phosModule-1];
 531 //       Int_t nCPVhits  = cpvHits->GetEntriesFast();
 532 //       for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
 533 //      cpvHit = (AliPHOSCPVHit*)cpvHits->UncheckedAt(ihit);
 534 //      xgen   = cpvHit->X();
 535 //      zgen   = cpvHit->Y();
 536 //      r2 = TMath::Power((xgen-xrec),2) + TMath::Power((zgen-zrec),2);
 537 //      if ( r2 < r2min ) {
 538 //        r2min = r2;
 539 //        dxmin = xgen - xrec;
 540 //        dzmin = zgen - zrec;
 541 //      }
 542 //       }
 543 //       hDx  ->Fill(dxmin);
 544 //       hDz  ->Fill(dzmin);
 545 //       hDr  ->Fill(TMath::Sqrt(r2min));
 546 //       hNrp ->Fill(rpMult);
 547 //       hNrpX->Fill(rpMultX);
 548 //       hNrpZ->Fill(rpMultZ);
 549 //     }
 550 //     delete [] hitsPerModule;
 551 //   }
 552 //   // Save histograms
 553
 554 //   Text_t outputname[80] ;
 555 //   sprintf(outputname,"%s.analyzed",fRootFile->GetName());
 556 //   TFile output(outputname,"RECREATE");
 557 //   output.cd();
 558
 559 //   hDx  ->Write() ;
 560 //   hDz  ->Write() ;
 561 //   hDr  ->Write() ;
 562 //   hNrp ->Write() ;
 563 //   hNrpX->Write() ;
 564 //   hNrpZ->Write() ;
 565
 566 //   // Plot histograms
 567
 568 //   TCanvas *cpvCanvas = new TCanvas("CPV","CPV analysis",20,20,800,400);
 569 //   gStyle->SetOptStat(111111);
 570 //   gStyle->SetOptFit(1);
 571 //   gStyle->SetOptDate(1);
 572 //   cpvCanvas->Divide(3,2);
 573
 574 //   cpvCanvas->cd(1);
 575 //   gPad->SetFillColor(10);
 576 //   hNrp->SetFillColor(16);
 577 //   hNrp->Draw();
 578
 579 //   cpvCanvas->cd(2);
 580 //   gPad->SetFillColor(10);
 581 //   hNrpX->SetFillColor(16);
 582 //   hNrpX->Draw();
 583
 584 //   cpvCanvas->cd(3);
 585 //   gPad->SetFillColor(10);
 586 //   hNrpZ->SetFillColor(16);
 587 //   hNrpZ->Draw();
 588
 589 //   cpvCanvas->cd(4);
 590 //   gPad->SetFillColor(10);
 591 //   hDx->SetFillColor(16);
 592 //   hDx->Fit("gaus");
 593 //   hDx->Draw();
 594
 595 //   cpvCanvas->cd(5);
 596 //   gPad->SetFillColor(10);
 597 //   hDz->SetFillColor(16);
 598 //   hDz->Fit("gaus");
 599 //   hDz->Draw();
 600
 601 //   cpvCanvas->cd(6);
 602 //   gPad->SetFillColor(10);
 603 //   hDr->SetFillColor(16);
 604 //   hDr->Draw();
 605
 606 //   cpvCanvas->Print("CPV.ps");
 607
 608 }
 609
 610 //____________________________________________________________________________
 611  void AliPHOSAnalyze::InvariantMass(Int_t Nevents )
 612 {
 613   // Calculates Real and Mixed invariant mass distributions
 614
 615   const Int_t knMixedEvents = 4 ; //# of events used for calculation of 'mixed' distribution
 616   Int_t mixedLoops = (Int_t )TMath::Ceil(Nevents/knMixedEvents) ;
 617
 618   //========== Booking Histograms
 619   TH2D * hRealEM   = new TH2D("hRealEM",   "Real for EM particles",      250,0.,1.,40,0.,4.) ;
 620   TH2D * hRealPhot = new TH2D("hRealPhot", "Real for kPhoton particles", 250,0.,1.,40,0.,4.) ;
 621   TH2D * hMixedEM  = new TH2D("hMixedEM",  "Mixed for EM particles",     250,0.,1.,40,0.,4.) ;
 622   TH2D * hMixedPhot= new TH2D("hMixedPhot","Mixed for kPhoton particles",250,0.,1.,40,0.,4.) ;
 623
 624   Int_t ievent;
 625   Int_t eventInMixedLoop ;
 626
 627   Int_t nRecParticles[4];//knMixedEvents] ;
 628
 629   AliPHOSRecParticle::RecParticlesList * allRecParticleList  = new TClonesArray("AliPHOSRecParticle", knMixedEvents*1000) ;
 630
 631   for(eventInMixedLoop = 0; eventInMixedLoop < mixedLoops; eventInMixedLoop++  ){
 632     Int_t iRecPhot = 0 ;
 633
 634     for ( ievent=0; ievent < knMixedEvents; ievent++){
 635
 636       Int_t absEventNumber = eventInMixedLoop*knMixedEvents + ievent ;
 637
 638       //=========== Connects the various Tree's for evt
 639       gAlice->GetEvent(absEventNumber);
 640
 641       //========== Creating branches ===================================
 642       fPHOS->SetTreeAddress() ;
 643
 644       gAlice->TreeD()->GetEvent(0) ;
 645       gAlice->TreeR()->GetEvent(0) ;
 646
 647       TClonesArray ** recParticleList  = fPHOS->RecParticles() ;
 648
 649
 650       AliPHOSRecParticle * recParticle ;
 651       Int_t iRecParticle ;
 652       for(iRecParticle = 0; iRecParticle < (*recParticleList)->GetEntries() ;iRecParticle++ )
 653         {
 654           recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(iRecParticle) ;
 655           if((recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)||
 656              (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM)){
 657             new( (*allRecParticleList)[iRecPhot] ) AliPHOSRecParticle(*recParticle) ;
 658             iRecPhot++;
 659           }
 660         }
 661
 662         nRecParticles[ievent] = iRecPhot-1 ;
 663     }
 664
 665     //Now calculate invariant mass:
 666     Int_t irp1,irp2 ;
 667     Int_t nCurEvent = 0 ;
 668
 669     for(irp1 = 0; irp1 < allRecParticleList->GetEntries()-1; irp1++){
 670       AliPHOSRecParticle * rp1 = (AliPHOSRecParticle *)allRecParticleList->At(irp1) ;
 671
 672       for(irp2 = irp1+1; irp2 < allRecParticleList->GetEntries(); irp2++){
 673         AliPHOSRecParticle * rp2 = (AliPHOSRecParticle *)allRecParticleList->At(irp2) ;
 674
 675         Double_t invMass ;
 676         invMass = (rp1->Energy()+rp2->Energy())*(rp1->Energy()+rp2->Energy())-
 677           (rp1->Px()+rp2->Px())*(rp1->Px()+rp2->Px())-
 678           (rp1->Py()+rp2->Py())*(rp1->Py()+rp2->Py())-
 679           (rp1->Pz()+rp2->Pz())*(rp1->Pz()+rp2->Pz()) ;
 680
 681         if(invMass> 0)
 682           invMass = TMath::Sqrt(invMass);
 683
 684         Double_t pt ;
 685         pt = TMath::Sqrt((rp1->Px()+rp2->Px() )*( rp1->Px()+rp2->Px() ) +(rp1->Py()+rp2->Py())*(rp1->Py()+rp2->Py()));
 686
 687         if(irp1 > nRecParticles[nCurEvent])
 688           nCurEvent++;
 689
 690         if(irp2 <= nRecParticles[nCurEvent]){ //'Real' event
 691           hRealEM->Fill(invMass,pt);
 692           if((rp1->GetType() == AliPHOSFastRecParticle::kGAMMA)&&(rp2->GetType() == AliPHOSFastRecParticle::kGAMMA))
 693             hRealPhot->Fill(invMass,pt);
 694         }
 695         else{
 696           hMixedEM->Fill(invMass,pt);
 697           if((rp1->GetType() == AliPHOSFastRecParticle::kGAMMA)&&(rp2->GetType() == AliPHOSFastRecParticle::kGAMMA))
 698             hMixedPhot->Fill(invMass,pt);
 699         } //real-mixed
 700
 701       } //loop over second rp
 702     }//loop over first rp
 703     allRecParticleList->Delete() ;
 704   } //Loop over events
 705
 706   delete allRecParticleList ;
 707
 708   //writing output
 709   TFile output("invmass.root","RECREATE");
 710   output.cd();
 711
 712   hRealEM->Write() ;
 713   hRealPhot->Write() ;
 714   hMixedEM->Write() ;
 715   hMixedPhot->Write() ;
 716
 717   output.Write();
 718   output.Close();
 719
 720 }
 721
 722 //____________________________________________________________________________
 723  void AliPHOSAnalyze::ReadAndPrintEMC(Int_t EvFirst, Int_t EvLast)
 724 {
 725 //   //
 726 //   // Read and print generated and reconstructed hits in EMC
 727 //   // for events from EvFirst to Nevent.
 728 //   // If only EvFirst is defined, print only this one event.
 729 //   // Author: Yuri Kharlov
 730 //   // 24 November 2000
 731 //   //
 732
 733 //   if (EvFirst!=0 && EvLast==0) EvLast=EvFirst;
 734 //   Int_t ievent;
 735 //   for (ievent=EvFirst; ievent<=EvLast; ievent++) {
 736
 737 //     //========== Event Number>
 738 //     cout << endl <<  "==== ReadAndPrintEMC ====> Event is " << ievent+1 << endl ;
 739
 740 //     //=========== Connects the various Tree's for evt
 741 //     Int_t ntracks = gAlice->GetEvent(ievent);
 742 //     fPHOS->SetTreeAddress() ;
 743
 744 //     gAlice->TreeD()->GetEvent(0) ;
 745 //     gAlice->TreeR()->GetEvent(0) ;
 746
 747 //     // Loop over reconstructed particles
 748
 749 //     TClonesArray ** recParticleList  = fPHOS->RecParticles() ;
 750 //     AliPHOSRecParticle * recParticle ;
 751 //     Int_t iRecParticle ;
 752 //     Int_t *primList;
 753 //     Int_t nPrimary;
 754 //     for(iRecParticle = 0; iRecParticle < (*recParticleList)->GetEntries() ;iRecParticle++ ) {
 755 //       recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(iRecParticle) ;
 756 //       Float_t recE = recParticle->Energy();
 757 //       primList     = recParticle->GetPrimaries(nPrimary);
 758 //       Int_t moduleNumberRec ;
 759 //       Double_t recX, recZ ;
 760 //       fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
 761 //       printf("Rec point: module %d, (X,Z) = (%8.4f,%8.4f) cm, E = %.3f GeV, primary = %d\n",
 762 //           moduleNumberRec,recX,recZ,recE,*primList);
 763 //     }
 764
 765 //     // Read and print EMC hits from EMCn branches
 766
 767 //     AliPHOSCPVModule emcModule;
 768 //     TClonesArray    *emcHits;
 769 //     Int_t           nEMChits;
 770 //     AliPHOSCPVHit   *emcHit;
 771 //     TLorentzVector   p;
 772 //     Float_t          xgen, zgen;
 773 //     Int_t            ipart, primary;
 774 //     Int_t            nGenHits = 0;
 775 //     for (Int_t itrack=0; itrack<ntracks; itrack++) {
 776 //       //=========== Get the Hits Tree for the Primary track itrack
 777 //       gAlice->ResetHits();
 778 //       gAlice->TreeH()->GetEvent(itrack);
 779 //       Int_t iModule = 0 ;
 780 //       for (iModule=0; iModule < fGeom->GetNModules(); iModule++) {
 781 //      emcModule = fPHOS->GetEMCModule(iModule);
 782 //      emcHits   = emcModule.Hits();
 783 //      nEMChits  = emcHits->GetEntriesFast();
 784 //      for (Int_t ihit=0; ihit<nEMChits; ihit++) {
 785 //        nGenHits++;
 786 //        emcHit = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
 787 //        p      = emcHit->GetMomentum();
 788 //        xgen   = emcHit->X();
 789 //        zgen   = emcHit->Y();
 790 //        ipart  = emcHit->GetIpart();
 791 //        primary= emcHit->GetTrack();
 792 //        printf("EMC hit A: module %d, ",iModule+1);
 793 //        printf("    p = (%f .4, %f .4, %f .4, %f .4) GeV,\n",
 794 //               p.Px(),p.Py(),p.Pz(),p.Energy());
 795 //        printf("                     (X,Z) = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n",
 796 //               xgen,zgen,ipart,primary);
 797 //      }
 798 //       }
 799 //     }
 800
 801 // //      // Read and print EMC hits from PHOS branch
 802
 803 // //      for (Int_t itrack=0; itrack<ntracks; itrack++) {
 804 // //        //=========== Get the Hits Tree for the Primary track itrack
 805 // //        gAlice->ResetHits();
 806 // //        gAlice->TreeH()->GetEvent(itrack);
 807 // //        TClonesArray *hits = fPHOS->Hits();
 808 // //        AliPHOSHit   *hit ;
 809 // //        Int_t ihit;
 810 // //        for ( ihit = 0 ; ihit < hits->GetEntries() ; ihit++ ) {
 811 // //   hit = (AliPHOSHit*)hits->At(ihit) ;
 812 // //   Float_t hitXYZ[3];
 813 // //   hitXYZ[0]   = hit->X();
 814 // //   hitXYZ[1]   = hit->Y();
 815 // //   hitXYZ[2]   = hit->Z();
 816 // //   ipart       = hit->GetPid();
 817 // //   primary     = hit->GetPrimary();
 818 // //   Int_t absId = hit->GetId();
 819 // //   Int_t relId[4];
 820 // //   fGeom->AbsToRelNumbering(absId, relId) ;
 821 // //   Int_t module = relId[0];
 822 // //   if (relId[1]==0 && !(hitXYZ[0]==0 && hitXYZ[2]==0))
 823 // //     printf("EMC hit B: module %d, (X,Z) = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n",
 824 // //            module,hitXYZ[0],hitXYZ[2],ipart,primary);
 825 // //        }
 826 // //      }
 827
 828 //   }
 829 }
 830
 831 //____________________________________________________________________________
 832  void AliPHOSAnalyze::AnalyzeEMC(Int_t Nevents)
 833 {
 834 //   //
 835 //   // Read generated and reconstructed hits in EMC for Nevents events.
 836 //   // Plots the coordinate and energy resolution histograms.
 837 //   // Coordinate resolution is a difference between the reconstructed
 838 //   // coordinate and the exact coordinate on the face of the PHOS
 839 //   // Author: Yuri Kharlov
 840 //   // 27 November 2000
 841 //   //
 842
 843 //   // Book histograms
 844
 845 //   TH1F *hDx1   = new TH1F("hDx1"  ,"EMC x-resolution", 100,-5. , 5.);
 846 //   TH1F *hDz1   = new TH1F("hDz1"  ,"EMC z-resolution", 100,-5. , 5.);
 847 //   TH1F *hDE1   = new TH1F("hDE1"  ,"EMC E-resolution", 100,-2. , 2.);
 848
 849 //   TH2F *hDx2   = new TH2F("hDx2"  ,"EMC x-resolution", 100, 0., 10., 100,-5. , 5.);
 850 //   TH2F *hDz2   = new TH2F("hDz2"  ,"EMC z-resolution", 100, 0., 10., 100,-5. , 5.);
 851 //   TH2F *hDE2   = new TH2F("hDE2"  ,"EMC E-resolution", 100, 0., 10., 100, 0. , 5.);
 852
 853 //   cout << "Start EMC Analysis"<< endl ;
 854 //   for (Int_t ievent=0; ievent<Nevents; ievent++) {
 855
 856 //     //========== Event Number>
 857 //     if ( (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
 858 //       cout << "==== AnalyzeEMC ====> Event is " << ievent+1 << endl ;
 859
 860 //     //=========== Connects the various Tree's for evt
 861 //     Int_t ntracks = gAlice->GetEvent(ievent);
 862
 863 //     fPHOS->SetTreeAddress() ;
 864
 865 //     gAlice->TreeD()->GetEvent(0) ;
 866 //     gAlice->TreeR()->GetEvent(0) ;
 867
 868 //     // Create and fill arrays of hits for each EMC module
 869
 870 //     Int_t nOfModules = fGeom->GetNModules();
 871 //     TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
 872 //     Int_t iModule;
 873 //     for (iModule=0; iModule < nOfModules; iModule++)
 874 //       hitsPerModule[iModule] = new TClonesArray("AliPHOSCPVHit",100);
 875
 876 //     AliPHOSCPVModule emcModule;
 877 //     TClonesArray    *emcHits;
 878 //     Int_t           nEMChits;
 879 //     AliPHOSCPVHit   *emcHit;
 880
 881 //     // First go through all primary tracks and fill the arrays
 882 //     // of hits per each EMC module
 883
 884 //     for (Int_t itrack=0; itrack<ntracks; itrack++) {
 885 //       // Get the Hits Tree for the Primary track itrack
 886 //       gAlice->ResetHits();
 887 //       gAlice->TreeH()->GetEvent(itrack);
 888 //       for (Int_t iModule=0; iModule < nOfModules; iModule++) {
 889 //      emcModule = fPHOS->GetEMCModule(iModule);
 890 //      emcHits   = emcModule.Hits();
 891 //      nEMChits  = emcHits->GetEntriesFast();
 892 //      for (Int_t ihit=0; ihit<nEMChits; ihit++) {
 893 //        emcHit   = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
 894 //        TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
 895 //        new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSCPVHit(*emcHit);
 896 //      }
 897 //      emcModule.Clear();
 898 //       }
 899 //     }
 900
 901 //     // Loop over reconstructed particles
 902
 903 //     TClonesArray ** recParticleList  = fPHOS->RecParticles() ;
 904 //     AliPHOSRecParticle * recParticle ;
 905 //     Int_t nEMCrecs = (*recParticleList)->GetEntries();
 906 //     if (nEMCrecs == 1) {
 907 //       recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(0) ;
 908 //       Float_t recE = recParticle->Energy();
 909 //       Int_t phosModule;
 910 //       Double_t recX, recZ ;
 911 //       fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), phosModule, recX, recZ) ;
 912
 913 //       // for this rec.point take the hit list in the same PHOS module
 914
 915 //       emcHits = hitsPerModule[phosModule-1];
 916 //       Int_t nEMChits  = emcHits->GetEntriesFast();
 917 //       if (nEMChits == 1) {
 918 //      Float_t genX, genZ, genE;
 919 //      for (Int_t ihit=0; ihit<nEMChits; ihit++) {
 920 //        emcHit = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
 921 //        genX   = emcHit->X();
 922 //        genZ   = emcHit->Y();
 923 //        genE   = emcHit->GetMomentum().E();
 924 //      }
 925 //      Float_t dx = recX - genX;
 926 //      Float_t dz = recZ - genZ;
 927 //      Float_t de = recE - genE;
 928 //      hDx1  ->Fill(dx);
 929 //      hDz1  ->Fill(dz);
 930 //      hDE1  ->Fill(de);
 931 //      hDx2  ->Fill(genE,dx);
 932 //      hDz2  ->Fill(genE,dz);
 933 //      hDE2  ->Fill(genE,recE);
 934 //       }
 935 //     }
 936 //     delete [] hitsPerModule;
 937 //   }
 938 //   // Save histograms
 939
 940 //   Text_t outputname[80] ;
 941 //   sprintf(outputname,"%s.analyzed",fRootFile->GetName());
 942 //   TFile output(outputname,"RECREATE");
 943 //   output.cd();
 944
 945 //   hDx1  ->Write() ;
 946 //   hDz1  ->Write() ;
 947 //   hDE1  ->Write() ;
 948 //   hDx2  ->Write() ;
 949 //   hDz2  ->Write() ;
 950 //   hDE2  ->Write() ;
 951
 952 //   // Plot histograms
 953
 954 //   TCanvas *emcCanvas = new TCanvas("EMC","EMC analysis",20,20,700,300);
 955 //   gStyle->SetOptStat(111111);
 956 //   gStyle->SetOptFit(1);
 957 //   gStyle->SetOptDate(1);
 958 //   emcCanvas->Divide(3,1);
 959
 960 //   emcCanvas->cd(1);
 961 //   gPad->SetFillColor(10);
 962 //   hDx1->SetFillColor(16);
 963 //   hDx1->Draw();
 964
 965 //   emcCanvas->cd(2);
 966 //   gPad->SetFillColor(10);
 967 //   hDz1->SetFillColor(16);
 968 //   hDz1->Draw();
 969
 970 //   emcCanvas->cd(3);
 971 //   gPad->SetFillColor(10);
 972 //   hDE1->SetFillColor(16);
 973 //   hDE1->Draw();
 974
 975 //   emcCanvas->Print("EMC.ps");
 976
 977 }
 978
 979 //____________________________________________________________________________
 980  void AliPHOSAnalyze::AnalyzeResolutions(Int_t Nevents )
 981 {
 982   // analyzes Nevents events and calculate Energy and Position resolution as well as
 983   // probaility of correct indentifiing of the incident particle
 984
 985   //========== Booking Histograms
 986   cout << "AnalyzeResolutions > " << "Booking Histograms" << endl ;
 987   BookResolutionHistograms();
 988
 989   Int_t counter[9][5] ;
 990   Int_t i1,i2,totalInd = 0 ;
 991   for(i1 = 0; i1<9; i1++)
 992     for(i2 = 0; i2<5; i2++)
 993       counter[i1][i2] = 0 ;
 994
 995   Int_t totalPrimary = 0 ;
 996   Int_t totalRecPart = 0 ;
 997   Int_t totalRPwithPrim = 0 ;
 998   Int_t ievent;
 999
1000   cout << "Start Analysing"<< endl ;
1001   for ( ievent=0; ievent<Nevents; ievent++)
1002     {
1003
1004       //========== Event Number>
1005       //      if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. )
1006         cout <<  "AnalyzeResolutions > " << "Event is " << ievent << endl ;
1007
1008       //=========== Connects the various Tree's for evt
1009       gAlice->GetEvent(ievent);
1010
1011       //=========== Gets the Kine TTree
1012       gAlice->TreeK()->GetEvent(0) ;
1013
1014       //=========== Gets the list of Primari Particles
1015
1016       TParticle * primary ;
1017       Int_t iPrimary ;
1018       for ( iPrimary = 0 ; iPrimary < gAlice->GetNtrack() ; iPrimary++)
1019         {
1020           primary = gAlice->Particle(iPrimary) ;
1021           Int_t primaryType = primary->GetPdgCode() ;
1022           if( primaryType == 22 ) {
1023             Int_t moduleNumber ;
1024             Double_t primX, primZ ;
1025             fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
1026             if(moduleNumber){
1027               fhPrimary->Fill(primary->Energy()) ;
1028               if(primary->Energy() > 0.3)
1029                 totalPrimary++ ;
1030             }
1031           }
1032         }
1033
1034       fPHOS->SetTreeAddress() ;
1035
1036       gAlice->TreeD()->GetEvent(0) ;
1037       gAlice->TreeR()->GetEvent(0) ;
1038
1039       TClonesArray ** recParticleList  = fPHOS->RecParticles() ;
1040
1041       AliPHOSRecParticle * recParticle ;
1042       Int_t iRecParticle ;
1043       for(iRecParticle = 0; iRecParticle < (*recParticleList)->GetEntries() ;iRecParticle++ )
1044         {
1045           recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(iRecParticle) ;
1046           fhAllRP->Fill(CorrectEnergy(recParticle->Energy())) ;
1047
1048           Int_t moduleNumberRec ;
1049           Double_t recX, recZ ;
1050           fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
1051
1052           Double_t minDistance  = 100. ;
1053           Int_t closestPrimary = -1 ;
1054
1055           Int_t numberofprimaries ;
1056           Int_t * listofprimaries  = recParticle->GetPrimaries(numberofprimaries)  ;
1057           Int_t index ;
1058           TParticle * primary ;
1059           Double_t distance = minDistance ;
1060           Double_t dX, dZ;
1061           Double_t dXmin = 0.;
1062           Double_t dZmin = 0. ;
1063           for ( index = 0 ; index < numberofprimaries ; index++){
1064             primary = gAlice->Particle(listofprimaries[index]) ;
1065             Int_t moduleNumber ;
1066             Double_t primX, primZ ;
1067             fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
1068             if(moduleNumberRec == moduleNumber) {
1069               dX = recX - primX;
1070               dZ = recZ - primZ;
1071               distance = TMath::Sqrt(dX*dX + dZ*dZ) ;
1072               if(minDistance > distance) {
1073                 minDistance = distance ;
1074                 dXmin = dX;
1075                 dZmin = dZ;
1076                 closestPrimary = listofprimaries[index] ;
1077               }
1078             }
1079           }
1080           totalRecPart++ ;
1081
1082           if(closestPrimary >=0 ){
1083             totalRPwithPrim++;
1084
1085             Int_t primaryType = gAlice->Particle(closestPrimary)->GetPdgCode() ;
1086 //          TParticlePDG* pDGparticle = gAlice->ParticleAt(closestPrimary)->GetPDG();
1087 //          Double_t charge =  PDGparticle->Charge() ;
1088 //          if(charge)
1089 //            cout <<"Primary " <<primaryType << " E " << ((TParticle *)primaryList->At(closestPrimary))->Energy() << endl ;
1090             Int_t primaryCode ;
1091             switch(primaryType)
1092               {
1093               case 22:
1094                 primaryCode = 0;  //Photon
1095                 fhAllEnergy   ->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy()) ;
1096                 fhAllPosition ->Fill(gAlice->Particle(closestPrimary)->Energy(), minDistance) ;
1097                 fhAllPositionX->Fill(dXmin);
1098                 fhAllPositionZ->Fill(dZmin);
1099                 break;
1100               case 11 :
1101                 primaryCode = 1;  //Electron
1102                 break;
1103               case -11 :
1104                 primaryCode = 1;  //positron
1105                 break;
1106               case 321 :
1107                 primaryCode = 4;  //K+
1108                 break;
1109               case -321 :
1110                 primaryCode = 4;  //K-
1111                 break;
1112               case 310 :
1113                 primaryCode = 4;  //K0s
1114                 break;
1115               case 130 :
1116                 primaryCode = 4;  //K0l
1117                 break;
1118               case 211 :
1119                 primaryCode = 2;  //K0l
1120                 break;
1121               case -211 :
1122                 primaryCode = 2;  //K0l
1123                 break;
1124               case 2212 :
1125                 primaryCode = 2;  //K0l
1126                 break;
1127               case -2212 :
1128                 primaryCode = 2;  //K0l
1129                 break;
1130               default:
1131                 primaryCode = 3; //ELSE
1132                 break;
1133               }
1134
1135             switch(recParticle->GetType())
1136               {
1137               case AliPHOSFastRecParticle::kGAMMA:
1138                 if(primaryType == 22){
1139                   fhPhotEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy() ) ;
1140                   fhEMEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy() ) ;
1141                   fhPPSDEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy() ) ;
1142
1143                   fhPhotPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance) ;
1144                   fhEMPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance) ;
1145                   fhPPSDPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance) ;
1146
1147                   fhPhotReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1148                   fhPhotEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1149                   fhPhotPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1150
1151                   fhPhotPhot->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1152                 }
1153                 if(primaryType == 2112){ //neutron
1154                   fhNReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1155                   fhNEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1156                   fhNPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1157                 }
1158
1159                 if(primaryType == -2112){ //neutron ~
1160                   fhNBarReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1161                   fhNBarEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1162                   fhNBarPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1163
1164                 }
1165                 if(primaryCode == 2){
1166                   fhChargedReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1167                   fhChargedEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1168                   fhChargedPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1169                 }
1170
1171                 fhAllReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1172                 fhAllEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1173                 fhAllPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1174                 fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1175                 fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1176                 fhPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1177                 counter[0][primaryCode]++;
1178                 break;
1179               case  AliPHOSFastRecParticle::kELECTRON:
1180                 if(primaryType == 22){
1181                   fhPhotElec->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1182                   fhEMEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy() ) ;
1183                   fhEMPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance) ;
1184                   fhPhotEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1185                   fhPhotPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1186                 }
1187                 if(primaryType == 2112){ //neutron
1188                   fhNEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1189                   fhNPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1190                 }
1191
1192                 if(primaryType == -2112){ //neutron ~
1193                   fhNBarEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1194                   fhNBarPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1195
1196                 }
1197                 if(primaryCode == 2){
1198                   fhChargedEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1199                   fhChargedPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1200                 }
1201
1202                 fhAllEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1203                 fhAllPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1204                 fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1205                 fhPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1206                 counter[1][primaryCode]++;
1207                 break;
1208               case  AliPHOSFastRecParticle::kNEUTRALHA:
1209                 if(primaryType == 22)
1210                   fhPhotNeuH->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1211
1212                 fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1213                 counter[2][primaryCode]++;
1214                 break ;
1215               case  AliPHOSFastRecParticle::kNEUTRALEM:
1216                 if(primaryType == 22){
1217                   fhEMEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(),recParticle->Energy() ) ;
1218                   fhEMPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance ) ;
1219
1220                   fhPhotNuEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1221                   fhPhotEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1222                 }
1223                 if(primaryType == 2112) //neutron
1224                   fhNEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1225
1226                 if(primaryType == -2112) //neutron ~
1227                   fhNBarEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1228
1229                 if(primaryCode == 2)
1230                   fhChargedEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1231
1232                 fhAllEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1233                 fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1234                 fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1235
1236                 counter[3][primaryCode]++;
1237                 break ;
1238               case  AliPHOSFastRecParticle::kCHARGEDHA:
1239                 if(primaryType == 22) //photon
1240                   fhPhotChHa->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1241
1242                 counter[4][primaryCode]++ ;
1243                 break ;
1244               case  AliPHOSFastRecParticle::kGAMMAHA:
1245                   if(primaryType == 22){ //photon
1246                     fhPhotGaHa->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1247                     fhPPSDEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy() ) ;
1248                     fhPPSDPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance) ;
1249                     fhPhotPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1250                   }
1251                   if(primaryType == 2112){ //neutron
1252                     fhNPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1253                   }
1254
1255                   if(primaryType == -2112){ //neutron ~
1256                     fhNBarPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1257                   }
1258                   if(primaryCode == 2){
1259                     fhChargedPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1260                   }
1261
1262                   fhAllPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1263                   fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1264                   fhPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1265                   counter[5][primaryCode]++ ;
1266                   break ;
1267               case  AliPHOSFastRecParticle::kABSURDEM:
1268                 counter[6][primaryCode]++ ;
1269                 fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1270                 break;
1271               case  AliPHOSFastRecParticle::kABSURDHA:
1272                 counter[7][primaryCode]++ ;
1273                 break;
1274               default:
1275                 counter[8][primaryCode]++ ;
1276                 break;
1277               }
1278           }
1279         }
1280     }   // endfor
1281   SaveHistograms();
1282   cout << "Resolutions: Analyzed " << Nevents << " event(s)" << endl ;
1283   cout << "Resolutions: Total primary       " << totalPrimary << endl ;
1284   cout << "Resoluitons: Total reconstracted " << totalRecPart << endl ;
1285   cout << "TotalReconstructed with Primarie " << totalRPwithPrim << endl ;
1286   cout << "                        Primary:   Photon   Electron   Ch. Hadr.  Neutr. Hadr  Kaons" << endl ;
1287   cout << "             Detected as photon       " << counter[0][0] << "          " << counter[0][1] << "          " << counter[0][2] << "          " <<counter[0][3] << "          " << counter[0][4] << endl ;
1288   cout << "           Detected as electron       " << counter[1][0] << "          " << counter[1][1] << "          " << counter[1][2] << "          " <<counter[1][3] << "          " << counter[1][4] << endl ;
1289   cout << "     Detected as neutral hadron       " << counter[2][0] << "          " << counter[2][1] << "          " << counter[2][2] << "          " <<counter[2][3] << "          " << counter[2][4] << endl ;
1290   cout << "         Detected as neutral EM       " << counter[3][0] << "          " << counter[3][1] << "          " << counter[3][2] << "          " <<counter[3][3] << "          " << counter[3][4] << endl ;
1291   cout << "     Detected as charged hadron       " << counter[4][0] << "          " << counter[4][1] << "          " << counter[4][2] << "          " <<counter[4][3] << "          " << counter[4][4] << endl ;
1292   cout << "       Detected as gamma-hadron       " << counter[5][0] << "          " << counter[5][1] << "          " << counter[5][2] << "          " <<counter[5][3] << "          " << counter[5][4] << endl ;
1293   cout << "          Detected as Absurd EM       " << counter[6][0] << "          " << counter[6][1] << "          " << counter[6][2] << "          " <<counter[6][3] << "          " << counter[6][4] << endl ;
1294   cout << "      Detected as absurd hadron       " << counter[7][0] << "          " << counter[7][1] << "          " << counter[7][2] << "          " <<counter[7][3] << "          " << counter[7][4] << endl ;
1295   cout << "          Detected as undefined       " << counter[8][0] << "          " << counter[8][1] << "          " << counter[8][2] << "          " <<counter[8][3] << "          " << counter[8][4] << endl ;
1296
1297       for(i1 = 0; i1<9; i1++)
1298         for(i2 = 0; i2<5; i2++)
1299           totalInd+=counter[i1][i2] ;
1300       cout << "Indentified particles            " << totalInd << endl ;
1301
1302 }           // endfunction
1303
1304
1305 //____________________________________________________________________________
1306 void  AliPHOSAnalyze::BookingHistograms()
1307 {
1308   // Books the histograms where the results of the analysis are stored (to be changed)
1309
1310   delete fhEmcDigit  ;
1311   delete fhVetoDigit  ;
1312   delete fhConvertorDigit   ;
1313   delete  fhEmcCluster   ;
1314   delete fhVetoCluster   ;
1315   delete fhConvertorCluster  ;
1316   delete fhConvertorEmc  ;
1317
1318   fhEmcDigit                = new TH1F("hEmcDigit",      "hEmcDigit",         1000,  0. ,  25.);
1319   fhVetoDigit               = new TH1F("hVetoDigit",     "hVetoDigit",         500,  0. ,  3.e-5);
1320   fhConvertorDigit          = new TH1F("hConvertorDigit","hConvertorDigit",    500,  0. ,  3.e-5);
1321   fhEmcCluster              = new TH1F("hEmcCluster",    "hEmcCluster",       1000,  0. ,  30.);
1322   fhVetoCluster             = new TH1F("hVetoCluster",   "hVetoCluster",       500,  0. ,  3.e-5);
1323   fhConvertorCluster        = new TH1F("hConvertorCluster","hConvertorCluster",500,  0. ,  3.e-5);
1324   fhConvertorEmc            = new TH2F("hConvertorEmc",  "hConvertorEmc",      200,  1. ,  3., 200, 0., 3.e-5);
1325
1326 }
1327 //____________________________________________________________________________
1328 void  AliPHOSAnalyze::BookResolutionHistograms()
1329 {
1330   // Books the histograms where the results of the Resolution analysis are stored
1331
1332 //   if(fhAllEnergy)
1333 //     delete fhAllEnergy ;
1334 //   if(fhPhotEnergy)
1335 //     delete fhPhotEnergy ;
1336 //   if(fhEMEnergy)
1337 //     delete fhEMEnergy ;
1338 //   if(fhPPSDEnergy)
1339 //     delete fhPPSDEnergy ;
1340
1341
1342   fhAllEnergy  = new TH2F("hAllEnergy",  "Energy of any RP with primary photon",100, 0., 5., 100, 0., 5.);
1343   fhPhotEnergy = new TH2F("hPhotEnergy", "Energy of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.);
1344   fhEMEnergy   = new TH2F("hEMEnergy",   "Energy of EM with primary photon",    100, 0., 5., 100, 0., 5.);
1345   fhPPSDEnergy = new TH2F("hPPSDEnergy", "Energy of PPSD with primary photon",  100, 0., 5., 100, 0., 5.);
1346
1347 //   if(fhAllPosition)
1348 //     delete fhAllPosition ;
1349 //   if(fhPhotPosition)
1350 //     delete fhPhotPosition ;
1351 //   if(fhEMPosition)
1352 //     delete fhEMPosition ;
1353 //   if(fhPPSDPosition)
1354 //     delete fhPPSDPosition ;
1355
1356
1357   fhAllPosition  = new TH2F("hAllPosition",  "Position of any RP with primary photon",100, 0., 5., 100, 0., 5.);
1358   fhPhotPosition = new TH2F("hPhotPosition", "Position of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.);
1359   fhEMPosition   = new TH2F("hEMPosition",   "Position of EM with primary photon",    100, 0., 5., 100, 0., 5.);
1360   fhPPSDPosition = new TH2F("hPPSDPosition", "Position of PPSD with primary photon",  100, 0., 5., 100, 0., 5.);
1361
1362   fhAllPositionX = new TH1F("hAllPositionX", "#Delta X of any RP with primary photon",100, -2., 2.);
1363   fhAllPositionZ = new TH1F("hAllPositionZ", "#Delta X of any RP with primary photon",100, -2., 2.);
1364
1365 //   if(fhAllReg)
1366 //     delete fhAllReg ;
1367 //   if(fhPhotReg)
1368 //     delete fhPhotReg ;
1369 //   if(fhNReg)
1370 //     delete fhNReg ;
1371 //   if(fhNBarReg)
1372 //     delete fhNBarReg ;
1373 //   if(fhChargedReg)
1374 //     delete fhChargedReg ;
1375
1376   fhAllReg    = new TH1F("hAllReg",    "All primaries registered as photon",  100, 0., 5.);
1377   fhPhotReg   = new TH1F("hPhotReg",   "Photon registered as photon",         100, 0., 5.);
1378   fhNReg      = new TH1F("hNReg",      "N registered as photon",              100, 0., 5.);
1379   fhNBarReg   = new TH1F("hNBarReg",   "NBar registered as photon",           100, 0., 5.);
1380   fhChargedReg= new TH1F("hChargedReg", "Charged hadron registered as photon",100, 0., 5.);
1381
1382 //   if(fhAllEM)
1383 //     delete fhAllEM ;
1384 //   if(fhPhotEM)
1385 //     delete fhPhotEM ;
1386 //   if(fhNEM)
1387 //     delete fhNEM ;
1388 //   if(fhNBarEM)
1389 //     delete fhNBarEM ;
1390 //   if(fhChargedEM)
1391 //     delete fhChargedEM ;
1392
1393   fhAllEM    = new TH1F("hAllEM",    "All primary registered as EM",100, 0., 5.);
1394   fhPhotEM   = new TH1F("hPhotEM",   "Photon registered as EM", 100, 0., 5.);
1395   fhNEM      = new TH1F("hNEM",      "N registered as EM",      100, 0., 5.);
1396   fhNBarEM   = new TH1F("hNBarEM",   "NBar registered as EM",   100, 0., 5.);
1397   fhChargedEM= new TH1F("hChargedEM","Charged registered as EM",100, 0., 5.);
1398
1399 //   if(fhAllPPSD)
1400 //     delete fhAllPPSD ;
1401 //   if(fhPhotPPSD)
1402 //     delete fhPhotPPSD ;
1403 //   if(fhNPPSD)
1404 //     delete fhNPPSD ;
1405 //   if(fhNBarPPSD)
1406 //     delete fhNBarPPSD ;
1407 //   if(fhChargedPPSD)
1408 //     delete fhChargedPPSD ;
1409
1410   fhAllPPSD    = new TH1F("hAllPPSD",    "All primary registered as PPSD",100, 0., 5.);
1411   fhPhotPPSD   = new TH1F("hPhotPPSD",   "Photon registered as PPSD", 100, 0., 5.);
1412   fhNPPSD      = new TH1F("hNPPSD",      "N registered as PPSD",      100, 0., 5.);
1413   fhNBarPPSD   = new TH1F("hNBarPPSD",   "NBar registered as PPSD",   100, 0., 5.);
1414   fhChargedPPSD= new TH1F("hChargedPPSD","Charged registered as PPSD",100, 0., 5.);
1415
1416 //   if(fhPrimary)
1417 //     delete fhPrimary ;
1418   fhPrimary= new TH1F("hPrimary", "hPrimary",  100, 0., 5.);
1419
1420 //   if(fhAllRP)
1421 //     delete fhAllRP ;
1422 //   if(fhVeto)
1423 //     delete fhVeto ;
1424 //   if(fhShape)
1425 //     delete fhShape ;
1426 //   if(fhPPSD)
1427 //     delete fhPPSD ;
1428
1429   fhAllRP = new TH1F("hAllRP","All Reconstructed particles",  100, 0., 5.);
1430   fhVeto  = new TH1F("hVeto", "All uncharged particles",      100, 0., 5.);
1431   fhShape = new TH1F("hShape","All particles with EM shaower",100, 0., 5.);
1432   fhPPSD  = new TH1F("hPPSD", "All PPSD photon particles",    100, 0., 5.);
1433
1434
1435 //   if(fhPhotPhot)
1436 //     delete fhPhotPhot ;
1437 //   if(fhPhotElec)
1438 //     delete fhPhotElec ;
1439 //   if(fhPhotNeuH)
1440 //     delete fhPhotNeuH ;
1441 //   if(fhPhotNuEM)
1442 //     delete fhPhotNuEM ;
1443 //   if(fhPhotChHa)
1444 //     delete fhPhotChHa ;
1445 //   if(fhPhotGaHa)
1446 //     delete fhPhotGaHa ;
1447
1448   fhPhotPhot = new TH1F("hPhotPhot","hPhotPhot", 100, 0., 5.);   //Photon registered as photon
1449   fhPhotElec = new TH1F("hPhotElec","hPhotElec", 100, 0., 5.);   //Photon registered as Electron
1450   fhPhotNeuH = new TH1F("hPhotNeuH","hPhotNeuH", 100, 0., 5.);   //Photon registered as Neutral Hadron
1451   fhPhotNuEM = new TH1F("hPhotNuEM","hPhotNuEM", 100, 0., 5.);   //Photon registered as Neutral EM
1452   fhPhotChHa = new TH1F("hPhotChHa","hPhotChHa", 100, 0., 5.);   //Photon registered as Charged Hadron
1453   fhPhotGaHa = new TH1F("hPhotGaHa","hPhotGaHa", 100, 0., 5.);   //Photon registered as Gamma-Hadron
1454 }
1455
1456 //____________________________________________________________________________
1457 Bool_t AliPHOSAnalyze::OpenRootFile(Text_t * name)
1458 {
1459   // Open the root file named "name"
1460
1461   fRootFile   = new TFile(name, "update") ;
1462   return  fRootFile->IsOpen() ;
1463 }
1464
1465 //____________________________________________________________________________
1466 void AliPHOSAnalyze::SaveHistograms()
1467 {
1468   // Saves the histograms in a root file named "name.analyzed"
1469
1470   Text_t outputname[80] ;
1471   sprintf(outputname,"%s.analyzed",fRootFile->GetName());
1472   TFile output(outputname,"RECREATE");
1473   output.cd();
1474
1475   if (fhAllEnergy)
1476     fhAllEnergy->Write() ;
1477   if (fhPhotEnergy)
1478     fhPhotEnergy->Write() ;
1479   if(fhEMEnergy)
1480     fhEMEnergy->Write()  ;
1481   if(fhPPSDEnergy)
1482     fhPPSDEnergy->Write() ;
1483   if(fhAllPosition)
1484     fhAllPosition->Write() ;
1485   if(fhAllPositionX)
1486     fhAllPositionX->Write() ;
1487   if(fhAllPositionZ)
1488     fhAllPositionZ->Write() ;
1489   if(fhPhotPosition)
1490     fhPhotPosition->Write() ;
1491   if(fhEMPosition)
1492     fhEMPosition->Write() ;
1493   if(fhPPSDPosition)
1494     fhPPSDPosition->Write() ;
1495   if (fhAllReg)
1496     fhAllReg->Write() ;
1497   if (fhPhotReg)
1498     fhPhotReg->Write() ;
1499   if(fhNReg)
1500     fhNReg->Write() ;
1501   if(fhNBarReg)
1502     fhNBarReg->Write() ;
1503   if(fhChargedReg)
1504     fhChargedReg->Write() ;
1505   if (fhAllEM)
1506     fhAllEM->Write() ;
1507   if (fhPhotEM)
1508     fhPhotEM->Write() ;
1509   if(fhNEM)
1510     fhNEM->Write() ;
1511   if(fhNBarEM)
1512     fhNBarEM->Write() ;
1513   if(fhChargedEM)
1514     fhChargedEM->Write() ;
1515   if (fhAllPPSD)
1516     fhAllPPSD->Write() ;
1517   if (fhPhotPPSD)
1518     fhPhotPPSD->Write() ;
1519   if(fhNPPSD)
1520     fhNPPSD->Write() ;
1521   if(fhNBarPPSD)
1522     fhNBarPPSD->Write() ;
1523   if(fhChargedPPSD)
1524     fhChargedPPSD->Write() ;
1525   if(fhPrimary)
1526     fhPrimary->Write() ;
1527   if(fhAllRP)
1528     fhAllRP->Write()  ;
1529   if(fhVeto)
1530     fhVeto->Write()  ;
1531   if(fhShape)
1532     fhShape->Write()  ;
1533   if(fhPPSD)
1534     fhPPSD->Write()  ;
1535   if(fhPhotPhot)
1536     fhPhotPhot->Write() ;
1537   if(fhPhotElec)
1538     fhPhotElec->Write() ;
1539   if(fhPhotNeuH)
1540     fhPhotNeuH->Write() ;
1541   if(fhPhotNuEM)
1542     fhPhotNuEM->Write() ;
1543   if(fhPhotNuEM)
1544     fhPhotNuEM->Write() ;
1545   if(fhPhotChHa)
1546     fhPhotChHa->Write() ;
1547   if(fhPhotGaHa)
1548     fhPhotGaHa->Write() ;
1549   if(fhEnergyCorrelations)
1550     fhEnergyCorrelations->Write() ;
1551
1552   output.Write();
1553   output.Close();
1554 }
1555 //____________________________________________________________________________
1556 Float_t AliPHOSAnalyze::CorrectEnergy(Float_t ERecPart)
1557 {
1558   return ERecPart/0.8783 ;
1559 }
1560
1561 //____________________________________________________________________________
1562 void AliPHOSAnalyze::ResetHistograms()
1563 {
1564    fhEnergyCorrelations = 0 ;     //Energy correlations between Eloss in Convertor and PPSD(2)
1565
1566    fhEmcDigit = 0 ;               // Histo of digit energies in the Emc
1567    fhVetoDigit = 0 ;              // Histo of digit energies in the Veto
1568    fhConvertorDigit = 0 ;         // Histo of digit energies in the Convertor
1569    fhEmcCluster = 0 ;             // Histo of Cluster energies in Emc
1570    fhVetoCluster = 0 ;            // Histo of Cluster energies in Veto
1571    fhConvertorCluster = 0 ;       // Histo of Cluster energies in Convertor
1572    fhConvertorEmc = 0 ;           // 2d Convertor versus Emc energies
1573
1574    fhAllEnergy = 0 ;
1575    fhPhotEnergy = 0 ;        // Total spectrum of detected photons
1576    fhEMEnergy = 0 ;         // Spectrum of detected electrons with electron primary
1577    fhPPSDEnergy = 0 ;
1578    fhAllPosition = 0 ;
1579    fhAllPositionX = 0 ;
1580    fhAllPositionZ = 0 ;
1581    fhPhotPosition = 0 ;
1582    fhEMPosition = 0 ;
1583    fhPPSDPosition = 0 ;
1584
1585    fhPhotReg = 0 ;
1586    fhAllReg = 0 ;
1587    fhNReg = 0 ;
1588    fhNBarReg = 0 ;
1589    fhChargedReg = 0 ;
1590    fhPhotEM = 0 ;
1591    fhAllEM = 0 ;
1592    fhNEM = 0 ;
1593    fhNBarEM = 0 ;
1594    fhChargedEM = 0 ;
1595    fhPhotPPSD = 0 ;
1596    fhAllPPSD = 0 ;
1597    fhNPPSD = 0 ;
1598    fhNBarPPSD = 0 ;
1599    fhChargedPPSD = 0 ;
1600
1601    fhPrimary = 0 ;
1602
1603    fhPhotPhot = 0 ;
1604    fhPhotElec = 0 ;
1605    fhPhotNeuH = 0 ;
1606    fhPhotNuEM = 0 ;
1607    fhPhotChHa = 0 ;
1608    fhPhotGaHa = 0 ;
1609
1610 }