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