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