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