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