Transition to NewIO
[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 "AliPHOSIhepAnalyze.h"
47 #include "AliPHOSDigit.h"
48 #include "AliPHOSRecParticle.h"
49 #include "AliPHOSLoader.h"
50 #include "AliPHOSHit.h"
51 #include "AliPHOSImpact.h"
52 #include "AliPHOSvImpacts.h"
53 #include "AliPHOSCpvRecPoint.h"
54 #include "AliRun.h"
55 #include "AliPHOSGeometry.h"
56 #include "AliPHOSEvalRecPoint.h"
57
58 ClassImp(AliPHOSIhepAnalyze)
59
60 //____________________________________________________________________________
61
62 AliPHOSIhepAnalyze::AliPHOSIhepAnalyze() 
63  {
64    fRunLoader = 0x0;
65  }
66
67 //____________________________________________________________________________
68
69 AliPHOSIhepAnalyze::AliPHOSIhepAnalyze(Text_t * name) : fFileName(name) 
70  {
71   fRunLoader = AliRunLoader::Open(fFileName);
72   if (fRunLoader == 0x0)
73    {
74      Fatal("AliPHOSIhepAnalyze","Can not load event from file %s",name);
75    }
76  }
77
78 //____________________________________________________________________________
79 void AliPHOSIhepAnalyze::AnalyzeCPV1(Int_t Nevents)
80 {
81   //
82   // Analyzes CPV characteristics: resolutions, cluster multiplicity,
83   // cluster lengths along Z and Phi.
84   // Author: Yuri Kharlov
85   // 9 October 2000
86   // Modified by Boris Polichtchouk, 3.07.2001
87   //
88
89   // Book histograms
90
91   TH1F *hDx   = new TH1F("hDx"  ,"CPV x-resolution@reconstruction",100,-5. , 5.);
92   TH1F *hDz   = new TH1F("hDz"  ,"CPV z-resolution@reconstruction",100,-5. , 5.);
93   TH1F *hChi2 = new TH1F("hChi2"  ,"Chi2/dof of one-gamma fit",30, 0. , 10.);
94   TH1S *hNrp  = new TH1S("hNrp" ,"CPV rec.point multiplicity",      21,-0.5,20.5);
95   TH1S *hNrpX = new TH1S("hNrpX","CPV rec.point Phi-length"  ,      21,-0.5,20.5);
96   TH1S *hNrpZ = new TH1S("hNrpZ","CPV rec.point Z-length"    ,      21,-0.5,20.5);
97
98   TH1F *hEg   = new TH1F("hEg"  ,"Energies of impacts",30,0.,6.);
99   TH1F *hEr   = new TH1F("hEr"  ,"Amplitudes of rec. points",50,0.,20.);
100   
101   TList * fCpvImpacts ;
102   TBranch * branchCPVimpacts;
103
104   
105   
106   AliPHOSLoader* please = dynamic_cast<AliPHOSLoader*>(fRunLoader->GetLoader("PHOSLoader"));
107   if ( please == 0 )
108    {
109      Error("AnalyzeCPV1","Could not obtain the Loader object !");
110      return ;
111    }
112
113   const AliPHOSGeometry *  fGeom  = please->PHOSGeometry();
114
115   Info("AnalyzeCPV1", "Start CPV Analysis-1. Resolutions, cluster multiplicity and lengths") ;
116   for ( Int_t ievent=0; ievent<Nevents; ievent++) {  
117     
118     Int_t nTotalGen = 0;
119     Int_t nChargedGen = 0;
120
121     Int_t ntracks = gAlice->GetEvent(ievent);
122     Info("AnalyzeCPV1", ">>>>>>>Event %d .<<<<<<<", ievent) ;
123     
124   /******************************************************************/
125       TTree* treeH = please->TreeH();
126       if (treeH == 0x0)
127        {
128          Error("AnalyzeCPV1","Can not get TreeH");
129          return;
130        }
131 /******************************************************************/     
132
133     // Get branch of CPV impacts
134     if (! (branchCPVimpacts =treeH->GetBranch("PHOSCpvImpacts")) ) {
135       Info("AnalyzeCPV1", "Couldn't find branch PHOSCpvImpacts. Exit.") ;
136       return;
137     }
138  
139     // Create and fill arrays of hits for each CPV module
140       
141     Int_t nOfModules = fGeom->GetNModules();
142     TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
143     Int_t iModule = 0;  
144     for (iModule=0; iModule < nOfModules; iModule++)
145       hitsPerModule[iModule] = new TClonesArray("AliPHOSImpact",100);
146
147     TClonesArray    *impacts;
148     AliPHOSImpact   *impact;
149     TLorentzVector   p;
150
151     // First go through all primary tracks and fill the arrays
152     // of hits per each CPV module
153
154     for (Int_t itrack=0; itrack<ntracks; itrack++) {
155       branchCPVimpacts ->SetAddress(&fCpvImpacts);
156       branchCPVimpacts ->GetEntry(itrack,0);
157
158       for (Int_t iModule=0; iModule < nOfModules; iModule++) {
159         impacts = (TClonesArray *)fCpvImpacts->At(iModule);
160         // Do loop over impacts in the module
161         for (Int_t iImpact=0; iImpact<impacts->GetEntries(); iImpact++) {
162           impact=(AliPHOSImpact*)impacts->At(iImpact);
163           hEg->Fill(impact->GetMomentum().E());
164           TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
165           if(IsCharged(impact->GetPid())) nChargedGen++;
166           new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSImpact(*impact);
167         }
168       }
169       fCpvImpacts->Clear();
170     }
171     for (iModule=0; iModule < nOfModules; iModule++) {
172       Int_t nsum = hitsPerModule[iModule]->GetEntriesFast();
173       printf("CPV module %d has %d impacts\n",iModule,nsum);
174       nTotalGen += nsum;
175     }
176
177     // Then go through reconstructed points and for each find
178     // the closeset hit
179     // The distance from the rec.point to the closest hit
180     // gives the coordinate resolution of the CPV
181
182     fRunLoader->GetEvent(ievent);
183     TIter nextRP(please->CpvRecPoints()) ;
184     AliPHOSCpvRecPoint *cpvRecPoint ;
185     Float_t xgen, ygen, zgen;
186     while( ( cpvRecPoint = (AliPHOSCpvRecPoint *)nextRP() ) ) {
187       
188       Float_t chi2dof = ((AliPHOSEvalRecPoint*)cpvRecPoint)->Chi2Dof();
189       hChi2->Fill(chi2dof);
190       hEr->Fill(cpvRecPoint->GetEnergy());
191
192       TVector3  locpos;
193       cpvRecPoint->GetLocalPosition(locpos);
194       Int_t phosModule = cpvRecPoint->GetPHOSMod();
195       Int_t rpMult     = cpvRecPoint->GetMultiplicity();
196       Int_t rpMultX, rpMultZ;
197       cpvRecPoint->GetClusterLengths(rpMultX,rpMultZ);
198       Float_t xrec  = locpos.X();
199       Float_t zrec  = locpos.Z();
200       Float_t dxmin = 1.e+10;
201       Float_t dzmin = 1.e+10;
202       Float_t r2min = 1.e+10;
203       Float_t r2;
204
205       Int_t nCPVhits  = (hitsPerModule[phosModule-1])->GetEntriesFast();
206       Float_t locImpX=1e10,locImpZ=1e10;         // local coords of closest impact
207       Float_t gImpX=1e10, gImpZ=1e10,gImpY=1e10; // global coords of closest impact
208       for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
209         impact = (AliPHOSImpact*)(hitsPerModule[phosModule-1])->UncheckedAt(ihit);
210         xgen   = impact->X();
211         zgen   = impact->Z();
212         ygen   = impact->Y();
213         
214         //Transform to the local ref.frame
215         const AliPHOSGeometry* geom = please->PHOSGeometry();
216         Float_t phig = geom->GetPHOSAngle(phosModule);
217         Float_t phi = TMath::Pi()/180*phig;
218         Float_t distanceIPtoCPV = geom->GetIPtoOuterCoverDistance() -
219                           (geom->GetFTPosition(1)+
220                           geom->GetFTPosition(2)+
221                           geom->GetCPVTextoliteThickness()
222                           )/2;
223         Float_t xoL,yoL,zoL ;
224 //      xoL = xgen*TMath::Cos(phig)+ygen*TMath::Sin(phig) ;
225 //      yoL = -xgen*TMath::Sin(phig)+ygen*TMath::Cos(phig) + distanceIPtoCPV;
226         xoL = xgen*TMath::Cos(phi)-ygen*TMath::Sin(phi) ;
227         yoL = xgen*TMath::Sin(phi)+ygen*TMath::Cos(phi) + distanceIPtoCPV;
228         zoL = zgen;
229
230         r2 = TMath::Power((xoL-xrec),2) + TMath::Power((zoL-zrec),2);
231         if ( r2 < r2min ) {
232           r2min = r2;
233           dxmin = xoL - xrec;
234           dzmin = zoL - zrec;
235           locImpX = xoL;
236           locImpZ = zoL;
237           gImpX = xgen;
238           gImpZ = zgen;
239           gImpY = ygen;
240         }
241       }
242       Info("AnalyzeCPV1", "Impact global (X,Z,Y) = %f %f %f", gImpX, gImpZ, gImpY);
243       Info("AnalyzeCPV1", "Impact local (X,Z) = %f %f", locImpX, locImpZ);
244       Info("AnalyzeCPV1", "Reconstructed (X,Z) = %f %f", xrec, zrec);
245       Info("AnalyzeCPV1", "dxmin %f dzmin %f", dxmin, dzmin);
246       hDx  ->Fill(dxmin);
247       hDz  ->Fill(dzmin);
248 //        hDr  ->Fill(TMath::Sqrt(r2min));
249       hNrp ->Fill(rpMult);
250       hNrpX->Fill(rpMultX);
251       hNrpZ->Fill(rpMultZ);
252     }
253     delete [] hitsPerModule;
254
255     Info("AnalyzeCPV1", "++++ Event %d : total %d impacts, %d charged impacts and %d  rec. points.", 
256           ievent, nTotalGen, nChargedGen, please->CpvRecPoints()->GetEntries()) ;
257   }
258   // Save histograms
259
260   Text_t outputname[80] ;
261   sprintf(outputname,"%s.analyzed",GetFileName().Data());
262   TFile output(outputname,"RECREATE");
263   output.cd();
264   
265   // Plot histograms
266
267   TCanvas *cpvCanvas = new TCanvas("Cpv1","CPV analysis-I",20,20,800,600);
268   gStyle->SetOptStat(111111);
269   gStyle->SetOptFit(1);
270   gStyle->SetOptDate(1);
271   cpvCanvas->Divide(3,3);
272
273   cpvCanvas->cd(1);
274   gPad->SetFillColor(10);
275   hNrp->SetFillColor(16);
276   hNrp->Draw();
277
278   cpvCanvas->cd(2);
279   gPad->SetFillColor(10);
280   hNrpX->SetFillColor(16);
281   hNrpX->Draw();
282
283   cpvCanvas->cd(3);
284   gPad->SetFillColor(10);
285   hNrpZ->SetFillColor(16);
286   hNrpZ->Draw();
287
288   cpvCanvas->cd(4);
289   gPad->SetFillColor(10);
290   hDx->SetFillColor(16);
291   hDx->Fit("gaus");
292   hDx->Draw();
293
294   cpvCanvas->cd(5);
295   gPad->SetFillColor(10);
296   hDz->SetFillColor(16);
297   hDz->Fit("gaus");
298   hDz->Draw();
299
300   cpvCanvas->cd(6);
301   gPad->SetFillColor(10);
302   hChi2->SetFillColor(16);
303   hChi2->Draw();
304
305   cpvCanvas->cd(7);
306   gPad->SetFillColor(10);
307   hEg->SetFillColor(16);
308   hEg->Draw();
309
310   cpvCanvas->cd(8);
311   gPad->SetFillColor(10);
312   hEr->SetFillColor(16);
313   hEr->Draw();
314
315   cpvCanvas->Write(0,kOverwrite);
316
317 }
318
319
320 void AliPHOSIhepAnalyze::AnalyzeEMC1(Int_t Nevents)
321 {
322   //
323   // Analyzes Emc characteristics: resolutions, cluster multiplicity,
324   // cluster lengths along Z and Phi.
325   // Author: Boris Polichtchouk, 24.08.2001
326   //
327
328   // Book histograms
329
330   TH1F *hDx   = new TH1F("hDx"  ,"EMC x-resolution@reconstruction",100,-5. , 5.);
331   TH1F *hDz   = new TH1F("hDz"  ,"EMC z-resolution@reconstruction",100,-5. , 5.);
332   TH1F *hChi2   = new TH1F("hChi2"  ,"Chi2/dof of one-gamma fit",30, 0. , 3.);
333   TH1S *hNrp  = new TH1S("hNrp" ,"EMC rec.point multiplicity",      21,-0.5,20.5);
334   TH1S *hNrpX = new TH1S("hNrpX","EMC rec.point Phi-length"  ,      21,-0.5,20.5);
335   TH1S *hNrpZ = new TH1S("hNrpZ","EMC rec.point Z-length"    ,      21,-0.5,20.5);
336
337   TList * fEmcImpacts ;
338   TBranch * branchEMCimpacts;
339
340   AliPHOSLoader* please = dynamic_cast<AliPHOSLoader*>(fRunLoader->GetLoader("PHOSLoader"));
341   if ( please == 0 )
342    {
343      Error("AnalyzeEMC1","Could not obtain the Loader object !");
344      return ;
345    }
346
347   const AliPHOSGeometry *  fGeom  = please->PHOSGeometry();
348
349   Info("AnalyzeCPV1", "Start EMC Analysis-1. Resolutions, cluster multiplicity and lengths");
350   for ( Int_t ievent=0; ievent<Nevents; ievent++) {  
351     
352     Int_t nTotalGen = 0;
353
354     Int_t ntracks = gAlice->GetEvent(ievent);
355
356     Info("AnalyzeCPV1", " >>>>>>>Event %d .<<<<<<<", ievent) ;
357
358     TTree* treeH = please->TreeH();
359     if (treeH == 0x0)
360      {
361        Error("AnalyzeEMC1","Can not get TreeH");
362        return;
363      }
364
365     
366     // Get branch of EMC impacts
367     if (! (branchEMCimpacts =treeH->GetBranch("PHOSEmcImpacts")) ) {
368       Info("AnalyzeCPV1", " Couldn't find branch PHOSEmcImpacts. Exit.");
369       return;
370     }
371  
372     // Create and fill arrays of hits for each EMC module
373       
374     Int_t nOfModules = fGeom->GetNModules();
375     TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
376     Int_t iModule = 0;  
377     for (iModule=0; iModule < nOfModules; iModule++)
378       hitsPerModule[iModule] = new TClonesArray("AliPHOSImpact",100);
379
380     TClonesArray    *impacts;
381     AliPHOSImpact   *impact;
382     TLorentzVector   p;
383
384     // First go through all primary tracks and fill the arrays
385     // of hits per each EMC module
386
387     for (Int_t itrack=0; itrack<ntracks; itrack++) {
388       branchEMCimpacts ->SetAddress(&fEmcImpacts);
389       branchEMCimpacts ->GetEntry(itrack,0);
390
391       for (Int_t iModule=0; iModule < nOfModules; iModule++) {
392         impacts = (TClonesArray *)fEmcImpacts->At(iModule);
393         // Do loop over impacts in the module
394         for (Int_t iImpact=0; iImpact<impacts->GetEntries(); iImpact++) {
395           impact=(AliPHOSImpact*)impacts->At(iImpact);
396           TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
397           new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSImpact(*impact);
398         }
399       }
400       fEmcImpacts->Clear();
401     }
402     for (iModule=0; iModule < nOfModules; iModule++) {
403       Int_t nsum = hitsPerModule[iModule]->GetEntriesFast();
404       printf("EMC module %d has %d hits\n",iModule,nsum);
405       nTotalGen += nsum;
406     }
407
408     // Then go through reconstructed points and for each find
409     // the closeset hit
410     // The distance from the rec.point to the closest hit
411     // gives the coordinate resolution of the EMC
412
413     fRunLoader->GetEvent(ievent);
414     TIter nextRP(please->EmcRecPoints()) ;
415     AliPHOSEmcRecPoint *emcRecPoint ;
416     Float_t xgen, ygen, zgen;
417     while( ( emcRecPoint = (AliPHOSEmcRecPoint *)nextRP() ) ) {
418       
419       Float_t chi2dof = ((AliPHOSEvalRecPoint*)emcRecPoint)->Chi2Dof();
420       hChi2->Fill(chi2dof);
421       
422       TVector3  locpos;
423       emcRecPoint->GetLocalPosition(locpos);
424       Int_t phosModule = emcRecPoint->GetPHOSMod();
425       Int_t rpMult     = emcRecPoint->GetMultiplicity();
426       Int_t rpMultX, rpMultZ;
427       ((AliPHOSEvalRecPoint*)emcRecPoint)->GetClusterLengths(rpMultX,rpMultZ);
428       Float_t xrec  = locpos.X();
429       Float_t zrec  = locpos.Z();
430       Float_t dxmin = 1.e+10;
431       Float_t dzmin = 1.e+10;
432       Float_t r2min = 1.e+10;
433       Float_t r2;
434
435       Int_t nEMChits  = (hitsPerModule[phosModule-1])->GetEntriesFast();
436       Float_t locImpX=1e10,locImpZ=1e10;         // local coords of closest impact
437       Float_t gImpX=1e10, gImpZ=1e10,gImpY=1e10; // global coords of closest impact
438       for (Int_t ihit=0; ihit<nEMChits; ihit++) {
439         impact = (AliPHOSImpact*)(hitsPerModule[phosModule-1])->UncheckedAt(ihit);
440         xgen   = impact->X();
441         zgen   = impact->Z();
442         ygen   = impact->Y();
443       
444         
445         //Transform to the local ref.frame
446         const AliPHOSGeometry* geom = please->PHOSGeometry();
447         Float_t phig = geom->GetPHOSAngle(phosModule);
448         Float_t phi = TMath::Pi()/180*phig;
449         Float_t distanceIPtoEMC = geom->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       Info("AnalyzeCPV1", " Impact global (X,Z,Y) = %f %f %f", gImpX, gImpZ, gImpY);
470       Info("AnalyzeCPV1", " Impact local (X,Z) = %f %f", locImpX, locImpZ);
471       Info("AnalyzeCPV1", " Reconstructed (X,Z) = %f %f", xrec, zrec);
472       Info("AnalyzeCPV1", " 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     Info("AnalyzeCPV1", "++++ 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* hDrij_cpv_r = new TH1F("Drij_cpv_r","Distance between reconstructed hits in CPV",140,0,50);
547   TH1F* hDrij_cpv_g = new TH1F("Drij_cpv_g","Distance between generated hits in CPV",140,0,50);
548   TH1F* hDrij_cpv_ratio = 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   hDrij_cpv_r->Sumw2();
553   hDrij_cpv_g->Sumw2();
554   hDrij_cpv_ratio->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      Error("AnalyzeCPV2","Could not obtain the Loader object !");
563      return ;
564    }
565   const AliPHOSGeometry *  fGeom  = please->PHOSGeometry();
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 nrec_cpv = 0; // Reconstructed points in event
575       Int_t ngen_cpv = 0; // Impacts in event
576
577       // Get branch of CPV impacts
578       TTree* treeH = please->TreeH();
579       if (treeH == 0x0)
580        {
581         Error("AnalyzeCPV2","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 = fGeom->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         Info("AnalyzeCPV1", " branchCPVimpacts ->SetAddress(&fCpvImpacts) OK.");
600         branchCPVimpacts ->GetEntry(itrack,0);
601
602         for (Int_t 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                 hDrij_cpv_g->Fill(dr);
632 //                      Info("AnalyzeCPV1", "(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       nrec_cpv =  cpvRecPoints->GetEntriesFast();
642
643       if(nrec_cpv)
644         {
645           AliPHOSCpvRecPoint* RecHit1;
646           AliPHOSCpvRecPoint* RecHit2;
647           TIter next_cpv_rec1(cpvRecPoints);
648           while(TObject* obj1 = next_cpv_rec1() )
649             {
650               TIter next_cpv_rec2(cpvRecPoints);
651               while (TObject* obj2 = next_cpv_rec2())
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                         hDrij_cpv_r->Fill(dr);
666                     }
667                 }
668             }   
669         }
670       
671       Info("AnalyzeCPV1", " Event %d . Total of %d hits, %d rec.points.", 
672            nev, ngen_cpv, nrec_cpv) ; 
673     
674       delete [] hitsPerModule;
675
676     } // End of loop over events.
677
678
679 //    hDrij_cpv_g->Draw();
680 //    hDrij_cpv_r->Draw();
681   hDrij_cpv_ratio->Divide(hDrij_cpv_r,hDrij_cpv_g);
682   hDrij_cpv_ratio->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      Error("CpvSingle","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         Info("AnalyzeCPV1", "+++++ 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(TClonesArray& hits, 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      Error("HitsCPV","Could not obtain the Loader object !");
792      return ;
793    }
794   const AliPHOSGeometry *  fGeom  = please->PHOSGeometry();
795
796      
797   printf("\n=================== Event %10d ===================\n",nev);
798   fRunLoader->GetEvent(nev);
799   Int_t ntracks = fRunLoader->GetHeader()->GetNtrack();
800     
801 //    Int_t nrec_cpv = 0; // Reconstructed points in event // 01.10.2001
802 //    Int_t ngen_cpv = 0; // Impacts in event
803
804   // Get branch of CPV impacts
805    TTree* treeH = please->TreeH();
806    if (treeH == 0x0)
807     {
808       Error("CPVSingle","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 = fGeom->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     Info("AnalyzeCPV1", " branchCPVimpacts ->SetAddress(&fCpvImpacts) OK.");
827     branchCPVimpacts ->GetEntry(itrack,0);
828
829     for (Int_t 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 //    Info("AnalyzeCPV1", " 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("AnalyzeCPV1", " branchCPVimpacts ->GetEntry(itrack,0) OK.");
865     
866 //      for (Int_t iModule=0; iModule < nOfModules; iModule++) {
867 //        impacts = (TClonesArray *)fCpvImpacts->At(iModule);
868 //        Info("AnalyzeCPV1", " 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("AnalyzeCPV1", " Add charged hit..";
876 //          new(hits[hits.GetEntriesFast()]) AliPHOSImpact(*impact);
877 //          Info("AnalyzeCPV1", "done.");
878 //        }
879 //        }
880 //      }
881 //      fCpvImpacts->Clear();
882 //    }
883
884 //    Info("AnalyzeCPV1", " 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("AnalyzeCPV1", " PHOS module "<<iModule<<": "<<hits->GetEntries()<<" charged CPV hits.");
908 //  }
909
910 Bool_t AliPHOSIhepAnalyze::IsCharged(Int_t pdg_code)
911 {
912   // For HIJING
913   Info("AnalyzeCPV1", "pdg_code %d", pdg_code);
914   if(pdg_code==211 || pdg_code==-211 || pdg_code==321 || pdg_code==-321 || pdg_code==11 || pdg_code==-11 || pdg_code==2212 || pdg_code==-2212) return kTRUE;
915   else
916     return kFALSE;
917 }
918 //---------------------------------------------------------------------------
919
920
921
922
923
924