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