Effective C++ corrections (T.Pocheptsov)
[u/mrichter/AliRoot.git] / PHOS / AliPHOSIhepAnalyze.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 //_________________________________________________________________________
17 // Algorythm class to analyze PHOSv1 events:
18 // Construct histograms and displays them.
19 // IHEP CPV/PHOS reconstruction algorithm used.
20 // Use the macro EditorBar.C for best access to the functionnalities
21 //*--
22 //*-- Author: B. Polichtchouk (IHEP)
23 //////////////////////////////////////////////////////////////////////////////
24
25 // --- ROOT system ---
26
27 #include "TFile.h"
28 #include "TH1.h"
29 #include "TPad.h"
30 #include "TH2.h"
31 #include "TParticle.h"
32 #include "TClonesArray.h"
33 #include "TTree.h"
34 #include "TMath.h"
35 #include "TCanvas.h" 
36 #include "TStyle.h" 
37
38 // --- Standard library ---
39
40 // --- AliRoot header files ---
41
42 #include "AliRunLoader.h"
43 #include "AliHeader.h"
44
45 // --- PHOS header files ---
46 #include "AliLog.h"
47 #include "AliPHOSIhepAnalyze.h"
48 #include "AliPHOSDigit.h"
49 #include "AliPHOSRecParticle.h"
50 #include "AliPHOSLoader.h"
51 #include "AliPHOSHit.h"
52 #include "AliPHOSImpact.h"
53 #include "AliPHOSvImpacts.h"
54 #include "AliPHOSCpvRecPoint.h"
55 #include "AliRun.h"
56 #include "AliPHOSGeometry.h"
57 #include "AliPHOSEvalRecPoint.h"
58
59 ClassImp(AliPHOSIhepAnalyze)
60
61 //____________________________________________________________________________
62 AliPHOSIhepAnalyze::AliPHOSIhepAnalyze():
63   fRunLoader(0),
64   fFileName()
65 {
66 }
67
68 //____________________________________________________________________________
69 AliPHOSIhepAnalyze::AliPHOSIhepAnalyze(Text_t * name) : 
70   fRunLoader(0),
71   fFileName(name)
72 {
73   // Constructor: open a header file
74   fRunLoader = AliRunLoader::Open(fFileName);
75   if (fRunLoader == 0x0)
76    {
77      AliFatal(Form("Can not load event from file %s",name));
78    }
79 }
80
81 //____________________________________________________________________________
82 void AliPHOSIhepAnalyze::AnalyzeCPV1(Int_t Nevents)
83 {
84   //
85   // Analyzes CPV characteristics: resolutions, cluster multiplicity,
86   // cluster lengths along Z and Phi.
87   // Author: Yuri Kharlov
88   // 9 October 2000
89   // Modified by Boris Polichtchouk, 3.07.2001
90   //
91
92   // Book histograms
93
94   TH1F *hDx   = new TH1F("hDx"  ,"CPV x-resolution@reconstruction",100,-5. , 5.);
95   TH1F *hDz   = new TH1F("hDz"  ,"CPV z-resolution@reconstruction",100,-5. , 5.);
96   TH1F *hChi2 = new TH1F("hChi2"  ,"Chi2/dof of one-gamma fit",30, 0. , 10.);
97   TH1S *hNrp  = new TH1S("hNrp" ,"CPV rec.point multiplicity",      21,-0.5,20.5);
98   TH1S *hNrpX = new TH1S("hNrpX","CPV rec.point Phi-length"  ,      21,-0.5,20.5);
99   TH1S *hNrpZ = new TH1S("hNrpZ","CPV rec.point Z-length"    ,      21,-0.5,20.5);
100
101   TH1F *hEg   = new TH1F("hEg"  ,"Energies of impacts",30,0.,6.);
102   TH1F *hEr   = new TH1F("hEr"  ,"Amplitudes of rec. points",50,0.,20.);
103   
104   TList * fCpvImpacts ;
105   TBranch * branchCPVimpacts;
106
107   
108   
109   AliPHOSLoader* please = dynamic_cast<AliPHOSLoader*>(fRunLoader->GetLoader("PHOSLoader"));
110   if ( please == 0 )
111    {
112      AliError(Form("Could not obtain the Loader object !"));
113      return ;
114    }
115
116   const AliPHOSGeometry *  fGeom  = please->PHOSGeometry();
117
118   AliInfo(Form("Start CPV Analysis-1. Resolutions, cluster multiplicity and lengths")) ;
119   for ( Int_t ievent=0; ievent<Nevents; ievent++) {  
120     
121     Int_t nTotalGen = 0;
122     Int_t nChargedGen = 0;
123
124     Int_t ntracks = gAlice->GetEvent(ievent);
125     AliInfo(Form(">>>>>>>Event %d .<<<<<<<", ievent)) ;
126     
127   /******************************************************************/
128       TTree* treeH = please->TreeH();
129       if (treeH == 0x0)
130        {
131         AliError(Form("Can not get TreeH"));
132          return;
133        }
134 /******************************************************************/     
135
136     // Get branch of CPV impacts
137     if (! (branchCPVimpacts =treeH->GetBranch("PHOSCpvImpacts")) ) {
138       AliWarning(Form("Couldn't find branch PHOSCpvImpacts. Exit.")) ;
139       return;
140     }
141  
142     // Create and fill arrays of hits for each CPV module
143       
144     Int_t nOfModules = fGeom->GetNModules();
145     TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
146     Int_t iModule = 0;  
147     for (iModule=0; iModule < nOfModules; iModule++)
148       hitsPerModule[iModule] = new TClonesArray("AliPHOSImpact",100);
149
150     TClonesArray    *impacts;
151     AliPHOSImpact   *impact;
152     TLorentzVector   p;
153
154     // First go through all primary tracks and fill the arrays
155     // of hits per each CPV module
156
157     for (Int_t itrack=0; itrack<ntracks; itrack++) {
158       branchCPVimpacts ->SetAddress(&fCpvImpacts);
159       branchCPVimpacts ->GetEntry(itrack,0);
160
161       for (Int_t iModule=0; iModule < nOfModules; iModule++) {
162         impacts = (TClonesArray *)fCpvImpacts->At(iModule);
163         // Do loop over impacts in the module
164         for (Int_t iImpact=0; iImpact<impacts->GetEntries(); iImpact++) {
165           impact=(AliPHOSImpact*)impacts->At(iImpact);
166           hEg->Fill(impact->GetMomentum().E());
167           TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
168           if(IsCharged(impact->GetPid())) nChargedGen++;
169           new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSImpact(*impact);
170         }
171       }
172       fCpvImpacts->Clear();
173     }
174     for (iModule=0; iModule < nOfModules; iModule++) {
175       Int_t nsum = hitsPerModule[iModule]->GetEntriesFast();
176       printf("CPV module %d has %d impacts\n",iModule,nsum);
177       nTotalGen += nsum;
178     }
179
180     // Then go through reconstructed points and for each find
181     // the closeset hit
182     // The distance from the rec.point to the closest hit
183     // gives the coordinate resolution of the CPV
184
185     fRunLoader->GetEvent(ievent);
186     TIter nextRP(please->CpvRecPoints()) ;
187     AliPHOSCpvRecPoint *cpvRecPoint ;
188     Float_t xgen, ygen, zgen;
189     while( ( cpvRecPoint = (AliPHOSCpvRecPoint *)nextRP() ) ) {
190       
191       Float_t chi2dof = ((AliPHOSEvalRecPoint*)cpvRecPoint)->Chi2Dof();
192       hChi2->Fill(chi2dof);
193       hEr->Fill(cpvRecPoint->GetEnergy());
194
195       TVector3  locpos;
196       cpvRecPoint->GetLocalPosition(locpos);
197       Int_t phosModule = cpvRecPoint->GetPHOSMod();
198       Int_t rpMult     = cpvRecPoint->GetMultiplicity();
199       Int_t rpMultX, rpMultZ;
200       cpvRecPoint->GetClusterLengths(rpMultX,rpMultZ);
201       Float_t xrec  = locpos.X();
202       Float_t zrec  = locpos.Z();
203       Float_t dxmin = 1.e+10;
204       Float_t dzmin = 1.e+10;
205       Float_t r2min = 1.e+10;
206       Float_t r2;
207
208       Int_t nCPVhits  = (hitsPerModule[phosModule-1])->GetEntriesFast();
209       Float_t locImpX=1e10,locImpZ=1e10;         // local coords of closest impact
210       Float_t gImpX=1e10, gImpZ=1e10,gImpY=1e10; // global coords of closest impact
211       for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
212         impact = (AliPHOSImpact*)(hitsPerModule[phosModule-1])->UncheckedAt(ihit);
213         xgen   = impact->X();
214         zgen   = impact->Z();
215         ygen   = impact->Y();
216         
217         //Transform to the local ref.frame
218         const AliPHOSGeometry* geom = please->PHOSGeometry();
219         Float_t phig = geom->GetPHOSAngle(phosModule);
220         Float_t phi = TMath::Pi()/180*phig;
221         Float_t distanceIPtoCPV = geom->GetIPtoOuterCoverDistance() -
222                           (geom->GetFTPosition(1)+
223                           geom->GetFTPosition(2)+
224                           geom->GetCPVTextoliteThickness()
225                           )/2;
226         Float_t xoL,yoL,zoL ;
227 //      xoL = xgen*TMath::Cos(phig)+ygen*TMath::Sin(phig) ;
228 //      yoL = -xgen*TMath::Sin(phig)+ygen*TMath::Cos(phig) + distanceIPtoCPV;
229         xoL = xgen*TMath::Cos(phi)-ygen*TMath::Sin(phi) ;
230         yoL = xgen*TMath::Sin(phi)+ygen*TMath::Cos(phi) + distanceIPtoCPV;
231         zoL = zgen;
232
233         r2 = TMath::Power((xoL-xrec),2) + TMath::Power((zoL-zrec),2);
234         if ( r2 < r2min ) {
235           r2min = r2;
236           dxmin = xoL - xrec;
237           dzmin = zoL - zrec;
238           locImpX = xoL;
239           locImpZ = zoL;
240           gImpX = xgen;
241           gImpZ = zgen;
242           gImpY = ygen;
243         }
244       }
245       AliInfo(Form("Impact global (X,Z,Y) = %f %f %f", gImpX, gImpZ, gImpY));
246       AliInfo(Form("Impact local (X,Z) = %f %f", locImpX, locImpZ));
247       AliInfo(Form("Reconstructed (X,Z) = %f %f", xrec, zrec));
248       AliInfo(Form("dxmin %f dzmin %f", dxmin, dzmin));
249       hDx  ->Fill(dxmin);
250       hDz  ->Fill(dzmin);
251 //        hDr  ->Fill(TMath::Sqrt(r2min));
252       hNrp ->Fill(rpMult);
253       hNrpX->Fill(rpMultX);
254       hNrpZ->Fill(rpMultZ);
255     }
256     delete [] hitsPerModule;
257
258     AliInfo(Form("++++ Event %d : total %d impacts, %d charged impacts and %d  rec. points.", 
259           ievent, nTotalGen, nChargedGen, please->CpvRecPoints()->GetEntries())) ;
260   }
261   // Save histograms
262
263   Text_t outputname[80] ;
264   sprintf(outputname,"%s.analyzed",GetFileName().Data());
265   TFile output(outputname,"RECREATE");
266   output.cd();
267   
268   // Plot histograms
269
270   TCanvas *cpvCanvas = new TCanvas("Cpv1","CPV analysis-I",20,20,800,600);
271   gStyle->SetOptStat(111111);
272   gStyle->SetOptFit(1);
273   gStyle->SetOptDate(1);
274   cpvCanvas->Divide(3,3);
275
276   cpvCanvas->cd(1);
277   gPad->SetFillColor(10);
278   hNrp->SetFillColor(16);
279   hNrp->Draw();
280
281   cpvCanvas->cd(2);
282   gPad->SetFillColor(10);
283   hNrpX->SetFillColor(16);
284   hNrpX->Draw();
285
286   cpvCanvas->cd(3);
287   gPad->SetFillColor(10);
288   hNrpZ->SetFillColor(16);
289   hNrpZ->Draw();
290
291   cpvCanvas->cd(4);
292   gPad->SetFillColor(10);
293   hDx->SetFillColor(16);
294   hDx->Fit("gaus");
295   hDx->Draw();
296
297   cpvCanvas->cd(5);
298   gPad->SetFillColor(10);
299   hDz->SetFillColor(16);
300   hDz->Fit("gaus");
301   hDz->Draw();
302
303   cpvCanvas->cd(6);
304   gPad->SetFillColor(10);
305   hChi2->SetFillColor(16);
306   hChi2->Draw();
307
308   cpvCanvas->cd(7);
309   gPad->SetFillColor(10);
310   hEg->SetFillColor(16);
311   hEg->Draw();
312
313   cpvCanvas->cd(8);
314   gPad->SetFillColor(10);
315   hEr->SetFillColor(16);
316   hEr->Draw();
317
318   cpvCanvas->Write(0,kOverwrite);
319
320 }
321
322
323 void AliPHOSIhepAnalyze::AnalyzeEMC1(Int_t Nevents)
324 {
325   //
326   // Analyzes Emc characteristics: resolutions, cluster multiplicity,
327   // cluster lengths along Z and Phi.
328   // Author: Boris Polichtchouk, 24.08.2001
329   //
330
331   // Book histograms
332
333   TH1F *hDx   = new TH1F("hDx"  ,"EMC x-resolution@reconstruction",100,-5. , 5.);
334   TH1F *hDz   = new TH1F("hDz"  ,"EMC z-resolution@reconstruction",100,-5. , 5.);
335   TH1F *hChi2   = new TH1F("hChi2"  ,"Chi2/dof of one-gamma fit",30, 0. , 3.);
336   TH1S *hNrp  = new TH1S("hNrp" ,"EMC rec.point multiplicity",      21,-0.5,20.5);
337   TH1S *hNrpX = new TH1S("hNrpX","EMC rec.point Phi-length"  ,      21,-0.5,20.5);
338   TH1S *hNrpZ = new TH1S("hNrpZ","EMC rec.point Z-length"    ,      21,-0.5,20.5);
339
340   TList * fEmcImpacts ;
341   TBranch * branchEMCimpacts;
342
343   AliPHOSLoader* please = dynamic_cast<AliPHOSLoader*>(fRunLoader->GetLoader("PHOSLoader"));
344   if ( please == 0 )
345    {
346      AliError(Form("Could not obtain the Loader object !"));
347      return ;
348    }
349
350   const AliPHOSGeometry *  fGeom  = please->PHOSGeometry();
351
352   AliInfo(Form("Start EMC Analysis-1. Resolutions, cluster multiplicity and lengths"));
353   for ( Int_t ievent=0; ievent<Nevents; ievent++) {  
354     
355     Int_t nTotalGen = 0;
356
357     Int_t ntracks = gAlice->GetEvent(ievent);
358
359     AliInfo(Form(" >>>>>>>Event %d .<<<<<<<", ievent)) ;
360
361     TTree* treeH = please->TreeH();
362     if (treeH == 0x0)
363      {
364       AliError(Form("Can not get TreeH"));
365        return;
366      }
367
368     
369     // Get branch of EMC impacts
370     if (! (branchEMCimpacts =treeH->GetBranch("PHOSEmcImpacts")) ) {
371       AliWarning(Form(" Couldn't find branch PHOSEmcImpacts. Exit."));
372       return;
373     }
374  
375     // Create and fill arrays of hits for each EMC module
376       
377     Int_t nOfModules = fGeom->GetNModules();
378     TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
379     Int_t iModule = 0;  
380     for (iModule=0; iModule < nOfModules; iModule++)
381       hitsPerModule[iModule] = new TClonesArray("AliPHOSImpact",100);
382
383     TClonesArray    *impacts;
384     AliPHOSImpact   *impact;
385     TLorentzVector   p;
386
387     // First go through all primary tracks and fill the arrays
388     // of hits per each EMC module
389
390     for (Int_t itrack=0; itrack<ntracks; itrack++) {
391       branchEMCimpacts ->SetAddress(&fEmcImpacts);
392       branchEMCimpacts ->GetEntry(itrack,0);
393
394       for (Int_t iModule=0; iModule < nOfModules; iModule++) {
395         impacts = (TClonesArray *)fEmcImpacts->At(iModule);
396         // Do loop over impacts in the module
397         for (Int_t iImpact=0; iImpact<impacts->GetEntries(); iImpact++) {
398           impact=(AliPHOSImpact*)impacts->At(iImpact);
399           TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
400           new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSImpact(*impact);
401         }
402       }
403       fEmcImpacts->Clear();
404     }
405     for (iModule=0; iModule < nOfModules; iModule++) {
406       Int_t nsum = hitsPerModule[iModule]->GetEntriesFast();
407       printf("EMC module %d has %d hits\n",iModule,nsum);
408       nTotalGen += nsum;
409     }
410
411     // Then go through reconstructed points and for each find
412     // the closeset hit
413     // The distance from the rec.point to the closest hit
414     // gives the coordinate resolution of the EMC
415
416     fRunLoader->GetEvent(ievent);
417     TIter nextRP(please->EmcRecPoints()) ;
418     AliPHOSEmcRecPoint *emcRecPoint ;
419     Float_t xgen, ygen, zgen;
420     while( ( emcRecPoint = (AliPHOSEmcRecPoint *)nextRP() ) ) {
421       
422       Float_t chi2dof = ((AliPHOSEvalRecPoint*)emcRecPoint)->Chi2Dof();
423       hChi2->Fill(chi2dof);
424       
425       TVector3  locpos;
426       emcRecPoint->GetLocalPosition(locpos);
427       Int_t phosModule = emcRecPoint->GetPHOSMod();
428       Int_t rpMult     = emcRecPoint->GetMultiplicity();
429       Int_t rpMultX, rpMultZ;
430       ((AliPHOSEvalRecPoint*)emcRecPoint)->GetClusterLengths(rpMultX,rpMultZ);
431       Float_t xrec  = locpos.X();
432       Float_t zrec  = locpos.Z();
433       Float_t dxmin = 1.e+10;
434       Float_t dzmin = 1.e+10;
435       Float_t r2min = 1.e+10;
436       Float_t r2;
437
438       Int_t nEMChits  = (hitsPerModule[phosModule-1])->GetEntriesFast();
439       Float_t locImpX=1e10,locImpZ=1e10;         // local coords of closest impact
440       Float_t gImpX=1e10, gImpZ=1e10,gImpY=1e10; // global coords of closest impact
441       for (Int_t ihit=0; ihit<nEMChits; ihit++) {
442         impact = (AliPHOSImpact*)(hitsPerModule[phosModule-1])->UncheckedAt(ihit);
443         xgen   = impact->X();
444         zgen   = impact->Z();
445         ygen   = impact->Y();
446       
447         
448         //Transform to the local ref.frame
449         const AliPHOSGeometry* geom = please->PHOSGeometry();
450         Float_t phig = geom->GetPHOSAngle(phosModule);
451         Float_t phi = TMath::Pi()/180*phig;
452         Float_t distanceIPtoEMC = geom->GetIPtoCrystalSurface();
453         Float_t xoL,yoL,zoL ;
454 //      xoL = xgen*TMath::Cos(phig)+ygen*TMath::Sin(phig) ;
455 //      yoL = -xgen*TMath::Sin(phig)+ygen*TMath::Cos(phig) + distanceIPtoEMC;
456         xoL = xgen*TMath::Cos(phi)-ygen*TMath::Sin(phi) ;
457         yoL = xgen*TMath::Sin(phi)+ygen*TMath::Cos(phi) + distanceIPtoEMC;
458         zoL = zgen;
459
460         r2 = TMath::Power((xoL-xrec),2) + TMath::Power((zoL-zrec),2);
461         if ( r2 < r2min ) {
462           r2min = r2;
463           dxmin = xoL - xrec;
464           dzmin = zoL - zrec;
465           locImpX = xoL;
466           locImpZ = zoL;
467           gImpX = xgen;
468           gImpZ = zgen;
469           gImpY = ygen;
470         }
471       }
472       AliInfo(Form(" Impact global (X,Z,Y) = %f %f %f", gImpX, gImpZ, gImpY));
473       AliInfo(Form(" Impact local (X,Z) = %f %f", locImpX, locImpZ));
474       AliInfo(Form(" Reconstructed (X,Z) = %f %f", xrec, zrec));
475       AliInfo(Form(" dxmin %f dzmin %f", dxmin, dzmin)) ;
476       hDx  ->Fill(dxmin);
477       hDz  ->Fill(dzmin);
478 //        hDr  ->Fill(TMath::Sqrt(r2min));
479       hNrp ->Fill(rpMult);
480       hNrpX->Fill(rpMultX);
481       hNrpZ->Fill(rpMultZ);
482     }
483     delete [] hitsPerModule;
484
485     AliInfo(Form("++++ Event %d : total  %d impacts,  %d Emc rec. points.", 
486          ievent, nTotalGen, please->EmcRecPoints()->GetEntriesFast())) ; 
487
488   }
489   // Save histograms
490
491   Text_t outputname[80] ;
492   sprintf(outputname,"%s.analyzed",GetFileName().Data());
493   TFile output(outputname,"update");
494   output.cd();
495   
496   // Plot histograms
497
498   TCanvas *emcCanvas = new TCanvas("Emc1","EMC analysis-I",20,20,800,400);
499   gStyle->SetOptStat(111111);
500   gStyle->SetOptFit(1);
501   gStyle->SetOptDate(1);
502   emcCanvas->Divide(3,2);
503
504   emcCanvas->cd(1);
505   gPad->SetFillColor(10);
506   hNrp->SetFillColor(16);
507   hNrp->Draw();
508
509   emcCanvas->cd(2);
510   gPad->SetFillColor(10);
511   hNrpX->SetFillColor(16);
512   hNrpX->Draw();
513
514   emcCanvas->cd(3);
515   gPad->SetFillColor(10);
516   hNrpZ->SetFillColor(16);
517   hNrpZ->Draw();
518
519   emcCanvas->cd(4);
520   gPad->SetFillColor(10);
521   hDx->SetFillColor(16);
522   hDx->Fit("gaus");
523   hDx->Draw();
524
525   emcCanvas->cd(5);
526   gPad->SetFillColor(10);
527   hDz->SetFillColor(16);
528   hDz->Fit("gaus");
529   hDz->Draw();
530
531   emcCanvas->cd(6);
532   gPad->SetFillColor(10);
533   hChi2->SetFillColor(16);
534   hChi2->Draw();
535
536   emcCanvas->Write(0,kOverwrite);
537 }
538
539 //____________________________________________________________________________
540 void AliPHOSIhepAnalyze::AnalyzeCPV2(Int_t Nevents)
541 {
542   // CPV analysis - part II.
543   // Ratio of combinatoric distances between generated
544   // and reconstructed hits.
545   // Author: Boris Polichtchouk (polishchuk@mx.ihep.su)
546   // 24 March 2001
547
548
549   TH1F* hDrijCPVr = new TH1F("Drij_cpv_r","Distance between reconstructed hits in CPV",140,0,50);
550   TH1F* hDrijCPVg = new TH1F("Drij_cpv_g","Distance between generated hits in CPV",140,0,50);
551   TH1F* hDrijCPVratio = new TH1F("Drij_cpv_ratio","R_{ij}^{rec}/R_{ij}^{gen} in CPV",140,0,50);
552
553 //    TH1F* hT0 = new TH1F("hT0","Type of entering particle",20000,-10000,10000);
554
555   hDrijCPVr->Sumw2();
556   hDrijCPVg->Sumw2();
557   hDrijCPVratio->Sumw2(); //correct treatment of errors
558
559   TList * fCpvImpacts = new TList();
560   TBranch * branchCPVimpacts;
561
562   AliPHOSLoader* please = dynamic_cast<AliPHOSLoader*>(fRunLoader->GetLoader("PHOSLoader"));
563   if ( please == 0 )
564    {
565      AliError(Form("Could not obtain the Loader object !"));
566      return ;
567    }
568   const AliPHOSGeometry *  fGeom  = please->PHOSGeometry();
569   fRunLoader->LoadHeader();
570
571   for (Int_t nev=0; nev<Nevents; nev++) 
572     { 
573       printf("\n=================== Event %10d ===================\n",nev);
574       fRunLoader->GetEvent(nev);
575       Int_t ntracks = fRunLoader->GetHeader()->GetNtrack();
576     
577       Int_t nRecCPV = 0; // Reconstructed points in event
578       Int_t nGenCPV = 0; // Impacts in event
579
580       // Get branch of CPV impacts
581       TTree* treeH = please->TreeH();
582       if (treeH == 0x0)
583        {
584          AliError(Form("Can not get TreeH"));
585         return;
586        }
587
588       if (! (branchCPVimpacts =treeH->GetBranch("PHOSCpvImpacts")) )  return;
589       
590       // Create and fill arrays of hits for each CPV module
591       Int_t nOfModules = fGeom->GetNModules();
592       TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
593       Int_t iModule = 0;        
594       for (iModule=0; iModule < nOfModules; iModule++)
595         hitsPerModule[iModule] = new TClonesArray("AliPHOSImpact",100);
596
597       TClonesArray    *impacts;
598       AliPHOSImpact   *impact;
599           
600       for (Int_t itrack=0; itrack<ntracks; itrack++) {
601         branchCPVimpacts ->SetAddress(&fCpvImpacts);
602         AliInfo(Form(" branchCPVimpacts ->SetAddress(&fCpvImpacts) OK."));
603         branchCPVimpacts ->GetEntry(itrack,0);
604
605         for (Int_t iModule=0; iModule < nOfModules; iModule++) {
606           impacts = (TClonesArray *)fCpvImpacts->At(iModule);
607           // Do loop over impacts in the module
608           for (Int_t iImpact=0; iImpact<impacts->GetEntries(); iImpact++) {
609             impact=(AliPHOSImpact*)impacts->At(iImpact);
610             TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
611             if(IsCharged(impact->GetPid()))
612               new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSImpact(*impact);
613           }
614         }
615         fCpvImpacts->Clear();
616       }
617
618       for (iModule=0; iModule < nOfModules; iModule++) {
619         Int_t nsum = hitsPerModule[iModule]->GetEntriesFast();
620         printf("CPV module %d has %d hits\n",iModule,nsum);
621
622         AliPHOSImpact* genHit1;
623         AliPHOSImpact* genHit2;
624         Int_t irp1,irp2;
625         for(irp1=0; irp1< nsum; irp1++)
626           {
627             genHit1 = (AliPHOSImpact*)((hitsPerModule[iModule])->At(irp1));
628             for(irp2 = irp1+1; irp2<nsum; irp2++)
629               {
630                 genHit2 = (AliPHOSImpact*)((hitsPerModule[iModule])->At(irp2));
631                 Float_t dx = genHit1->X() - genHit2->X();
632                 Float_t dz = genHit1->Z() - genHit2->Z();
633                 Float_t dr = TMath::Sqrt(dx*dx + dz*dz);
634                 hDrijCPVg->Fill(dr);
635 //                      AliInfo(Form("(dx dz dr): %f %f", dx, dz));
636               }
637           }
638       }
639
640  
641   //--------- Combinatoric distance between rec. hits in CPV
642
643       TObjArray* cpvRecPoints = please->CpvRecPoints();
644       nRecCPV =  cpvRecPoints->GetEntriesFast();
645
646       if(nRecCPV)
647         {
648           AliPHOSCpvRecPoint* recHit1;
649           AliPHOSCpvRecPoint* recHit2;
650           TIter nextCPVrec1(cpvRecPoints);
651           while(TObject* obj1 = nextCPVrec1() )
652             {
653               TIter nextCPVrec2(cpvRecPoints);
654               while (TObject* obj2 = nextCPVrec2())
655                 {
656                   if(!obj2->IsEqual(obj1))
657                     {
658                       recHit1 = (AliPHOSCpvRecPoint*)obj1;
659                       recHit2 = (AliPHOSCpvRecPoint*)obj2;
660                       TVector3 locpos1;
661                       TVector3 locpos2;
662                       recHit1->GetLocalPosition(locpos1);
663                       recHit2->GetLocalPosition(locpos2);
664                       Float_t dx = locpos1.X() - locpos2.X();
665                       Float_t dz = locpos1.Z() - locpos2.Z();                 
666                       Float_t dr = TMath::Sqrt(dx*dx + dz*dz);
667                       if(recHit1->GetPHOSMod() == recHit2->GetPHOSMod())
668                         hDrijCPVr->Fill(dr);
669                     }
670                 }
671             }   
672         }
673       
674       AliInfo(Form(" Event %d . Total of %d hits, %d rec.points.", 
675            nev, nGenCPV, nRecCPV)) ; 
676     
677       delete [] hitsPerModule;
678
679     } // End of loop over events.
680
681
682 //    hDrijCPVg->Draw();
683 //    hDrijCPVr->Draw();
684   hDrijCPVratio->Divide(hDrijCPVr,hDrijCPVg);
685   hDrijCPVratio->Draw();
686
687 //    hT0->Draw();
688
689 }
690
691
692 void AliPHOSIhepAnalyze::CpvSingle(Int_t nevents)
693 {
694   // Distributions of coordinates and cluster lengths of rec. points
695   // in the case of single initial particle.
696
697   TH1F* hZr = new TH1F("Zrec","Reconstructed Z distribution",150,-5,5);
698   TH1F* hXr = new TH1F("Xrec","Reconstructed X distribution",150,-14,-2);
699   TH1F *hChi2   = new TH1F("hChi2"  ,"Chi2/dof of one-gamma fit",100, 0. , 20.);
700
701   TH1S *hNrp  = new TH1S("hNrp" ,"CPV rec.point multiplicity",21,-0.5,20.5);
702   TH1S *hNrpX = new TH1S("hNrpX","CPV rec.point Phi-length"  ,21,-0.5,20.5);
703   TH1S *hNrpZ = new TH1S("hNrpZ","CPV rec.point Z-length"    ,21,-0.5,20.5);
704  
705   AliPHOSLoader* gime = dynamic_cast<AliPHOSLoader*>(fRunLoader->GetLoader("PHOSLoader"));
706   if ( gime == 0 )
707    {
708      AliError(Form("Could not obtain the Loader object !"));
709      return ;
710    }
711   
712   for(Int_t ievent=0; ievent<nevents; ievent++)
713     {
714       fRunLoader->GetEvent(ievent);
715       if(gime->CpvRecPoints()->GetEntriesFast()>1) continue;
716
717       AliPHOSCpvRecPoint* pt = (AliPHOSCpvRecPoint*)(gime->CpvRecPoints())->At(0);
718       if(pt) {
719         TVector3 lpos;
720         pt->GetLocalPosition(lpos);
721         hXr->Fill(lpos.X());
722         hZr->Fill(lpos.Z());
723
724         Int_t rpMult = pt->GetMultiplicity();
725         hNrp->Fill(rpMult);
726         Int_t rpMultX, rpMultZ;
727         pt->GetClusterLengths(rpMultX,rpMultZ);
728         hNrpX->Fill(rpMultX);
729         hNrpZ->Fill(rpMultZ);
730         hChi2->Fill(((AliPHOSEvalRecPoint*)pt)->Chi2Dof());
731         AliInfo(Form("+++++ Event %d . (Mult,MultX,MultZ) = %d %d %d +++++", 
732              ievent, rpMult, rpMultX, rpMultZ)) ;
733
734       }
735
736     }
737         
738   Text_t outputname[80] ;
739   sprintf(outputname,"%s.analyzed.single",GetFileName().Data());
740   TFile output(outputname,"RECREATE");
741   output.cd();
742     
743   // Plot histograms
744   TCanvas *cpvCanvas = new TCanvas("SingleParticle","Single particle events",20,20,800,400);
745   gStyle->SetOptStat(111111);
746   gStyle->SetOptFit(1);
747   gStyle->SetOptDate(1);
748   cpvCanvas->Divide(3,2);
749
750   cpvCanvas->cd(1);
751   gPad->SetFillColor(10);
752   hXr->SetFillColor(16);
753   hXr->Draw();
754
755   cpvCanvas->cd(2);
756   gPad->SetFillColor(10);
757   hZr->SetFillColor(16);
758   hZr->Draw();
759
760   cpvCanvas->cd(3);
761   gPad->SetFillColor(10);
762   hChi2->SetFillColor(16);
763   hChi2->Draw();
764
765   cpvCanvas->cd(4);
766   gPad->SetFillColor(10);
767   hNrp->SetFillColor(16);
768   hNrp->Draw();
769
770   cpvCanvas->cd(5);
771   gPad->SetFillColor(10);
772   hNrpX->SetFillColor(16);
773   hNrpX->Draw();
774
775   cpvCanvas->cd(6);
776   gPad->SetFillColor(10);
777   hNrpZ->SetFillColor(16);
778   hNrpZ->Draw();
779
780   cpvCanvas->Write(0,kOverwrite);
781   
782 }
783
784 void AliPHOSIhepAnalyze::HitsCPV(Int_t nev)
785 {
786   // Cumulative list of charged CPV impacts in event nev.
787
788   TList * fCpvImpacts ;
789   TBranch * branchCPVimpacts;
790
791   AliPHOSLoader* please = dynamic_cast<AliPHOSLoader*>(fRunLoader->GetLoader("PHOSLoader"));
792   if ( please == 0 )
793    {
794      AliError(Form("Could not obtain the Loader object !"));
795      return ;
796    }
797   const AliPHOSGeometry *  fGeom  = please->PHOSGeometry();
798
799      
800   printf("\n=================== Event %10d ===================\n",nev);
801   fRunLoader->GetEvent(nev);
802   Int_t ntracks = fRunLoader->GetHeader()->GetNtrack();
803     
804 //    Int_t nRecCPV = 0; // Reconstructed points in event // 01.10.2001
805 //    Int_t nGenCPV = 0; // Impacts in event
806
807   // Get branch of CPV impacts
808    TTree* treeH = please->TreeH();
809    if (treeH == 0x0)
810     {
811       AliError(Form("Can not get TreeH"));
812       return;
813     }
814
815   if (! (branchCPVimpacts =treeH->GetBranch("PHOSCpvImpacts")) )  return;
816       
817   // Create and fill arrays of hits for each CPV module
818   Int_t nOfModules = fGeom->GetNModules();
819   TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
820   Int_t iModule = 0;    
821   for (iModule=0; iModule < nOfModules; iModule++)
822     hitsPerModule[iModule] = new TClonesArray("AliPHOSImpact",100);
823   
824   TClonesArray    *impacts;
825   AliPHOSImpact   *impact;
826           
827   for (Int_t itrack=0; itrack<ntracks; itrack++) {
828     branchCPVimpacts ->SetAddress(&fCpvImpacts);
829     AliInfo(Form(" branchCPVimpacts ->SetAddress(&fCpvImpacts) OK."));
830     branchCPVimpacts ->GetEntry(itrack,0);
831
832     for (Int_t iModule=0; iModule < nOfModules; iModule++) {
833       impacts = (TClonesArray *)fCpvImpacts->At(iModule);
834       // Do loop over impacts in the module
835       for (Int_t iImpact=0; iImpact<impacts->GetEntries(); iImpact++) {
836         impact=(AliPHOSImpact*)impacts->At(iImpact);
837         TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
838         new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSImpact(*impact);
839       }
840     }
841     fCpvImpacts->Clear();
842   }
843
844   for (iModule=0; iModule < nOfModules; iModule++) {
845     Int_t nsum = hitsPerModule[iModule]->GetEntriesFast();
846     printf("CPV module %d has %d hits\n",iModule,nsum);
847   }
848
849 //    TList * fCpvImpacts ;
850 //    TBranch * branchCPVimpacts;
851 //    AliPHOSImpact* impact;
852 //    TClonesArray* impacts;
853
854 //    AliPHOSLoader * please = AliPHOSLoader::GetInstance(GetFileName().Data(),"PHOS");
855 //    const AliPHOSGeometry *  fGeom  = please->PHOSGeometry();
856
857 //    Int_t ntracks = gAlice->GetEvent(ievent);
858 //    Int_t nOfModules = fGeom->GetNModules();
859 //    AliInfo(Form(" Tracks: "<<ntracks<<"  Modules: "<<nOfModules));
860
861 //    if (! (branchCPVimpacts =gAlice->TreeH()->GetBranch("PHOSCpvImpacts")) )  return;
862
863 //    for (Int_t itrack=0; itrack<ntracks; itrack++) {
864 //      branchCPVimpacts ->SetAddress(&fCpvImpacts);
865 //      Info("AnalyzeCPV1", " branchCPVimpacts ->SetAddress(&fCpvImpacts) OK.");
866 //      branchCPVimpacts ->GetEntry(itrack,0);
867 //      Info(Form(" branchCPVimpacts ->GetEntry(itrack,0) OK."));
868     
869 //      for (Int_t iModule=0; iModule < nOfModules; iModule++) {
870 //        impacts = (TClonesArray *)fCpvImpacts->At(iModule);
871 //        Info(Form(" fCpvImpacts->At(iModule) OK."));
872 //        // Do loop over impacts in the module
873 //        for (Int_t iImpact=0; iImpact<impacts->GetEntries(); iImpact++) {
874 //      impact=(AliPHOSImpact*)impacts->At(iImpact);
875 //      impact->Print();
876 //      if(IsCharged(impact->GetPid()))
877 //        {
878 //          Info(Form(" Add charged hit.."));
879 //          new(hits[hits.GetEntriesFast()]) AliPHOSImpact(*impact);
880 //          Info(Form("done."));
881 //        }
882 //        }
883 //      }
884 //      fCpvImpacts->Clear();
885 //    }
886
887 //    Info(Form(" PHOS event "<<ievent<<": "<<hits.GetEntries()<<" charged CPV hits."));
888
889 }
890
891
892 //  void AliPHOSIhepAnalyze::ChargedHitsCPV(TClonesArray* hits, Int_t ievent, Int_t iModule)
893 //  {
894 //    // Cumulative list of charged CPV hits in event ievent 
895 //    // in PHOS module iModule.
896
897 //    HitsCPV(hits,ievent,iModule);
898 //    TIter next(hits);
899
900 //    while(AliPHOSCPVHit* cpvhit = (AliPHOSCPVHit*)next())
901 //      {
902 //        if(!IsCharged(cpvhit->GetIpart()))
903 //      {
904 //        hits->Remove(cpvhit);
905 //        delete cpvhit;
906 //        hits->Compress();
907 //      }
908 //      }
909
910 //    Info(Form(" PHOS module "<<iModule<<": "<<hits->GetEntries()<<" charged CPV hits."));
911 //  }
912
913 Bool_t AliPHOSIhepAnalyze::IsCharged(Int_t pdgCode)
914 {
915   // For HIJING
916   AliInfo(Form("pdgCode %d", pdgCode));
917   if(pdgCode==211 || pdgCode==-211 || pdgCode==321 || pdgCode==-321 || pdgCode==11 || pdgCode==-11 || pdgCode==2212 || pdgCode==-2212) return kTRUE;
918   else
919     return kFALSE;
920 }
921 //---------------------------------------------------------------------------
922
923
924
925
926
927