PHOS/AliPHOSAnalyze.cxx

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