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