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