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 #include "TStyle.h"
  38
  39 // --- Standard library ---
  40
  41 #include <iostream.h>
  42 #include <stdio.h>
  43
  44 // --- AliRoot header files ---
  45
  46 #include "AliRun.h"
  47 #include "AliPHOSAnalyze.h"
  48 #include "AliPHOSClusterizerv1.h"
  49 #include "AliPHOSTrackSegmentMakerv1.h"
  50 #include "AliPHOSPIDv1.h"
  51 #include "AliPHOSReconstructioner.h"
  52 #include "AliPHOSDigit.h"
  53 #include "AliPHOSTrackSegment.h"
  54 #include "AliPHOSRecParticle.h"
  55 #include "AliPHOSIndexToObject.h"
  56 #include "AliPHOSCPV.h"
  57
  58 ClassImp(AliPHOSAnalyze)
  59
  60
  61 //____________________________________________________________________________
  62   AliPHOSAnalyze::AliPHOSAnalyze()
  63 {
  64   // default ctor (useless)
  65
  66   fRootFile = 0 ;
  67 }
  68
  69 //____________________________________________________________________________
  70 AliPHOSAnalyze::AliPHOSAnalyze(Text_t * name)
  71 {
  72   // ctor: analyze events from root file "name"
  73
  74   Bool_t ok = OpenRootFile(name)  ;
  75   if ( !ok ) {
  76     cout << " AliPHOSAnalyze > Error opening " << name << endl ;
  77   }
  78   else {
  79       //========== Get AliRun object from file
  80       gAlice = (AliRun*) fRootFile->Get("gAlice") ;
  81
  82       //=========== Get the PHOS object and associated geometry from the file
  83       fPHOS  = (AliPHOSv1 *)gAlice->GetDetector("PHOS") ;
  84       fGeom  = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() );
  85
  86       //========== Initializes the Index to Object converter
  87       fObjGetter = AliPHOSIndexToObject::GetInstance(fPHOS) ;
  88       //========== Current event number
  89       fEvt = -999 ;
  90
  91   }
  92   fDebugLevel = 0;
  93   ResetHistograms() ;
  94 }
  95
  96 //____________________________________________________________________________
  97 AliPHOSAnalyze::AliPHOSAnalyze(const AliPHOSAnalyze & ana)
  98 {
  99   // copy ctor
 100   ( (AliPHOSAnalyze &)ana ).Copy(*this) ;
 101 }
 102
 103 //____________________________________________________________________________
 104 void AliPHOSAnalyze::Copy(TObject & obj)
 105 {
 106   // copy an analysis into an other one
 107   TObject::Copy(obj) ;
 108   // I do nothing more because the copy is silly but the Code checkers requires one
 109 }
 110
 111 //____________________________________________________________________________
 112 AliPHOSAnalyze::~AliPHOSAnalyze()
 113 {
 114   // dtor
 115
 116   if (fRootFile->IsOpen() )
 117     fRootFile->Close() ;
 118   if(fRootFile) {delete fRootFile ; fRootFile=0 ;}
 119   if(fPHOS)     {delete fPHOS     ; fPHOS    =0 ;}
 120   if(fClu)      {delete fClu      ; fClu     =0 ;}
 121   if(fPID)      {delete fPID      ; fPID     =0 ;}
 122   if(fRec)      {delete fRec      ; fRec     =0 ;}
 123   if(fTrs)      {delete fTrs      ; fTrs     =0 ;}
 124
 125 }
 126
 127 //____________________________________________________________________________
 128 void AliPHOSAnalyze::ActivePPSD(Int_t Nevents=1){
 129
 130   fhEnergyCorrelations  = new TH2F("hEnergyCorrelations","hEnergyCorrelations",40,  0., 0.15, 30, 0., 3.e-5);
 131   //========== Create the Clusterizer
 132   fClu = new AliPHOSClusterizerv1() ;
 133   fClu->SetEmcEnergyThreshold(0.05) ;
 134   fClu->SetEmcClusteringThreshold(0.20) ;
 135   fClu->SetPpsdEnergyThreshold    (0.0000002) ;
 136   fClu->SetPpsdClusteringThreshold(0.0000001) ;
 137   fClu->SetLocalMaxCut(0.03) ;
 138   fClu->SetCalibrationParameters(0., 0.00000001) ;
 139
 140   Int_t ievent;
 141
 142   for ( ievent=0; ievent<Nevents; ievent++)
 143     {
 144
 145       //========== Event Number>
 146       if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. )
 147         cout <<  "AnalyzeResolutions > " << "Event is " << ievent << endl ;
 148
 149       //=========== Connects the various Tree's for evt
 150       gAlice->GetEvent(ievent);
 151
 152       //=========== Gets the Kine TTree
 153       gAlice->TreeK()->GetEvent(0) ;
 154
 155       //=========== Get the Digit Tree
 156       gAlice->TreeD()->GetEvent(0) ;
 157
 158       //========== Creating branches ===================================
 159       AliPHOSRecPoint::RecPointsList ** EmcRecPoints =  fPHOS->EmcRecPoints() ;
 160       gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP", EmcRecPoints ) ;
 161
 162       AliPHOSRecPoint::RecPointsList ** PpsdRecPoints = fPHOS->PpsdRecPoints() ;
 163       gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", PpsdRecPoints ) ;
 164
 165       AliPHOSTrackSegment::TrackSegmentsList **  TrackSegmentsList = fPHOS->TrackSegments() ;
 166       if( (*TrackSegmentsList) )
 167         (*TrackSegmentsList)->Clear() ;
 168       gAlice->TreeR()->SetBranchAddress( "PHOSTS", TrackSegmentsList ) ;
 169
 170       AliPHOSRecParticle::RecParticlesList ** RecParticleList  = fPHOS->RecParticles() ;
 171       if( (*RecParticleList) )
 172         (*RecParticleList)->Clear() ;
 173       gAlice->TreeR()->SetBranchAddress( "PHOSRP", RecParticleList ) ;
 174
 175
 176       //=========== Gets the Reconstraction TTree
 177       gAlice->TreeR()->GetEvent(0) ;
 178
 179       AliPHOSPpsdRecPoint * RecPoint ;
 180       Int_t relid[4] ;
 181       TIter nextRP(*fPHOS->PpsdRecPoints() ) ;
 182       while( ( RecPoint = (AliPHOSPpsdRecPoint *)nextRP() ) )
 183         {
 184           if(!(RecPoint->GetUp()) ) {
 185             AliPHOSDigit *digit ;
 186             Int_t iDigit ;
 187             for(iDigit = 0; iDigit < fPHOS->Digits()->GetEntries(); iDigit++)
 188               {
 189                 digit = (AliPHOSDigit *) fPHOS->Digits()->At(iDigit) ;
 190                 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
 191                 if((relid[2]==1)&&(relid[3]==1)&&(relid[0]==RecPoint->GetPHOSMod())){
 192                   Float_t ConvertorEnergy = fClu->Calibrate(digit->GetAmp()) ;
 193                   fhEnergyCorrelations->Fill(ConvertorEnergy,RecPoint->GetTotalEnergy() );
 194                   break ;
 195                 }
 196               }
 197             break ;
 198           }
 199         }
 200     }
 201   SaveHistograms() ;
 202   fhEnergyCorrelations->Draw("BOX") ;
 203 }
 204
 205 //____________________________________________________________________________
 206 void AliPHOSAnalyze::AnalyzeManyEvents(Int_t Nevents, Int_t module)
 207 {
 208   // analyzes Nevents events in a single PHOS module
 209   // Events should be reconstructed by Reconstruct()
 210
 211   if ( fRootFile == 0 )
 212     cout << "AnalyzeManyEvents > " << "Root File not openned" << endl ;
 213   else
 214     {
 215       //========== Booking Histograms
 216       cout << "AnalyzeManyEvents > " << "Booking Histograms" << endl ;
 217       BookingHistograms();
 218
 219       Int_t ievent;
 220       Int_t relid[4] ;
 221       AliPHOSDigit * digit ;
 222       AliPHOSEmcRecPoint * emc ;
 223       AliPHOSPpsdRecPoint * ppsd ;
 224       //      AliPHOSTrackSegment * tracksegment ;
 225       AliPHOSRecParticle * recparticle;
 226
 227       for ( ievent=0; ievent<Nevents; ievent++)
 228         {
 229           //========== Event Number>
 230           if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. )
 231             cout <<  "AnalyzeManyEvents > " << "Event is " << ievent << endl ;
 232
 233           //=========== Connects the various Tree's for evt
 234           gAlice->GetEvent(ievent);
 235
 236           //=========== Gets the Digit TTree
 237           gAlice->TreeD()->GetEvent(0) ;
 238
 239           //=========== Gets the number of entries in the Digits array
 240           TIter nextdigit(fPHOS->Digits()) ;
 241           while( ( digit = (AliPHOSDigit *)nextdigit() ) )
 242             {
 243               fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
 244               if (fClu->IsInEmc(digit)) fhEmcDigit->Fill(fClu->Calibrate(digit->GetAmp())) ;
 245               else
 246                 {
 247                   if (relid[1]<17) fhVetoDigit->Fill(fClu->Calibrate(digit->GetAmp()));
 248                   if (relid[1]>16) fhConvertorDigit->Fill(fClu->Calibrate(digit->GetAmp()));
 249                 }
 250             }
 251
 252
 253           //=========== Cluster in module
 254           TIter nextEmc(*fPHOS->EmcRecPoints()  ) ;
 255           while((emc = (AliPHOSEmcRecPoint *)nextEmc()))
 256             {
 257               if ( emc->GetPHOSMod() == module )
 258                 {
 259                   fhEmcCluster->Fill(  emc->GetTotalEnergy()  );
 260                   TIter nextPpsd( *fPHOS->PpsdRecPoints()) ;
 261                   while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd()))
 262                     {
 263                       if ( ppsd->GetPHOSMod() == module )
 264                         {
 265                           if (!ppsd->GetUp()) fhConvertorEmc->Fill(emc->GetTotalEnergy(),ppsd->GetTotalEnergy()) ;
 266                         }
 267                     }
 268                 }
 269             }
 270
 271           //=========== Cluster in module PPSD Down
 272           TIter nextPpsd(*fPHOS->PpsdRecPoints() ) ;
 273           while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd()))
 274             {
 275               if ( ppsd->GetPHOSMod() == module )
 276                 {
 277                   if (!ppsd->GetUp()) fhConvertorCluster->Fill(ppsd->GetTotalEnergy()) ;
 278                   if (ppsd->GetUp())  fhVetoCluster     ->Fill(ppsd->GetTotalEnergy()) ;
 279                 }
 280             }
 281
 282           //========== TRackSegments in the event
 283           TIter nextRecParticle(*fPHOS->RecParticles() ) ;
 284           while((recparticle = (AliPHOSRecParticle *)nextRecParticle()))
 285             {
 286               if ( recparticle->GetPHOSTrackSegment()->GetPHOSMod() == module )
 287                 {
 288                   cout << "Particle type is " << recparticle->GetType() << endl ;
 289                   Int_t numberofprimaries = 0 ;
 290                   Int_t * listofprimaries = recparticle->GetPrimaries(numberofprimaries) ;
 291                   cout << "Number of primaries = " << numberofprimaries << endl ;
 292                   Int_t index ;
 293                   for ( index = 0 ; index < numberofprimaries ; index++)
 294                     cout << "    primary # " << index << " =  " << listofprimaries[index] << endl ;
 295                 }
 296             }
 297         }   // endfor
 298       SaveHistograms();
 299     }       // endif
 300 }           // endfunction
 301
 302 //____________________________________________________________________________
 303 void AliPHOSAnalyze::ReconstructCPV(Int_t Nevents )
 304 {
 305
 306   // Perform reconstruction of EMC and CPV (GPS2 or IHEP) for <Nevents> events
 307   // Yuri Kharlov. 19 October 2000
 308
 309   Int_t ievent ;
 310   for ( ievent=0; ievent<Nevents; ievent++) {
 311     if (ievent==0) {
 312       cout << "Analyze > Starting Reconstructing " << endl ;
 313       //========== Create the Clusterizer
 314       fClu = new AliPHOSClusterizerv1() ;
 315       fClu->SetEmcEnergyThreshold(0.05) ;
 316       fClu->SetEmcClusteringThreshold(0.20) ;
 317       fClu->SetLocalMaxCut(0.03) ;
 318       if      (strcmp(fGeom->GetName(),"GPS2") == 0) {
 319         fClu->SetPpsdEnergyThreshold    (0.0000002) ;
 320         fClu->SetPpsdClusteringThreshold(0.0000001) ;
 321       }
 322       else if (strcmp(fGeom->GetName(),"IHEP") == 0) {
 323         fClu->SetLocalMaxCutCPV(0.03) ;
 324         fClu->SetLogWeightCutCPV(4.0) ;
 325         fClu->SetPpsdEnergyThreshold    (0.09) ;
 326       }
 327       fClu->SetCalibrationParameters(0., 0.00000001) ;
 328
 329       //========== Creates the track segment maker
 330       fTrs = new AliPHOSTrackSegmentMakerv1()  ;
 331
 332       //========== Creates the particle identifier for GPS2 only
 333       if      (strcmp(fGeom->GetName(),"GPS2") == 0) {
 334         fPID = new AliPHOSPIDv1() ;
 335         fPID->SetShowerProfileCuts(0.3, 1.8, 0.3, 1.8 ) ;
 336       }
 337
 338       //========== Creates the Reconstructioner
 339       fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;
 340       if (fDebugLevel != 0) fRec -> SetDebugReconstruction(kTRUE);
 341     }
 342
 343     if (fDebugLevel != 0 ||
 344         (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
 345       cout <<  "======= Analyze ======> Event " << ievent+1 << endl ;
 346
 347     //=========== Connects the various Tree's for evt
 348     gAlice->GetEvent(ievent);
 349
 350     //=========== Gets the Digit TTree
 351     gAlice->TreeD()->GetEvent(0) ;
 352
 353     //=========== Do the reconstruction
 354     fPHOS->Reconstruction(fRec);
 355   }
 356
 357   if(fClu)      {delete fClu      ; fClu     =0 ;}
 358   if(fPID)      {delete fPID      ; fPID     =0 ;}
 359   if(fRec)      {delete fRec      ; fRec     =0 ;}
 360   if(fTrs)      {delete fTrs      ; fTrs     =0 ;}
 361
 362 }
 363 //-------------------------------------------------------------------------------------
 364
 365 void AliPHOSAnalyze::Reconstruct(Int_t Nevents )
 366 {
 367   Int_t ievent ;
 368   for ( ievent=0; ievent<Nevents; ievent++)
 369     {
 370       if (ievent==0)
 371         {
 372           cout << "Analyze > Starting Reconstructing " << endl ;
 373           //========== Create the Clusterizer
 374           fClu = new AliPHOSClusterizerv1() ;
 375           fClu->SetEmcEnergyThreshold(0.05) ;
 376           fClu->SetEmcClusteringThreshold(0.20) ;
 377           fClu->SetPpsdEnergyThreshold    (0.0000002) ;
 378           fClu->SetPpsdClusteringThreshold(0.0000001) ;
 379           fClu->SetLocalMaxCut(0.03) ;
 380           fClu->SetCalibrationParameters(0., 0.00000001) ;
 381
 382           //========== Creates the track segment maker
 383           fTrs = new AliPHOSTrackSegmentMakerv1()  ;
 384           //      fTrs->UnsetUnfoldFlag() ;
 385
 386           //========== Creates the particle identifier
 387           fPID = new AliPHOSPIDv1() ;
 388           fPID->SetShowerProfileCuts(0.3, 1.8, 0.3, 1.8 ) ;
 389
 390           //========== Creates the Reconstructioner
 391           fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;
 392 //        fRec -> SetDebugReconstruction(kTRUE);
 393
 394         }
 395
 396       //========== Event Number>
 397     if ((ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
 398       cout <<  "======= Analyze ======> Event " << ievent+1 << endl ;
 399
 400       //=========== Connects the various Tree's for evt
 401       gAlice->GetEvent(ievent);
 402
 403       //=========== Gets the Digit TTree
 404       gAlice->TreeD()->GetEvent(0) ;
 405
 406       //=========== Do the reconstruction
 407       fPHOS->Reconstruction(fRec);
 408     }
 409
 410   if(fClu)      {delete fClu      ; fClu     =0 ;}
 411   if(fPID)      {delete fPID      ; fPID     =0 ;}
 412   if(fRec)      {delete fRec      ; fRec     =0 ;}
 413   if(fTrs)      {delete fTrs      ; fTrs     =0 ;}
 414
 415 }
 416 //-------------------------------------------------------------------------------------
 417
 418 //   TClonesArray AllDigitArray = TClonesArray("AliPHOSDigit",1000) ;
 419 //   TClonesArray * PhotonsList ;
 420 //   TClonesArray * FalsDigitsList ;
 421 //   TClonesArray AllPrimary = TClonesArray("TParticle",5000) ;
 422 //   TFile * file2 = new TFile("ph100.root") ; // file with added photons
 423 //   gAlice = (AliRun*) file2->Get("gAlice") ;
 424 //   Int_t ievent;
 425 //   Int_t NDigits[Nevents+1] ;
 426 //   NDigits[0]=0 ;
 427 //   Int_t NAllDigits = 0;
 428 //   Int_t NprimPerEvent = 20 ;
 429 //   for (ievent=0; ievent <Nevents; ievent++)
 430 //     {
 431 //       PhotonsList  = gAlice->Particles();  //Primary
 432 //       FalsDigitsList  = ((AliPHOSv1 *)gAlice->GetDetector("PHOS"))->Digits();  //Digits
 433 //       gAlice->GetEvent(ievent) ;
 434 //       gAlice->TreeD()->GetEvent(0) ;
 435 //       gAlice->TreeK()->GetEvent(0) ;
 436 //       //Copy Primary
 437 //       Int_t Nprim ;
 438 //       for(Nprim = 0 ;Nprim < NprimPerEvent ; Nprim++)
 439 //      new (AllPrimary[Nprim+ievent*NprimPerEvent])  TParticle(*((TParticle *) PhotonsList->At(Nprim))) ;
 440
 441 //       //Copy Digits
 442 //       TIter nextDigit(FalsDigitsList) ;
 443 //       AliPHOSDigit * FalseDigit ;
 444 //       NDigits[ievent+1] = NDigits[ievent]+ FalsDigitsList->GetEntriesFast() ;
 445 //       while( (FalseDigit = (AliPHOSDigit *) nextDigit()))
 446 //      {
 447 //        new (AllDigitArray[NAllDigits])  AliPHOSDigit(FalseDigit->GetPrimary(1),FalseDigit->GetId(),FalseDigit->GetAmp()) ;
 448 //        NAllDigits++ ;
 449 //      }
 450 //     }
 451 //   file2->Close() ;
 452
 453
 454
 455 //       //Add primary particles
 456 //       cout << "# of Primaries before add " << PrimaryList->GetEntriesFast() << endl;
 457 //      Int_t NTruePrimary = 0 ;  //PrimaryList->GetEntriesFast() ;
 458 //       Int_t Nprim ;
 459 //       for(Nprim = 0; Nprim < NprimPerEvent; Nprim++)
 460 //      new ((*PrimaryList)[NTruePrimary+Nprim])  TParticle(*((TParticle *) AllPrimary.At(Nprim+ievent*NprimPerEvent))) ;
 461
 462 //       cout << "# of Primaries after add " << PrimaryList->GetEntriesFast() <<endl;
 463
 464 //       cout << "Digits before add " << DigitsList->GetEntries() << endl ;
 465 //       cout << "Digits to add " <<  NDigits[ievent+1]-  NDigits[ievent]<< endl ;
 466
 467       //=========== Add fals digits ==============================
 468 //       TIter nextDigit(DigitsList) ;
 469 //       AliPHOSDigit * FalseDigit ;
 470 //       AliPHOSDigit * RealDigit ;
 471 //       Int_t NTrueDigits = DigitsList->GetEntriesFast() ;
 472 //       Int_t Ndigit ;
 473 //       for(Ndigit=NDigits[ievent];Ndigit<NDigits[ievent+1];Ndigit++)
 474 //      {
 475 //        FalseDigit = (AliPHOSDigit*) AllDigitArray.At(Ndigit) ;
 476 //        Bool_t Add = kTRUE ;
 477 //        AliPHOSDigit tmpDigit=AliPHOSDigit(FalseDigit->GetPrimary(1)+NTruePrimary,FalseDigit->GetId(),FalseDigit->GetAmp()) ;
 478
 479 //        while( (RealDigit = (AliPHOSDigit *) nextDigit()) && Add)
 480 //          {
 481 //            if((*RealDigit) == (tmpDigit))
 482 //              {
 483 //                *RealDigit=*RealDigit+tmpDigit ;
 484 //                Add = kFALSE ;
 485 //              }
 486 //          }
 487 //        if(Add)
 488 //          {
 489 //            new ((*DigitsList)[NTrueDigits])  AliPHOSDigit(FalseDigit->GetPrimary(1)+NTruePrimary,FalseDigit->GetId(),FalseDigit->GetAmp()) ;
 490 //            ((AliPHOSDigit *)DigitsList->At(NTrueDigits))->SetIndexInList(NTrueDigits) ;
 491 //            NTrueDigits++ ;
 492 //          }
 493 //      }
 494 //       cout << "Digits after add " << DigitsList->GetEntries() << endl ;
 495
 496
 497 //____________________________________________________________________________
 498 void AliPHOSAnalyze::ReadAndPrintCPV(Int_t Nevents)
 499 {
 500   //
 501   // Read and print generated and reconstructed hits in CPV
 502   // Author: Yuri Kharlov
 503   // 12 October 2000
 504   //
 505
 506   cout << "Start CPV Analysis"<< endl ;
 507   for ( Int_t ievent=0; ievent<Nevents; ievent++) {
 508
 509     //========== Event Number>
 510     cout << endl <<  "==== ReadAndPrintCPV ====> Event is " << ievent+1 << endl ;
 511
 512     //=========== Connects the various Tree's for evt
 513     gAlice->GetEvent(ievent);
 514
 515     //=========== Get the Hits Tree
 516     gAlice->ResetHits();
 517     gAlice->TreeH()->GetEvent(0);
 518
 519     //========== Creating branches ===================================
 520     AliPHOSRecPoint::RecPointsList ** EmcRecPoints = fPHOS->EmcRecPoints() ;
 521     gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP" , EmcRecPoints  ) ;
 522
 523     AliPHOSRecPoint::RecPointsList ** CpvRecPoints = fPHOS->PpsdRecPoints() ;
 524     gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", CpvRecPoints ) ;
 525
 526     //=========== Gets the Reconstruction TTree
 527     gAlice->TreeR()->GetEvent(0) ;
 528
 529     // Read and print CPV hits
 530
 531     TClonesArray *CPVhits;
 532     for (Int_t iModule=0; iModule < fGeom->GetNModules(); iModule++) {
 533       CPVModule   cpvModule = fPHOS->GetCPVModule(iModule);
 534       CPVhits   = cpvModule.Hits();
 535       Int_t nCPVhits  = CPVhits->GetEntriesFast();
 536       for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
 537         CPVHit        *cpvHit = (CPVHit*)CPVhits->UncheckedAt(ihit);
 538         TLorentzVector p      = cpvHit->GetMomentum();
 539         Float_t        xgen   = cpvHit->GetX();
 540         Float_t        zgen   = cpvHit->GetY();
 541         Int_t          ipart  = cpvHit->GetIpart();
 542         printf("CPV hit in module %d: ",iModule+1);
 543         printf(" p = (%f, %f, %f, %f) GeV,\n",
 544                p.Px(),p.Py(),p.Pz(),p.Energy());
 545         printf("               xy = (%8.4f, %8.4f) cm, ipart = %d\n",
 546                xgen,zgen,ipart);
 547       }
 548     }
 549
 550     // Read and print CPV reconstructed points
 551
 552     TIter nextRP(*fPHOS->PpsdRecPoints() ) ;
 553     AliPHOSPpsdRecPoint *cpvRecPoint ;
 554     while( ( cpvRecPoint = (AliPHOSPpsdRecPoint *)nextRP() ) ) {
 555       TVector3  locpos;
 556       cpvRecPoint->GetLocalPosition(locpos);
 557       Int_t PHOSModule = cpvRecPoint->GetPHOSMod();
 558       printf("CPV recpoint in module %d: (X,Y,Z) = (%f,%f,%f) cm\n",
 559              PHOSModule,locpos.X(),locpos.Y(),locpos.Z());
 560     }
 561   }
 562 }
 563
 564 //____________________________________________________________________________
 565 void AliPHOSAnalyze::AnalyzeCPV(Int_t Nevents)
 566 {
 567   //
 568   // Analyzes CPV characteristics
 569   // Author: Yuri Kharlov
 570   // 9 October 2000
 571   //
 572
 573   // Book histograms
 574
 575   TH1F *hDx   = new TH1F("hDx"  ,"CPV x-resolution@reconstruction",100,-5. , 5.);
 576   TH1F *hDz   = new TH1F("hDz"  ,"CPV z-resolution@reconstruction",100,-5. , 5.);
 577   TH1S *hNrp  = new TH1S("hNrp" ,"CPV rec.point multiplicity",      21,-0.5,20.5);
 578
 579   cout << "Start CPV Analysis"<< endl ;
 580   for ( Int_t ievent=0; ievent<Nevents; ievent++) {
 581
 582     //========== Event Number>
 583     if ( (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
 584       cout << endl <<  "==== AnalyzeCPV ====> Event is " << ievent+1 << endl ;
 585
 586     //=========== Connects the various Tree's for evt
 587     gAlice->GetEvent(ievent);
 588
 589     //=========== Get the Hits Tree
 590     gAlice->ResetHits();
 591     gAlice->TreeH()->GetEvent(0);
 592
 593     //========== Creating branches ===================================
 594     AliPHOSRecPoint::RecPointsList ** EmcRecPoints = fPHOS->EmcRecPoints() ;
 595     gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP" , EmcRecPoints  ) ;
 596
 597     AliPHOSRecPoint::RecPointsList ** CpvRecPoints = fPHOS->PpsdRecPoints() ;
 598     gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", CpvRecPoints ) ;
 599
 600     //=========== Gets the Reconstruction TTree
 601     gAlice->TreeR()->GetEvent(0) ;
 602
 603     TIter nextRP(*fPHOS->PpsdRecPoints() ) ;
 604     AliPHOSEmcRecPoint *cpvRecPoint ;
 605     CPVModule           cpvModule;
 606     TClonesArray        *CPVhits;
 607     while( ( cpvRecPoint = (AliPHOSEmcRecPoint *)nextRP() ) ) {
 608       TVector3  locpos;
 609       cpvRecPoint->GetLocalPosition(locpos);
 610       Int_t PHOSModule = cpvRecPoint->GetPHOSMod();
 611       Int_t rpMult     = cpvRecPoint->GetDigitsMultiplicity();
 612       Float_t xrec  = locpos.X();
 613       Float_t zrec  = locpos.Z();
 614       Float_t dxmin = 1.e+10;
 615       Float_t dzmin = 1.e+10;
 616
 617       cpvModule = fPHOS->GetCPVModule(PHOSModule-1);
 618       CPVhits   = cpvModule.Hits();
 619       Int_t nCPVhits  = CPVhits->GetEntriesFast();
 620       for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
 621         CPVHit        *cpvHit = (CPVHit*)CPVhits->UncheckedAt(ihit);
 622         Float_t        xgen   = cpvHit->GetX();
 623         Float_t        zgen   = cpvHit->GetY();
 624         if ( TMath::Abs(xgen-xrec) < TMath::Abs(dxmin) ) dxmin = xgen-xrec;
 625         if ( TMath::Abs(zgen-zrec) < TMath::Abs(dzmin) ) dzmin = zgen-zrec;
 626       }
 627       cpvModule.Clear();
 628       hDx  ->Fill(dxmin);
 629       hDz  ->Fill(dzmin);
 630       hNrp ->Fill(rpMult);
 631     }
 632
 633   }
 634
 635   // Save histograms
 636
 637   Text_t outputname[80] ;
 638   sprintf(outputname,"%s.analyzed",fRootFile->GetName());
 639   TFile output(outputname,"RECREATE");
 640   output.cd();
 641
 642   hDx  ->Write() ;
 643   hDz  ->Write() ;
 644   hNrp ->Write() ;
 645
 646   // Plot histograms
 647
 648   TCanvas *CPVcanvas = new TCanvas("CPV","CPV analysis",20,20,300,900);
 649   gStyle->SetOptStat(111111);
 650   gStyle->SetOptFit(1);
 651   gStyle->SetOptDate(1);
 652   CPVcanvas->Divide(3,1);
 653
 654   CPVcanvas->cd(1);
 655   gPad->SetFillColor(10);
 656   hNrp->SetFillColor(16);
 657   hNrp->Draw();
 658
 659   CPVcanvas->cd(2);
 660   gPad->SetFillColor(10);
 661   hDx->SetFillColor(16);
 662   hDx->Fit("gaus");
 663   hDx->Draw();
 664
 665   CPVcanvas->cd(3);
 666   gPad->SetFillColor(10);
 667   hDz->SetFillColor(16);
 668   hDz->Fit("gaus");
 669   hDz->Draw();
 670
 671   CPVcanvas->Print("CPV.ps");
 672
 673 }
 674
 675 //____________________________________________________________________________
 676  void AliPHOSAnalyze::AnalyzeResolutions(Int_t Nevents )
 677 {
 678   // analyzes Nevents events and calculate Energy and Position resolution as well as
 679   // probaility of correct indentifiing of the incident particle
 680
 681   //========== Booking Histograms
 682   cout << "AnalyzeResolutions > " << "Booking Histograms" << endl ;
 683   BookResolutionHistograms();
 684
 685   Int_t Counter[9][5] ;
 686   Int_t i1,i2,TotalInd = 0 ;
 687   for(i1 = 0; i1<9; i1++)
 688     for(i2 = 0; i2<5; i2++)
 689       Counter[i1][i2] = 0 ;
 690
 691   Int_t TotalPrimary = 0 ;
 692   Int_t TotalRecPart = 0 ;
 693   Int_t TotalRPwithPrim = 0 ;
 694   Int_t ievent;
 695
 696   cout << "Start Analysing"<< endl ;
 697   for ( ievent=0; ievent<Nevents; ievent++)
 698     {
 699
 700       //========== Event Number>
 701       if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. )
 702         cout <<  "AnalyzeResolutions > " << "Event is " << ievent << endl ;
 703
 704       //=========== Connects the various Tree's for evt
 705       gAlice->GetEvent(ievent);
 706
 707
 708
 709             //=========== Gets the Kine TTree
 710       gAlice->TreeK()->GetEvent(0) ;
 711
 712       //=========== Gets the list of Primari Particles
 713       TClonesArray * PrimaryList  = gAlice->Particles();
 714
 715       TParticle * Primary ;
 716       Int_t iPrimary ;
 717       for ( iPrimary = 0 ; iPrimary < PrimaryList->GetEntries() ; iPrimary++)
 718             {
 719               Primary = (TParticle*)PrimaryList->UncheckedAt(iPrimary) ;
 720               Int_t PrimaryType = Primary->GetPdgCode() ;
 721               if( PrimaryType == 22 )
 722                 fhPrimary->Fill(Primary->Energy()) ;
 723             }
 724
 725       //=========== Get the Digit Tree
 726       gAlice->TreeD()->GetEvent(0) ;
 727
 728       //========== Creating branches ===================================
 729       AliPHOSRecPoint::RecPointsList ** EmcRecPoints =  fPHOS->EmcRecPoints() ;
 730       gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP", EmcRecPoints ) ;
 731
 732       AliPHOSRecPoint::RecPointsList ** PpsdRecPoints = fPHOS->PpsdRecPoints() ;
 733       gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", PpsdRecPoints ) ;
 734
 735       AliPHOSTrackSegment::TrackSegmentsList **  TrackSegmentsList = fPHOS->TrackSegments() ;
 736       if( (*TrackSegmentsList) )
 737         (*TrackSegmentsList)->Clear() ;
 738       gAlice->TreeR()->SetBranchAddress( "PHOSTS", TrackSegmentsList ) ;
 739
 740       AliPHOSRecParticle::RecParticlesList ** RecParticleList  = fPHOS->RecParticles() ;
 741       if( (*RecParticleList) )
 742         (*RecParticleList)->Clear() ;
 743       gAlice->TreeR()->SetBranchAddress( "PHOSRP", RecParticleList ) ;
 744
 745
 746       //=========== Gets the Reconstraction TTree
 747       gAlice->TreeR()->GetEvent(0) ;
 748
 749       cout << ievent << " " << (*EmcRecPoints) << " " <<(*PpsdRecPoints) <<fPHOS->Digits()<< endl ;
 750       cout << "    " << " " << (*EmcRecPoints)->GetEntries() << " " <<(*PpsdRecPoints)->GetEntries() <<fPHOS->Digits()->GetEntries()<< endl ;
 751
 752       AliPHOSRecParticle * RecParticle ;
 753       Int_t iRecParticle ;
 754       for(iRecParticle = 0; iRecParticle < (*RecParticleList)->GetEntries() ;iRecParticle++ )
 755         {
 756           RecParticle = (AliPHOSRecParticle *) (*RecParticleList)->At(iRecParticle) ;
 757
 758           Int_t ModuleNumberRec ;
 759           Double_t RecX, RecZ ;
 760           fGeom->ImpactOnEmc(RecParticle->Theta(), RecParticle->Phi(), ModuleNumberRec, RecX, RecZ) ;
 761
 762           Double_t MinDistance = 10000 ;
 763           Int_t ClosestPrimary = -1 ;
 764
 765           Int_t numberofprimaries ;
 766           Int_t * listofprimaries  = RecParticle->GetPrimaries(numberofprimaries)  ;
 767           Int_t index ;
 768           TParticle * Primary ;
 769           Double_t Distance = MinDistance ;
 770           for ( index = 0 ; index < numberofprimaries ; index++)
 771             {
 772               Primary = (TParticle*)PrimaryList->UncheckedAt(listofprimaries[index]) ;
 773               Int_t ModuleNumber ;
 774               Double_t PrimX, PrimZ ;
 775               fGeom->ImpactOnEmc(Primary->Theta(), Primary->Phi(), ModuleNumber, PrimX, PrimZ) ;
 776               if(ModuleNumberRec == ModuleNumber)
 777                 Distance = TMath::Sqrt((RecX-PrimX)*(RecX-PrimX)+(RecZ-PrimZ)*(RecZ-PrimZ) ) ;
 778               if(MinDistance > Distance)
 779                 {
 780                   MinDistance = Distance ;
 781                   ClosestPrimary = listofprimaries[index] ;
 782                 }
 783             }
 784           TotalRecPart++ ;
 785           if(ClosestPrimary >=0 )
 786             {
 787               fhPhotonAllEnergy->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(), RecParticle->Energy() ) ;
 788               fhPhotonAllPosition->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(),Distance) ;
 789               TotalRPwithPrim++;
 790               Int_t PrimaryType = ((TParticle *)PrimaryList->At(ClosestPrimary))->GetPdgCode() ;
 791               TParticlePDG* PDGparticle = ((TParticle *)PrimaryList->At(ClosestPrimary))->GetPDG();
 792               Double_t charge =  PDGparticle->Charge() ;
 793               Int_t PrimaryCode ;
 794               switch(PrimaryType)
 795                 {
 796                 case 22:
 797                   PrimaryCode = 0;  //Photon
 798                   break;
 799                 case 11 :
 800                   PrimaryCode = 1;  //Electron
 801                   break;
 802                 case -11 :
 803                   PrimaryCode = 1;  //positron
 804                   break;
 805                 case 321 :
 806                   PrimaryCode = 4;  //K+
 807                   break;
 808                 case -321 :
 809                   PrimaryCode = 4;  //K-
 810                   break;
 811                 case 310 :
 812                   PrimaryCode = 4;  //K0s
 813                   break;
 814                 case 130 :
 815                   PrimaryCode = 4;  //K0l
 816                   break;
 817                 default:
 818                   if(charge)
 819                     PrimaryCode = 2; //Charged hadron
 820                   else
 821                     PrimaryCode = 3; //Neutral hadron
 822                   break;
 823                 }
 824
 825               switch(RecParticle->GetType())
 826                 {
 827                 case AliPHOSFastRecParticle::kGAMMA:
 828                   if(PrimaryType == 22){
 829                     fhPhotonEnergy->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(), RecParticle->Energy() ) ;
 830                     fhPhotonPosition->Fill(((TParticle *) PrimaryList->At(ClosestPrimary))->Energy(),Distance) ;
 831                     fhPhotonReg->Fill(RecParticle->Energy() ) ;
 832                     fhPhotonEM->Fill(RecParticle->Energy() ) ;
 833                     fhPhotPhot->Fill(RecParticle->Energy() ) ;
 834                   }
 835                   if(PrimaryType == 2112){ //neutron
 836                     fhNReg->Fill(RecParticle->Energy() ) ;
 837                     fhNEM->Fill(RecParticle->Energy() ) ;
 838                   }
 839
 840                   if(PrimaryType == -2112){ //neutron ~
 841                     fhNBarReg->Fill(RecParticle->Energy() ) ;
 842                     fhNBarEM->Fill(RecParticle->Energy() ) ;
 843
 844                   }
 845                   if(PrimaryCode == 2){
 846                     fhChargedReg->Fill(RecParticle->Energy() ) ;
 847                     fhChargedEM->Fill(RecParticle->Energy() ) ;
 848                   }
 849
 850                   fhAllReg->Fill(RecParticle->Energy() ) ;
 851                   fhAllEM->Fill(RecParticle->Energy() ) ;
 852                   Counter[0][PrimaryCode]++;
 853                   break;
 854                 case  AliPHOSFastRecParticle::kELECTRON:
 855                   if(PrimaryType == 11 || PrimaryType == -11){
 856                     fhElectronEnergy->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(), RecParticle->Energy() ) ;
 857                     fhElectronPosition->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(),Distance ) ;
 858                   }
 859                   if(PrimaryType == 22)
 860                     fhPhotElec->Fill(RecParticle->Energy() ) ;
 861
 862                   Counter[1][PrimaryCode]++;
 863                   break;
 864                 case  AliPHOSFastRecParticle::kNEUTRALHA:
 865                   if(PrimaryType == 22)
 866                     fhPhotNeuH->Fill(RecParticle->Energy() ) ;
 867
 868                   fhNeutralHadronEnergy->Fill( ((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(), RecParticle->Energy()) ;
 869                   fhNeutralHadronPosition->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy() ,Distance  ) ;
 870                   Counter[2][PrimaryCode]++;
 871                   break ;
 872                 case  AliPHOSFastRecParticle::kNEUTRALEM:
 873                   if(PrimaryType == 22 || PrimaryType == 11 || PrimaryType == -11){
 874                     fhNeutralEMEnergy->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(),RecParticle->Energy() ) ;
 875                     fhNeutralEMPosition->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(),Distance ) ;
 876                   }
 877
 878                   if(PrimaryType == 22){ //photon
 879                     fhPhotNuEM->Fill(RecParticle->Energy() ) ;
 880                     fhPhotonEM->Fill(RecParticle->Energy() ) ;
 881                   }
 882                   if(PrimaryType == 2112) //neutron
 883                     fhNEM->Fill(RecParticle->Energy() ) ;
 884
 885                   if(PrimaryType == -2112) //neutron ~
 886                     fhNBarEM->Fill(RecParticle->Energy() ) ;
 887
 888                   if(PrimaryCode == 2)
 889                     fhChargedEM->Fill(RecParticle->Energy() ) ;
 890
 891                   fhAllEM->Fill(RecParticle->Energy() ) ;
 892
 893                   Counter[3][PrimaryCode]++;
 894                   break ;
 895                 case  AliPHOSFastRecParticle::kCHARGEDHA:
 896                   if(PrimaryType == 22) //photon
 897                     fhPhotChHa->Fill(RecParticle->Energy() ) ;
 898
 899                   fhChargedHadronEnergy->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(),RecParticle->Energy() ) ;
 900                   fhChargedHadronPosition->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(),Distance ) ;
 901                   Counter[4][PrimaryCode]++ ;
 902                   break ;
 903                 case  AliPHOSFastRecParticle::kGAMMAHA:
 904                   if(PrimaryType == 22) //photon
 905                     fhPhotGaHa->Fill(RecParticle->Energy() ) ;
 906                   fhPhotonHadronEnergy->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(), RecParticle->Energy()) ;
 907                   fhPhotonHadronPosition->Fill(((TParticle *)PrimaryList->At(ClosestPrimary))->Energy(),Distance ) ;
 908                   Counter[5][PrimaryCode]++ ;
 909                   break ;
 910                 case  AliPHOSFastRecParticle::kABSURDEM:
 911                   Counter[6][PrimaryCode]++ ;
 912                   break;
 913                 case  AliPHOSFastRecParticle::kABSURDHA:
 914                   Counter[7][PrimaryCode]++ ;
 915                   break;
 916                 default:
 917                   Counter[8][PrimaryCode]++ ;
 918                   break;
 919                 }
 920             }
 921         }
 922     }   // endfor
 923   SaveHistograms();
 924   cout << "Resolutions: Analyzed " << Nevents << " event(s)" << endl ;
 925   cout << "Resolutions: Total primary       " << TotalPrimary << endl ;
 926   cout << "Resoluitons: Total reconstracted " << TotalRecPart << endl ;
 927   cout << "TotalReconstructed with Primarie " << TotalRPwithPrim << endl ;
 928   cout << "                        Primary:   Photon   Electron   Ch. Hadr.  Neutr. Hadr  Kaons" << endl ;
 929   cout << "             Detected as photon       " << Counter[0][0] << "          " << Counter[0][1] << "          " << Counter[0][2] << "          " <<Counter[0][3] << "          " << Counter[0][4] << endl ;
 930   cout << "           Detected as electron       " << Counter[1][0] << "          " << Counter[1][1] << "          " << Counter[1][2] << "          " <<Counter[1][3] << "          " << Counter[1][4] << endl ;
 931   cout << "     Detected as neutral hadron       " << Counter[2][0] << "          " << Counter[2][1] << "          " << Counter[2][2] << "          " <<Counter[2][3] << "          " << Counter[2][4] << endl ;
 932   cout << "         Detected as neutral EM       " << Counter[3][0] << "          " << Counter[3][1] << "          " << Counter[3][2] << "          " <<Counter[3][3] << "          " << Counter[3][4] << endl ;
 933   cout << "     Detected as charged hadron       " << Counter[4][0] << "          " << Counter[4][1] << "          " << Counter[4][2] << "          " <<Counter[4][3] << "          " << Counter[4][4] << endl ;
 934   cout << "       Detected as gamma-hadron       " << Counter[5][0] << "          " << Counter[5][1] << "          " << Counter[5][2] << "          " <<Counter[5][3] << "          " << Counter[5][4] << endl ;
 935   cout << "          Detected as Absurd EM       " << Counter[6][0] << "          " << Counter[6][1] << "          " << Counter[6][2] << "          " <<Counter[6][3] << "          " << Counter[6][4] << endl ;
 936   cout << "      Detected as absurd hadron       " << Counter[7][0] << "          " << Counter[7][1] << "          " << Counter[7][2] << "          " <<Counter[7][3] << "          " << Counter[7][4] << endl ;
 937   cout << "          Detected as undefined       " << Counter[8][0] << "          " << Counter[8][1] << "          " << Counter[8][2] << "          " <<Counter[8][3] << "          " << Counter[8][4] << endl ;
 938
 939       for(i1 = 0; i1<9; i1++)
 940         for(i2 = 0; i2<5; i2++)
 941           TotalInd+=Counter[i1][i2] ;
 942       cout << "Indentified particles            " << TotalInd << endl ;
 943
 944 }           // endfunction
 945
 946
 947 //____________________________________________________________________________
 948 void  AliPHOSAnalyze::BookingHistograms()
 949 {
 950   // Books the histograms where the results of the analysis are stored (to be changed)
 951
 952   delete fhEmcDigit  ;
 953   delete fhVetoDigit  ;
 954   delete fhConvertorDigit   ;
 955   delete  fhEmcCluster   ;
 956   delete fhVetoCluster   ;
 957   delete fhConvertorCluster  ;
 958   delete fhConvertorEmc  ;
 959
 960   fhEmcDigit                = new TH1F("hEmcDigit",      "hEmcDigit",         1000,  0. ,  25.);
 961   fhVetoDigit               = new TH1F("hVetoDigit",     "hVetoDigit",         500,  0. ,  3.e-5);
 962   fhConvertorDigit          = new TH1F("hConvertorDigit","hConvertorDigit",    500,  0. ,  3.e-5);
 963   fhEmcCluster              = new TH1F("hEmcCluster",    "hEmcCluster",       1000,  0. ,  30.);
 964   fhVetoCluster             = new TH1F("hVetoCluster",   "hVetoCluster",       500,  0. ,  3.e-5);
 965   fhConvertorCluster        = new TH1F("hConvertorCluster","hConvertorCluster",500,  0. ,  3.e-5);
 966   fhConvertorEmc            = new TH2F("hConvertorEmc",  "hConvertorEmc",      200,  1. ,  3., 200, 0., 3.e-5);
 967
 968 }
 969 //____________________________________________________________________________
 970 void  AliPHOSAnalyze::BookResolutionHistograms()
 971 {
 972   // Books the histograms where the results of the Resolution analysis are stored
 973
 974   if(fhPhotonEnergy)
 975     delete fhPhotonEnergy ;
 976   if(fhPhotonAllEnergy)
 977     delete fhPhotonAllEnergy ;
 978   if(fhElectronEnergy)
 979     delete fhElectronEnergy ;
 980   if(fhElectronAllEnergy)
 981     delete fhElectronAllEnergy ;
 982   if(fhNeutralHadronEnergy)
 983     delete fhNeutralHadronEnergy ;
 984   if(fhNeutralEMEnergy)
 985     delete fhNeutralEMEnergy ;
 986   if(fhNeutralEMAllEnergy)
 987     delete fhNeutralEMAllEnergy ;
 988   if(fhChargedHadronEnergy)
 989     delete fhChargedHadronEnergy ;
 990   if(fhPhotonHadronEnergy)
 991     delete fhPhotonHadronEnergy ;
 992   if(fhPhotonPosition)
 993     delete fhPhotonPosition ;
 994   if(fhPhotonAllPosition)
 995     delete fhPhotonAllPosition ;
 996   if(fhElectronPosition)
 997     delete fhElectronPosition ;
 998   if(fhElectronAllPosition)
 999     delete fhElectronAllPosition ;
1000   if(fhNeutralHadronPosition)
1001     delete fhNeutralHadronPosition ;
1002   if(fhNeutralEMPosition)
1003     delete fhNeutralEMPosition ;
1004   if(fhNeutralEMAllPosition)
1005     delete fhNeutralEMAllPosition ;
1006   if(fhChargedHadronPosition)
1007     delete fhChargedHadronPosition ;
1008   if(fhPhotonHadronPosition)
1009     delete fhPhotonHadronPosition ;
1010
1011   fhPhotonEnergy            = new TH2F("hPhotonEnergy",  "hPhotonEnergy",              100, 0., 5., 100, 0., 5.);
1012   fhPhotonAllEnergy         = new TH2F("hPhotonAllEnergy",  "hPhotonAllEnergy",        100, 0., 5., 100, 0., 5.);
1013   fhElectronEnergy          = new TH2F("hElectronEnergy","hElectronEnergy",            100, 0., 5., 100, 0., 5.);
1014   fhElectronAllEnergy       = new TH2F("hElectronAllEnergy","hElectronAllEnergy",      100, 0., 5., 100, 0., 5.);
1015   fhNeutralHadronEnergy     = new TH2F("hNeutralHadronEnergy", "hNeutralHadronEnergy", 100, 0., 5., 100, 0., 5.);
1016   fhNeutralEMEnergy         = new TH2F("hNeutralEMEnergy", "hNeutralEMEnergy",         100, 0., 5., 100, 0., 5.);
1017   fhNeutralEMAllEnergy      = new TH2F("hNeutralEMAllEnergy", "hNeutralEMAllEnergy",   100, 0., 5., 100, 0., 5.);
1018   fhChargedHadronEnergy     = new TH2F("hChargedHadronEnergy", "hChargedHadronEnergy", 100, 0., 5., 100, 0., 5.);
1019   fhPhotonHadronEnergy      = new TH2F("hPhotonHadronEnergy","hPhotonHadronEnergy",    100, 0., 5., 100, 0., 5.);
1020   fhPhotonPosition          = new TH2F("hPhotonPosition","hPhotonPosition",                20, 0., 5., 100, 0., 5.);
1021   fhPhotonAllPosition       = new TH2F("hPhotonAllPosition","hPhotonAllPosition",          20, 0., 5., 100, 0., 5.);
1022   fhElectronPosition        = new TH2F("hElectronPosition","hElectronPosition",            20, 0., 5., 100, 0., 5.);
1023   fhElectronAllPosition     = new TH2F("hElectronAllPosition","hElectronAllPosition",      20, 0., 5., 100, 0., 5.);
1024   fhNeutralHadronPosition   = new TH2F("hNeutralHadronPosition","hNeutralHadronPosition",  20, 0., 5., 100, 0., 5.);
1025   fhNeutralEMPosition       = new TH2F("hNeutralEMPosition","hNeutralEMPosition",          20, 0., 5., 100, 0., 5.);
1026   fhNeutralEMAllPosition    = new TH2F("hNeutralEMAllPosition","hNeutralEMAllPosition",    20, 0., 5., 100, 0., 5.);
1027   fhChargedHadronPosition   = new TH2F("hChargedHadronPosition","hChargedHadronPosition",  20, 0., 5., 100, 0., 5.);
1028   fhPhotonHadronPosition    = new TH2F("hPhotonHadronPosition","hPhotonHadronPosition",    20, 0., 5., 100, 0., 5.);
1029
1030   if(fhPhotonReg)
1031     delete fhPhotonReg ;
1032   if(fhAllReg)
1033     delete fhAllReg ;
1034   if(fhNReg)
1035     delete fhNReg ;
1036   if(fhNReg)
1037     delete fhNReg ;
1038   if(fhNReg)
1039     delete fhNReg ;
1040
1041   fhPhotonReg = new TH1F("hPhotonReg","hPhotonReg", 20, 0., 5.);
1042   fhAllReg    = new TH1F("hAllReg", "hAllReg",  20, 0., 5.);
1043   fhNReg      = new TH1F("hNReg", "hNReg",  20, 0., 5.);
1044   fhNBarReg   = new TH1F("hNBarReg", "hNBarReg",  20, 0., 5.);
1045   fhChargedReg= new TH1F("hChargedReg", "hChargedReg",  20, 0., 5.);
1046
1047   if(fhPhotonEM)
1048     delete fhPhotonEM ;
1049   if(fhAllEM)
1050     delete fhAllEM ;
1051   if(fhNEM)
1052     delete fhNEM ;
1053   if(fhNBarEM)
1054     delete fhNBarEM ;
1055   if(fhChargedEM)
1056     delete fhChargedEM ;
1057
1058   fhPhotonEM = new TH1F("hPhotonEM","hPhotonEM", 20, 0., 5.);
1059   fhAllEM    = new TH1F("hAllEM", "hAllEM",  20, 0., 5.);
1060   fhNEM      = new TH1F("hNEM", "hNEM",  20, 0., 5.);
1061   fhNBarEM   = new TH1F("hNBarEM", "hNBarEM",  20, 0., 5.);
1062   fhChargedEM= new TH1F("hChargedEM", "hChargedEM",  20, 0., 5.);
1063
1064   if(fhPrimary)
1065     delete fhPrimary ;
1066   fhPrimary= new TH1F("hPrimary", "hPrimary",  20, 0., 5.);
1067
1068   if(fhPhotPhot)
1069     delete fhPhotPhot ;
1070   if(fhPhotElec)
1071     delete fhPhotElec ;
1072   if(fhPhotNeuH)
1073     delete fhPhotNeuH ;
1074   if(fhPhotNuEM)
1075     delete fhPhotNuEM ;
1076   if(fhPhotChHa)
1077     delete fhPhotChHa ;
1078   if(fhPhotGaHa)
1079     delete fhPhotGaHa ;
1080
1081   fhPhotPhot = new TH1F("hPhotPhot","hPhotPhot", 20, 0., 5.);   //Photon registered as photon
1082   fhPhotElec = new TH1F("hPhotElec","hPhotElec", 20, 0., 5.);   //Photon registered as Electron
1083   fhPhotNeuH = new TH1F("hPhotNeuH","hPhotNeuH", 20, 0., 5.);   //Photon registered as Neutral Hadron
1084   fhPhotNuEM = new TH1F("hPhotNuEM","hPhotNuEM", 20, 0., 5.);   //Photon registered as Neutral EM
1085   fhPhotChHa = new TH1F("hPhotChHa","hPhotChHa", 20, 0., 5.);   //Photon registered as Charged Hadron
1086   fhPhotGaHa = new TH1F("hPhotGaHa","hPhotGaHa", 20, 0., 5.);   //Photon registered as Gamma-Hadron
1087
1088
1089 }
1090 //____________________________________________________________________________
1091 Bool_t AliPHOSAnalyze::Init(Int_t evt)
1092 {
1093   // Do a few initializations: open the root file
1094   //                           get the AliRun object
1095   //                           defines the clusterizer, tracksegment maker and particle identifier
1096   //                           sets the associated parameters
1097
1098   Bool_t ok = kTRUE ;
1099
1100    //========== Open galice root file
1101
1102   if ( fRootFile == 0 ) {
1103     Text_t * name  = new Text_t[80] ;
1104     cout << "AnalyzeOneEvent > Enter file root file name : " ;
1105     cin >> name ;
1106     Bool_t ok = OpenRootFile(name) ;
1107     if ( !ok )
1108       cout << " AliPHOSAnalyze > Error opening " << name << endl ;
1109     else {
1110       //========== Get AliRun object from file
1111
1112       gAlice = (AliRun*) fRootFile->Get("gAlice") ;
1113
1114       //=========== Get the PHOS object and associated geometry from the file
1115
1116       fPHOS  = (AliPHOSv1 *)gAlice->GetDetector("PHOS") ;
1117       fGeom = fPHOS->GetGeometry();
1118       //      fGeom  = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() );
1119
1120     } // else !ok
1121   } // if fRootFile
1122
1123   if ( ok ) {
1124
1125     //========== Create the Clusterizer
1126
1127     fClu =  new AliPHOSClusterizerv1() ;
1128     fClu->SetEmcEnergyThreshold(0.030) ;
1129     fClu->SetEmcClusteringThreshold(0.20) ;
1130     fClu->SetPpsdEnergyThreshold    (0.0000002) ;
1131     fClu->SetPpsdClusteringThreshold(0.0000001) ;
1132     fClu->SetLocalMaxCut(0.03) ;
1133     fClu->SetCalibrationParameters(0., 0.00000001) ;
1134     cout <<  "AnalyzeOneEvent > using clusterizer " << fClu->GetName() << endl ;
1135     fClu->PrintParameters() ;
1136
1137     //========== Creates the track segment maker
1138
1139     fTrs = new AliPHOSTrackSegmentMakerv1() ;
1140     cout <<  "AnalyzeOneEvent > using tack segment maker " << fTrs->GetName() << endl ;
1141     //   fTrs->UnsetUnfoldFlag() ;
1142
1143     //========== Creates the particle identifier
1144
1145     fPID = new AliPHOSPIDv1() ;
1146     cout <<  "AnalyzeOneEvent > using particle identifier " << fPID->GetName() << endl ;
1147     //fPID->SetShowerProfileCuts(Float_t l1m, Float_t l1M, Float_t l2m, Float_t l2M) ;
1148     fPID->SetShowerProfileCuts(0.7, 2.0 , 0.6 , 1.5) ;
1149
1150     //========== Creates the Reconstructioner
1151
1152     fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;
1153 //      fRec -> SetDebugReconstruction(kFALSE);
1154     fRec -> SetDebugReconstruction(kTRUE);
1155
1156     //=========== Connect the various Tree's for evt
1157
1158     if ( evt == -999 ) {
1159       cout <<  "AnalyzeOneEvent > Enter event number : " ;
1160       cin >> evt ;
1161       cout << evt << endl ;
1162     }
1163     fEvt = evt ;
1164     gAlice->GetEvent(evt);
1165
1166     //=========== Get the Digit TTree
1167
1168     gAlice->TreeD()->GetEvent(0) ;
1169
1170   } // ok
1171
1172   return ok ;
1173 }
1174
1175
1176 //____________________________________________________________________________
1177 void AliPHOSAnalyze::DisplayKineEvent(Int_t evt)
1178 {
1179   // Display particles from the Kine Tree in global Alice (theta, phi) coordinates.
1180   // One PHOS module at the time.
1181   // The particle type can be selected.
1182
1183   if (evt == -999)
1184     evt = fEvt ;
1185
1186   Int_t module ;
1187   cout <<  "DisplayKineEvent > which module (1-5,  -1: all) ? " ;
1188   cin >> module ; cout << module << endl ;
1189
1190   Int_t testparticle ;
1191   cout << " 22      : PHOTON " << endl
1192        << " (-)11   : (POSITRON)ELECTRON " << endl
1193        << " (-)2112 : (ANTI)NEUTRON " << endl
1194        << " -999    : Everything else " << endl ;
1195   cout  <<  "DisplayKineEvent > enter PDG particle code to display " ;
1196   cin >> testparticle ; cout << testparticle << endl ;
1197
1198   Text_t histoname[80] ;
1199   sprintf(histoname,"Event %d: Incident particles in module %d", evt, module) ;
1200
1201   Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
1202   fGeom->EmcModuleCoverage(module, tm, tM, pm, pM, AliPHOSGeometry::Degre() ) ;
1203
1204   Double_t theta, phi ;
1205   fGeom->EmcXtalCoverage(theta, phi, AliPHOSGeometry::Degre() ) ;
1206
1207   Int_t tdim = (Int_t)( (tM - tm) / theta ) ;
1208   Int_t pdim = (Int_t)( (pM - pm) / phi ) ;
1209
1210   tm -= theta ;
1211   tM += theta ;
1212   pm -= phi ;
1213   pM += phi ;
1214
1215   TH2F * histoparticle = new TH2F("histoparticle",  histoname,
1216                                           pdim, pm, pM, tdim, tm, tM) ;
1217   histoparticle->SetStats(kFALSE) ;
1218
1219   // Get pointers to Alice Particle TClonesArray
1220
1221   TParticle * particle;
1222   TClonesArray * particlearray  = gAlice->Particles();
1223
1224   Text_t canvasname[80];
1225   sprintf(canvasname,"Particles incident in PHOS/EMC module # %d",module) ;
1226   TCanvas * kinecanvas = new TCanvas("kinecanvas", canvasname, 650, 500) ;
1227
1228   // get the KINE Tree
1229
1230   TTree * kine =  gAlice->TreeK() ;
1231   Stat_t nParticles =  kine->GetEntries() ;
1232   cout << "DisplayKineEvent > events in kine " << nParticles << endl ;
1233
1234   // loop over particles
1235
1236   Double_t kRADDEG = 180. / TMath::Pi() ;
1237   Int_t index1 ;
1238   Int_t nparticlein = 0 ;
1239   for (index1 = 0 ; index1 < nParticles ; index1++){
1240     Int_t nparticle = particlearray->GetEntriesFast() ;
1241     Int_t index2 ;
1242     for( index2 = 0 ; index2 < nparticle ; index2++) {
1243       particle            = (TParticle*)particlearray->UncheckedAt(index2) ;
1244       Int_t  particletype = particle->GetPdgCode() ;
1245       if (testparticle == -999 || testparticle == particletype) {
1246         Double_t phi        = particle->Phi() ;
1247         Double_t theta      = particle->Theta() ;
1248         Int_t mod ;
1249         Double_t x, z ;
1250         fGeom->ImpactOnEmc(theta, phi, mod, z, x) ;
1251         if ( mod == module ) {
1252           nparticlein++ ;
1253           if (particle->Energy() >  fClu->GetEmcClusteringThreshold()  )
1254             histoparticle->Fill(phi*kRADDEG, theta*kRADDEG, particle->Energy() ) ;
1255         }
1256       }
1257     }
1258   }
1259   kinecanvas->Draw() ;
1260   histoparticle->Draw("color") ;
1261   TPaveText *  pavetext = new TPaveText(294, 100, 300, 101);
1262   Text_t text[40] ;
1263   sprintf(text, "Particles: %d ", nparticlein) ;
1264   pavetext->AddText(text) ;
1265   pavetext->Draw() ;
1266   kinecanvas->Update();
1267
1268 }
1269 //____________________________________________________________________________
1270 void AliPHOSAnalyze::DisplayRecParticles()
1271 {
1272   // Display reconstructed particles in global Alice(theta, phi) coordinates.
1273   // One PHOS module at the time.
1274   // Click on symbols indicate the reconstructed particle type.
1275
1276   if (fEvt == -999) {
1277     cout << "DisplayRecParticles > Analyze an event first ... (y/n) " ;
1278     Text_t answer[1] ;
1279     cin >> answer ; cout << answer ;
1280 //     if ( answer == "y" )
1281 //       AnalyzeOneEvent() ;
1282   }
1283     if (fEvt != -999) {
1284
1285       Int_t module ;
1286       cout <<  "DisplayRecParticles > which module (1-5,  -1: all) ? " ;
1287       cin >> module ; cout << module << endl ;
1288       Text_t histoname[80] ;
1289       sprintf(histoname,"Event %d: Reconstructed particles in module %d", fEvt, module) ;
1290       Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
1291       fGeom->EmcModuleCoverage(module, tm, tM, pm, pM, AliPHOSGeometry::Degre() ) ;
1292       Double_t theta, phi ;
1293       fGeom->EmcXtalCoverage(theta, phi, AliPHOSGeometry::Degre() ) ;
1294       Int_t tdim = (Int_t)( (tM - tm) / theta ) ;
1295       Int_t pdim = (Int_t)( (pM - pm) / phi ) ;
1296       tm -= theta ;
1297       tM += theta ;
1298       pm -= phi ;
1299       TH2F * histoRparticle = new TH2F("histoRparticle",  histoname,
1300                                        pdim, pm, pM, tdim, tm, tM) ;
1301       histoRparticle->SetStats(kFALSE) ;
1302       Text_t canvasname[80] ;
1303       sprintf(canvasname, "Reconstructed particles in PHOSmodule # %d", module) ;
1304       TCanvas * rparticlecanvas = new TCanvas("RparticleCanvas", canvasname, 650, 500) ;
1305       AliPHOSRecParticle::RecParticlesList * rpl = *fPHOS->RecParticles() ;
1306       Int_t nRecParticles = rpl->GetEntries() ;
1307       Int_t nRecParticlesInModule = 0 ;
1308       TIter nextRecPart(rpl) ;
1309       AliPHOSRecParticle * rp ;
1310       cout << "DisplayRecParticles > " << nRecParticles << " reconstructed particles " << endl ;
1311       Double_t kRADDEG = 180. / TMath::Pi() ;
1312       while ( (rp = (AliPHOSRecParticle *)nextRecPart() ) ) {
1313         AliPHOSTrackSegment * ts = rp->GetPHOSTrackSegment() ;
1314         if ( ts->GetPHOSMod() == module ) {
1315           Int_t numberofprimaries = 0 ;
1316           Int_t * listofprimaries = 0;
1317           rp->GetPrimaries(numberofprimaries) ;
1318           cout << "Number of primaries = " << numberofprimaries << endl ;
1319           Int_t index ;
1320           for ( index = 0 ; index < numberofprimaries ; index++)
1321             cout << "    primary # " << index << " =  " << listofprimaries[index] << endl ;
1322
1323           nRecParticlesInModule++ ;
1324           Double_t theta = rp->Theta() * kRADDEG ;
1325           Double_t phi   = rp->Phi() * kRADDEG ;
1326           Double_t energy = rp->Energy() ;
1327           histoRparticle->Fill(phi, theta, energy) ;
1328         }
1329       }
1330       histoRparticle->Draw("color") ;
1331
1332       nextRecPart.Reset() ;
1333       while ( (rp = (AliPHOSRecParticle *)nextRecPart() ) ) {
1334         AliPHOSTrackSegment * ts = rp->GetPHOSTrackSegment() ;
1335         if ( ts->GetPHOSMod() == module )
1336           rp->Draw("P") ;
1337       }
1338
1339       Text_t text[80] ;
1340       sprintf(text, "reconstructed particles: %d", nRecParticlesInModule) ;
1341       TPaveText *  pavetext = new TPaveText(292, 100, 300, 101);
1342       pavetext->AddText(text) ;
1343       pavetext->Draw() ;
1344       rparticlecanvas->Update() ;
1345     }
1346 }
1347
1348 //____________________________________________________________________________
1349 void AliPHOSAnalyze::DisplayRecPoints()
1350 {
1351   // Display reconstructed points in local PHOS-module (x, z) coordinates.
1352   // One PHOS module at the time.
1353   // Click on symbols displays the EMC cluster, or PPSD information.
1354
1355   if (fEvt == -999) {
1356     cout << "DisplayRecPoints > Analyze an event first ... (y/n) " ;
1357     Text_t answer[1] ;
1358     cin >> answer ; cout << answer ;
1359 //     if ( answer == "y" )
1360 //       AnalyzeOneEvent() ;
1361   }
1362     if (fEvt != -999) {
1363
1364       Int_t module ;
1365       cout <<  "DisplayRecPoints > which module (1-5,  -1: all) ? " ;
1366       cin >> module ; cout << module << endl ;
1367
1368       Text_t canvasname[80];
1369       sprintf(canvasname,"Digits in PHOS/EMC module # %d",module) ;
1370       TCanvas * modulecanvas = new TCanvas("module", canvasname, 650, 500) ;
1371       modulecanvas->Draw() ;
1372
1373       //=========== Creating 2d-histogram of the PHOS module
1374       // a little bit junkie but is used to test Geom functinalities
1375
1376       Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
1377
1378       fGeom->EmcModuleCoverage(module, tm, tM, pm, pM);
1379       // convert angles into coordinates local to the EMC module of interest
1380
1381       Int_t emcModuleNumber ;
1382       Double_t emcModulexm, emcModulezm ; // minimum local coordinate in a given EMCA module
1383       Double_t emcModulexM, emcModulezM ; // maximum local coordinate in a given EMCA module
1384       fGeom->ImpactOnEmc(tm, pm, emcModuleNumber, emcModulezm, emcModulexm) ;
1385       fGeom->ImpactOnEmc(tM, pM, emcModuleNumber, emcModulezM, emcModulexM) ;
1386       Int_t xdim = (Int_t)( ( emcModulexM - emcModulexm ) / fGeom->GetCrystalSize(0) ) ;
1387       Int_t zdim = (Int_t)( ( emcModulezM - emcModulezm ) / fGeom->GetCrystalSize(2) ) ;
1388       Float_t xmin = emcModulexm - fGeom->GetCrystalSize(0) ;
1389       Float_t xMax = emcModulexM + fGeom->GetCrystalSize(0) ;
1390       Float_t zmin = emcModulezm - fGeom->GetCrystalSize(2) ;
1391       Float_t zMax = emcModulezM + fGeom->GetCrystalSize(2) ;
1392       Text_t histoname[80];
1393       sprintf(histoname,"Event %d: Digits and RecPoints in module %d", fEvt, module) ;
1394       TH2F * hModule = new TH2F("HistoReconstructed", histoname,
1395                                 xdim, xmin, xMax, zdim, zmin, zMax) ;
1396       hModule->SetMaximum(2.0);
1397       hModule->SetMinimum(0.0);
1398       hModule->SetStats(kFALSE);
1399
1400       TIter next(fPHOS->Digits()) ;
1401       Float_t energy, y, z;
1402       Float_t etot=0.;
1403       Int_t relid[4]; Int_t nDigits = 0 ;
1404       AliPHOSDigit * digit ;
1405
1406       // Making 2D histogram of the EMC module
1407       while((digit = (AliPHOSDigit *)next()))
1408         {
1409           fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
1410           if (relid[0] == module && relid[1] == 0)
1411             {
1412               energy = fClu->Calibrate(digit->GetAmp()) ;
1413               cout << "Energy is " << energy << " and threshold is " << fClu->GetEmcEnergyThreshold() << endl;
1414               if (energy >  fClu->GetEmcEnergyThreshold()  ){
1415                 nDigits++ ;
1416                 etot += energy ;
1417                 fGeom->RelPosInModule(relid,y,z) ;
1418                 hModule->Fill(y, z, energy) ;
1419               }
1420             }
1421         }
1422       cout <<"DrawRecPoints >  Found in module "
1423            << module << " " << nDigits << "  digits with total energy " << etot << endl ;
1424       hModule->Draw("col2") ;
1425
1426       //=========== Cluster in module
1427
1428       //      TClonesArray * emcRP = fPHOS->EmcClusters() ;
1429       TObjArray * emcRP = *(fPHOS->EmcRecPoints()) ;
1430
1431       etot = 0.;
1432       Int_t totalnClusters = 0 ;
1433       Int_t nClusters = 0 ;
1434       TIter nextemc(emcRP) ;
1435       AliPHOSEmcRecPoint * emc ;
1436       while((emc = (AliPHOSEmcRecPoint *)nextemc()))
1437         {
1438           //      Int_t numberofprimaries ;
1439           //      Int_t * primariesarray = new Int_t[10] ;
1440           //      emc->GetPrimaries(numberofprimaries, primariesarray) ;
1441           totalnClusters++ ;
1442           if ( emc->GetPHOSMod() == module )
1443             {
1444               nClusters++ ;
1445               energy = emc->GetTotalEnergy() ;
1446               etot+= energy ;
1447               emc->Draw("M") ;
1448             }
1449         }
1450       cout << "DrawRecPoints > Found " << totalnClusters << " EMC Clusters in PHOS" << endl ;
1451       cout << "DrawRecPoints > Found in module " << module << "  " << nClusters << " EMC Clusters " << endl ;
1452       cout << "DrawRecPoints > total energy  " << etot << endl ;
1453
1454       TPaveText *  pavetext = new TPaveText(22, 80, 83, 90);
1455       Text_t text[40] ;
1456       sprintf(text, "digits: %d;  clusters: %d", nDigits, nClusters) ;
1457       pavetext->AddText(text) ;
1458       pavetext->Draw() ;
1459       modulecanvas->Update();
1460
1461       //=========== Cluster in module PPSD Down
1462
1463       //      TClonesArray * ppsdRP = fPHOS->PpsdClusters() ;
1464       TObjArray * ppsdRP = *(fPHOS->PpsdRecPoints() );
1465
1466       etot = 0.;
1467       TIter nextPpsd(ppsdRP) ;
1468       AliPHOSPpsdRecPoint * ppsd ;
1469       while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd()))
1470         {
1471           totalnClusters++ ;
1472           if ( ppsd->GetPHOSMod() == module )
1473             {
1474               nClusters++ ;
1475               energy = ppsd->GetEnergy() ;
1476               etot+=energy ;
1477               if (!ppsd->GetUp()) ppsd->Draw("P") ;
1478             }
1479         }
1480       cout << "DrawRecPoints > Found " << totalnClusters << " Ppsd Down Clusters in PHOS" << endl ;
1481       cout << "DrawRecPoints > Found in module " << module << "  " << nClusters << " Ppsd Down Clusters " << endl ;
1482       cout << "DrawRecPoints > total energy  " << etot << endl ;
1483
1484       //=========== Cluster in module PPSD Up
1485
1486       ppsdRP = *(fPHOS->PpsdRecPoints()) ;
1487
1488       etot = 0.;
1489       TIter nextPpsdUp(ppsdRP) ;
1490       while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsdUp()))
1491         {
1492           totalnClusters++ ;
1493           if ( ppsd->GetPHOSMod() == module )
1494             {
1495               nClusters++ ;
1496               energy = ppsd->GetEnergy() ;
1497               etot+=energy ;
1498               if (ppsd->GetUp()) ppsd->Draw("P") ;
1499             }
1500         }
1501   cout << "DrawRecPoints > Found " << totalnClusters << " Ppsd Up Clusters in PHOS" << endl ;
1502   cout << "DrawRecPoints > Found in module " << module << "  " << nClusters << " Ppsd Up Clusters " << endl ;
1503   cout << "DrawRecPoints > total energy  " << etot << endl ;
1504
1505     } // if !-999
1506 }
1507
1508 //____________________________________________________________________________
1509 void AliPHOSAnalyze::DisplayTrackSegments()
1510 {
1511   // Display track segments in local PHOS-module (x, z) coordinates.
1512   // One PHOS module at the time.
1513   // One symbol per PHOS subsystem: EMC, upper PPSD, lower PPSD.
1514
1515   if (fEvt == -999) {
1516     cout << "DisplayTrackSegments > Analyze an event first ... (y/n) " ;
1517     Text_t answer[1] ;
1518     cin >> answer ; cout << answer ;
1519 //     if ( answer == "y" )
1520 //       AnalyzeOneEvent() ;
1521   }
1522     if (fEvt != -999) {
1523
1524       Int_t module ;
1525       cout <<  "DisplayTrackSegments > which module (1-5,  -1: all) ? " ;
1526       cin >> module ; cout << module << endl ;
1527       //=========== Creating 2d-histogram of the PHOS module
1528       // a little bit junkie but is used to test Geom functinalities
1529
1530       Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module
1531
1532       fGeom->EmcModuleCoverage(module, tm, tM, pm, pM);
1533       // convert angles into coordinates local to the EMC module of interest
1534
1535       Int_t emcModuleNumber ;
1536       Double_t emcModulexm, emcModulezm ; // minimum local coordinate in a given EMCA module
1537       Double_t emcModulexM, emcModulezM ; // maximum local coordinate in a given EMCA module
1538       fGeom->ImpactOnEmc(tm, pm, emcModuleNumber, emcModulezm, emcModulexm) ;
1539       fGeom->ImpactOnEmc(tM, pM, emcModuleNumber, emcModulezM, emcModulexM) ;
1540       Int_t xdim = (Int_t)( ( emcModulexM - emcModulexm ) / fGeom->GetCrystalSize(0) ) ;
1541       Int_t zdim = (Int_t)( ( emcModulezM - emcModulezm ) / fGeom->GetCrystalSize(2) ) ;
1542       Float_t xmin = emcModulexm - fGeom->GetCrystalSize(0) ;
1543       Float_t xMax = emcModulexM + fGeom->GetCrystalSize(0) ;
1544       Float_t zmin = emcModulezm - fGeom->GetCrystalSize(2) ;
1545       Float_t zMax = emcModulezM + fGeom->GetCrystalSize(2) ;
1546       Text_t histoname[80];
1547       sprintf(histoname,"Event %d: Track Segments in module %d", fEvt, module) ;
1548       TH2F * histotrack = new TH2F("histotrack",  histoname,
1549                                    xdim, xmin, xMax, zdim, zmin, zMax) ;
1550       histotrack->SetStats(kFALSE);
1551       Text_t canvasname[80];
1552       sprintf(canvasname,"Track segments in PHOS/EMC-PPSD module # %d", module) ;
1553       TCanvas * trackcanvas = new TCanvas("TrackSegmentCanvas", canvasname, 650, 500) ;
1554       histotrack->Draw() ;
1555
1556       AliPHOSTrackSegment::TrackSegmentsList * trsegl = *(fPHOS->TrackSegments()) ;
1557       AliPHOSTrackSegment * trseg ;
1558
1559       Int_t nTrackSegments = trsegl->GetEntries() ;
1560       Int_t index ;
1561       Float_t etot = 0 ;
1562       Int_t nTrackSegmentsInModule = 0 ;
1563       for(index = 0; index < nTrackSegments ; index++){
1564         trseg = (AliPHOSTrackSegment * )trsegl->At(index) ;
1565         etot+= trseg->GetEnergy() ;
1566         if ( trseg->GetPHOSMod() == module ) {
1567           nTrackSegmentsInModule++ ;
1568           trseg->Draw("P");
1569         }
1570       }
1571       Text_t text[80] ;
1572       sprintf(text, "track segments: %d", nTrackSegmentsInModule) ;
1573       TPaveText *  pavetext = new TPaveText(22, 80, 83, 90);
1574       pavetext->AddText(text) ;
1575       pavetext->Draw() ;
1576       trackcanvas->Update() ;
1577       cout << "DisplayTrackSegments > Found " << trsegl->GetEntries() << " Track segments with total energy "<< etot << endl ;
1578
1579    }
1580 }
1581 //____________________________________________________________________________
1582 Bool_t AliPHOSAnalyze::OpenRootFile(Text_t * name)
1583 {
1584   // Open the root file named "name"
1585
1586   fRootFile   = new TFile(name, "update") ;
1587   return  fRootFile->IsOpen() ;
1588 }
1589 //____________________________________________________________________________
1590 // void AliPHOSAnalyze::SavingHistograms()
1591 // {
1592 //   // Saves the histograms in a root file named "name.analyzed"
1593
1594 //   Text_t outputname[80] ;
1595 //   sprintf(outputname,"%s.analyzed",fRootFile->GetName());
1596 //   TFile output(outputname,"RECREATE");
1597 //   output.cd();
1598 //   if (fhEmcDigit )
1599 //     fhEmcDigit->Write()  ;
1600 //   if (fhVetoDigit )
1601 //     fhVetoDigit->Write()  ;
1602 //   if (fhConvertorDigit )
1603 //     fhConvertorDigit->Write()   ;
1604 //   if (fhEmcCluster   )
1605 //     fhEmcCluster->Write()   ;
1606 //   if (fhVetoCluster )
1607 //     fhVetoCluster->Write()   ;
1608 //   if (fhConvertorCluster )
1609 //     fhConvertorCluster->Write()  ;
1610 //   if (fhConvertorEmc )
1611 //     fhConvertorEmc->Write()  ;
1612 //   if (fhPhotonEnergy)
1613 //     fhPhotonEnergy->Write() ;
1614 //   if (fhPhotonPositionX)
1615 //     fhPhotonPositionX->Write() ;
1616 //   if (fhPhotonPositionY)
1617 //     fhPhotonPositionX->Write() ;
1618 //   if (fhElectronEnergy)
1619 //     fhElectronEnergy->Write() ;
1620 //   if (fhElectronPositionX)
1621 //     fhElectronPositionX->Write() ;
1622 //   if (fhElectronPositionY)
1623 //     fhElectronPositionX->Write() ;
1624 //   if (fhNeutralHadronEnergy)
1625 //     fhNeutralHadronEnergy->Write() ;
1626 //   if (fhNeutralHadronPositionX)
1627 //     fhNeutralHadronPositionX->Write() ;
1628 //   if (fhNeutralHadronPositionY)
1629 //     fhNeutralHadronPositionX->Write() ;
1630 //   if (fhNeutralEMEnergy)
1631 //     fhNeutralEMEnergy->Write() ;
1632 //   if (fhNeutralEMPositionX)
1633 //     fhNeutralEMPositionX->Write() ;
1634 //   if (fhNeutralEMPositionY)
1635 //     fhNeutralEMPositionX->Write() ;
1636 //   if (fhChargedHadronEnergy)
1637 //     fhChargedHadronEnergy->Write() ;
1638 //   if (fhChargedHadronPositionX)
1639 //     fhChargedHadronPositionX->Write() ;
1640 //   if (fhChargedHadronPositionY)
1641 //     fhChargedHadronPositionX->Write() ;
1642 //   if (fhPhotonHadronEnergy)
1643 //     fhPhotonHadronEnergy->Write() ;
1644 //   if (fhPhotonHadronPositionX)
1645 //     fhPhotonHadronPositionX->Write() ;
1646 //   if (fhPhotonHadronPositionY)
1647 //     fhPhotonHadronPositionX->Write() ;
1648
1649 //   output.Write();
1650 //   output.Close();
1651 // }
1652 //____________________________________________________________________________
1653 void AliPHOSAnalyze::SaveHistograms()
1654 {
1655   // Saves the histograms in a root file named "name.analyzed"
1656
1657   Text_t outputname[80] ;
1658   sprintf(outputname,"%s.analyzed",fRootFile->GetName());
1659   TFile output(outputname,"RECREATE");
1660   output.cd();
1661
1662   if (fhPhotonEnergy)
1663     fhPhotonEnergy->Write() ;
1664   if (fhPhotonAllEnergy)
1665     fhPhotonAllEnergy->Write() ;
1666   if (fhPhotonPosition)
1667     fhPhotonPosition->Write() ;
1668   if (fhPhotonAllPosition)
1669     fhPhotonAllPosition->Write() ;
1670   if (fhElectronEnergy)
1671     fhElectronEnergy->Write() ;
1672   if (fhElectronAllEnergy)
1673     fhElectronAllEnergy->Write() ;
1674   if (fhElectronPosition)
1675     fhElectronPosition->Write() ;
1676   if (fhElectronAllPosition)
1677     fhElectronAllPosition->Write() ;
1678   if (fhNeutralHadronEnergy)
1679     fhNeutralHadronEnergy->Write() ;
1680   if (fhNeutralHadronPosition)
1681     fhNeutralHadronPosition->Write() ;
1682   if (fhNeutralEMEnergy)
1683     fhNeutralEMEnergy->Write() ;
1684   if (fhNeutralEMAllEnergy)
1685     fhNeutralEMAllEnergy->Write() ;
1686   if (fhNeutralEMPosition)
1687     fhNeutralEMPosition->Write() ;
1688   if (fhNeutralEMAllPosition)
1689     fhNeutralEMAllPosition->Write() ;
1690   if (fhChargedHadronEnergy)
1691     fhChargedHadronEnergy->Write() ;
1692   if (fhChargedHadronPosition)
1693     fhChargedHadronPosition->Write() ;
1694   if (fhPhotonHadronEnergy)
1695     fhPhotonHadronEnergy->Write() ;
1696   if (fhPhotonHadronPosition)
1697     fhPhotonHadronPosition->Write() ;
1698   if (fhPhotonReg)
1699     fhPhotonReg->Write() ;
1700   if (fhAllReg)
1701     fhAllReg->Write() ;
1702   if(fhNReg)
1703     fhNReg->Write() ;
1704   if(fhNBarReg)
1705     fhNBarReg->Write() ;
1706   if(fhChargedReg)
1707     fhChargedReg->Write() ;
1708   if (fhPhotonEM)
1709     fhPhotonEM->Write() ;
1710   if (fhAllEM)
1711     fhAllEM->Write() ;
1712   if(fhNEM)
1713     fhNEM->Write() ;
1714   if(fhNBarEM)
1715     fhNBarEM->Write() ;
1716   if(fhChargedEM)
1717     fhChargedEM->Write() ;
1718   if(fhPrimary)
1719     fhPrimary->Write() ;
1720   if(fhPhotPhot)
1721     fhPhotPhot->Write() ;
1722   if(fhPhotElec)
1723     fhPhotElec->Write() ;
1724   if(fhPhotNeuH)
1725     fhPhotNeuH->Write() ;
1726   if(fhPhotNuEM)
1727     fhPhotNuEM->Write() ;
1728   if(fhPhotNuEM)
1729     fhPhotNuEM->Write() ;
1730   if(fhPhotChHa)
1731     fhPhotChHa->Write() ;
1732   if(fhPhotGaHa)
1733     fhPhotGaHa->Write() ;
1734   if(fhEnergyCorrelations)
1735     fhEnergyCorrelations->Write() ;
1736
1737   output.Write();
1738   output.Close();
1739 }
1740 //____________________________________________________________________________
1741 void AliPHOSAnalyze::ResetHistograms()
1742 {
1743    fhEnergyCorrelations = 0 ;     //Energy correlations between Eloss in Convertor and PPSD(2)
1744
1745    fhEmcDigit = 0 ;               // Histo of digit energies in the Emc
1746    fhVetoDigit = 0 ;              // Histo of digit energies in the Veto
1747    fhConvertorDigit = 0 ;         // Histo of digit energies in the Convertor
1748    fhEmcCluster = 0 ;             // Histo of Cluster energies in Emc
1749    fhVetoCluster = 0 ;            // Histo of Cluster energies in Veto
1750    fhConvertorCluster = 0 ;       // Histo of Cluster energies in Convertor
1751    fhConvertorEmc = 0 ;           // 2d Convertor versus Emc energies
1752
1753    fhPhotonEnergy = 0 ;           // Spectrum of detected photons with photon primary
1754    fhPhotonAllEnergy = 0 ;        // Total spectrum of detected photons
1755    fhElectronEnergy = 0 ;         // Spectrum of detected electrons with electron primary
1756    fhElectronAllEnergy = 0 ;      // Total spectrum of detected electrons
1757    fhNeutralHadronEnergy = 0 ;    // Spectrum of detected neutral hadron
1758    fhNeutralEMEnergy = 0 ;        // Spectrum of detected neutral EM with EM primary
1759    fhNeutralEMAllEnergy = 0 ;     // Spectrum of detected neutral EM
1760    fhChargedHadronEnergy = 0 ;    // Spectrum of detected charged
1761    fhPhotonHadronEnergy = 0 ;     // Spectrum of detected Photon-Hadron
1762    fhPhotonPosition = 0 ;        // Position Resolution of  photons with photon primary
1763    fhPhotonAllPosition = 0 ;     // Position Resolution of  photons
1764    fhElectronPosition = 0 ;      // Position Resolution of electrons with electron primary
1765    fhElectronAllPosition = 0 ;   // Position Resolution of electrons
1766    fhNeutralHadronPosition = 0 ; // Position Resolution of neutral hadron
1767    fhNeutralEMPosition = 0 ;     // Position Resolution of neutral EM with EM primary
1768    fhNeutralEMAllPosition = 0 ;  // Position Resolution of neutral EM
1769    fhChargedHadronPosition = 0 ; // Position Resolution of charged
1770    fhPhotonHadronPosition = 0 ;  // Position Resolution of Photon-Hadron
1771    fhPhotonPositionY = 0 ;        // Y distribution of detected photons
1772    fhElectronPositionY = 0 ;      // Y distribution of detected electrons
1773    fhNeutralHadronPositionY = 0 ; // Y distribution of detected neutral hadron
1774    fhNeutralEMPositionY = 0 ;     // Y distribution of detected neutral EM
1775    fhChargedHadronPositionY = 0 ; // Y distribution of detected charged
1776    fhPhotonHadronPositionY = 0 ;  // Y distribution of detected Photon-Hadron
1777    fhPhotonReg = 0 ;
1778    fhAllReg = 0 ;
1779    fhNReg = 0 ;
1780    fhNBarReg = 0 ;
1781    fhChargedReg = 0 ;
1782    fhPhotonEM = 0 ;
1783    fhAllEM = 0 ;
1784    fhNEM = 0 ;
1785    fhNBarEM = 0 ;
1786    fhChargedEM = 0 ;
1787    fhPrimary = 0 ;
1788
1789    fhPhotPhot = 0 ;
1790    fhPhotElec = 0 ;
1791    fhPhotNeuH = 0 ;
1792    fhPhotNuEM = 0 ;
1793    fhPhotChHa = 0 ;
1794    fhPhotGaHa = 0 ;
1795
1796
1797 }
1798
1799