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