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