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