]> 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
GetRandom returning impact parameter bin and flag for hard process 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
cbd576a6 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
54ClassImp(AliPHOSIhepAnalyze)
55
56//____________________________________________________________________________
57
58 AliPHOSIhepAnalyze::AliPHOSIhepAnalyze() {}
59
60//____________________________________________________________________________
61
62AliPHOSIhepAnalyze::AliPHOSIhepAnalyze(Text_t * name) : fFileName(name) {}
63
64//____________________________________________________________________________
65void 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
21cd0c07 93 Info("AnalyzeCPV1", "Start CPV Analysis-1. Resolutions, cluster multiplicity and lengths") ;
cbd576a6 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);
21cd0c07 100 Info("AnalyzeCPV1", ">>>>>>>Event %d .<<<<<<<", ievent) ;
cbd576a6 101
102 // Get branch of CPV impacts
103 if (! (branchCPVimpacts =gAlice->TreeH()->GetBranch("PHOSCpvImpacts")) ) {
21cd0c07 104 Info("AnalyzeCPV1", "Couldn't find branch PHOSCpvImpacts. Exit.") ;
105 return;
cbd576a6 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 }
21cd0c07 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);
cbd576a6 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
21cd0c07 224 Info("AnalyzeCPV1", "++++ Event %d : total %d impacts, %d charged impacts and %d rec. points.",
225 ievent, nTotalGen, nChargedGen, please->CpvRecPoints()->GetEntries()) ;
cbd576a6 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
289void 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
21cd0c07 312 Info("AnalyzeCPV1", "Start EMC Analysis-1. Resolutions, cluster multiplicity and lengths");
cbd576a6 313 for ( Int_t ievent=0; ievent<Nevents; ievent++) {
314
315 Int_t nTotalGen = 0;
316
317 Int_t ntracks = gAlice->GetEvent(ievent);
21cd0c07 318 Info("AnalyzeCPV1", " >>>>>>>Event %d .<<<<<<<", ievent) ;
cbd576a6 319
320 // Get branch of EMC impacts
321 if (! (branchEMCimpacts =gAlice->TreeH()->GetBranch("PHOSEmcImpacts")) ) {
21cd0c07 322 Info("AnalyzeCPV1", " Couldn't find branch PHOSEmcImpacts. Exit.");
cbd576a6 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 }
21cd0c07 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) ;
cbd576a6 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
21cd0c07 436 Info("AnalyzeCPV1", "++++ Event %d : total %d impacts, %d Emc rec. points.",
437 ievent, nTotalGen, please->EmcRecPoints()->GetEntriesFast()) ;
cbd576a6 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//____________________________________________________________________________
491void 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
0bc3b8ed 500 TH1F* hDrijCpvR = new TH1F("DrijCpvR","Distance between reconstructed hits in CPV",140,0,50);
501 TH1F* hDrijCpvG = new TH1F("Drij_cpv_g","Distance between generated hits in CPV",140,0,50);
502 TH1F* hDrijCpvRatio = new TH1F("DrijCpvRatio","R_{ij}^{rec}/R_{ij}^{gen} in CPV",140,0,50);
cbd576a6 503
504// TH1F* hT0 = new TH1F("hT0","Type of entering particle",20000,-10000,10000);
505
0bc3b8ed 506 hDrijCpvR->Sumw2();
507 hDrijCpvG->Sumw2();
508 hDrijCpvRatio->Sumw2(); //correct treatment of errors
cbd576a6 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
0bc3b8ed 522 Int_t nrecCPV = 0; // Reconstructed points in event
523 Int_t ngenCPV = 0; // Impacts in event
cbd576a6 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);
21cd0c07 540 Info("AnalyzeCPV1", " branchCPVimpacts ->SetAddress(&fCpvImpacts) OK.");
cbd576a6 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
0bc3b8ed 560 AliPHOSImpact* genHit1;
561 AliPHOSImpact* genHit2;
cbd576a6 562 Int_t irp1,irp2;
563 for(irp1=0; irp1< nsum; irp1++)
564 {
0bc3b8ed 565 genHit1 = (AliPHOSImpact*)((hitsPerModule[iModule])->At(irp1));
cbd576a6 566 for(irp2 = irp1+1; irp2<nsum; irp2++)
567 {
0bc3b8ed 568 genHit2 = (AliPHOSImpact*)((hitsPerModule[iModule])->At(irp2));
569 Float_t dx = genHit1->X() - genHit2->X();
570 Float_t dz = genHit1->Z() - genHit2->Z();
cbd576a6 571 Float_t dr = TMath::Sqrt(dx*dx + dz*dz);
0bc3b8ed 572 hDrijCpvG->Fill(dr);
21cd0c07 573// Info("AnalyzeCPV1", "(dx dz dr): %f %f", dx, dz);
cbd576a6 574 }
575 }
576 }
577
578
579 //--------- Combinatoric distance between rec. hits in CPV
580
581 TObjArray* cpvRecPoints = please->CpvRecPoints();
0bc3b8ed 582 nrecCPV = cpvRecPoints->GetEntriesFast();
cbd576a6 583
0bc3b8ed 584 if(nrecCPV)
cbd576a6 585 {
0bc3b8ed 586 AliPHOSCpvRecPoint* recHit1;
587 AliPHOSCpvRecPoint* recHit2;
588 TIter nextCPVrec1(cpvRecPoints);
589 while(TObject* obj1 = nextCPVrec1() )
cbd576a6 590 {
0bc3b8ed 591 TIter nextCPVrec2(cpvRecPoints);
592 while (TObject* obj2 = nextCPVrec2())
cbd576a6 593 {
594 if(!obj2->IsEqual(obj1))
595 {
0bc3b8ed 596 recHit1 = (AliPHOSCpvRecPoint*)obj1;
597 recHit2 = (AliPHOSCpvRecPoint*)obj2;
cbd576a6 598 TVector3 locpos1;
599 TVector3 locpos2;
0bc3b8ed 600 recHit1->GetLocalPosition(locpos1);
601 recHit2->GetLocalPosition(locpos2);
cbd576a6 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);
0bc3b8ed 605 if(recHit1->GetPHOSMod() == recHit2->GetPHOSMod())
606 hDrijCpvR->Fill(dr);
cbd576a6 607 }
608 }
609 }
610 }
611
21cd0c07 612 Info("AnalyzeCPV1", " Event %d . Total of %d hits, %d rec.points.",
0bc3b8ed 613 nev, ngenCPV, nrecCPV) ;
cbd576a6 614
615 delete [] hitsPerModule;
616
617 } // End of loop over events.
618
619
0bc3b8ed 620// hDrijCpvG->Draw();
621// hDrijCpvR->Draw();
622 hDrijCpvRatio->Divide(hDrijCpvR,hDrijCpvG);
623 hDrijCpvRatio->Draw();
cbd576a6 624
625// hT0->Draw();
626
627}
628
629
630void 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());
21cd0c07 664 Info("AnalyzeCPV1", "+++++ Event %d . (Mult,MultX,MultZ) = %d %d %d +++++",
665 ievent, rpMult, rpMultX, rpMultZ) ;
cbd576a6 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
717void 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
0bc3b8ed 732// Int_t nrecCPV = 0; // Reconstructed points in event // 01.10.2001
733// Int_t ngenCPV = 0; // Impacts in event
cbd576a6 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);
21cd0c07 750 Info("AnalyzeCPV1", " branchCPVimpacts ->SetAddress(&fCpvImpacts) OK.");
cbd576a6 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();
21cd0c07 780// Info("AnalyzeCPV1", " Tracks: "<<ntracks<<" Modules: "<<nOfModules);
cbd576a6 781
782// if (! (branchCPVimpacts =gAlice->TreeH()->GetBranch("PHOSCpvImpacts")) ) return;
783
784// for (Int_t itrack=0; itrack<ntracks; itrack++) {
785// branchCPVimpacts ->SetAddress(&fCpvImpacts);
21cd0c07 786// Info("AnalyzeCPV1", " branchCPVimpacts ->SetAddress(&fCpvImpacts) OK.");
cbd576a6 787// branchCPVimpacts ->GetEntry(itrack,0);
21cd0c07 788// Info("AnalyzeCPV1", " branchCPVimpacts ->GetEntry(itrack,0) OK.");
cbd576a6 789
790// for (Int_t iModule=0; iModule < nOfModules; iModule++) {
791// impacts = (TClonesArray *)fCpvImpacts->At(iModule);
21cd0c07 792// Info("AnalyzeCPV1", " fCpvImpacts->At(iModule) OK.");
cbd576a6 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// {
21cd0c07 799// Info("AnalyzeCPV1", " Add charged hit..";
cbd576a6 800// new(hits[hits.GetEntriesFast()]) AliPHOSImpact(*impact);
21cd0c07 801// Info("AnalyzeCPV1", "done.");
cbd576a6 802// }
803// }
804// }
805// fCpvImpacts->Clear();
806// }
807
21cd0c07 808// Info("AnalyzeCPV1", " PHOS event "<<ievent<<": "<<hits.GetEntries()<<" charged CPV hits.");
cbd576a6 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
21cd0c07 831// Info("AnalyzeCPV1", " PHOS module "<<iModule<<": "<<hits->GetEntries()<<" charged CPV hits.");
cbd576a6 832// }
833
0bc3b8ed 834Bool_t AliPHOSIhepAnalyze::IsCharged(Int_t pdgCode)
cbd576a6 835{
836 // For HIJING
0bc3b8ed 837 Info("AnalyzeCPV1", "pdgCode %d", pdgCode);
838 if(pdgCode==211 || pdgCode==-211 || pdgCode==321 || pdgCode==-321 || pdgCode==11 || pdgCode==-11 || pdgCode==2212 || pdgCode==-2212) return kTRUE;
cbd576a6 839 else
840 return kFALSE;
841}
842//---------------------------------------------------------------------------
843
844
845
846
847
848