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 "AliPHOSParticleGuesserv1.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 fPag ;
  94   fPag = 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       for ( ievent=0; ievent<Nevents; ievent++)
 153         {
 154           if (ievent==0)  cout << "AnalyzeManyEvents > " << "Starting Analyzing " << endl ;
 155           //========== Create the Clusterizer
 156           fClu = new AliPHOSClusterizerv1() ;
 157           fClu->SetEmcEnergyThreshold(0.025) ;
 158           fClu->SetEmcClusteringThreshold(0.75) ;
 159           fClu->SetPpsdEnergyThreshold    (0.0000002) ;
 160           fClu->SetPpsdClusteringThreshold(0.0000001) ;
 161           fClu->SetLocalMaxCut(0.03) ;
 162           fClu->SetCalibrationParameters(0., 0.00000001) ;
 163           //========== Creates the track segment maker
 164           fTrs = new AliPHOSTrackSegmentMakerv1()  ;
 165           //========== Creates the particle guesser
 166           fPag = new AliPHOSParticleGuesserv1() ;
 167           //========== Creates the Reconstructioner
 168           fRec = new AliPHOSReconstructioner(fClu, fTrs, fPag) ;
 169           //========== Event Number
 170           if ( ( log10(ievent+1) - (Int_t)(log10(ievent+1)) ) == 0. ) cout <<  "AnalyzeManyEvents > " << "Event is " << ievent << endl ;
 171           //=========== Connects the various Tree's for evt
 172           gAlice->GetEvent(ievent);
 173           //=========== Gets the Digit TTree
 174           gAlice->TreeD()->GetEvent(0) ;
 175           //=========== Gets the number of entries in the Digits array
 176           TIter nextdigit(fPHOS->Digits()) ;
 177           while( ( digit = (AliPHOSDigit *)nextdigit() ) )
 178             {
 179               fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
 180               if (fClu->IsInEmc(digit)) fhEmcDigit->Fill(fClu->Calibrate(digit->GetAmp())) ;
 181               else
 182                 {
 183                   if (relid[1]<17) fhVetoDigit->Fill(fClu->Calibrate(digit->GetAmp()));
 184                   if (relid[1]>16) fhConvertorDigit->Fill(fClu->Calibrate(digit->GetAmp()));
 185                 }
 186             }
 187           //=========== Do the reconstruction
 188           fPHOS->Reconstruction(fRec);
 189           //=========== Cluster in module
 190           TIter nextEmc(fPHOS->EmcClusters()  ) ;
 191           while((emc = (AliPHOSEmcRecPoint *)nextEmc()))
 192             {
 193               if ( emc->GetPHOSMod() == module )
 194                 {
 195                   fhEmcCluster->Fill(  emc->GetTotalEnergy()  );
 196                   TIter nextPpsd( fPHOS->PpsdClusters()) ;
 197                   while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd()))
 198                     {
 199                       if ( ppsd->GetPHOSMod() == module )
 200                         {
 201                           if (!ppsd->GetUp()) fhConvertorEmc->Fill(emc->GetTotalEnergy(),ppsd->GetTotalEnergy()) ;
 202                         }
 203                     }
 204                 }
 205             }
 206           //=========== Cluster in module PPSD Down
 207           TIter nextPpsd(fPHOS->PpsdClusters() ) ;
 208           while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd()))
 209             {
 210               if ( ppsd->GetPHOSMod() == module )
 211                 {
 212                   if (!ppsd->GetUp()) fhConvertorCluster->Fill(ppsd->GetTotalEnergy()) ;
 213                   if (ppsd->GetUp())  fhVetoCluster     ->Fill(ppsd->GetTotalEnergy()) ;
 214                 }
 215             }
 216           //========== TRackSegments in the event
 217           TIter nextTrackSegment(fPHOS->TrackSegments() ) ;
 218           while((tracksegment = (AliPHOSTrackSegment *)nextTrackSegment()))
 219             {
 220               if ( tracksegment->GetPHOSMod() == module )
 221                 {
 222                   switch(tracksegment->GetPartType())
 223                     {
 224                     case 0:
 225                       fhPhotonEnergy->Fill(tracksegment->GetEnergy() ) ;
 226                       //                      fhPhotonPositionX->Fill(tracksegment-> ) ;
 227                       //fhPhotonPositionY->Fill(tracksegment-> ) ;
 228                       //cout << "PHOTON" << endl;
 229                       break;
 230                     case 1 :
 231                       fhElectronEnergy->Fill(tracksegment->GetEnergy() ) ;
 232                       //fhElectronPositionX->Fill(tracksegment-> ) ;
 233                       //fhElectronPositionY->Fill(tracksegment-> ) ;
 234                       //cout << "ELECTRON" << endl;
 235                       break;
 236                     case 2 :
 237                       fhNeutralEnergy->Fill(tracksegment->GetEnergy() ) ;
 238                       //fhNeutralPositionX->Fill(tracksegment-> ) ;
 239                       //fhNeutralPositionY->Fill(tracksegment-> ) ;
 240                       //cout << "NEUTRAL" << endl;
 241                       break ;
 242                     case 3 :
 243                       fhChargedEnergy->Fill(tracksegment->GetEnergy() ) ;
 244                       //fhChargedPositionX->Fill(tracksegment-> ) ;
 245                       //fhChargedPositionY->Fill(tracksegment-> ) ;
 246                       //cout << "CHARGED" << endl;
 247                       break ;
 248
 249                     }
 250                 }
 251             }
 252           // Deleting fClu, fTrs, fPag et fRec
 253           fClu->Delete();
 254           fTrs->Delete();
 255           fPag->Delete();
 256           fRec->Delete();
 257
 258         }   // endfor
 259       SavingHistograms();
 260     }       // endif
 261 }           // endfunction
 262
 263
 264 //____________________________________________________________________________
 265 void  AliPHOSAnalyze::BookingHistograms()
 266 {
 267   if (fhEmcDigit )         delete fhEmcDigit  ;
 268   if (fhVetoDigit )        delete fhVetoDigit  ;
 269   if (fhConvertorDigit )   delete fhConvertorDigit   ;
 270   if (fhEmcCluster   )     delete  fhEmcCluster   ;
 271   if (fhVetoCluster )      delete fhVetoCluster   ;
 272   if (fhConvertorCluster ) delete fhConvertorCluster  ;
 273   if (fhConvertorEmc )     delete fhConvertorEmc  ;
 274   fhEmcDigit         = new TH1F("hEmcDigit",      "hEmcDigit",          500,  0. ,  5.);
 275   fhVetoDigit        = new TH1F("hVetoDigit",     "hVetoDigit",         500,  0. ,  3.e-5);
 276   fhConvertorDigit   = new TH1F("hConvertorDigit","hConvertorDigit",    500,  0. ,  3.e-5);
 277   fhEmcCluster       = new TH1F("hEmcCluster",    "hEmcCluster",        500, -5. ,  5.);
 278   fhVetoCluster      = new TH1F("hVetoCluster",   "hVetoCluster",       500,  0. ,  3.e-5);
 279   fhConvertorCluster = new TH1F("hConvertorCluster","hConvertorCluster",500,  0. ,  3.e-5);
 280   fhConvertorEmc     = new TH2F("hConvertorEmc",  "hConvertorEmc",      200,  1. ,  3., 200, 0., 3.e-5);
 281   fhPhotonEnergy     = new TH1F("hPhotonEnergy",  "hPhotonEnergy",      500,  0. ,  5.);
 282   fhElectronEnergy   = new TH1F("hElectronEnergy","hElectronEnergy",    500,  0. ,  5.);
 283   fhNeutralEnergy    = new TH1F("hNeutralEnergy", "hNeutralEnergy",     500,  0. ,  5.);
 284   fhChargedEnergy    = new TH1F("hChargedEnergy", "hChargedEnergy",     500,  0. ,  5.);
 285   fhPhotonPositionX  = new TH1F("hPhotonPositionX","hPhotonPositionX",   500,-80. , 80.);
 286   fhElectronPositionX= new TH1F("hElectronPositionX","hElectronPositionX",500,-80. , 80.);
 287   fhNeutralPositionX  = new TH1F("hNeutralPositionX","hNeutralPositionX",500,-80. , 80.);
 288   fhChargedPositionX  = new TH1F("hChargedPositionX","hChargedPositionX",500,-80. , 80.);
 289   fhPhotonPositionY  = new TH1F("hPhotonPositionY","hPhotonPositionY",   500,-80. , 80.);
 290   fhElectronPositionY= new TH1F("hElectronPositionY","hElectronPositionY",500,-80. , 80.);
 291   fhNeutralPositionY  = new TH1F("hNeutralPositionY","hNeutralPositionY",500,-80. , 80.);
 292   fhChargedPositionY  = new TH1F("hChargedPositionY","hChargedPositionY",500,-80. , 80.);
 293
 294 }
 295 //____________________________________________________________________________
 296 Bool_t AliPHOSAnalyze::Init(Int_t evt)
 297 {
 298
 299   Bool_t ok = kTRUE ;
 300
 301    //========== Open galice root file
 302
 303   if ( fRootFile == 0 ) {
 304     Text_t * name  = new Text_t[80] ;
 305     cout << "AnalyzeOneEvent > Enter file root file name : " ;
 306     cin >> name ;
 307     Bool_t ok = OpenRootFile(name) ;
 308     if ( !ok )
 309       cout << " AliPHOSAnalyze > Error opening " << name << endl ;
 310     else {
 311       //========== Get AliRun object from file
 312
 313       gAlice = (AliRun*) fRootFile->Get("gAlice") ;
 314
 315       //=========== Get the PHOS object and associated geometry from the file
 316
 317       fPHOS  = (AliPHOSv0 *)gAlice->GetDetector("PHOS") ;
 318       fGeom  = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() );
 319     } // else !ok
 320   } // if fRootFile
 321
 322   if ( ok ) {
 323
 324     //========== Create the Clusterizer
 325
 326     fClu = new AliPHOSClusterizerv1() ;
 327     fClu->SetEmcEnergyThreshold(0.025) ;
 328     fClu->SetEmcClusteringThreshold(0.75) ;
 329     fClu->SetPpsdEnergyThreshold    (0.0000002) ;
 330     fClu->SetPpsdClusteringThreshold(0.0000001) ;
 331     fClu->SetLocalMaxCut(0.03) ;
 332     fClu->SetCalibrationParameters(0., 0.00000001) ;
 333     cout <<  "AnalyzeOneEvent > using clusterizer " << fClu->GetName() << endl ;
 334     fClu->PrintParameters() ;
 335
 336     //========== Creates the track segment maker
 337
 338     fTrs = new AliPHOSTrackSegmentMakerv1() ;
 339     cout <<  "AnalyzeOneEvent > using tack segment maker " << fTrs->GetName() << endl ;
 340
 341     //========== Creates the particle guesser
 342
 343     fPag = new AliPHOSParticleGuesserv1() ;
 344     cout <<  "AnalyzeOneEvent > using particle guess " << fPag->GetName() << endl ;
 345
 346     //========== Creates the Reconstructioner
 347
 348     fRec = new AliPHOSReconstructioner(fClu, fTrs, fPag) ;
 349
 350     //=========== Connect the various Tree's for evt
 351
 352     if ( evt == -999 ) {
 353       cout <<  "AnalyzeOneEvent > Enter event number : " ;
 354       cin >> evt ;
 355       cout << evt << endl ;
 356     }
 357     fEvt = evt ;
 358     gAlice->GetEvent(evt);
 359
 360     //=========== Get the Digit TTree
 361
 362     gAlice->TreeD()->GetEvent(0) ;
 363
 364   } // ok
 365
 366   return ok ;
 367 }
 368
 369
 370 //____________________________________________________________________________
 371 void AliPHOSAnalyze::DisplayKineEvent(Int_t evt)
 372 {
 373   if (evt == -999)
 374     evt = fEvt ;
 375
 376   Int_t module ;
 377   cout <<  "DisplayKineEvent > which module (1-5,  -1: all) ? " ;
 378   cin >> module ; cout << module << endl ;
 379
 380   Int_t testparticle ;
 381   cout << " 22      : PHOTON " << endl
 382        << " (-)11   : (POSITRON)ELECTRON " << endl
 383        << " (-)2112 : (ANTI)NEUTRON " << endl
 384        << " -999    : Everything else " << endl ;
 385   cout  <<  "DisplayKineEvent > enter PDG particle code to display " ;
 386   cin >> testparticle ; cout << testparticle << endl ;
 387
 388   Text_t histoname[80] ;
 389   sprintf(histoname,"Event %d: Incident particles in module %d", evt, module) ;
 390
 391   Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
 392   fGeom->EmcModuleCoverage(module, tm, tM, pm, pM, kDegre) ;
 393
 394   Double_t theta, phi ;
 395   fGeom->EmcXtalCoverage(theta, phi, kDegre) ;
 396
 397   Int_t tdim = (Int_t)( (tM - tm) / theta ) ;
 398   Int_t pdim = (Int_t)( (pM - pm) / phi ) ;
 399
 400   tm -= theta ;
 401   tM += theta ;
 402   pm -= phi ;
 403   pM += phi ;
 404
 405   TH2F * histoparticle = new TH2F("histoparticle",  histoname,
 406                                           pdim, pm, pM, tdim, tm, tM) ;
 407   histoparticle->SetStats(kFALSE) ;
 408
 409   // Get pointers to Alice Particle TClonesArray
 410
 411   TParticle * particle;
 412   TClonesArray * particlearray  = gAlice->Particles();
 413
 414   Text_t canvasname[80];
 415   sprintf(canvasname,"Particles incident in PHOS/EMC module # %d",module) ;
 416   TCanvas * kinecanvas = new TCanvas("kinecanvas", canvasname, 650, 500) ;
 417
 418   // get the KINE Tree
 419
 420   TTree * kine =  gAlice->TreeK() ;
 421   Stat_t nParticles =  kine->GetEntries() ;
 422   cout << "DisplayKineEvent > events in kine " << nParticles << endl ;
 423
 424   // loop over particles
 425
 426   Double_t kRADDEG = 180. / TMath::Pi() ;
 427   Int_t index1 ;
 428   Int_t nparticlein = 0 ;
 429   for (index1 = 0 ; index1 < nParticles ; index1++){
 430     Int_t nparticle = particlearray->GetEntriesFast() ;
 431     Int_t index2 ;
 432     for( index2 = 0 ; index2 < nparticle ; index2++) {
 433       particle            = (TParticle*)particlearray->UncheckedAt(index2) ;
 434       Int_t  particletype = particle->GetPdgCode() ;
 435       if (testparticle == -999 || testparticle == particletype) {
 436         Double_t phi        = particle->Phi() ;
 437         Double_t theta      = particle->Theta() ;
 438         Int_t mod ;
 439         Double_t x, z ;
 440         fGeom->ImpactOnEmc(theta, phi, mod, z, x) ;
 441         if ( mod == module ) {
 442           nparticlein++ ;
 443           histoparticle->Fill(phi*kRADDEG, theta*kRADDEG, particle->Energy() ) ;
 444         }
 445       }
 446     }
 447   }
 448   kinecanvas->Draw() ;
 449   histoparticle->Draw("color") ;
 450   TPaveText *  pavetext = new TPaveText(294, 100, 300, 101);
 451   Text_t text[40] ;
 452   sprintf(text, "Particles: %d ", nparticlein) ;
 453   pavetext->AddText(text) ;
 454   pavetext->Draw() ;
 455   kinecanvas->Update();
 456
 457 }
 458 //____________________________________________________________________________
 459 void AliPHOSAnalyze::DisplayRecParticles()
 460 {
 461   if (fEvt == -999) {
 462     cout << "DisplayRecPoints > Analyze an event first ... (y/n) " ;
 463     Text_t answer[1] ;
 464     cin >> answer ; cout << answer ;
 465     if ( answer == "y" )
 466       AnalyzeOneEvent() ;
 467   }
 468     if (fEvt != -999) {
 469
 470       Int_t module ;
 471       cout <<  "DisplayRecPoints > which module (1-5,  -1: all) ? " ;
 472       cin >> module ; cout << module << endl ;
 473       Text_t histoname[80] ;
 474       sprintf(histoname,"Event %d: Reconstructed particles in module %d", fEvt, module) ;
 475       Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
 476       fGeom->EmcModuleCoverage(module, tm, tM, pm, pM, kDegre) ;
 477       Double_t theta, phi ;
 478       fGeom->EmcXtalCoverage(theta, phi, kDegre) ;
 479       Int_t tdim = (Int_t)( (tM - tm) / theta ) ;
 480       Int_t pdim = (Int_t)( (pM - pm) / phi ) ;
 481       tm -= theta ;
 482       tM += theta ;
 483       pm -= phi ;
 484       TH2F * histoRparticle = new TH2F("histoRparticle",  histoname,
 485                                        pdim, pm, pM, tdim, tm, tM) ;
 486       histoRparticle->SetStats(kFALSE) ;
 487       Text_t canvasname[80] ;
 488       sprintf(canvasname, "Reconstructed particles in PHOSmodule # %d", module) ;
 489       TCanvas * rparticlecanvas = new TCanvas("RparticleCanvas", canvasname, 650, 500) ;
 490       RecParticlesList * rpl = fPHOS->RecParticles() ;
 491       Int_t nRecParticles = rpl->GetEntries() ;
 492       Int_t nRecParticlesInModule = 0 ;
 493       TIter nextRecPart(rpl) ;
 494       AliPHOSRecParticle * rp ;
 495       cout << "DisplayRecParticles > " << nRecParticles << " reconstructed particles " << endl ;
 496       Double_t kRADDEG = 180. / TMath::Pi() ;
 497       while ( (rp = (AliPHOSRecParticle *)nextRecPart() ) ) {
 498         AliPHOSTrackSegment * ts = rp->GetPHOSTrackSegment() ;
 499         if ( ts->GetPHOSMod() == module ) {
 500           nRecParticlesInModule++ ;
 501           Double_t theta = rp->Theta() * kRADDEG ;
 502           Double_t phi   = rp->Phi() * kRADDEG ;
 503           Double_t energy = rp->Energy() ;
 504           histoRparticle->Fill(phi, theta, energy) ;
 505         }
 506       }
 507       histoRparticle->Draw("color") ;
 508       Text_t text[80] ;
 509       sprintf(text, "reconstructed particles: %d", nRecParticlesInModule) ;
 510       TPaveText *  pavetext = new TPaveText(292, 100, 300, 101);
 511       pavetext->AddText(text) ;
 512       pavetext->Draw() ;
 513       rparticlecanvas->Update() ;
 514     }
 515 }
 516
 517 //____________________________________________________________________________
 518 void AliPHOSAnalyze::DisplayRecPoints()
 519 {
 520   if (fEvt == -999) {
 521     cout << "DisplayRecPoints > Analyze an event first ... (y/n) " ;
 522     Text_t answer[1] ;
 523     cin >> answer ; cout << answer ;
 524     if ( answer == "y" )
 525       AnalyzeOneEvent() ;
 526   }
 527     if (fEvt != -999) {
 528
 529       Int_t module ;
 530       cout <<  "DisplayRecPoints > which module (1-5,  -1: all) ? " ;
 531       cin >> module ; cout << module << endl ;
 532
 533       Text_t canvasname[80];
 534       sprintf(canvasname,"Digits in PHOS/EMC module # %d",module) ;
 535       TCanvas * modulecanvas = new TCanvas("module", canvasname, 650, 500) ;
 536       modulecanvas->Draw() ;
 537
 538       //=========== Creating 2d-histogram of the PHOS module
 539       // a little bit junkie but is used to test Geom functinalities
 540
 541       Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
 542
 543       fGeom->EmcModuleCoverage(module, tm, tM, pm, pM);
 544       // convert angles into coordinates local to the EMC module of interest
 545
 546       Int_t emcModuleNumber ;
 547       Double_t emcModulexm, emcModulezm ; // minimum local coordinate in a given EMCA module
 548       Double_t emcModulexM, emcModulezM ; // maximum local coordinate in a given EMCA module
 549       fGeom->ImpactOnEmc(tm, pm, emcModuleNumber, emcModulezm, emcModulexm) ;
 550       fGeom->ImpactOnEmc(tM, pM, emcModuleNumber, emcModulezM, emcModulexM) ;
 551       Int_t xdim = (Int_t)( ( emcModulexM - emcModulexm ) / fGeom->GetCrystalSize(0) ) ;
 552       Int_t zdim = (Int_t)( ( emcModulezM - emcModulezm ) / fGeom->GetCrystalSize(2) ) ;
 553       Float_t xmin = emcModulexm - fGeom->GetCrystalSize(0) ;
 554       Float_t xMax = emcModulexM + fGeom->GetCrystalSize(0) ;
 555       Float_t zmin = emcModulezm - fGeom->GetCrystalSize(2) ;
 556       Float_t zMax = emcModulezM + fGeom->GetCrystalSize(2) ;
 557       Text_t histoname[80];
 558       sprintf(histoname,"Event %d: Digits and RecPoints in module %d", fEvt, module) ;
 559       TH2F * hModule = new TH2F("HistoReconstructed", histoname,
 560                                 xdim, xmin, xMax, zdim, zmin, zMax) ;
 561       hModule->SetMaximum(2.0);
 562       hModule->SetMinimum(0.0);
 563       hModule->SetStats(kFALSE);
 564
 565       TIter next(fPHOS->Digits()) ;
 566       Float_t energy, y, z;
 567       Float_t etot=0.;
 568       Int_t relid[4]; Int_t nDigits = 0 ;
 569       AliPHOSDigit * digit ;
 570       while((digit = (AliPHOSDigit *)next()))
 571         {
 572           fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
 573           if (relid[0] == module)
 574             {
 575               nDigits++ ;
 576               energy = fClu->Calibrate(digit->GetAmp()) ;
 577               etot += energy ;
 578               fGeom->RelPosInModule(relid,y,z) ;
 579               if (energy > 0.01 )
 580                 hModule->Fill(y, z, energy) ;
 581             }
 582         }
 583       cout <<"DrawRecPoints >  Found in module "
 584            << module << " " << nDigits << "  digits with total energy " << etot << endl ;
 585       hModule->Draw("col2") ;
 586
 587       //=========== Cluster in module
 588
 589       TClonesArray * emcRP = fPHOS->EmcClusters() ;
 590       etot = 0.;
 591       Int_t totalnClusters = 0 ;
 592       Int_t nClusters = 0 ;
 593       TIter nextemc(emcRP) ;
 594       AliPHOSEmcRecPoint * emc ;
 595       while((emc = (AliPHOSEmcRecPoint *)nextemc()))
 596         {
 597           totalnClusters++ ;
 598           if ( emc->GetPHOSMod() == module )
 599             {
 600               nClusters++ ;
 601               energy = emc->GetTotalEnergy() ;
 602               etot+= energy ;
 603               emc->Draw("P") ;
 604             }
 605         }
 606       cout << "DrawRecPoints > Found " << totalnClusters << " EMC Clusters in PHOS" << endl ;
 607       cout << "DrawRecPoints > Found in module " << module << "  " << nClusters << " EMC Clusters " << endl ;
 608       cout << "DrawRecPoints > total energy  " << etot << endl ;
 609
 610       TPaveText *  pavetext = new TPaveText(22, 80, 83, 90);
 611       Text_t text[40] ;
 612       sprintf(text, "digits: %d;  clusters: %d", nDigits, nClusters) ;
 613       pavetext->AddText(text) ;
 614       pavetext->Draw() ;
 615       modulecanvas->Update();
 616
 617       //=========== Cluster in module PPSD Down
 618
 619       TClonesArray * ppsdRP = fPHOS->PpsdClusters() ;
 620       etot = 0.;
 621       TIter nextPpsd(ppsdRP) ;
 622       AliPHOSPpsdRecPoint * ppsd ;
 623       while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd()))
 624         {
 625           totalnClusters++ ;
 626           if ( ppsd->GetPHOSMod() == module )
 627             {
 628               nClusters++ ;
 629               energy = ppsd->GetEnergy() ;
 630               etot+=energy ;
 631               if (!ppsd->GetUp()) ppsd->Draw("P") ;
 632             }
 633         }
 634       cout << "DrawRecPoints > Found " << totalnClusters << " Ppsd Down Clusters in PHOS" << endl ;
 635       cout << "DrawRecPoints > Found in module " << module << "  " << nClusters << " Ppsd Down Clusters " << endl ;
 636       cout << "DrawRecPoints > total energy  " << etot << endl ;
 637
 638       //=========== Cluster in module PPSD Up
 639
 640       ppsdRP = fPHOS->PpsdClusters() ;
 641       etot = 0.;
 642       TIter nextPpsdUp(ppsdRP) ;
 643       while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsdUp()))
 644         {
 645           totalnClusters++ ;
 646           if ( ppsd->GetPHOSMod() == module )
 647             {
 648               nClusters++ ;
 649               energy = ppsd->GetEnergy() ;
 650               etot+=energy ;
 651               if (ppsd->GetUp()) ppsd->Draw("P") ;
 652             }
 653         }
 654   cout << "DrawRecPoints > Found " << totalnClusters << " Ppsd Up Clusters in PHOS" << endl ;
 655   cout << "DrawRecPoints > Found in module " << module << "  " << nClusters << " Ppsd Up Clusters " << endl ;
 656   cout << "DrawRecPoints > total energy  " << etot << endl ;
 657
 658     } // if !-999
 659 }
 660
 661 //____________________________________________________________________________
 662 void AliPHOSAnalyze::DisplayTrackSegments()
 663 {
 664   if (fEvt == -999) {
 665     cout << "DisplayTrackSegments > Analyze an event first ... (y/n) " ;
 666     Text_t answer[1] ;
 667     cin >> answer ; cout << answer ;
 668     if ( answer == "y" )
 669       AnalyzeOneEvent() ;
 670   }
 671     if (fEvt != -999) {
 672
 673       Int_t module ;
 674       cout <<  "DisplayTrackSegments > which module (1-5,  -1: all) ? " ;
 675       cin >> module ; cout << module << endl ;
 676       //=========== Creating 2d-histogram of the PHOS module
 677       // a little bit junkie but is used to test Geom functinalities
 678
 679       Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
 680
 681       fGeom->EmcModuleCoverage(module, tm, tM, pm, pM);
 682       // convert angles into coordinates local to the EMC module of interest
 683
 684       Int_t emcModuleNumber ;
 685       Double_t emcModulexm, emcModulezm ; // minimum local coordinate in a given EMCA module
 686       Double_t emcModulexM, emcModulezM ; // maximum local coordinate in a given EMCA module
 687       fGeom->ImpactOnEmc(tm, pm, emcModuleNumber, emcModulezm, emcModulexm) ;
 688       fGeom->ImpactOnEmc(tM, pM, emcModuleNumber, emcModulezM, emcModulexM) ;
 689       Int_t xdim = (Int_t)( ( emcModulexM - emcModulexm ) / fGeom->GetCrystalSize(0) ) ;
 690       Int_t zdim = (Int_t)( ( emcModulezM - emcModulezm ) / fGeom->GetCrystalSize(2) ) ;
 691       Float_t xmin = emcModulexm - fGeom->GetCrystalSize(0) ;
 692       Float_t xMax = emcModulexM + fGeom->GetCrystalSize(0) ;
 693       Float_t zmin = emcModulezm - fGeom->GetCrystalSize(2) ;
 694       Float_t zMax = emcModulezM + fGeom->GetCrystalSize(2) ;
 695       Text_t histoname[80];
 696       sprintf(histoname,"Event %d: Track Segments in module %d", fEvt, module) ;
 697       TH2F * histotrack = new TH2F("histotrack",  histoname,
 698                                    xdim, xmin, xMax, zdim, zmin, zMax) ;
 699       histotrack->SetStats(kFALSE);
 700       Text_t canvasname[80];
 701       sprintf(canvasname,"Track segments in PHOS/EMC-PPSD module # %d", module) ;
 702       TCanvas * trackcanvas = new TCanvas("TrackSegmentCanvas", canvasname, 650, 500) ;
 703       histotrack->Draw() ;
 704
 705       TrackSegmentsList * trsegl = fPHOS->TrackSegments() ;
 706       AliPHOSTrackSegment * trseg ;
 707
 708       Int_t nTrackSegments = trsegl->GetEntries() ;
 709       Int_t index ;
 710       Float_t etot = 0 ;
 711       Int_t nTrackSegmentsInModule = 0 ;
 712       for(index = 0; index < nTrackSegments ; index++){
 713         trseg = (AliPHOSTrackSegment * )trsegl->At(index) ;
 714         etot+= trseg->GetEnergy() ;
 715         if ( trseg->GetPHOSMod() == module ) {
 716           nTrackSegmentsInModule++ ;
 717           trseg->Draw("P");
 718         }
 719       }
 720       Text_t text[80] ;
 721       sprintf(text, "track segments: %d", nTrackSegmentsInModule) ;
 722       TPaveText *  pavetext = new TPaveText(22, 80, 83, 90);
 723       pavetext->AddText(text) ;
 724       pavetext->Draw() ;
 725       trackcanvas->Update() ;
 726       cout << "DisplayTrackSegments > Found " << trsegl->GetEntries() << " Track segments with total energy "<< etot << endl ;
 727
 728    }
 729 }
 730 //____________________________________________________________________________
 731 Bool_t AliPHOSAnalyze::OpenRootFile(Text_t * name)
 732 {
 733   fRootFile   = new TFile(name) ;
 734   return  fRootFile->IsOpen() ;
 735 }
 736 //____________________________________________________________________________
 737 void AliPHOSAnalyze::SavingHistograms()
 738 {
 739   Text_t outputname[80] ;// = fRootFile->GetName();
 740   sprintf(outputname,"%s.analyzed",fRootFile->GetName());
 741   TFile output(outputname,"RECREATE");
 742   output.cd();
 743   if (fhEmcDigit )         fhEmcDigit->Write()  ;
 744   if (fhVetoDigit )        fhVetoDigit->Write()  ;
 745   if (fhConvertorDigit )   fhConvertorDigit->Write()   ;
 746   if (fhEmcCluster   )     fhEmcCluster->Write()   ;
 747   if (fhVetoCluster )      fhVetoCluster->Write()   ;
 748   if (fhConvertorCluster ) fhConvertorCluster->Write()  ;
 749   if (fhConvertorEmc )     fhConvertorEmc->Write()  ;
 750   if (fhPhotonEnergy)      fhPhotonEnergy->Write() ;
 751   if (fhPhotonPositionX)   fhPhotonPositionX->Write() ;
 752   if (fhPhotonPositionY)   fhPhotonPositionX->Write() ;
 753   if (fhElectronEnergy)      fhElectronEnergy->Write() ;
 754   if (fhElectronPositionX) fhElectronPositionX->Write() ;
 755   if (fhElectronPositionY) fhElectronPositionX->Write() ;
 756   if (fhNeutralEnergy)     fhNeutralEnergy->Write() ;
 757   if (fhNeutralPositionX)  fhNeutralPositionX->Write() ;
 758   if (fhNeutralPositionY)  fhNeutralPositionX->Write() ;
 759   if (fhChargedEnergy)     fhChargedEnergy->Write() ;
 760   if (fhChargedPositionX)  fhChargedPositionX->Write() ;
 761   if (fhChargedPositionY)  fhChargedPositionX->Write() ;
 762
 763   output.Write();
 764   output.Close();
 765 }