PHOS/ReconstructionTest.cxx

   1
   2 /**************************************************************************
   3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   4  *                                                                        *
   5  * Author: The ALICE Off-line Project.                                    *
   6  * Contributors are mentioned in the code where appropriate.              *
   7  *                                                                        *
   8  * Permission to use, copy, modify and distribute this software and its   *
   9  * documentation strictly for non-commercial purposes is hereby granted   *
  10  * without fee, provided that the above copyright notice appears in all   *
  11  * copies and that both the copyright notice and this permission notice   *
  12  * appear in the supporting documentation. The authors make no claims     *
  13  * about the suitability of this software for any purpose. It is          *
  14  * provided "as is" without express or implied warranty.                  *
  15  **************************************************************************/
  16
  17 //_________________________________________________________________________
  18 // Cheking Reconstruction procedure of PHOS
  19 //*-- Author : Gines MARTINEZ  SUBATECH january 2000
  20 //////////////////////////////////////////////////////////////////////////////
  21
  22 // --- ROOT system ---
  23
  24 #include "TFile.h"
  25 #include "TH1.h"
  26 #include "TPad.h"
  27 #include "TTree.h"
  28 #include "TCanvas.h"
  29 #include "TParticle.h"
  30
  31 // --- Standard library ---
  32
  33 #include <iostream>
  34 #include <cstdio>
  35
  36
  37 // --- AliRoot header files ---
  38 #include "AliRun.h"
  39 #include "TFile.h"
  40 #include "AliPHOSGeometry.h"
  41 #include "AliPHOSv0.h"
  42 #include "AliPHOSDigit.h"
  43 #include "AliPHOSRecPoint.h"
  44 #include "AliPHOSEmcRecPoint.h"
  45 #include "AliPHOSPpsdRecPoint.h"
  46 #include "AliPHOSClusterizerv1.h"
  47 #include "AliPHOSReconstructioner.h"
  48 #include "AliPHOSTrackSegment.h"
  49 #include "AliPHOSTrackSegmentMakerv1.h"
  50 #include "ReconstructionTest.h"
  51
  52
  53 void ReconstructionTest(Text_t* infile,Int_t evt, Int_t Module)
  54 {
  55   //========== Opening galice.root file
  56   TFile * file   = new TFile(infile);
  57   //========== Get AliRun object from file
  58   gAlice = (AliRun*) file->Get("gAlice");
  59   //=========== Gets the PHOS object and associated geometry from the file
  60   AliPHOSv0 * PHOS  = (AliPHOSv0 *)gAlice->GetDetector("PHOS");
  61   AliPHOSGeometry * Geom = AliPHOSGeometry::GetInstance(PHOS->GetGeometry()->GetName(),PHOS->GetGeometry()->GetTitle());
  62   //========== Creates the Clusterizer
  63   AliPHOSClusterizerv1 clusterizer;
  64   clusterizer.SetEmcEnergyThreshold(0.01) ;
  65   clusterizer.SetEmcClusteringThreshold(0.1) ;
  66   clusterizer.SetPpsdEnergyThreshold(0.0000001) ;
  67   clusterizer.SetPpsdClusteringThreshold(0.0000002) ;
  68   clusterizer.SetLocalMaxCut(0.03) ;
  69   clusterizer.SetCalibrationParameters(0., 0.0000001) ;
  70   //========== Creates the track segment maker
  71   AliPHOSTrackSegmentMakerv1 tracksegmentmaker ;
  72   //========== Creates the Reconstructioner
  73   AliPHOSReconstructioner Reconstructioner(clusterizer,tracksegmentmaker);
  74   //=========== Connects the various Tree's for evt
  75   gAlice->GetEvent(evt);
  76   //=========== Gets the Digit TTree
  77   gAlice->TreeD()->GetEvent(0) ;
  78   //=========== Gets the number of entries in the Digits array
  79   Int_t nId = PHOS->Digits()->GetEntries();
  80   printf("ReconstructionTest> Number of entries in the Digit array is %d \n",nId);
  81   //=========== Do the reconstruction
  82   AliPHOSDigit * digit ;
  83   TIter next(PHOS->Digits()) ;
  84   Float_t Etot=0 ;
  85   while( ( digit = (AliPHOSDigit *)next() ) )
  86     Etot+=clusterizer.Calibrate(digit->GetAmp()) ;
  87   cout <<"ReconstructionTest> Found  " << nId << "  digits in PHOS with total energy " << Etot << endl ;
  88   PHOS->Reconstruction(Reconstructioner);
  89   //================Make checks===========================
  90
  91   //=========== Creating Canvas
  92   TCanvas * ModuleCanvas = new TCanvas("Module","Events in a single PHOS Module", 650, 500) ;
  93   ModuleCanvas->Draw() ;
  94
  95   //=========== Creating 2d-histogram of the PHOS Module
  96   // a little bit junkie but is used to test Geom functinalities
  97
  98   Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by Module 1
  99
 100   Geom->EmcModuleCoverage(1, tm, tM, pm, pM);
 101   // convert angles into coordinates local to the EMC module of interest
 102
 103   Int_t EmcModuleNumber ;
 104   Double_t EmcModulexm, EmcModulezm ; // minimum local coordinate in a given EMCA module
 105   Double_t EmcModulexM, EmcModulezM ; // maximum local coordinate in a given EMCA module
 106   Geom->ImpactOnEmc(tm, pm, EmcModuleNumber, EmcModulezm, EmcModulexm) ;
 107   Geom->ImpactOnEmc(tM, pM, EmcModuleNumber, EmcModulezM, EmcModulexM) ;
 108   Int_t xdim = (Int_t)( ( EmcModulexM - EmcModulexm ) / Geom->GetCrystalSize(0) ) ;
 109   Int_t zdim = (Int_t)( ( EmcModulezM - EmcModulezm ) / Geom->GetCrystalSize(2) ) ;
 110   Float_t xmin = EmcModulexm - Geom->GetCrystalSize(0) ;
 111   Float_t xMax = EmcModulexM + Geom->GetCrystalSize(0) ;
 112   Float_t zmin = EmcModulezm - Geom->GetCrystalSize(2) ;
 113   Float_t zMax = EmcModulezM + Geom->GetCrystalSize(2) ;
 114   // histogram of reconstructed events
 115   Text_t HistoName[80];
 116   sprintf(HistoName,"Event %d: Reconstructed particles in module %d", evt, Module) ;
 117   TH2F * hModule = new TH2F("HistoReconstructed", HistoName,
 118                             xdim, xmin, xMax, zdim, zmin, zMax) ;
 119   hModule->SetMaximum(2.0);
 120   hModule->SetMinimum(0.0);
 121   hModule->SetStats(kFALSE);
 122   // histogram of generated  particles
 123   sprintf(HistoName,"Event %d: Incident particles in module %d", evt, Module) ;
 124   TH2F * HistoParticle = new TH2F("HistoParticle",  HistoName,
 125                                   xdim, xmin, xMax, zdim, zmin, zMax) ;
 126   HistoParticle->SetStats(kFALSE) ;
 127
 128   //=========== Digits of Module
 129   TIter next2(PHOS->Digits()) ;
 130   Float_t Energy, y, z;
 131   Int_t RelId[4]; Int_t NumberOfDigits = 0 ;
 132   while((digit = (AliPHOSDigit *)next2()))
 133     {
 134       Geom->AbsToRelNumbering(digit->GetId(), RelId) ;
 135       if (RelId[0] == Module)
 136         {
 137           NumberOfDigits++ ;
 138           Energy = clusterizer.Calibrate(digit->GetAmp()) ;
 139           Etot+=Energy ;
 140           Geom->RelPosInModule(RelId,y,z) ;
 141           if (Energy>0.01 )
 142             hModule->Fill(y,z,Energy) ;
 143         }
 144   }
 145   cout <<"TestRec> Found in Module " << Module << " " << NumberOfDigits << "  digits with total energy " << Etot << endl ;
 146   hModule->Draw("col2") ;
 147
 148   //=========== Cluster in Module
 149   TClonesArray * EmcRP = PHOS->EmcClusters() ;
 150   Etot = 0.;
 151   Int_t TotalNumberOfClusters = 0 ;
 152   Int_t NumberOfClusters = 0 ;
 153   TIter nextemc(EmcRP) ;
 154   AliPHOSEmcRecPoint * emc ;
 155   while((emc = (AliPHOSEmcRecPoint *)nextemc()))
 156     {
 157       TotalNumberOfClusters++ ;
 158       if ( emc->GetPHOSMod() == Module )
 159         {
 160           NumberOfClusters++ ;
 161           Energy = emc->GetTotalEnergy() ;
 162           Etot+=Energy ;
 163           emc->Draw("P") ;
 164     }
 165     }
 166   cout << "TestRec> Found " << TotalNumberOfClusters << " EMC Clusters in PHOS" << endl ;
 167   cout << "TestRec> Found in Module " << Module << "  " << NumberOfClusters << " EMC Clusters " << endl ;
 168   cout << "TestRec> Total energy  " <<Etot << endl ;
 169
 170   TPaveText *  PaveText = new TPaveText(22, 80, 83, 90);
 171   Text_t text[40] ;
 172   sprintf(text, "digits: %d;  clusters: %d", NumberOfDigits, NumberOfClusters) ;
 173   PaveText->AddText(text) ;
 174   PaveText->Draw() ;
 175   ModuleCanvas->Update();
 176
 177   //=========== Cluster in Module PPSD Down
 178   TClonesArray * PpsdRP = PHOS->PpsdClusters() ;
 179   Etot = 0.;
 180   TIter nextPpsd(PpsdRP) ;
 181   AliPHOSPpsdRecPoint * Ppsd ;
 182   while((Ppsd = (AliPHOSPpsdRecPoint *)nextPpsd()))
 183   {
 184     TotalNumberOfClusters++ ;
 185     if ( Ppsd->GetPHOSMod() == Module )
 186     {
 187       NumberOfClusters++ ;
 188       Energy = Ppsd->GetEnergy() ;
 189       Etot+=Energy ;
 190       if (!Ppsd->GetUp()) Ppsd->Draw("P") ;
 191     }
 192   }
 193   cout << "TestRec> Found " << TotalNumberOfClusters << " Ppsd Down Clusters in PHOS" << endl ;
 194   cout << "TestRec> Found in Module " << Module << "  " << NumberOfClusters << " Ppsd Down Clusters " << endl ;
 195   cout << "TestRec> Total energy  " <<Etot << endl ;
 196
 197 //=========== Cluster in Module PPSD Up
 198   PpsdRP = PHOS->PpsdClusters() ;
 199   Etot = 0.;
 200   TIter nextPpsdUp(PpsdRP) ;
 201   while((Ppsd = (AliPHOSPpsdRecPoint *)nextPpsdUp()))
 202   {
 203     TotalNumberOfClusters++ ;
 204     if ( Ppsd->GetPHOSMod() == Module )
 205     {
 206       NumberOfClusters++ ;
 207       Energy = Ppsd->GetEnergy() ;
 208       Etot+=Energy ;
 209       if (Ppsd->GetUp()) Ppsd->Draw("P") ;
 210     }
 211   }
 212   cout << "TestRec> Found " << TotalNumberOfClusters << " Ppsd Up Clusters in PHOS" << endl ;
 213   cout << "TestRec> Found in Module " << Module << "  " << NumberOfClusters << " Ppsd Up Clusters " << endl ;
 214   cout << "TestRec> Total energy  " <<Etot << endl ;
 215
 216    // Get pointers to Alice Particle TClonesArray
 217   TParticle * Particle;
 218   TClonesArray * ArrayOfParticles  = gAlice->Particles();
 219   TCanvas * KineCanvas = new TCanvas("KineCnvas", "Incident particles", 650, 500) ;
 220   // get the KINE Tree
 221   TTree * Kine =  gAlice->TreeK() ;
 222   Stat_t NumberOfParticles =  Kine->GetEntries() ;
 223   cout << "events in Kine " << NumberOfParticles << endl ;
 224
 225   // loop over particles
 226   Int_t index1 ;
 227   Int_t nparticlein = 0 ;
 228   for (index1 = 0 ; index1 < NumberOfParticles ; index1++){
 229     Int_t nparticle = ArrayOfParticles->GetEntriesFast() ;
 230     cout << nparticle << endl ;
 231     Int_t index2 ;
 232     for( index2 = 0 ; index2 < nparticle; index2++) {
 233       Particle            = (TParticle*)ArrayOfParticles->UncheckedAt(index2) ;
 234       Int_t  ParticleType = Particle->GetPdgCode() ;
 235       Double_t Phi        = Particle->Phi() ;
 236       Double_t Theta      = Particle->Theta() ;
 237       Int_t mod ;
 238       Double_t x, z ;
 239       Geom->ImpactOnEmc(Theta, Phi, mod, z, x) ;
 240       if ( mod == Module ) {
 241         nparticlein++ ;
 242         HistoParticle->Fill( x, -z, Particle->Energy() ) ; //-z don't know why, but that is how it works
 243       }
 244     }
 245   }
 246   KineCanvas->Draw() ;
 247   HistoParticle->Draw("color") ;
 248   TPaveText *  PaveText2 = new TPaveText(22, 80, 83, 90);
 249   sprintf(text, "Particles: %d ", nparticlein) ;
 250   PaveText2->AddText(text) ;
 251   PaveText2->Draw() ;
 252   KineCanvas->Update();
 253
 254 //  TObjArray * trsegl = PHOS->TrackSegments() ;
 255 //    AliPHOSTrackSegment trseg ;
 256
 257 //    Int_t NTrackSegments = trsegl->GetEntries() ;
 258 //    Int_t index ;
 259 //    Etot = 0 ;
 260 //    for(index = 0; index < NTrackSegments ; index++){
 261 //      trseg = (AliPHOSTrackSegment * )trsegl->At(index) ;
 262 //      Etot+= trseg->GetEnergy() ;
 263 //      if ( trseg->GetPHOSMod() == Module )
 264 //      {
 265 // //       trseg->Draw("P");
 266 //        trseg->Print() ;
 267 //      }
 268 //    }
 269 //    cout << "Found " << trsegl->GetEntries() << " Track segments with total energy "<< Etot << endl ;
 270 //
 271
 272
 273
 274 }
 275
 276