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     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  = (AliPHOSv1 *)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
 236 //        //=========== Cluster in module
 237 //        TIter nextEmc(fPHOS->EmcRecPoints()  ) ;
 238 //        while((emc = (AliPHOSEmcRecPoint *)nextEmc()))
 239 //          {
 240 //            if ( emc->GetPHOSMod() == module )
 241 //              {
 242 //                fhEmcCluster->Fill(  emc->GetTotalEnergy()  );
 243 //                TIter nextPpsd( fPHOS->PpsdRecPoints()) ;
 244 //                while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd()))
 245 //                  {
 246 //                    if ( ppsd->GetPHOSMod() == module )
 247 //                      {
 248 //                        if (!ppsd->GetUp()) fhConvertorEmc->Fill(emc->GetTotalEnergy(),ppsd->GetTotalEnergy()) ;
 249 //                      }
 250 //                  }
 251 //              }
 252 //          }
 253 //        //=========== Cluster in module PPSD Down
 254 //        TIter nextPpsd(fPHOS->PpsdRecPoints() ) ;
 255 //        while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd()))
 256 //          {
 257 //            if ( ppsd->GetPHOSMod() == module )
 258 //              {
 259 //                if (!ppsd->GetUp()) fhConvertorCluster->Fill(ppsd->GetTotalEnergy()) ;
 260 //                if (ppsd->GetUp())  fhVetoCluster     ->Fill(ppsd->GetTotalEnergy()) ;
 261 //              }
 262 //          }
 263           //========== TRackSegments in the event
 264           TIter nextRecParticle(fPHOS->RecParticles() ) ;
 265           while((recparticle = (AliPHOSRecParticle *)nextRecParticle()))
 266             {
 267               if ( recparticle->GetPHOSTrackSegment()->GetPHOSMod() == module )
 268                 {
 269                   cout << "Particle type is " << recparticle->GetType() << endl ;
 270                   Int_t numberofprimaries = 0 ;
 271                   Int_t * listofprimaries = recparticle->GetPrimaries(numberofprimaries) ;
 272                   cout << "Number of primaries = " << numberofprimaries << endl ;
 273                   Int_t index ;
 274                   for ( index = 0 ; index < numberofprimaries ; index++)
 275                     cout << "    primary # " << index << " =  " << listofprimaries[index] << endl ;
 276 //                switch(recparticle->GetType())
 277 //                  {
 278 //                  case AliPHOSFastRecParticle::kGAMMA:
 279 //                    fhPhotonEnergy->Fill(recparticle->Energy() ) ;
 280 //                    //fhPhotonPositionX->Fill(recpart. ) ;
 281 //                    //fhPhotonPositionY->Fill(recpart. ) ;
 282 //                    cout << "PHOTON" << endl;
 283 //                    break;
 284 //                  case  AliPHOSFastRecParticle::kELECTRON:
 285 //                    fhElectronEnergy->Fill(recparticle->Energy() ) ;
 286 //                    //fhElectronPositionX->Fill(recpart. ) ;
 287 //                    //fhElectronPositionY->Fill(recpart. ) ;
 288 //                    cout << "ELECTRON" << endl;
 289 //                    break;
 290 //                  case  AliPHOSFastRecParticle::kNEUTRALHA:
 291 //                    fhNeutralHadronEnergy->Fill(recparticle->Energy() ) ;
 292 //                    //fhNeutralHadronPositionX->Fill(recpart. ) ;
 293 //                    //fhNeutralHadronPositionY->Fill(recpart. ) ;
 294 //                    cout << "NEUTRAl HADRON" << endl;
 295 //                    break ;
 296 //                  case  AliPHOSFastRecParticle::kNEUTRALEM:
 297 //                    fhNeutralEMEnergy->Fill(recparticle->Energy() ) ;
 298 //                    //fhNeutralEMPositionX->Fill(recpart. ) ;
 299 //                    //fhNeutralEMPositionY->Fill(recpart. ) ;
 300 //                    //cout << "NEUTRAL EM" << endl;
 301 //                    break ;
 302 //                  case  AliPHOSFastRecParticle::kCHARGEDHA:
 303 //                    fhChargedHadronEnergy->Fill(recparticle->Energy() ) ;
 304 //                    //fhChargedHadronPositionX->Fill(recpart. ) ;
 305 //                    //fhChargedHadronPositionY->Fill(recpart. ) ;
 306 //                    cout << "CHARGED HADRON" << endl;
 307 //                    break ;
 308 //                  case  AliPHOSFastRecParticle::kGAMMAHA:
 309 //                    fhPhotonHadronEnergy->Fill(recparticle->Energy() ) ;
 310 //                    //fhPhotonHadronPositionX->Fill(recpart. ) ;
 311 //                    //fhPhotonHadronPositionY->Fill(recpart. ) ;
 312 //                    cout << "PHOTON HADRON" << endl;
 313 //                    break ;
 314 //                  }
 315                 }
 316             }
 317           // Deleting fClu, fTrs, fPID et fRec
 318           fClu->Delete();
 319           fTrs->Delete();
 320           fPID->Delete();
 321           fRec->Delete();
 322
 323         }   // endfor
 324       SavingHistograms();
 325     }       // endif
 326 }           // endfunction
 327
 328 //____________________________________________________________________________
 329  void AliPHOSAnalyze::AnalyzeResolutions(Int_t Nevents )
 330 {
 331   // analyzes Nevents events and calculate Energy and Position resolution as well as
 332   // probaility of correct indentifiing of the incident particle
 333
 334   if ( fRootFile == 0 )
 335     cout << "AnalyzeResolutions > " << "Root File not openned" << endl ;
 336   else
 337     {
 338       //========== Get AliRun object from file
 339       gAlice = (AliRun*) fRootFile->Get("gAlice") ;
 340       //=========== Get the PHOS object and associated geometry from the file
 341       fPHOS  = (AliPHOSv1 *)gAlice->GetDetector("PHOS") ;
 342       fGeom  = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() );
 343
 344       //========== Initializes the Index to Object converter
 345       fObjGetter = AliPHOSIndexToObject::GetInstance(fPHOS) ;
 346
 347       //========== Booking Histograms
 348       cout << "AnalyzeResolutions > " << "Booking Histograms" << endl ;
 349       BookResolutionHistograms();
 350
 351       Int_t ievent;
 352       for ( ievent=0; ievent<Nevents; ievent++)
 353         {
 354           if (ievent==0)  cout << "AnalyzeResolutions > " << "Starting Analyzing " << endl ;
 355
 356           //========== Create the Clusterizer
 357           fClu = new AliPHOSClusterizerv1() ;
 358           fClu->SetEmcEnergyThreshold(0.05) ;
 359           fClu->SetEmcClusteringThreshold(0.50) ;
 360           fClu->SetPpsdEnergyThreshold    (0.0000002) ;
 361           fClu->SetPpsdClusteringThreshold(0.0000001) ;
 362           fClu->SetLocalMaxCut(0.03) ;
 363           fClu->SetCalibrationParameters(0., 0.00000001) ;
 364
 365           //========== Creates the track segment maker
 366           fTrs = new AliPHOSTrackSegmentMakerv1()  ;
 367           //      fTrs->UnsetUnfoldFlag() ;
 368
 369           //========== Creates the particle identifier
 370           fPID = new AliPHOSPIDv1() ;
 371           fPID->SetShowerProfileCuts(0.3, 1.8, 0.3, 1.8 ) ;
 372
 373           //========== Creates the Reconstructioner
 374           fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;
 375
 376           //========== Event Number>
 377           if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. )
 378             cout <<  "AnalyzeResolutions > " << "Event is " << ievent << endl ;
 379
 380           //=========== Connects the various Tree's for evt
 381           gAlice->GetEvent(ievent);
 382
 383           //=========== Get the Digit Tree
 384           gAlice->TreeD()->GetEvent(0) ;
 385
 386           //=========== Do the reconstruction
 387           fPHOS->Reconstruction(fRec);
 388
 389           //==========
 390           TClonesArray * RecParticleList     = new  TClonesArray("AliPHOSRecParticle", 2000) ;
 391           TBranch * RecBranch = gAlice->TreeR()->GetBranch("PHOSRP");
 392           RecBranch-> SetAddress(&RecParticleList);
 393
 394           //=========== Gets the Reconstraction TTree
 395           gAlice->TreeR()->GetEvent(0) ;
 396
 397           //=========== Gets the Kine TTree
 398           gAlice->TreeK()->GetEvent(0) ;
 399
 400           //=========== Gets the list of Primari Particles
 401           TClonesArray * PrimaryList  = gAlice->Particles();
 402
 403           //========== TRackSegments in the event
 404           TIter nextRecParticle(RecParticleList) ;
 405           AliPHOSRecParticle * RecParticle ;
 406
 407           while( (RecParticle = (AliPHOSRecParticle *) nextRecParticle()))
 408             {
 409               Int_t ModuleNumberRec ;
 410               Double_t RecX, RecZ ;
 411               fGeom->ImpactOnEmc(RecParticle->Theta(), RecParticle->Phi(), ModuleNumberRec, RecX, RecZ) ;
 412
 413               Double_t MinDistance = 10000 ;
 414               Int_t ClosestPrimary = -1 ;
 415
 416               Int_t numberofprimaries ;
 417               Int_t * listofprimaries  = RecParticle->GetPrimaries(numberofprimaries)  ;
 418
 419               Int_t index ;
 420               TParticle * Primary ;
 421               Double_t Distance = MinDistance ;
 422               for ( index = 0 ; index < numberofprimaries ; index++)
 423                 {
 424                   Primary = (TParticle*)PrimaryList->UncheckedAt(listofprimaries[index]) ;
 425
 426                   Int_t ModuleNumber ;
 427                   Double_t PrimX, PrimZ ;
 428                   fGeom->ImpactOnEmc(Primary->Theta(), Primary->Phi(), ModuleNumber, PrimX, PrimZ) ;
 429                   if(ModuleNumberRec == ModuleNumber)
 430                     Distance = TMath::Sqrt((RecX-PrimX)*(RecX-PrimX)+(RecZ-PrimZ)*(RecZ-PrimZ) ) ;
 431                   if(MinDistance > Distance)
 432                     {
 433                       MinDistance = Distance ;
 434                       ClosestPrimary = listofprimaries[index] ;
 435                     }
 436                 }
 437
 438               if(ClosestPrimary >=0 )
 439                 {
 440                   switch(RecParticle->GetType())
 441                     {
 442                     case AliPHOSFastRecParticle::kGAMMA:
 443                       fhPhotonEnergy->Fill(RecParticle->Energy(),((TParticle *) PrimaryList->At(ClosestPrimary))->Energy() ) ;
 444                       fhPhotonPosition->Fill(RecParticle->Energy(),Distance) ;
 445                       break;
 446                     case  AliPHOSFastRecParticle::kELECTRON:
 447                       fhElectronEnergy->Fill(RecParticle->Energy(),((TParticle *)PrimaryList->At(ClosestPrimary))->Energy() ) ;
 448                       fhElectronPosition->Fill(RecParticle->Energy(),Distance ) ;
 449                       break;
 450                     case  AliPHOSFastRecParticle::kNEUTRALHA:
 451                       fhNeutralHadronEnergy->Fill( RecParticle->Energy(),((TParticle *)PrimaryList->At(ClosestPrimary))->Energy()) ;
 452                       fhNeutralHadronPosition->Fill(RecParticle->Energy(),Distance  ) ;
 453                       break ;
 454                     case  AliPHOSFastRecParticle::kNEUTRALEM:
 455                       fhNeutralEMEnergy->Fill( RecParticle->Energy(),((TParticle *)PrimaryList->At(ClosestPrimary))->Energy() ) ;
 456                       fhNeutralEMPosition->Fill(RecParticle->Energy(),Distance ) ;
 457                       break ;
 458                     case  AliPHOSFastRecParticle::kCHARGEDHA:
 459                       fhChargedHadronEnergy->Fill(RecParticle->Energy(),((TParticle *)PrimaryList->At(ClosestPrimary))->Energy() ) ;
 460                       fhChargedHadronPosition->Fill(RecParticle->Energy(),Distance ) ;
 461                       break ;
 462                     case  AliPHOSFastRecParticle::kGAMMAHA:
 463                       fhPhotonHadronEnergy->Fill( RecParticle->Energy(),((TParticle *)PrimaryList->At(ClosestPrimary))->Energy()) ;
 464                       fhPhotonHadronPosition->Fill(RecParticle->Energy(),Distance ) ;
 465                       break ;
 466                     }
 467                 }
 468             }
 469
 470           // Deleting fClu, fTrs, fPID et fRec
 471           fClu->Delete();
 472           fTrs->Delete();
 473           fPID->Delete();
 474           fRec->Delete();
 475           PrimaryList->Delete() ;
 476
 477         }   // endfor
 478       SaveResolutionHistograms();
 479     }       // endif
 480 }           // endfunction
 481
 482
 483 //____________________________________________________________________________
 484 void  AliPHOSAnalyze::BookingHistograms()
 485 {
 486   // Books the histograms where the results of the analysis are stored (to be changed)
 487
 488   if (fhEmcDigit )
 489     delete fhEmcDigit  ;
 490   if (fhVetoDigit )
 491     delete fhVetoDigit  ;
 492   if (fhConvertorDigit )
 493     delete fhConvertorDigit   ;
 494   if (fhEmcCluster   )
 495     delete  fhEmcCluster   ;
 496   if (fhVetoCluster )
 497     delete fhVetoCluster   ;
 498   if (fhConvertorCluster )
 499     delete fhConvertorCluster  ;
 500   if (fhConvertorEmc )
 501     delete fhConvertorEmc  ;
 502
 503 //   fhEmcDigit                = new TH1F("hEmcDigit",      "hEmcDigit",         1000,  0. ,  25.);
 504 //   fhVetoDigit               = new TH1F("hVetoDigit",     "hVetoDigit",         500,  0. ,  3.e-5);
 505 //   fhConvertorDigit          = new TH1F("hConvertorDigit","hConvertorDigit",    500,  0. ,  3.e-5);
 506 //   fhEmcCluster              = new TH1F("hEmcCluster",    "hEmcCluster",       1000,  0. ,  30.);
 507 //   fhVetoCluster             = new TH1F("hVetoCluster",   "hVetoCluster",       500,  0. ,  3.e-5);
 508 //   fhConvertorCluster        = new TH1F("hConvertorCluster","hConvertorCluster",500,  0. ,  3.e-5);
 509 //   fhConvertorEmc            = new TH2F("hConvertorEmc",  "hConvertorEmc",      200,  1. ,  3., 200, 0., 3.e-5);
 510 //   fhPhotonEnergy            = new TH1F("hPhotonEnergy",  "hPhotonEnergy",     1000,  0. ,  30.);
 511 //   fhElectronEnergy          = new TH1F("hElectronEnergy","hElectronEnergy",   1000,  0. ,  30.);
 512 //   fhNeutralHadronEnergy     = new TH1F("hNeutralHadronEnergy", "hNeutralHadronEnergy",    1000,  0. ,  30.);
 513 //   fhNeutralEMEnergy         = new TH1F("hNeutralEMEnergy", "hNeutralEMEnergy",    1000,  0. ,  30.);
 514 //   fhChargedHadronEnergy     = new TH1F("hChargedHadronEnergy", "hChargedHadronEnergy",    1000,  0. ,  30.);
 515 //   fhPhotonHadronEnergy      = new TH1F("hPhotonHadronEnergy","hPhotonHadronEnergy",500,-80. , 80.);
 516 //   fhPhotonPositionX         = new TH1F("hPhotonPositionX","hPhotonPositionX",   500,-80. , 80.);
 517 //   fhElectronPositionX       = new TH1F("hElectronPositionX","hElectronPositionX",500,-80. , 80.);
 518 //   fhNeutralHadronPositionX  = new TH1F("hNeutralHadronPositionX","hNeutralHadronPositionX",500,-80. , 80.);
 519 //   fhNeutralEMPositionX      = new TH1F("hNeutralEMPositionX","hNeutralEMPositionX",500,-80. , 80.);
 520 //   fhChargedHadronPositionX  = new TH1F("hChargedHadronPositionX","hChargedHadronPositionX",500,-80. , 80.);
 521 //   fhPhotonHadronPositionX   = new TH1F("hPhotonHadronPositionX","hPhotonHadronPositionX",500,-80. , 80.);
 522 //   fhPhotonPositionY         = new TH1F("hPhotonPositionY","hPhotonPositionY",   500,-80. , 80.);
 523 //   fhElectronPositionY       = new TH1F("hElectronPositionY","hElectronPositionY",500,-80. , 80.);
 524 //   fhNeutralHadronPositionY  = new TH1F("hNeutralHadronPositionY","hNeutralHadronPositionY",500,-80. , 80.);
 525 //   fhNeutralEMPositionY      = new TH1F("hNeutralEMPositionY","hNeutralEMPositionY",500,-80. , 80.);
 526 //   fhChargedHadronPositionY  = new TH1F("hChargedHadronPositionY","hChargedHadronPositionY",500,-80. , 80.);
 527 //   fhPhotonHadronPositionY   = new TH1F("hPhotonHadronPositionY","hPhotonHadronPositionY",500,-80. , 80.);
 528
 529
 530 }
 531 //____________________________________________________________________________
 532 void  AliPHOSAnalyze::BookResolutionHistograms()
 533 {
 534   // Books the histograms where the results of the Resolution analysis are stored
 535
 536   if(fhPhotonEnergy)
 537     delete fhPhotonEnergy ;
 538   if(fhElectronEnergy)
 539     delete fhElectronEnergy ;
 540   if(fhNeutralHadronEnergy)
 541     delete fhNeutralHadronEnergy ;
 542   if(fhNeutralEMEnergy)
 543     delete fhNeutralEMEnergy ;
 544   if(fhChargedHadronEnergy)
 545     delete fhChargedHadronEnergy ;
 546   if(fhPhotonHadronEnergy)
 547     delete fhPhotonHadronEnergy ;
 548   if(fhPhotonPosition)
 549     delete fhPhotonPosition ;
 550   if(fhElectronPosition)
 551     delete fhElectronPosition ;
 552   if(fhNeutralHadronPosition)
 553     delete fhNeutralHadronPosition ;
 554   if(fhNeutralEMPosition)
 555     delete fhNeutralEMPosition ;
 556   if(fhChargedHadronPosition)
 557     delete fhChargedHadronPosition ;
 558   if(fhPhotonHadronPosition)
 559     delete fhPhotonHadronPosition ;
 560
 561   fhPhotonEnergy            = new TH2F("hPhotonEnergy",  "hPhotonEnergy",              100, 0., 5., 100, 0., 5.);
 562   fhElectronEnergy          = new TH2F("hElectronEnergy","hElectronEnergy",            100, 0., 5., 100, 0., 5.);
 563   fhNeutralHadronEnergy     = new TH2F("hNeutralHadronEnergy", "hNeutralHadronEnergy", 100, 0., 5., 100, 0., 5.);
 564   fhNeutralEMEnergy         = new TH2F("hNeutralEMEnergy", "hNeutralEMEnergy",         100, 0., 5., 100, 0., 5.);
 565   fhChargedHadronEnergy     = new TH2F("hChargedHadronEnergy", "hChargedHadronEnergy", 100, 0., 5., 100, 0., 5.);
 566   fhPhotonHadronEnergy      = new TH2F("hPhotonHadronEnergy","hPhotonHadronEnergy",    100, 0., 5., 100, 0., 5.);
 567   fhPhotonPosition          = new TH2F("hPhotonPosition","hPhotonPosition",                100, 0., 5., 100, 0., 5.);
 568   fhElectronPosition        = new TH2F("hElectronPosition","hElectronPosition",            100, 0., 5., 100, 0., 5.);
 569   fhNeutralHadronPosition   = new TH2F("hNeutralHadronPosition","hNeutralHadronPosition",  100, 0., 5., 100, 0., 5.);
 570   fhNeutralEMPosition       = new TH2F("hNeutralEMPosition","hNeutralEMPosition",          100, 0., 5., 100, 0., 5.);
 571   fhChargedHadronPosition   = new TH2F("hChargedHadronPosition","hChargedHadronPosition",  100, 0., 5., 100, 0., 5.);
 572   fhPhotonHadronPosition    = new TH2F("hPhotonHadronPosition","hPhotonHadronPosition",    100, 0., 5., 100, 0., 5.);
 573
 574 }
 575 //____________________________________________________________________________
 576 Bool_t AliPHOSAnalyze::Init(Int_t evt)
 577 {
 578   // Do a few initializations: open the root file
 579   //                           get the AliRun object
 580   //                           defines the clusterizer, tracksegment maker and particle identifier
 581   //                           sets the associated parameters
 582
 583   Bool_t ok = kTRUE ;
 584
 585    //========== Open galice root file
 586
 587   if ( fRootFile == 0 ) {
 588     Text_t * name  = new Text_t[80] ;
 589     cout << "AnalyzeOneEvent > Enter file root file name : " ;
 590     cin >> name ;
 591     Bool_t ok = OpenRootFile(name) ;
 592     if ( !ok )
 593       cout << " AliPHOSAnalyze > Error opening " << name << endl ;
 594     else {
 595       //========== Get AliRun object from file
 596
 597       gAlice = (AliRun*) fRootFile->Get("gAlice") ;
 598
 599       //=========== Get the PHOS object and associated geometry from the file
 600
 601       fPHOS  = (AliPHOSv1 *)gAlice->GetDetector("PHOS") ;
 602       fGeom = fPHOS->GetGeometry();
 603       //      fGeom  = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() );
 604
 605     } // else !ok
 606   } // if fRootFile
 607
 608   if ( ok ) {
 609
 610
 611     //========== Initializes the Index to Object converter
 612
 613     fObjGetter = AliPHOSIndexToObject::GetInstance(fPHOS) ;
 614
 615     //========== Create the Clusterizer
 616
 617     fClu =  new AliPHOSClusterizerv1() ;
 618     fClu->SetEmcEnergyThreshold(0.030) ;
 619     fClu->SetEmcClusteringThreshold(0.20) ;
 620     fClu->SetPpsdEnergyThreshold    (0.0000002) ;
 621     fClu->SetPpsdClusteringThreshold(0.0000001) ;
 622     fClu->SetLocalMaxCut(0.03) ;
 623     fClu->SetCalibrationParameters(0., 0.00000001) ;
 624     cout <<  "AnalyzeOneEvent > using clusterizer " << fClu->GetName() << endl ;
 625     fClu->PrintParameters() ;
 626
 627     //========== Creates the track segment maker
 628
 629     fTrs = new AliPHOSTrackSegmentMakerv1() ;
 630     cout <<  "AnalyzeOneEvent > using tack segment maker " << fTrs->GetName() << endl ;
 631     //   fTrs->UnsetUnfoldFlag() ;
 632
 633     //========== Creates the particle identifier
 634
 635     fPID = new AliPHOSPIDv1() ;
 636     cout <<  "AnalyzeOneEvent > using particle identifier " << fPID->GetName() << endl ;
 637     //fPID->SetShowerProfileCuts(Float_t l1m, Float_t l1M, Float_t l2m, Float_t l2M) ;
 638     fPID->SetShowerProfileCuts(0.7, 2.0 , 0.6 , 1.5) ;
 639
 640     //========== Creates the Reconstructioner
 641
 642     fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;
 643     fRec -> SetDebugReconstruction(kFALSE);
 644
 645     //=========== Connect the various Tree's for evt
 646
 647     if ( evt == -999 ) {
 648       cout <<  "AnalyzeOneEvent > Enter event number : " ;
 649       cin >> evt ;
 650       cout << evt << endl ;
 651     }
 652     fEvt = evt ;
 653     gAlice->GetEvent(evt);
 654
 655     //=========== Get the Digit TTree
 656
 657     gAlice->TreeD()->GetEvent(0) ;
 658
 659   } // ok
 660
 661   return ok ;
 662 }
 663
 664
 665 //____________________________________________________________________________
 666 void AliPHOSAnalyze::DisplayKineEvent(Int_t evt)
 667 {
 668   // Display particles from the Kine Tree in global Alice (theta, phi) coordinates.
 669   // One PHOS module at the time.
 670   // The particle type can be selected.
 671
 672   if (evt == -999)
 673     evt = fEvt ;
 674
 675   Int_t module ;
 676   cout <<  "DisplayKineEvent > which module (1-5,  -1: all) ? " ;
 677   cin >> module ; cout << module << endl ;
 678
 679   Int_t testparticle ;
 680   cout << " 22      : PHOTON " << endl
 681        << " (-)11   : (POSITRON)ELECTRON " << endl
 682        << " (-)2112 : (ANTI)NEUTRON " << endl
 683        << " -999    : Everything else " << endl ;
 684   cout  <<  "DisplayKineEvent > enter PDG particle code to display " ;
 685   cin >> testparticle ; cout << testparticle << endl ;
 686
 687   Text_t histoname[80] ;
 688   sprintf(histoname,"Event %d: Incident particles in module %d", evt, module) ;
 689
 690   Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
 691   fGeom->EmcModuleCoverage(module, tm, tM, pm, pM, AliPHOSGeometry::Degre() ) ;
 692
 693   Double_t theta, phi ;
 694   fGeom->EmcXtalCoverage(theta, phi, AliPHOSGeometry::Degre() ) ;
 695
 696   Int_t tdim = (Int_t)( (tM - tm) / theta ) ;
 697   Int_t pdim = (Int_t)( (pM - pm) / phi ) ;
 698
 699   tm -= theta ;
 700   tM += theta ;
 701   pm -= phi ;
 702   pM += phi ;
 703
 704   TH2F * histoparticle = new TH2F("histoparticle",  histoname,
 705                                           pdim, pm, pM, tdim, tm, tM) ;
 706   histoparticle->SetStats(kFALSE) ;
 707
 708   // Get pointers to Alice Particle TClonesArray
 709
 710   TParticle * particle;
 711   TClonesArray * particlearray  = gAlice->Particles();
 712
 713   Text_t canvasname[80];
 714   sprintf(canvasname,"Particles incident in PHOS/EMC module # %d",module) ;
 715   TCanvas * kinecanvas = new TCanvas("kinecanvas", canvasname, 650, 500) ;
 716
 717   // get the KINE Tree
 718
 719   TTree * kine =  gAlice->TreeK() ;
 720   Stat_t nParticles =  kine->GetEntries() ;
 721   cout << "DisplayKineEvent > events in kine " << nParticles << endl ;
 722
 723   // loop over particles
 724
 725   Double_t kRADDEG = 180. / TMath::Pi() ;
 726   Int_t index1 ;
 727   Int_t nparticlein = 0 ;
 728   for (index1 = 0 ; index1 < nParticles ; index1++){
 729     Int_t nparticle = particlearray->GetEntriesFast() ;
 730     Int_t index2 ;
 731     for( index2 = 0 ; index2 < nparticle ; index2++) {
 732       particle            = (TParticle*)particlearray->UncheckedAt(index2) ;
 733       Int_t  particletype = particle->GetPdgCode() ;
 734       if (testparticle == -999 || testparticle == particletype) {
 735         Double_t phi        = particle->Phi() ;
 736         Double_t theta      = particle->Theta() ;
 737         Int_t mod ;
 738         Double_t x, z ;
 739         fGeom->ImpactOnEmc(theta, phi, mod, z, x) ;
 740         if ( mod == module ) {
 741           nparticlein++ ;
 742           if (particle->Energy() >  fClu->GetEmcClusteringThreshold()  )
 743             histoparticle->Fill(phi*kRADDEG, theta*kRADDEG, particle->Energy() ) ;
 744         }
 745       }
 746     }
 747   }
 748   kinecanvas->Draw() ;
 749   histoparticle->Draw("color") ;
 750   TPaveText *  pavetext = new TPaveText(294, 100, 300, 101);
 751   Text_t text[40] ;
 752   sprintf(text, "Particles: %d ", nparticlein) ;
 753   pavetext->AddText(text) ;
 754   pavetext->Draw() ;
 755   kinecanvas->Update();
 756
 757 }
 758 //____________________________________________________________________________
 759 void AliPHOSAnalyze::DisplayRecParticles()
 760 {
 761   // Display reconstructed particles in global Alice(theta, phi) coordinates.
 762   // One PHOS module at the time.
 763   // Click on symbols indicate the reconstructed particle type.
 764
 765   if (fEvt == -999) {
 766     cout << "DisplayRecParticles > Analyze an event first ... (y/n) " ;
 767     Text_t answer[1] ;
 768     cin >> answer ; cout << answer ;
 769     if ( answer == "y" )
 770       AnalyzeOneEvent() ;
 771   }
 772     if (fEvt != -999) {
 773
 774       Int_t module ;
 775       cout <<  "DisplayRecParticles > which module (1-5,  -1: all) ? " ;
 776       cin >> module ; cout << module << endl ;
 777       Text_t histoname[80] ;
 778       sprintf(histoname,"Event %d: Reconstructed particles in module %d", fEvt, module) ;
 779       Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
 780       fGeom->EmcModuleCoverage(module, tm, tM, pm, pM, AliPHOSGeometry::Degre() ) ;
 781       Double_t theta, phi ;
 782       fGeom->EmcXtalCoverage(theta, phi, AliPHOSGeometry::Degre() ) ;
 783       Int_t tdim = (Int_t)( (tM - tm) / theta ) ;
 784       Int_t pdim = (Int_t)( (pM - pm) / phi ) ;
 785       tm -= theta ;
 786       tM += theta ;
 787       pm -= phi ;
 788       TH2F * histoRparticle = new TH2F("histoRparticle",  histoname,
 789                                        pdim, pm, pM, tdim, tm, tM) ;
 790       histoRparticle->SetStats(kFALSE) ;
 791       Text_t canvasname[80] ;
 792       sprintf(canvasname, "Reconstructed particles in PHOSmodule # %d", module) ;
 793       TCanvas * rparticlecanvas = new TCanvas("RparticleCanvas", canvasname, 650, 500) ;
 794       AliPHOSRecParticle::RecParticlesList * rpl = fPHOS->RecParticles() ;
 795       Int_t nRecParticles = rpl->GetEntries() ;
 796       Int_t nRecParticlesInModule = 0 ;
 797       TIter nextRecPart(rpl) ;
 798       AliPHOSRecParticle * rp ;
 799       cout << "DisplayRecParticles > " << nRecParticles << " reconstructed particles " << endl ;
 800       Double_t kRADDEG = 180. / TMath::Pi() ;
 801       while ( (rp = (AliPHOSRecParticle *)nextRecPart() ) ) {
 802         AliPHOSTrackSegment * ts = rp->GetPHOSTrackSegment() ;
 803         if ( ts->GetPHOSMod() == module ) {
 804           Int_t numberofprimaries = 0 ;
 805           Int_t * listofprimaries = 0;
 806           rp->GetPrimaries(numberofprimaries) ;
 807           cout << "Number of primaries = " << numberofprimaries << endl ;
 808           Int_t index ;
 809           for ( index = 0 ; index < numberofprimaries ; index++)
 810             cout << "    primary # " << index << " =  " << listofprimaries[index] << endl ;
 811
 812           nRecParticlesInModule++ ;
 813           Double_t theta = rp->Theta() * kRADDEG ;
 814           Double_t phi   = rp->Phi() * kRADDEG ;
 815           Double_t energy = rp->Energy() ;
 816           histoRparticle->Fill(phi, theta, energy) ;
 817         }
 818       }
 819       histoRparticle->Draw("color") ;
 820
 821       nextRecPart.Reset() ;
 822       while ( (rp = (AliPHOSRecParticle *)nextRecPart() ) ) {
 823         AliPHOSTrackSegment * ts = rp->GetPHOSTrackSegment() ;
 824         if ( ts->GetPHOSMod() == module )
 825           rp->Draw("P") ;
 826       }
 827
 828       Text_t text[80] ;
 829       sprintf(text, "reconstructed particles: %d", nRecParticlesInModule) ;
 830       TPaveText *  pavetext = new TPaveText(292, 100, 300, 101);
 831       pavetext->AddText(text) ;
 832       pavetext->Draw() ;
 833       rparticlecanvas->Update() ;
 834     }
 835 }
 836
 837 //____________________________________________________________________________
 838 void AliPHOSAnalyze::DisplayRecPoints()
 839 {
 840   // Display reconstructed points in local PHOS-module (x, z) coordinates.
 841   // One PHOS module at the time.
 842   // Click on symbols displays the EMC cluster, or PPSD information.
 843
 844   if (fEvt == -999) {
 845     cout << "DisplayRecPoints > Analyze an event first ... (y/n) " ;
 846     Text_t answer[1] ;
 847     cin >> answer ; cout << answer ;
 848     if ( answer == "y" )
 849       AnalyzeOneEvent() ;
 850   }
 851     if (fEvt != -999) {
 852
 853       Int_t module ;
 854       cout <<  "DisplayRecPoints > which module (1-5,  -1: all) ? " ;
 855       cin >> module ; cout << module << endl ;
 856
 857       Text_t canvasname[80];
 858       sprintf(canvasname,"Digits in PHOS/EMC module # %d",module) ;
 859       TCanvas * modulecanvas = new TCanvas("module", canvasname, 650, 500) ;
 860       modulecanvas->Draw() ;
 861
 862       //=========== Creating 2d-histogram of the PHOS module
 863       // a little bit junkie but is used to test Geom functinalities
 864
 865       Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
 866
 867       fGeom->EmcModuleCoverage(module, tm, tM, pm, pM);
 868       // convert angles into coordinates local to the EMC module of interest
 869
 870       Int_t emcModuleNumber ;
 871       Double_t emcModulexm, emcModulezm ; // minimum local coordinate in a given EMCA module
 872       Double_t emcModulexM, emcModulezM ; // maximum local coordinate in a given EMCA module
 873       fGeom->ImpactOnEmc(tm, pm, emcModuleNumber, emcModulezm, emcModulexm) ;
 874       fGeom->ImpactOnEmc(tM, pM, emcModuleNumber, emcModulezM, emcModulexM) ;
 875       Int_t xdim = (Int_t)( ( emcModulexM - emcModulexm ) / fGeom->GetCrystalSize(0) ) ;
 876       Int_t zdim = (Int_t)( ( emcModulezM - emcModulezm ) / fGeom->GetCrystalSize(2) ) ;
 877       Float_t xmin = emcModulexm - fGeom->GetCrystalSize(0) ;
 878       Float_t xMax = emcModulexM + fGeom->GetCrystalSize(0) ;
 879       Float_t zmin = emcModulezm - fGeom->GetCrystalSize(2) ;
 880       Float_t zMax = emcModulezM + fGeom->GetCrystalSize(2) ;
 881       Text_t histoname[80];
 882       sprintf(histoname,"Event %d: Digits and RecPoints in module %d", fEvt, module) ;
 883       TH2F * hModule = new TH2F("HistoReconstructed", histoname,
 884                                 xdim, xmin, xMax, zdim, zmin, zMax) ;
 885       hModule->SetMaximum(2.0);
 886       hModule->SetMinimum(0.0);
 887       hModule->SetStats(kFALSE);
 888
 889       TIter next(fPHOS->Digits()) ;
 890       Float_t energy, y, z;
 891       Float_t etot=0.;
 892       Int_t relid[4]; Int_t nDigits = 0 ;
 893       AliPHOSDigit * digit ;
 894
 895       // Making 2D histogram of the EMC module
 896       while((digit = (AliPHOSDigit *)next()))
 897         {
 898           fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
 899           if (relid[0] == module && relid[1] == 0)
 900             {
 901               energy = fClu->Calibrate(digit->GetAmp()) ;
 902               cout << "Energy is " << energy << " and threshold is " << fClu->GetEmcEnergyThreshold() << endl;
 903               if (energy >  fClu->GetEmcEnergyThreshold()  ){
 904                 nDigits++ ;
 905                 etot += energy ;
 906                 fGeom->RelPosInModule(relid,y,z) ;
 907                 hModule->Fill(y, z, energy) ;
 908               }
 909             }
 910         }
 911       cout <<"DrawRecPoints >  Found in module "
 912            << module << " " << nDigits << "  digits with total energy " << etot << endl ;
 913       hModule->Draw("col2") ;
 914
 915       //=========== Cluster in module
 916
 917       //      TClonesArray * emcRP = fPHOS->EmcClusters() ;
 918       TObjArray * emcRP = fPHOS->EmcRecPoints() ;
 919
 920       etot = 0.;
 921       Int_t totalnClusters = 0 ;
 922       Int_t nClusters = 0 ;
 923       TIter nextemc(emcRP) ;
 924       AliPHOSEmcRecPoint * emc ;
 925       while((emc = (AliPHOSEmcRecPoint *)nextemc()))
 926         {
 927           //      Int_t numberofprimaries ;
 928           //      Int_t * primariesarray = new Int_t[10] ;
 929           //      emc->GetPrimaries(numberofprimaries, primariesarray) ;
 930           totalnClusters++ ;
 931           if ( emc->GetPHOSMod() == module )
 932             {
 933               nClusters++ ;
 934               energy = emc->GetTotalEnergy() ;
 935               etot+= energy ;
 936               emc->Draw("M") ;
 937             }
 938         }
 939       cout << "DrawRecPoints > Found " << totalnClusters << " EMC Clusters in PHOS" << endl ;
 940       cout << "DrawRecPoints > Found in module " << module << "  " << nClusters << " EMC Clusters " << endl ;
 941       cout << "DrawRecPoints > total energy  " << etot << endl ;
 942
 943       TPaveText *  pavetext = new TPaveText(22, 80, 83, 90);
 944       Text_t text[40] ;
 945       sprintf(text, "digits: %d;  clusters: %d", nDigits, nClusters) ;
 946       pavetext->AddText(text) ;
 947       pavetext->Draw() ;
 948       modulecanvas->Update();
 949
 950       //=========== Cluster in module PPSD Down
 951
 952       //      TClonesArray * ppsdRP = fPHOS->PpsdClusters() ;
 953       TObjArray * ppsdRP = fPHOS->PpsdRecPoints() ;
 954
 955       etot = 0.;
 956       TIter nextPpsd(ppsdRP) ;
 957       AliPHOSPpsdRecPoint * ppsd ;
 958       while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd()))
 959         {
 960           totalnClusters++ ;
 961           if ( ppsd->GetPHOSMod() == module )
 962             {
 963               nClusters++ ;
 964               energy = ppsd->GetEnergy() ;
 965               etot+=energy ;
 966               if (!ppsd->GetUp()) ppsd->Draw("P") ;
 967             }
 968         }
 969       cout << "DrawRecPoints > Found " << totalnClusters << " Ppsd Down Clusters in PHOS" << endl ;
 970       cout << "DrawRecPoints > Found in module " << module << "  " << nClusters << " Ppsd Down Clusters " << endl ;
 971       cout << "DrawRecPoints > total energy  " << etot << endl ;
 972
 973       //=========== Cluster in module PPSD Up
 974
 975       ppsdRP = fPHOS->PpsdRecPoints() ;
 976
 977       etot = 0.;
 978       TIter nextPpsdUp(ppsdRP) ;
 979       while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsdUp()))
 980         {
 981           totalnClusters++ ;
 982           if ( ppsd->GetPHOSMod() == module )
 983             {
 984               nClusters++ ;
 985               energy = ppsd->GetEnergy() ;
 986               etot+=energy ;
 987               if (ppsd->GetUp()) ppsd->Draw("P") ;
 988             }
 989         }
 990   cout << "DrawRecPoints > Found " << totalnClusters << " Ppsd Up Clusters in PHOS" << endl ;
 991   cout << "DrawRecPoints > Found in module " << module << "  " << nClusters << " Ppsd Up Clusters " << endl ;
 992   cout << "DrawRecPoints > total energy  " << etot << endl ;
 993
 994     } // if !-999
 995 }
 996
 997 //____________________________________________________________________________
 998 void AliPHOSAnalyze::DisplayTrackSegments()
 999 {
1000   // Display track segments in local PHOS-module (x, z) coordinates.
1001   // One PHOS module at the time.
1002   // One symbol per PHOS subsystem: EMC, upper PPSD, lower PPSD.
1003
1004   if (fEvt == -999) {
1005     cout << "DisplayTrackSegments > Analyze an event first ... (y/n) " ;
1006     Text_t answer[1] ;
1007     cin >> answer ; cout << answer ;
1008     if ( answer == "y" )
1009       AnalyzeOneEvent() ;
1010   }
1011     if (fEvt != -999) {
1012
1013       Int_t module ;
1014       cout <<  "DisplayTrackSegments > which module (1-5,  -1: all) ? " ;
1015       cin >> module ; cout << module << endl ;
1016       //=========== Creating 2d-histogram of the PHOS module
1017       // a little bit junkie but is used to test Geom functinalities
1018
1019       Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
1020
1021       fGeom->EmcModuleCoverage(module, tm, tM, pm, pM);
1022       // convert angles into coordinates local to the EMC module of interest
1023
1024       Int_t emcModuleNumber ;
1025       Double_t emcModulexm, emcModulezm ; // minimum local coordinate in a given EMCA module
1026       Double_t emcModulexM, emcModulezM ; // maximum local coordinate in a given EMCA module
1027       fGeom->ImpactOnEmc(tm, pm, emcModuleNumber, emcModulezm, emcModulexm) ;
1028       fGeom->ImpactOnEmc(tM, pM, emcModuleNumber, emcModulezM, emcModulexM) ;
1029       Int_t xdim = (Int_t)( ( emcModulexM - emcModulexm ) / fGeom->GetCrystalSize(0) ) ;
1030       Int_t zdim = (Int_t)( ( emcModulezM - emcModulezm ) / fGeom->GetCrystalSize(2) ) ;
1031       Float_t xmin = emcModulexm - fGeom->GetCrystalSize(0) ;
1032       Float_t xMax = emcModulexM + fGeom->GetCrystalSize(0) ;
1033       Float_t zmin = emcModulezm - fGeom->GetCrystalSize(2) ;
1034       Float_t zMax = emcModulezM + fGeom->GetCrystalSize(2) ;
1035       Text_t histoname[80];
1036       sprintf(histoname,"Event %d: Track Segments in module %d", fEvt, module) ;
1037       TH2F * histotrack = new TH2F("histotrack",  histoname,
1038                                    xdim, xmin, xMax, zdim, zmin, zMax) ;
1039       histotrack->SetStats(kFALSE);
1040       Text_t canvasname[80];
1041       sprintf(canvasname,"Track segments in PHOS/EMC-PPSD module # %d", module) ;
1042       TCanvas * trackcanvas = new TCanvas("TrackSegmentCanvas", canvasname, 650, 500) ;
1043       histotrack->Draw() ;
1044
1045       AliPHOSTrackSegment::TrackSegmentsList * trsegl = fPHOS->TrackSegments() ;
1046       AliPHOSTrackSegment * trseg ;
1047
1048       Int_t nTrackSegments = trsegl->GetEntries() ;
1049       Int_t index ;
1050       Float_t etot = 0 ;
1051       Int_t nTrackSegmentsInModule = 0 ;
1052       for(index = 0; index < nTrackSegments ; index++){
1053         trseg = (AliPHOSTrackSegment * )trsegl->At(index) ;
1054         etot+= trseg->GetEnergy() ;
1055         if ( trseg->GetPHOSMod() == module ) {
1056           nTrackSegmentsInModule++ ;
1057           trseg->Draw("P");
1058         }
1059       }
1060       Text_t text[80] ;
1061       sprintf(text, "track segments: %d", nTrackSegmentsInModule) ;
1062       TPaveText *  pavetext = new TPaveText(22, 80, 83, 90);
1063       pavetext->AddText(text) ;
1064       pavetext->Draw() ;
1065       trackcanvas->Update() ;
1066       cout << "DisplayTrackSegments > Found " << trsegl->GetEntries() << " Track segments with total energy "<< etot << endl ;
1067
1068    }
1069 }
1070 //____________________________________________________________________________
1071 Bool_t AliPHOSAnalyze::OpenRootFile(Text_t * name)
1072 {
1073   // Open the root file named "name"
1074
1075   fRootFile   = new TFile(name, "update") ;
1076   return  fRootFile->IsOpen() ;
1077 }
1078 //____________________________________________________________________________
1079 void AliPHOSAnalyze::SavingHistograms()
1080 {
1081   // Saves the histograms in a root file named "name.analyzed"
1082
1083   Text_t outputname[80] ;
1084   sprintf(outputname,"%s.analyzed",fRootFile->GetName());
1085   TFile output(outputname,"RECREATE");
1086   output.cd();
1087 //   if (fhEmcDigit )
1088 //     fhEmcDigit->Write()  ;
1089 //   if (fhVetoDigit )
1090 //     fhVetoDigit->Write()  ;
1091 //   if (fhConvertorDigit )
1092 //     fhConvertorDigit->Write()   ;
1093 //   if (fhEmcCluster   )
1094 //     fhEmcCluster->Write()   ;
1095 //   if (fhVetoCluster )
1096 //     fhVetoCluster->Write()   ;
1097 //   if (fhConvertorCluster )
1098 //     fhConvertorCluster->Write()  ;
1099 //   if (fhConvertorEmc )
1100 //     fhConvertorEmc->Write()  ;
1101 //   if (fhPhotonEnergy)
1102 //     fhPhotonEnergy->Write() ;
1103 //   if (fhPhotonPositionX)
1104 //     fhPhotonPositionX->Write() ;
1105 //   if (fhPhotonPositionY)
1106 //     fhPhotonPositionX->Write() ;
1107 //   if (fhElectronEnergy)
1108 //     fhElectronEnergy->Write() ;
1109 //   if (fhElectronPositionX)
1110 //     fhElectronPositionX->Write() ;
1111 //   if (fhElectronPositionY)
1112 //     fhElectronPositionX->Write() ;
1113 //   if (fhNeutralHadronEnergy)
1114 //     fhNeutralHadronEnergy->Write() ;
1115 //   if (fhNeutralHadronPositionX)
1116 //     fhNeutralHadronPositionX->Write() ;
1117 //   if (fhNeutralHadronPositionY)
1118 //     fhNeutralHadronPositionX->Write() ;
1119 //   if (fhNeutralEMEnergy)
1120 //     fhNeutralEMEnergy->Write() ;
1121 //   if (fhNeutralEMPositionX)
1122 //     fhNeutralEMPositionX->Write() ;
1123 //   if (fhNeutralEMPositionY)
1124 //     fhNeutralEMPositionX->Write() ;
1125 //   if (fhChargedHadronEnergy)
1126 //     fhChargedHadronEnergy->Write() ;
1127 //   if (fhChargedHadronPositionX)
1128 //     fhChargedHadronPositionX->Write() ;
1129 //   if (fhChargedHadronPositionY)
1130 //     fhChargedHadronPositionX->Write() ;
1131 //   if (fhPhotonHadronEnergy)
1132 //     fhPhotonHadronEnergy->Write() ;
1133 //   if (fhPhotonHadronPositionX)
1134 //     fhPhotonHadronPositionX->Write() ;
1135 //   if (fhPhotonHadronPositionY)
1136 //     fhPhotonHadronPositionX->Write() ;
1137
1138   output.Write();
1139   output.Close();
1140 }
1141 //____________________________________________________________________________
1142 void AliPHOSAnalyze::SaveResolutionHistograms()
1143 {
1144   // Saves the histograms in a root file named "name.analyzed"
1145
1146   Text_t outputname[80] ;
1147   sprintf(outputname,"%s.analyzed",fRootFile->GetName());
1148   TFile output(outputname,"RECREATE");
1149   output.cd();
1150   if (fhPhotonEnergy)
1151     fhPhotonEnergy->Write() ;
1152   if (fhPhotonPosition)
1153     fhPhotonPosition->Write() ;
1154   if (fhElectronEnergy)
1155     fhElectronEnergy->Write() ;
1156   if (fhElectronPosition)
1157     fhElectronPosition->Write() ;
1158   if (fhNeutralHadronEnergy)
1159     fhNeutralHadronEnergy->Write() ;
1160   if (fhNeutralHadronPosition)
1161     fhNeutralHadronPosition->Write() ;
1162   if (fhNeutralEMEnergy)
1163     fhNeutralEMEnergy->Write() ;
1164   if (fhNeutralEMPosition)
1165     fhNeutralEMPosition->Write() ;
1166   if (fhChargedHadronEnergy)
1167     fhChargedHadronEnergy->Write() ;
1168   if (fhChargedHadronPosition)
1169     fhChargedHadronPosition->Write() ;
1170   if (fhPhotonHadronEnergy)
1171     fhPhotonHadronEnergy->Write() ;
1172   if (fhPhotonHadronPosition)
1173     fhPhotonHadronPosition->Write() ;
1174
1175   output.Write();
1176   output.Close();
1177 }
1178
1179
1180