]> 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
Adding ruleckecher files into the distribbution
[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
6c8cd883 116 AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance() ;
cbd576a6 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
6c8cd883 144 Int_t nOfModules = phosgeom->GetNModules();
cbd576a6 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
02d7f1c9 161 for (iModule=0; iModule < nOfModules; iModule++) {
cbd576a6 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
6c8cd883 218 Float_t phig = phosgeom->GetPHOSAngle(phosModule);
cbd576a6 219 Float_t phi = TMath::Pi()/180*phig;
6c8cd883 220 Float_t distanceIPtoCPV = phosgeom->GetIPtoOuterCoverDistance() -
221 (phosgeom->GetFTPosition(1)+
222 phosgeom->GetFTPosition(2)+
223 phosgeom->GetCPVTextoliteThickness())/2;
cbd576a6 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] ;
63b20aec 262 snprintf(outputname,80,"%s.analyzed",GetFileName().Data());
cbd576a6 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
6c8cd883 348 AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance() ;
cbd576a6 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
6c8cd883 375 Int_t nOfModules = phosgeom->GetNModules();
cbd576a6 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
02d7f1c9 392 for (iModule=0; iModule < nOfModules; iModule++) {
cbd576a6 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
6c8cd883 447 Float_t phig = phosgeom->GetPHOSAngle(phosModule);
cbd576a6 448 Float_t phi = TMath::Pi()/180*phig;
6c8cd883 449 Float_t distanceIPtoEMC = phosgeom->GetIPtoCrystalSurface();
cbd576a6 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 }
351dd634 469 AliInfo(Form(" Impact global (X,Z,Y) = %f %f %f", gImpX, gImpZ, gImpY));
470 AliInfo(Form(" Impact local (X,Z) = %f %f", locImpX, locImpZ));
471 AliInfo(Form(" Reconstructed (X,Z) = %f %f", xrec, zrec));
472 AliInfo(Form(" 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
351dd634 482 AliInfo(Form("++++ 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] ;
63b20aec 489 snprintf(outputname,80,"%s.analyzed",GetFileName().Data());
cbd576a6 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 {
351dd634 562 AliError(Form("Could not obtain the Loader object !"));
88cb7938 563 return ;
564 }
6c8cd883 565 AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance() ;
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 {
351dd634 581 AliError(Form("Can not get TreeH"));
88cb7938 582 return;
583 }
584
585 if (! (branchCPVimpacts =treeH->GetBranch("PHOSCpvImpacts")) ) return;
cbd576a6 586
587 // Create and fill arrays of hits for each CPV module
6c8cd883 588 Int_t nOfModules = phosgeom->GetNModules();
cbd576a6 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);
351dd634 599 AliInfo(Form(" branchCPVimpacts ->SetAddress(&fCpvImpacts) OK."));
cbd576a6 600 branchCPVimpacts ->GetEntry(itrack,0);
601
02d7f1c9 602 for (iModule=0; iModule < nOfModules; iModule++) {
cbd576a6 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);
351dd634 632// AliInfo(Form("(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
351dd634 671 AliInfo(Form(" Event %d . Total of %d hits, %d rec.points.",
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 {
351dd634 705 AliError(Form("Could not obtain the Loader object !"));
88cb7938 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());
351dd634 728 AliInfo(Form("+++++ Event %d . (Mult,MultX,MultZ) = %d %d %d +++++",
729 ievent, rpMult, rpMultX, rpMultZ)) ;
cbd576a6 730
731 }
732
733 }
734
735 Text_t outputname[80] ;
63b20aec 736 snprintf(outputname,80,"%s.analyzed.single",GetFileName().Data());
cbd576a6 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
88cb7938 785 AliPHOSLoader* please = dynamic_cast<AliPHOSLoader*>(fRunLoader->GetLoader("PHOSLoader"));
786 if ( please == 0 )
787 {
351dd634 788 AliError(Form("Could not obtain the Loader object !"));
88cb7938 789 return ;
790 }
cbd576a6 791
792
793 printf("\n=================== Event %10d ===================\n",nev);
63b20aec 794 //16.03.2011: DP. Code below seems to be obsollete
795 //Comment it to sutisfy Coverity
796/*
215730f8 797 TList * fCpvImpacts ;
798 TBranch * branchCPVimpacts;
799
8789edd0 800 AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance() ;
63b20aec 801
88cb7938 802 fRunLoader->GetEvent(nev);
803 Int_t ntracks = fRunLoader->GetHeader()->GetNtrack();
cbd576a6 804
6a2b5f60 805// Int_t nRecCPV = 0; // Reconstructed points in event // 01.10.2001
806// Int_t nGenCPV = 0; // Impacts in event
cbd576a6 807
808 // Get branch of CPV impacts
88cb7938 809 TTree* treeH = please->TreeH();
810 if (treeH == 0x0)
811 {
351dd634 812 AliError(Form("Can not get TreeH"));
88cb7938 813 return;
814 }
815
816 if (! (branchCPVimpacts =treeH->GetBranch("PHOSCpvImpacts")) ) return;
cbd576a6 817
818 // Create and fill arrays of hits for each CPV module
6c8cd883 819 Int_t nOfModules = phosgeom->GetNModules();
cbd576a6 820 TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
821 Int_t iModule = 0;
822 for (iModule=0; iModule < nOfModules; iModule++)
823 hitsPerModule[iModule] = new TClonesArray("AliPHOSImpact",100);
824
825 TClonesArray *impacts;
826 AliPHOSImpact *impact;
827
828 for (Int_t itrack=0; itrack<ntracks; itrack++) {
829 branchCPVimpacts ->SetAddress(&fCpvImpacts);
351dd634 830 AliInfo(Form(" branchCPVimpacts ->SetAddress(&fCpvImpacts) OK."));
cbd576a6 831 branchCPVimpacts ->GetEntry(itrack,0);
832
02d7f1c9 833 for (iModule=0; iModule < nOfModules; iModule++) {
cbd576a6 834 impacts = (TClonesArray *)fCpvImpacts->At(iModule);
835 // Do loop over impacts in the module
836 for (Int_t iImpact=0; iImpact<impacts->GetEntries(); iImpact++) {
837 impact=(AliPHOSImpact*)impacts->At(iImpact);
838 TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
839 new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSImpact(*impact);
840 }
841 }
842 fCpvImpacts->Clear();
843 }
844
845 for (iModule=0; iModule < nOfModules; iModule++) {
846 Int_t nsum = hitsPerModule[iModule]->GetEntriesFast();
847 printf("CPV module %d has %d hits\n",iModule,nsum);
848 }
849
63b20aec 850*/
851
852
cbd576a6 853// TList * fCpvImpacts ;
854// TBranch * branchCPVimpacts;
855// AliPHOSImpact* impact;
856// TClonesArray* impacts;
857
88cb7938 858// AliPHOSLoader * please = AliPHOSLoader::GetInstance(GetFileName().Data(),"PHOS");
cbd576a6 859// const AliPHOSGeometry * fGeom = please->PHOSGeometry();
860
861// Int_t ntracks = gAlice->GetEvent(ievent);
862// Int_t nOfModules = fGeom->GetNModules();
351dd634 863// AliInfo(Form(" Tracks: "<<ntracks<<" Modules: "<<nOfModules));
cbd576a6 864
865// if (! (branchCPVimpacts =gAlice->TreeH()->GetBranch("PHOSCpvImpacts")) ) return;
866
867// for (Int_t itrack=0; itrack<ntracks; itrack++) {
868// branchCPVimpacts ->SetAddress(&fCpvImpacts);
21cd0c07 869// Info("AnalyzeCPV1", " branchCPVimpacts ->SetAddress(&fCpvImpacts) OK.");
cbd576a6 870// branchCPVimpacts ->GetEntry(itrack,0);
351dd634 871// Info(Form(" branchCPVimpacts ->GetEntry(itrack,0) OK."));
cbd576a6 872
873// for (Int_t iModule=0; iModule < nOfModules; iModule++) {
874// impacts = (TClonesArray *)fCpvImpacts->At(iModule);
351dd634 875// Info(Form(" fCpvImpacts->At(iModule) OK."));
cbd576a6 876// // Do loop over impacts in the module
877// for (Int_t iImpact=0; iImpact<impacts->GetEntries(); iImpact++) {
878// impact=(AliPHOSImpact*)impacts->At(iImpact);
879// impact->Print();
880// if(IsCharged(impact->GetPid()))
881// {
351dd634 882// Info(Form(" Add charged hit.."));
cbd576a6 883// new(hits[hits.GetEntriesFast()]) AliPHOSImpact(*impact);
351dd634 884// Info(Form("done."));
cbd576a6 885// }
886// }
887// }
888// fCpvImpacts->Clear();
889// }
890
351dd634 891// Info(Form(" PHOS event "<<ievent<<": "<<hits.GetEntries()<<" charged CPV hits."));
cbd576a6 892
893}
894
895
896// void AliPHOSIhepAnalyze::ChargedHitsCPV(TClonesArray* hits, Int_t ievent, Int_t iModule)
897// {
898// // Cumulative list of charged CPV hits in event ievent
899// // in PHOS module iModule.
900
901// HitsCPV(hits,ievent,iModule);
902// TIter next(hits);
903
904// while(AliPHOSCPVHit* cpvhit = (AliPHOSCPVHit*)next())
905// {
906// if(!IsCharged(cpvhit->GetIpart()))
907// {
908// hits->Remove(cpvhit);
909// delete cpvhit;
910// hits->Compress();
911// }
912// }
913
351dd634 914// Info(Form(" PHOS module "<<iModule<<": "<<hits->GetEntries()<<" charged CPV hits."));
cbd576a6 915// }
916
6a2b5f60 917Bool_t AliPHOSIhepAnalyze::IsCharged(Int_t pdgCode)
cbd576a6 918{
919 // For HIJING
351dd634 920 AliInfo(Form("pdgCode %d", pdgCode));
6a2b5f60 921 if(pdgCode==211 || pdgCode==-211 || pdgCode==321 || pdgCode==-321 || pdgCode==11 || pdgCode==-11 || pdgCode==2212 || pdgCode==-2212) return kTRUE;
cbd576a6 922 else
923 return kFALSE;
924}
925//---------------------------------------------------------------------------
926
927
928
929
930
931