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