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