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