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