PHOS/AliPHOSAnalyze.cxx

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