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 //_________________________________________________________________________
  17 // Algorythm class to analyze PHOS events
  18 //*-- Y. Schutz :   SUBATECH
  19 //////////////////////////////////////////////////////////////////////////////
  20
  21 // --- ROOT system ---
  22
  23 #include "TFile.h"
  24 #include "TH1.h"
  25 #include "TPad.h"
  26 #include "TH2.h"
  27 #include "TParticle.h"
  28 #include "TClonesArray.h"
  29 #include "TTree.h"
  30 #include "TMath.h"
  31 #include "TCanvas.h"
  32
  33 // --- Standard library ---
  34
  35 #include <iostream>
  36 #include <cstdio>
  37
  38 // --- AliRoot header files ---
  39
  40 #include "AliRun.h"
  41 #include "AliPHOSAnalyze.h"
  42 #include "AliPHOSClusterizerv1.h"
  43 #include "AliPHOSTrackSegmentMakerv1.h"
  44 #include "AliPHOSPIDv1.h"
  45 #include "AliPHOSReconstructioner.h"
  46 #include "AliPHOSDigit.h"
  47 #include "AliPHOSTrackSegment.h"
  48 #include "AliPHOSRecParticle.h"
  49
  50 ClassImp(AliPHOSAnalyze)
  51
  52
  53 //____________________________________________________________________________
  54   AliPHOSAnalyze::AliPHOSAnalyze()
  55 {
  56   // ctor
  57
  58   fRootFile = 0 ;
  59 }
  60
  61 //____________________________________________________________________________
  62 AliPHOSAnalyze::AliPHOSAnalyze(Text_t * name)
  63 {
  64   // ctor
  65
  66   Bool_t ok = OpenRootFile(name)  ;
  67   if ( !ok ) {
  68     cout << " AliPHOSAnalyze > Error opening " << name << endl ;
  69   }
  70   else {
  71     gAlice = (AliRun*) fRootFile->Get("gAlice");
  72     fPHOS  = (AliPHOSv0 *)gAlice->GetDetector("PHOS") ;
  73     fGeom  = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() ) ;
  74     fEvt = -999 ;
  75   }
  76 }
  77
  78 //____________________________________________________________________________
  79 AliPHOSAnalyze::~AliPHOSAnalyze()
  80 {
  81   // dtor
  82
  83   fRootFile->Close() ;
  84   delete fRootFile ;
  85   fRootFile = 0 ;
  86
  87   delete fPHOS ;
  88   fPHOS = 0 ;
  89
  90   delete fClu ;
  91   fClu = 0 ;
  92
  93   delete fPID ;
  94   fPID = 0 ;
  95
  96   delete fRec ;
  97   fRec = 0 ;
  98
  99   delete fTrs ;
 100   fTrs = 0 ;
 101
 102 }
 103
 104 //____________________________________________________________________________
 105 void AliPHOSAnalyze::AnalyzeOneEvent(Int_t evt)
 106 {
 107   Bool_t ok = Init(evt) ;
 108
 109   if ( ok ) {
 110     //=========== Get the number of entries in the Digits array
 111
 112     Int_t nId = fPHOS->Digits()->GetEntries();
 113     printf("AnalyzeOneEvent > Number of entries in the Digit array is %d \n",nId);
 114
 115     //=========== Do the reconstruction
 116
 117     cout << "AnalyzeOneEvent > Found  " << nId << "  digits in PHOS"   << endl ;
 118
 119     fPHOS->Reconstruction(fRec);
 120
 121     // =========== End of reconstruction
 122
 123     cout << "AnalyzeOneEvent > event # " << fEvt << " processed" << endl ;
 124   } // ok
 125   else
 126     cout << "AnalyzeOneEvent > filed to process event # " << evt << endl ;
 127
 128 }
 129
 130 //____________________________________________________________________________
 131  void AliPHOSAnalyze::AnalyzeManyEvents(Int_t Nevents, Int_t module)   // analyzes many events
 132 {
 133
 134   if ( fRootFile == 0 )
 135     cout << "AnalyzeManyEvents > " << "Root File not openned" << endl ;
 136   else
 137     {
 138       //========== Get AliRun object from file
 139       gAlice = (AliRun*) fRootFile->Get("gAlice") ;
 140       //=========== Get the PHOS object and associated geometry from the file
 141       fPHOS  = (AliPHOSv0 *)gAlice->GetDetector("PHOS") ;
 142       fGeom  = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() );
 143       //========== Booking Histograms
 144       cout << "AnalyzeManyEvents > " << "Booking Histograms" << endl ;
 145       BookingHistograms();
 146       Int_t ievent;
 147       Int_t relid[4] ;
 148       AliPHOSDigit * digit ;
 149       AliPHOSEmcRecPoint * emc ;
 150       AliPHOSPpsdRecPoint * ppsd ;
 151       //      AliPHOSTrackSegment * tracksegment ;
 152       AliPHOSRecParticle * recparticle;
 153       for ( ievent=0; ievent<Nevents; ievent++)
 154         {
 155           if (ievent==0)  cout << "AnalyzeManyEvents > " << "Starting Analyzing " << endl ;
 156           //========== Create the Clusterizer
 157           fClu = new AliPHOSClusterizerv1() ;
 158           fClu->SetEmcEnergyThreshold(0.025) ;
 159           fClu->SetEmcClusteringThreshold(0.75) ;
 160           fClu->SetPpsdEnergyThreshold    (0.0000002) ;
 161           fClu->SetPpsdClusteringThreshold(0.0000001) ;
 162           fClu->SetLocalMaxCut(0.03) ;
 163           fClu->SetCalibrationParameters(0., 0.00000001) ;
 164           //========== Creates the track segment maker
 165           fTrs = new AliPHOSTrackSegmentMakerv1()  ;
 166           fTrs->SetUnfoldFlag() ;
 167           //========== Creates the particle identifier
 168           fPID = new AliPHOSPIDv1() ;
 169           fPID->SetShowerProfileCuts(0.5, 1.5, 0.5, 1.5 ) ;
 170           fPID->Print() ;
 171           //========== Creates the Reconstructioner
 172           fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;
 173           //========== Event Number
 174           if ( ( log10(ievent+1) - (Int_t)(log10(ievent+1)) ) == 0. ) cout <<  "AnalyzeManyEvents > " << "Event is " << ievent << endl ;
 175           //=========== Connects the various Tree's for evt
 176           gAlice->GetEvent(ievent);
 177           //=========== Gets the Digit TTree
 178           gAlice->TreeD()->GetEvent(0) ;
 179           //=========== Gets the number of entries in the Digits array
 180           TIter nextdigit(fPHOS->Digits()) ;
 181           while( ( digit = (AliPHOSDigit *)nextdigit() ) )
 182             {
 183               fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
 184               if (fClu->IsInEmc(digit)) fhEmcDigit->Fill(fClu->Calibrate(digit->GetAmp())) ;
 185               else
 186                 {
 187                   if (relid[1]<17) fhVetoDigit->Fill(fClu->Calibrate(digit->GetAmp()));
 188                   if (relid[1]>16) fhConvertorDigit->Fill(fClu->Calibrate(digit->GetAmp()));
 189                 }
 190             }
 191           //=========== Do the reconstruction
 192           fPHOS->Reconstruction(fRec);
 193           //=========== Cluster in module
 194           TIter nextEmc(fPHOS->EmcClusters()  ) ;
 195           while((emc = (AliPHOSEmcRecPoint *)nextEmc()))
 196             {
 197               if ( emc->GetPHOSMod() == module )
 198                 {
 199                   fhEmcCluster->Fill(  emc->GetTotalEnergy()  );
 200                   TIter nextPpsd( fPHOS->PpsdClusters()) ;
 201                   while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd()))
 202                     {
 203                       if ( ppsd->GetPHOSMod() == module )
 204                         {
 205                           if (!ppsd->GetUp()) fhConvertorEmc->Fill(emc->GetTotalEnergy(),ppsd->GetTotalEnergy()) ;
 206                         }
 207                     }
 208                 }
 209             }
 210           //=========== Cluster in module PPSD Down
 211           TIter nextPpsd(fPHOS->PpsdClusters() ) ;
 212           while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd()))
 213             {
 214               if ( ppsd->GetPHOSMod() == module )
 215                 {
 216                   if (!ppsd->GetUp()) fhConvertorCluster->Fill(ppsd->GetTotalEnergy()) ;
 217                   if (ppsd->GetUp())  fhVetoCluster     ->Fill(ppsd->GetTotalEnergy()) ;
 218                 }
 219             }
 220           //========== TRackSegments in the event
 221           TIter nextRecParticle(fPHOS->RecParticles() ) ;
 222           while((recparticle = (AliPHOSRecParticle *)nextRecParticle()))
 223             {
 224               if ( recparticle->GetPHOSTrackSegment()->GetPHOSMod() == module )
 225                 {
 226                   cout << "Particle type is " << recparticle->GetType() << endl ;
 227                   switch(recparticle->GetType())
 228                     {
 229                     case kGAMMA:
 230                       fhPhotonEnergy->Fill(recparticle->Energy() ) ;
 231                       //fhPhotonPositionX->Fill(recpart. ) ;
 232                       //fhPhotonPositionY->Fill(recpart. ) ;
 233                       cout << "PHOTON" << endl;
 234                       break;
 235                     case kELECTRON:
 236                       fhElectronEnergy->Fill(recparticle->Energy() ) ;
 237                       //fhElectronPositionX->Fill(recpart. ) ;
 238                       //fhElectronPositionY->Fill(recpart. ) ;
 239                       cout << "ELECTRON" << endl;
 240                       break;
 241                     case kNEUTRALHADRON:
 242                       fhNeutralHadronEnergy->Fill(recparticle->Energy() ) ;
 243                       //fhNeutralHadronPositionX->Fill(recpart. ) ;
 244                       //fhNeutralHadronPositionY->Fill(recpart. ) ;
 245                       cout << "NEUTRAl HADRON" << endl;
 246                       break ;
 247                     case kNEUTRALEM:
 248                       fhNeutralEMEnergy->Fill(recparticle->Energy() ) ;
 249                       //fhNeutralEMPositionX->Fill(recpart. ) ;
 250                       //fhNeutralEMPositionY->Fill(recpart. ) ;
 251                       //cout << "NEUTRAL EM" << endl;
 252                       break ;
 253                     case kCHARGEDHADRON:
 254                       fhChargedHadronEnergy->Fill(recparticle->Energy() ) ;
 255                       //fhChargedHadronPositionX->Fill(recpart. ) ;
 256                       //fhChargedHadronPositionY->Fill(recpart. ) ;
 257                       cout << "CHARGED HADRON" << endl;
 258                       break ;
 259
 260                     }
 261                 }
 262             }
 263           // Deleting fClu, fTrs, fPID et fRec
 264           fClu->Delete();
 265           fTrs->Delete();
 266           fPID->Delete();
 267           fRec->Delete();
 268
 269         }   // endfor
 270       SavingHistograms();
 271     }       // endif
 272 }           // endfunction
 273
 274
 275 //____________________________________________________________________________
 276 void  AliPHOSAnalyze::BookingHistograms()
 277 {
 278   if (fhEmcDigit )
 279     delete fhEmcDigit  ;
 280   if (fhVetoDigit )
 281     delete fhVetoDigit  ;
 282   if (fhConvertorDigit )
 283     delete fhConvertorDigit   ;
 284   if (fhEmcCluster   )
 285     delete  fhEmcCluster   ;
 286   if (fhVetoCluster )
 287     delete fhVetoCluster   ;
 288   if (fhConvertorCluster )
 289     delete fhConvertorCluster  ;
 290   if (fhConvertorEmc )
 291     delete fhConvertorEmc  ;
 292
 293   fhEmcDigit                = new TH1F("hEmcDigit",      "hEmcDigit",         1000,  0. ,  25.);
 294   fhVetoDigit               = new TH1F("hVetoDigit",     "hVetoDigit",         500,  0. ,  3.e-5);
 295   fhConvertorDigit          = new TH1F("hConvertorDigit","hConvertorDigit",    500,  0. ,  3.e-5);
 296   fhEmcCluster              = new TH1F("hEmcCluster",    "hEmcCluster",       1000,  0. ,  30.);
 297   fhVetoCluster             = new TH1F("hVetoCluster",   "hVetoCluster",       500,  0. ,  3.e-5);
 298   fhConvertorCluster        = new TH1F("hConvertorCluster","hConvertorCluster",500,  0. ,  3.e-5);
 299   fhConvertorEmc            = new TH2F("hConvertorEmc",  "hConvertorEmc",      200,  1. ,  3., 200, 0., 3.e-5);
 300   fhPhotonEnergy            = new TH1F("hPhotonEnergy",  "hPhotonEnergy",     1000,  0. ,  30.);
 301   fhElectronEnergy          = new TH1F("hElectronEnergy","hElectronEnergy",   1000,  0. ,  30.);
 302   fhNeutralHadronEnergy     = new TH1F("hNeutralHadronEnergy", "hNeutralHadronEnergy",    1000,  0. ,  30.);
 303   fhNeutralEMEnergy         = new TH1F("hNeutralEMEnergy", "hNeutralEMEnergy",    1000,  0. ,  30.);
 304   fhChargedHadronEnergy     = new TH1F("hChargedHadronEnergy", "hChargedHadronEnergy",    1000,  0. ,  30.);
 305   fhPhotonPositionX         = new TH1F("hPhotonPositionX","hPhotonPositionX",   500,-80. , 80.);
 306   fhElectronPositionX       = new TH1F("hElectronPositionX","hElectronPositionX",500,-80. , 80.);
 307   fhNeutralHadronPositionX  = new TH1F("hNeutralHadronPositionX","hNeutralHadronPositionX",500,-80. , 80.);
 308   fhNeutralEMPositionX      = new TH1F("hNeutralEMPositionX","hNeutralEMPositionX",500,-80. , 80.);
 309   fhChargedHadronPositionX  = new TH1F("hChargedHadronPositionX","hChargedHadronPositionX",500,-80. , 80.);
 310   fhPhotonPositionY         = new TH1F("hPhotonPositionY","hPhotonPositionY",   500,-80. , 80.);
 311   fhElectronPositionY       = new TH1F("hElectronPositionY","hElectronPositionY",500,-80. , 80.);
 312   fhNeutralHadronPositionY  = new TH1F("hNeutralHadronPositionY","hNeutralHadronPositionY",500,-80. , 80.);
 313   fhNeutralEMPositionY      = new TH1F("hNeutralEMPositionY","hNeutralEMPositionY",500,-80. , 80.);
 314   fhChargedHadronPositionY  = new TH1F("hChargedHadronPositionY","hChargedHadronPositionY",500,-80. , 80.);
 315
 316 }
 317 //____________________________________________________________________________
 318 Bool_t AliPHOSAnalyze::Init(Int_t evt)
 319 {
 320
 321   Bool_t ok = kTRUE ;
 322
 323    //========== Open galice root file
 324
 325   if ( fRootFile == 0 ) {
 326     Text_t * name  = new Text_t[80] ;
 327     cout << "AnalyzeOneEvent > Enter file root file name : " ;
 328     cin >> name ;
 329     Bool_t ok = OpenRootFile(name) ;
 330     if ( !ok )
 331       cout << " AliPHOSAnalyze > Error opening " << name << endl ;
 332     else {
 333       //========== Get AliRun object from file
 334
 335       gAlice = (AliRun*) fRootFile->Get("gAlice") ;
 336
 337       //=========== Get the PHOS object and associated geometry from the file
 338
 339       fPHOS  = (AliPHOSv0 *)gAlice->GetDetector("PHOS") ;
 340       fGeom  = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() );
 341     } // else !ok
 342   } // if fRootFile
 343
 344   if ( ok ) {
 345
 346     //========== Create the Clusterizer
 347
 348     fClu = new AliPHOSClusterizerv1() ;
 349     fClu->SetEmcEnergyThreshold(0.025) ;
 350     fClu->SetEmcClusteringThreshold(0.75) ;
 351     fClu->SetPpsdEnergyThreshold    (0.0000002) ;
 352     fClu->SetPpsdClusteringThreshold(0.0000001) ;
 353     fClu->SetLocalMaxCut(0.03) ;
 354     fClu->SetCalibrationParameters(0., 0.00000001) ;
 355     cout <<  "AnalyzeOneEvent > using clusterizer " << fClu->GetName() << endl ;
 356     fClu->PrintParameters() ;
 357
 358     //========== Creates the track segment maker
 359
 360     fTrs = new AliPHOSTrackSegmentMakerv1() ;
 361     cout <<  "AnalyzeOneEvent > using tack segment maker " << fTrs->GetName() << endl ;
 362
 363     //========== Creates the particle identifier
 364
 365     fPID = new AliPHOSPIDv1() ;
 366     cout <<  "AnalyzeOneEvent > using particle identifier " << fPID->GetName() << endl ;
 367
 368     //========== Creates the Reconstructioner
 369
 370     fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;
 371
 372     //=========== Connect the various Tree's for evt
 373
 374     if ( evt == -999 ) {
 375       cout <<  "AnalyzeOneEvent > Enter event number : " ;
 376       cin >> evt ;
 377       cout << evt << endl ;
 378     }
 379     fEvt = evt ;
 380     gAlice->GetEvent(evt);
 381
 382     //=========== Get the Digit TTree
 383
 384     gAlice->TreeD()->GetEvent(0) ;
 385
 386   } // ok
 387
 388   return ok ;
 389 }
 390
 391
 392 //____________________________________________________________________________
 393 void AliPHOSAnalyze::DisplayKineEvent(Int_t evt)
 394 {
 395   if (evt == -999)
 396     evt = fEvt ;
 397
 398   Int_t module ;
 399   cout <<  "DisplayKineEvent > which module (1-5,  -1: all) ? " ;
 400   cin >> module ; cout << module << endl ;
 401
 402   Int_t testparticle ;
 403   cout << " 22      : PHOTON " << endl
 404        << " (-)11   : (POSITRON)ELECTRON " << endl
 405        << " (-)2112 : (ANTI)NEUTRON " << endl
 406        << " -999    : Everything else " << endl ;
 407   cout  <<  "DisplayKineEvent > enter PDG particle code to display " ;
 408   cin >> testparticle ; cout << testparticle << endl ;
 409
 410   Text_t histoname[80] ;
 411   sprintf(histoname,"Event %d: Incident particles in module %d", evt, module) ;
 412
 413   Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
 414   fGeom->EmcModuleCoverage(module, tm, tM, pm, pM, kDegre) ;
 415
 416   Double_t theta, phi ;
 417   fGeom->EmcXtalCoverage(theta, phi, kDegre) ;
 418
 419   Int_t tdim = (Int_t)( (tM - tm) / theta ) ;
 420   Int_t pdim = (Int_t)( (pM - pm) / phi ) ;
 421
 422   tm -= theta ;
 423   tM += theta ;
 424   pm -= phi ;
 425   pM += phi ;
 426
 427   TH2F * histoparticle = new TH2F("histoparticle",  histoname,
 428                                           pdim, pm, pM, tdim, tm, tM) ;
 429   histoparticle->SetStats(kFALSE) ;
 430
 431   // Get pointers to Alice Particle TClonesArray
 432
 433   TParticle * particle;
 434   TClonesArray * particlearray  = gAlice->Particles();
 435
 436   Text_t canvasname[80];
 437   sprintf(canvasname,"Particles incident in PHOS/EMC module # %d",module) ;
 438   TCanvas * kinecanvas = new TCanvas("kinecanvas", canvasname, 650, 500) ;
 439
 440   // get the KINE Tree
 441
 442   TTree * kine =  gAlice->TreeK() ;
 443   Stat_t nParticles =  kine->GetEntries() ;
 444   cout << "DisplayKineEvent > events in kine " << nParticles << endl ;
 445
 446   // loop over particles
 447
 448   Double_t kRADDEG = 180. / TMath::Pi() ;
 449   Int_t index1 ;
 450   Int_t nparticlein = 0 ;
 451   for (index1 = 0 ; index1 < nParticles ; index1++){
 452     Int_t nparticle = particlearray->GetEntriesFast() ;
 453     Int_t index2 ;
 454     for( index2 = 0 ; index2 < nparticle ; index2++) {
 455       particle            = (TParticle*)particlearray->UncheckedAt(index2) ;
 456       Int_t  particletype = particle->GetPdgCode() ;
 457       if (testparticle == -999 || testparticle == particletype) {
 458         Double_t phi        = particle->Phi() ;
 459         Double_t theta      = particle->Theta() ;
 460         Int_t mod ;
 461         Double_t x, z ;
 462         fGeom->ImpactOnEmc(theta, phi, mod, z, x) ;
 463         if ( mod == module ) {
 464           nparticlein++ ;
 465           histoparticle->Fill(phi*kRADDEG, theta*kRADDEG, particle->Energy() ) ;
 466         }
 467       }
 468     }
 469   }
 470   kinecanvas->Draw() ;
 471   histoparticle->Draw("color") ;
 472   TPaveText *  pavetext = new TPaveText(294, 100, 300, 101);
 473   Text_t text[40] ;
 474   sprintf(text, "Particles: %d ", nparticlein) ;
 475   pavetext->AddText(text) ;
 476   pavetext->Draw() ;
 477   kinecanvas->Update();
 478
 479 }
 480 //____________________________________________________________________________
 481 void AliPHOSAnalyze::DisplayRecParticles()
 482 {
 483   if (fEvt == -999) {
 484     cout << "DisplayRecPoints > Analyze an event first ... (y/n) " ;
 485     Text_t answer[1] ;
 486     cin >> answer ; cout << answer ;
 487     if ( answer == "y" )
 488       AnalyzeOneEvent() ;
 489   }
 490     if (fEvt != -999) {
 491
 492       Int_t module ;
 493       cout <<  "DisplayRecPoints > which module (1-5,  -1: all) ? " ;
 494       cin >> module ; cout << module << endl ;
 495       Text_t histoname[80] ;
 496       sprintf(histoname,"Event %d: Reconstructed particles in module %d", fEvt, module) ;
 497       Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
 498       fGeom->EmcModuleCoverage(module, tm, tM, pm, pM, kDegre) ;
 499       Double_t theta, phi ;
 500       fGeom->EmcXtalCoverage(theta, phi, kDegre) ;
 501       Int_t tdim = (Int_t)( (tM - tm) / theta ) ;
 502       Int_t pdim = (Int_t)( (pM - pm) / phi ) ;
 503       tm -= theta ;
 504       tM += theta ;
 505       pm -= phi ;
 506       TH2F * histoRparticle = new TH2F("histoRparticle",  histoname,
 507                                        pdim, pm, pM, tdim, tm, tM) ;
 508       histoRparticle->SetStats(kFALSE) ;
 509       Text_t canvasname[80] ;
 510       sprintf(canvasname, "Reconstructed particles in PHOSmodule # %d", module) ;
 511       TCanvas * rparticlecanvas = new TCanvas("RparticleCanvas", canvasname, 650, 500) ;
 512       RecParticlesList * rpl = fPHOS->RecParticles() ;
 513       Int_t nRecParticles = rpl->GetEntries() ;
 514       Int_t nRecParticlesInModule = 0 ;
 515       TIter nextRecPart(rpl) ;
 516       AliPHOSRecParticle * rp ;
 517       cout << "DisplayRecParticles > " << nRecParticles << " reconstructed particles " << endl ;
 518       Double_t kRADDEG = 180. / TMath::Pi() ;
 519       while ( (rp = (AliPHOSRecParticle *)nextRecPart() ) ) {
 520         AliPHOSTrackSegment * ts = rp->GetPHOSTrackSegment() ;
 521         if ( ts->GetPHOSMod() == module ) {
 522           nRecParticlesInModule++ ;
 523           Double_t theta = rp->Theta() * kRADDEG ;
 524           Double_t phi   = rp->Phi() * kRADDEG ;
 525           Double_t energy = rp->Energy() ;
 526           histoRparticle->Fill(phi, theta, energy) ;
 527         }
 528       }
 529       histoRparticle->Draw("color") ;
 530       Text_t text[80] ;
 531       sprintf(text, "reconstructed particles: %d", nRecParticlesInModule) ;
 532       TPaveText *  pavetext = new TPaveText(292, 100, 300, 101);
 533       pavetext->AddText(text) ;
 534       pavetext->Draw() ;
 535       rparticlecanvas->Update() ;
 536     }
 537 }
 538
 539 //____________________________________________________________________________
 540 void AliPHOSAnalyze::DisplayRecPoints()
 541 {
 542   if (fEvt == -999) {
 543     cout << "DisplayRecPoints > Analyze an event first ... (y/n) " ;
 544     Text_t answer[1] ;
 545     cin >> answer ; cout << answer ;
 546     if ( answer == "y" )
 547       AnalyzeOneEvent() ;
 548   }
 549     if (fEvt != -999) {
 550
 551       Int_t module ;
 552       cout <<  "DisplayRecPoints > which module (1-5,  -1: all) ? " ;
 553       cin >> module ; cout << module << endl ;
 554
 555       Text_t canvasname[80];
 556       sprintf(canvasname,"Digits in PHOS/EMC module # %d",module) ;
 557       TCanvas * modulecanvas = new TCanvas("module", canvasname, 650, 500) ;
 558       modulecanvas->Draw() ;
 559
 560       //=========== Creating 2d-histogram of the PHOS module
 561       // a little bit junkie but is used to test Geom functinalities
 562
 563       Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
 564
 565       fGeom->EmcModuleCoverage(module, tm, tM, pm, pM);
 566       // convert angles into coordinates local to the EMC module of interest
 567
 568       Int_t emcModuleNumber ;
 569       Double_t emcModulexm, emcModulezm ; // minimum local coordinate in a given EMCA module
 570       Double_t emcModulexM, emcModulezM ; // maximum local coordinate in a given EMCA module
 571       fGeom->ImpactOnEmc(tm, pm, emcModuleNumber, emcModulezm, emcModulexm) ;
 572       fGeom->ImpactOnEmc(tM, pM, emcModuleNumber, emcModulezM, emcModulexM) ;
 573       Int_t xdim = (Int_t)( ( emcModulexM - emcModulexm ) / fGeom->GetCrystalSize(0) ) ;
 574       Int_t zdim = (Int_t)( ( emcModulezM - emcModulezm ) / fGeom->GetCrystalSize(2) ) ;
 575       Float_t xmin = emcModulexm - fGeom->GetCrystalSize(0) ;
 576       Float_t xMax = emcModulexM + fGeom->GetCrystalSize(0) ;
 577       Float_t zmin = emcModulezm - fGeom->GetCrystalSize(2) ;
 578       Float_t zMax = emcModulezM + fGeom->GetCrystalSize(2) ;
 579       Text_t histoname[80];
 580       sprintf(histoname,"Event %d: Digits and RecPoints in module %d", fEvt, module) ;
 581       TH2F * hModule = new TH2F("HistoReconstructed", histoname,
 582                                 xdim, xmin, xMax, zdim, zmin, zMax) ;
 583       hModule->SetMaximum(2.0);
 584       hModule->SetMinimum(0.0);
 585       hModule->SetStats(kFALSE);
 586
 587       TIter next(fPHOS->Digits()) ;
 588       Float_t energy, y, z;
 589       Float_t etot=0.;
 590       Int_t relid[4]; Int_t nDigits = 0 ;
 591       AliPHOSDigit * digit ;
 592       while((digit = (AliPHOSDigit *)next()))
 593         {
 594           fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
 595           if (relid[0] == module)
 596             {
 597               nDigits++ ;
 598               energy = fClu->Calibrate(digit->GetAmp()) ;
 599               etot += energy ;
 600               fGeom->RelPosInModule(relid,y,z) ;
 601               if (energy > 0.01 )
 602                 hModule->Fill(y, z, energy) ;
 603             }
 604         }
 605       cout <<"DrawRecPoints >  Found in module "
 606            << module << " " << nDigits << "  digits with total energy " << etot << endl ;
 607       hModule->Draw("col2") ;
 608
 609       //=========== Cluster in module
 610
 611       TClonesArray * emcRP = fPHOS->EmcClusters() ;
 612       etot = 0.;
 613       Int_t totalnClusters = 0 ;
 614       Int_t nClusters = 0 ;
 615       TIter nextemc(emcRP) ;
 616       AliPHOSEmcRecPoint * emc ;
 617       while((emc = (AliPHOSEmcRecPoint *)nextemc()))
 618         {
 619           Int_t numberofprimaries ;
 620           Int_t * primariesarray = new Int_t[10] ;
 621           emc->GetPrimaries(numberofprimaries, primariesarray) ;
 622           totalnClusters++ ;
 623           if ( emc->GetPHOSMod() == module )
 624             {
 625               nClusters++ ;
 626               energy = emc->GetTotalEnergy() ;
 627               etot+= energy ;
 628               emc->Draw("M") ;
 629             }
 630         }
 631       cout << "DrawRecPoints > Found " << totalnClusters << " EMC Clusters in PHOS" << endl ;
 632       cout << "DrawRecPoints > Found in module " << module << "  " << nClusters << " EMC Clusters " << endl ;
 633       cout << "DrawRecPoints > total energy  " << etot << endl ;
 634
 635       TPaveText *  pavetext = new TPaveText(22, 80, 83, 90);
 636       Text_t text[40] ;
 637       sprintf(text, "digits: %d;  clusters: %d", nDigits, nClusters) ;
 638       pavetext->AddText(text) ;
 639       pavetext->Draw() ;
 640       modulecanvas->Update();
 641
 642       //=========== Cluster in module PPSD Down
 643
 644       TClonesArray * ppsdRP = fPHOS->PpsdClusters() ;
 645       etot = 0.;
 646       TIter nextPpsd(ppsdRP) ;
 647       AliPHOSPpsdRecPoint * ppsd ;
 648       while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd()))
 649         {
 650           totalnClusters++ ;
 651           if ( ppsd->GetPHOSMod() == module )
 652             {
 653               nClusters++ ;
 654               energy = ppsd->GetEnergy() ;
 655               etot+=energy ;
 656               if (!ppsd->GetUp()) ppsd->Draw("P") ;
 657             }
 658         }
 659       cout << "DrawRecPoints > Found " << totalnClusters << " Ppsd Down Clusters in PHOS" << endl ;
 660       cout << "DrawRecPoints > Found in module " << module << "  " << nClusters << " Ppsd Down Clusters " << endl ;
 661       cout << "DrawRecPoints > total energy  " << etot << endl ;
 662
 663       //=========== Cluster in module PPSD Up
 664
 665       ppsdRP = fPHOS->PpsdClusters() ;
 666       etot = 0.;
 667       TIter nextPpsdUp(ppsdRP) ;
 668       while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsdUp()))
 669         {
 670           totalnClusters++ ;
 671           if ( ppsd->GetPHOSMod() == module )
 672             {
 673               nClusters++ ;
 674               energy = ppsd->GetEnergy() ;
 675               etot+=energy ;
 676               if (ppsd->GetUp()) ppsd->Draw("P") ;
 677             }
 678         }
 679   cout << "DrawRecPoints > Found " << totalnClusters << " Ppsd Up Clusters in PHOS" << endl ;
 680   cout << "DrawRecPoints > Found in module " << module << "  " << nClusters << " Ppsd Up Clusters " << endl ;
 681   cout << "DrawRecPoints > total energy  " << etot << endl ;
 682
 683     } // if !-999
 684 }
 685
 686 //____________________________________________________________________________
 687 void AliPHOSAnalyze::DisplayTrackSegments()
 688 {
 689   if (fEvt == -999) {
 690     cout << "DisplayTrackSegments > Analyze an event first ... (y/n) " ;
 691     Text_t answer[1] ;
 692     cin >> answer ; cout << answer ;
 693     if ( answer == "y" )
 694       AnalyzeOneEvent() ;
 695   }
 696     if (fEvt != -999) {
 697
 698       Int_t module ;
 699       cout <<  "DisplayTrackSegments > which module (1-5,  -1: all) ? " ;
 700       cin >> module ; cout << module << endl ;
 701       //=========== Creating 2d-histogram of the PHOS module
 702       // a little bit junkie but is used to test Geom functinalities
 703
 704       Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
 705
 706       fGeom->EmcModuleCoverage(module, tm, tM, pm, pM);
 707       // convert angles into coordinates local to the EMC module of interest
 708
 709       Int_t emcModuleNumber ;
 710       Double_t emcModulexm, emcModulezm ; // minimum local coordinate in a given EMCA module
 711       Double_t emcModulexM, emcModulezM ; // maximum local coordinate in a given EMCA module
 712       fGeom->ImpactOnEmc(tm, pm, emcModuleNumber, emcModulezm, emcModulexm) ;
 713       fGeom->ImpactOnEmc(tM, pM, emcModuleNumber, emcModulezM, emcModulexM) ;
 714       Int_t xdim = (Int_t)( ( emcModulexM - emcModulexm ) / fGeom->GetCrystalSize(0) ) ;
 715       Int_t zdim = (Int_t)( ( emcModulezM - emcModulezm ) / fGeom->GetCrystalSize(2) ) ;
 716       Float_t xmin = emcModulexm - fGeom->GetCrystalSize(0) ;
 717       Float_t xMax = emcModulexM + fGeom->GetCrystalSize(0) ;
 718       Float_t zmin = emcModulezm - fGeom->GetCrystalSize(2) ;
 719       Float_t zMax = emcModulezM + fGeom->GetCrystalSize(2) ;
 720       Text_t histoname[80];
 721       sprintf(histoname,"Event %d: Track Segments in module %d", fEvt, module) ;
 722       TH2F * histotrack = new TH2F("histotrack",  histoname,
 723                                    xdim, xmin, xMax, zdim, zmin, zMax) ;
 724       histotrack->SetStats(kFALSE);
 725       Text_t canvasname[80];
 726       sprintf(canvasname,"Track segments in PHOS/EMC-PPSD module # %d", module) ;
 727       TCanvas * trackcanvas = new TCanvas("TrackSegmentCanvas", canvasname, 650, 500) ;
 728       histotrack->Draw() ;
 729
 730       TrackSegmentsList * trsegl = fPHOS->TrackSegments() ;
 731       AliPHOSTrackSegment * trseg ;
 732
 733       Int_t nTrackSegments = trsegl->GetEntries() ;
 734       Int_t index ;
 735       Float_t etot = 0 ;
 736       Int_t nTrackSegmentsInModule = 0 ;
 737       for(index = 0; index < nTrackSegments ; index++){
 738         trseg = (AliPHOSTrackSegment * )trsegl->At(index) ;
 739         etot+= trseg->GetEnergy() ;
 740         if ( trseg->GetPHOSMod() == module ) {
 741           nTrackSegmentsInModule++ ;
 742           trseg->Draw("P");
 743         }
 744       }
 745       Text_t text[80] ;
 746       sprintf(text, "track segments: %d", nTrackSegmentsInModule) ;
 747       TPaveText *  pavetext = new TPaveText(22, 80, 83, 90);
 748       pavetext->AddText(text) ;
 749       pavetext->Draw() ;
 750       trackcanvas->Update() ;
 751       cout << "DisplayTrackSegments > Found " << trsegl->GetEntries() << " Track segments with total energy "<< etot << endl ;
 752
 753    }
 754 }
 755 //____________________________________________________________________________
 756 Bool_t AliPHOSAnalyze::OpenRootFile(Text_t * name)
 757 {
 758   fRootFile   = new TFile(name) ;
 759   return  fRootFile->IsOpen() ;
 760 }
 761 //____________________________________________________________________________
 762 void AliPHOSAnalyze::SavingHistograms()
 763 {
 764   Text_t outputname[80] ;// = fRootFile->GetName();
 765   sprintf(outputname,"%s.analyzed",fRootFile->GetName());
 766   TFile output(outputname,"RECREATE");
 767   output.cd();
 768   if (fhEmcDigit )
 769     fhEmcDigit->Write()  ;
 770   if (fhVetoDigit )
 771     fhVetoDigit->Write()  ;
 772   if (fhConvertorDigit )
 773     fhConvertorDigit->Write()   ;
 774   if (fhEmcCluster   )
 775     fhEmcCluster->Write()   ;
 776   if (fhVetoCluster )
 777     fhVetoCluster->Write()   ;
 778   if (fhConvertorCluster )
 779     fhConvertorCluster->Write()  ;
 780   if (fhConvertorEmc )
 781     fhConvertorEmc->Write()  ;
 782   if (fhPhotonEnergy)
 783     fhPhotonEnergy->Write() ;
 784   if (fhPhotonPositionX)
 785     fhPhotonPositionX->Write() ;
 786   if (fhPhotonPositionY)
 787     fhPhotonPositionX->Write() ;
 788   if (fhElectronEnergy)
 789     fhElectronEnergy->Write() ;
 790   if (fhElectronPositionX)
 791     fhElectronPositionX->Write() ;
 792   if (fhElectronPositionY)
 793     fhElectronPositionX->Write() ;
 794   if (fhNeutralHadronEnergy)
 795     fhNeutralHadronEnergy->Write() ;
 796   if (fhNeutralHadronPositionX)
 797     fhNeutralHadronPositionX->Write() ;
 798   if (fhNeutralHadronPositionY)
 799     fhNeutralHadronPositionX->Write() ;
 800   if (fhNeutralEMEnergy)
 801     fhNeutralEMEnergy->Write() ;
 802   if (fhNeutralEMPositionX)
 803     fhNeutralEMPositionX->Write() ;
 804   if (fhNeutralEMPositionY)
 805     fhNeutralEMPositionX->Write() ;
 806   if (fhChargedHadronEnergy)
 807     fhChargedHadronEnergy->Write() ;
 808   if (fhChargedHadronPositionX)
 809     fhChargedHadronPositionX->Write() ;
 810   if (fhChargedHadronPositionY)
 811     fhChargedHadronPositionX->Write() ;
 812
 813   output.Write();
 814   output.Close();
 815 }