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