Clean up to correct for the mess introduced by my eratic branching !
[u/mrichter/AliRoot.git] / 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