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