PHOS/AliPHOSAnalyze.cxx

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