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