]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PHOS/AliPHOS.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Medium array is now private for each module
[u/mrichter/AliRoot.git] / PHOS / AliPHOS.cxx
CommitLineData
fe4da5cc 1////////////////////////////////////////////////
2// Manager and hits classes for set:PHOS //
3////////////////////////////////////////////////
4
5// --- ROOT system ---
6#include "TH1.h"
7#include "TRandom.h"
8#include "TFile.h"
9#include "TTree.h"
10#include "TBRIK.h"
11#include "TNode.h"
bc9ab547 12#include "TMath.h"
fe4da5cc 13
14// --- Standard library ---
15#include <stdio.h>
16#include <string.h>
17#include <stdlib.h>
18
19// --- galice header files ---
20#include "AliPHOS.h"
21#include "AliRun.h"
22#include "AliMC.h"
23#include "TGeant3.h"
24
25//______________________________________________________________________________
26
27
28ClassImp(AliPHOS)
29
30//______________________________________________________________________________
31
32AliPHOS::~AliPHOS(void)
33{
bc9ab547 34 delete fHits; // 28.12.1998
35 delete fTreePHOS; // 28.12.1998
fe4da5cc 36 fCradles->Delete();
37 delete fCradles;
38}
39
40//______________________________________________________________________________
41
42AliPHOS::AliPHOS() :
43 fDebugLevel (0),
44 fTreePHOS (NULL),
45 fBranchNameOfCradles ("AliPHOSCradles"),
46 fTreeName ("PHOS")
47{
48 fIshunt = 0;
49
50 if( NULL==(fCradles=new TObjArray) )
51 {
52 Error("AliPHOS","Can not create fCradles");
53 exit(1);
54 }
55 DefPars();
56}
57
58//______________________________________________________________________________
59
60AliPHOS::AliPHOS(const char *name, const char *title)
61 : AliDetector (name,title),
62 fDebugLevel (0),
63 fTreePHOS (NULL),
64 fBranchNameOfCradles ("AliPHOSCradles"),
65 fTreeName ("PHOS")
66{
67//Begin_Html
68/*
1439f98e 69<img src="picts/aliphos.gif">
fe4da5cc 70*/
71//End_Html
72
73 fHits = new TClonesArray("AliPHOShit", 405);
74
75 fIshunt = 0;
76
77 SetMarkerColor(kGreen);
78 SetMarkerStyle(2);
79 SetMarkerSize(0.4);
80
81 if( NULL==(fCradles=new TObjArray) ) {
82 Error("AliPHOS","Can not create fCradles");
83 exit(1);
84 }
85 DefPars();
86}
87
88//______________________________________________________________________________
89
90void AliPHOS::DefPars()
91{
92 PHOSflags[0]=0;
93 PHOSflags[1]=1;
94 PHOSflags[2]=0;
95 PHOSflags[3]=0;
96 PHOSflags[4]=0;
97 PHOSflags[5]=0;
98 PHOSflags[6]=0;
99 PHOSflags[7]=0;
100 PHOSflags[8]=0;
101 PHOScell[0]=2.2;
102 PHOScell[1]=18.;
103 PHOScell[2]=0.01;
104 PHOScell[3]=0.01;
105 PHOScell[4]=1.0;
106 PHOScell[5]=0.1;
107 PHOScell[6]=0.;
108 PHOScell[7]=0.;
109 PHOScell[8]=0.;
110 PHOSradius=460.;
111 PHOSsize[0]=104;
112 PHOSsize[1]=88;
113 PHOSsize[2]=4;
114 PHOScradlesA=0.;
115 PHOSCPV[0]=1.;
116 PHOSCPV[1]=2.;
117 PHOSCPV[2]=0.;
118 PHOSCPV[3]=0.;
119 PHOSCPV[4]=0.;
120 PHOSCPV[5]=0.;
121 PHOSCPV[6]=0.;
122 PHOSCPV[7]=0.;
123 PHOSCPV[8]=0.;
124 PHOSextra[0]=0.001;
125 PHOSextra[1]=6.95;
126 PHOSextra[2]=4.;
127 PHOSextra[3]=5.;
128 PHOSextra[4]=2.;
129 PHOSextra[5]=0.06;
130 PHOSextra[6]=10.;
131 PHOSextra[7]=3.;
132 PHOSextra[8]=1.;
133 PHOSTXW[0]=209.;
134 PHOSTXW[1]=71.;
135 PHOSTXW[2]=250.;
136 PHOSAIR[0]=206.;
137 PHOSAIR[1]=66.;
138 PHOSAIR[2]=244.;
139 PHOSFTI[0]=214.6;
140 PHOSFTI[1]=80.;
141 PHOSFTI[2]=260.;
142 PHOSFTI[3]=467.;
143}
144//______________________________________________________________________________
145
146void AliPHOS::AddHit(Int_t track, Int_t *vol, Float_t *hits)
147{
148 TClonesArray &lhits = *fHits;
149 new(lhits[fNhits++]) AliPHOShit(fIshunt,track,vol,hits);
150}
151
152//___________________________________________
153void AliPHOS::BuildGeometry()
154{
155
156 TNode *Node, *Top;
157
158 const int kColorPHOS = kRed;
159 //
160 Top=gAlice->GetGeometry()->GetNode("alice");
161
162
163 // PHOS
164 Float_t pphi=12.9399462;
165 new TRotMatrix("rot988","rot988",90,-3*pphi,90,90-3*pphi,0,0);
166 new TRotMatrix("rot989","rot989",90,- pphi,90,90- pphi,0,0);
167 new TRotMatrix("rot990","rot990",90, pphi,90,90+ pphi,0,0);
168 new TRotMatrix("rot991","rot991",90, 3*pphi,90,90+3*pphi,0,0);
169 new TBRIK("S_PHOS","PHOS box","void",107.3,40,130);
170 Top->cd();
171 Node = new TNode("PHOS1","PHOS1","S_PHOS",-317.824921,-395.014343,0,"rot988");
172 Node->SetLineColor(kColorPHOS);
173 fNodes->Add(Node);
174 Top->cd();
175 Node = new TNode("PHOS2","PHOS2","S_PHOS",-113.532333,-494.124908,0,"rot989");
176 fNodes->Add(Node);
177 Node->SetLineColor(kColorPHOS);
178 Top->cd();
179 Node = new TNode("PHOS3","PHOS3","S_PHOS", 113.532333,-494.124908,0,"rot990");
180 Node->SetLineColor(kColorPHOS);
181 fNodes->Add(Node);
182 Top->cd();
183 Node = new TNode("PHOS4","PHOS4","S_PHOS", 317.824921,-395.014343,0,"rot991");
184 Node->SetLineColor(kColorPHOS);
185 fNodes->Add(Node);
186}
187
188//___________________________________________
189void AliPHOS::CreateMaterials()
190{
191// *** DEFINITION OF AVAILABLE PHOS MATERIALS ***
192
193// CALLED BY : PHOS_MEDIA
194// ORIGIN : NICK VAN EIJNDHOVEN
195
196
197 AliMC* pMC = AliMC::GetMC();
198
199 Int_t ISXFLD = gAlice->Field()->Integ();
200 Float_t SXMGMX = gAlice->Field()->Max();
201
202// --- The PbWO4 crystals ---
203 Float_t ax[3] = { 207.19,183.85,16. };
204 Float_t zx[3] = { 82.,74.,8. };
205 Float_t wx[3] = { 1.,1.,4. };
206 Float_t dx = 8.28;
207// --- Stainless Steel ---
208 Float_t as[5] = { 55.847,12.011,51.9961,58.69,28.0855 };
209 Float_t zs[5] = { 26.,6.,24.,28.,14. };
210 Float_t ws[5] = { .6392,8e-4,.2,.14,.02 };
211 Float_t ds = 8.;
212// --- The polysterene scintillator (CH) ---
213 Float_t ap[2] = { 12.011,1.00794 };
214 Float_t zp[2] = { 6.,1. };
215 Float_t wp[2] = { 1.,1. };
216 Float_t dp = 1.032;
217// --- Tyvek (CnH2n)
218 Float_t at[2] = { 12.011,1.00794 };
219 Float_t zt[2] = { 6.,1. };
220 Float_t wt[2] = { 1.,2. };
221 Float_t dt = .331;
222// --- Polystyrene foam ---
223 Float_t af[2] = { 12.011,1.00794 };
224 Float_t zf[2] = { 6.,1. };
225 Float_t wf[2] = { 1.,1. };
226 Float_t df = .12;
227//--- Foam thermo insulation (actual chemical composition unknown yet!) ---
228 Float_t ati[2] = { 12.011,1.00794 };
229 Float_t zti[2] = { 6.,1. };
230 Float_t wti[2] = { 1.,1. };
231 Float_t dti = .1;
232// --- Textolit (actual chemical composition unknown yet!) ---
233 Float_t atx[2] = { 12.011,1.00794 };
234 Float_t ztx[2] = { 6.,1. };
235 Float_t wtx[2] = { 1.,1. };
236 Float_t dtx = 1.83;
237
ad51aeb0 238 Int_t *idtmed = fIdtmed->GetArray()-699;
fe4da5cc 239
240 AliMixture( 0, "PbWO4$", ax, zx, dx, -3, wx);
241 AliMixture( 1, "Polystyrene$", ap, zp, dp, -2, wp);
242 AliMaterial( 2, "Al$", 26.98, 13., 2.7, 8.9, 999);
243// --- Absorption length^ is ignored ---
244 AliMixture( 3, "Tyvek$", at, zt, dt, -2, wt);
245 AliMixture( 4, "Foam$", af, zf, df, -2, wf);
246 AliMixture( 5, "Stainless Steel$", as, zs, ds, 5, ws);
247 AliMaterial( 6, "Si$", 28.09, 14., 2.33, 9.36, 42.3);
248 AliMixture( 7, "Thermo Insul.$", ati, zti, dti, -2, wti);
249 AliMixture( 8, "Textolit$", atx, ztx, dtx, -2, wtx);
250 AliMaterial(99, "Air$", 14.61, 7.3, .001205, 30420., 67500);
251
ad51aeb0 252 AliMedium(0, "PHOS Xtal $", 0, 1, ISXFLD, SXMGMX, 10., .1, .1, .1, .1);
253 AliMedium(1, "CPV scint. $", 1, 1, ISXFLD, SXMGMX, 10., .1, .1, .1, .1);
254 AliMedium(2, "Al parts $", 2, 0, ISXFLD, SXMGMX, 10., .1, .1, .001, .001);
255 AliMedium(3, "Tyvek wrapper$", 3, 0, ISXFLD, SXMGMX, 10., .1, .1, .001, .001);
256 AliMedium(4, "Polyst. foam $", 4, 0, ISXFLD, SXMGMX, 10., .1, .1, .1, .1);
257 AliMedium(5, "Steel cover $", 5, 0, ISXFLD, SXMGMX, 10., .1, .1, 1e-4, 1e-4);
258 AliMedium(6, "Si PIN $", 6, 0, ISXFLD, SXMGMX, 10., .1, .1, .01, .01);
259 AliMedium(7, "Thermo Insul.$", 7, 0, ISXFLD, SXMGMX, 10., .1, .1, .1, .1);
260 AliMedium(8, "Textolit $", 8, 0, ISXFLD, SXMGMX, 10., .1, .1, .1, .1);
261 AliMedium(99, "Air $",99, 0, ISXFLD, SXMGMX, 10., 1., .1, .1, 10);
fe4da5cc 262
263// --- Generate explicitly delta rays in the steel cover ---
264 pMC->Gstpar(idtmed[704], "LOSS", 3.);
265 pMC->Gstpar(idtmed[704], "DRAY", 1.);
266// --- and in aluminium parts ---
267 pMC->Gstpar(idtmed[701], "LOSS", 3.);
268 pMC->Gstpar(idtmed[701], "DRAY", 1.);
269}
270
271//______________________________________________________________________________
272
273void AliPHOS::AddPHOSCradles()
274{
275 Int_t i;
276 for(i=0;i<GetCradlesAmount();i++) {
277
278 int n = fCradles->GetEntries();
279 fCradles->Add(new AliPHOSCradle( IsVersion(), // geometry.
280 GetCrystalSideSize (),
281 GetCrystalLength (),
282 GetWrapThickness (),
283 GetAirThickness (),
284 GetPIN_SideSize (),
285 GetPIN_Length (),
286 GetRadius (),
287 GetCPV_Thickness (),
288 GetCPV_PHOS_Distance (),
289 GetNz (),
290 GetNphi (),
291 GetCradleAngle (i)));
292
293 if( n+1 != fCradles->GetEntries() || NULL == fCradles->At(n) )
294 {
295 cout << " Can not create or add AliPHOSCradle.\n";
296 exit(1);
297 }
298 }
299}
300
301//______________________________________________________________________________
302
303Int_t AliPHOS::DistancetoPrimitive(Int_t , Int_t )
304{
305 return 9999;
306}
307
308//___________________________________________
309void AliPHOS::Init()
310{
311 Int_t i;
312 //
313 printf("\n");
314 for(i=0;i<35;i++) printf("*");
315 printf(" PHOS_INIT ");
316 for(i=0;i<35;i++) printf("*");
317 printf("\n");
318 //
319 // Here the ABSO initialisation code (if any!)
320 for(i=0;i<80;i++) printf("*");
321 printf("\n");
322}
323
324//______________________________________________________________________________
325
326void AliPHOS::MakeBranch(Option_t *)
327{
328// ROOT output initialization to ROOT file.
329//
330// AliDetector::MakeBranch() is always called.
331//
332// There will be also special tree "PHOS" with one branch "AliPHOSCradles"
333// if it was set next flag in the galice card file:
334// * PHOSflags: YES: X<>0 NO: X=0
335// * PHOSflags(1) : -----X. Create branch for TObjArray of AliPHOSCradle
336// Examples:
337// PHOSflags 1.
338// PHOSflags 636301.
339// In that case special bit CradlesBranch_Bit will be set for AliPHOS
340
341 AliDetector::MakeBranch();
342
343 int i;
344 float t = GetPHOS_flag(0)/10;
345 i = (int) t;
346 i = (int) ((t-i)*10);
347 if( !i )
348 return;
349
350 SetBit(CradlesBranch_Bit);
351
352 if( NULL==(fTreePHOS=new TTree(fTreeName.Data(),"PHOS events tree")) )
353 {
354 Error("MakeBranch","Can not create TTree");
355 exit(1);
356 }
357
358 if( NULL==fTreePHOS->GetCurrentFile() )
359 {
360 Error("MakeBranch","There is no opened ROOT file");
361 exit(1);
362 }
363
364 // Create a new branch in the current Root Tree.
365
366 if( NULL==fTreePHOS->Branch(fBranchNameOfCradles.Data(),"TObjArray",&fCradles,4000,0) )
367 {
368 Error("MakeBranch","Can not create branch");
369 exit(1);
370 }
371
372 printf("The branch %s has been created\n",fBranchNameOfCradles.Data());
373}
374
375//______________________________________________________________________________
376
377void AliPHOS::SetTreeAddress(void)
378{
379// ROOT input initialization.
380//
381// AliDetector::SetTreeAddress() is always called.
382//
383// If CradlesBranch_Bit is set (see AliPHOS::MakeBranch) than fTreePHOS is
384// initilized.
385
386 AliDetector::SetTreeAddress();
387
388 if( !TestBit(CradlesBranch_Bit) )
389 return;
390
391 if( NULL==(fTreePHOS=(TTree*)gDirectory->Get((char*)(fTreeName.Data())) ) )
392 {
393 Error("Can not find Tree \"%s\"\n",fTreeName.Data());
394 exit(1);
395 }
396
397 TBranch *branch = fTreePHOS->GetBranch(fBranchNameOfCradles.Data());
398 if( NULL==branch )
399 {
400 Error("SetTreeAddress","Can not find branch %s in TTree:%s",fBranchNameOfCradles.Data(),fTreeName.Data());
401 exit(1);
402 }
403
404 branch->SetAddress(&fCradles);
405}
406
407//______________________________________________________________________________
408
bc9ab547 409AliPHOSCradle *AliPHOS::GetCradleOfTheParticle(const TVector3 &p,const TVector3 &v) const
fe4da5cc 410{
411// For a given direction 'p' and source point 'v' returns pointer to AliPHOSCradle
412// in that direction or NULL if AliPHOSCradle was not found.
413
414 for( int m=0; m<fCradles->GetEntries(); m++ )
415 {
416 AliPHOS *PHOS = (AliPHOS *)this; // Removing 'const'...
417 AliPHOSCradle *cradle = (AliPHOSCradle *)PHOS->fCradles->operator[](m);
418
419 float x,y,l;
420 const float d = cradle->GetRadius()-cradle->GetCPV_PHOS_Distance()-cradle->GetCPV_Thikness();
421 cradle->GetXY(p,v,d,x,y,l);
422
bc9ab547 423 if( l>0 && TMath::Abs(x)<cradle->GetNz ()*cradle->GetCellSideSize()/2
424 && TMath::Abs(y)<cradle->GetNphi()*cradle->GetCellSideSize()/2 )
fe4da5cc 425 return cradle;
426 }
427
428 return NULL;
429}
430
431//______________________________________________________________________________
432
433void AliPHOS::Reconstruction(Float_t signal_step, UInt_t min_signal_reject)
434{
435// Call AliPHOSCradle::Reconstruction(Float_t signal_step, UInt_t min_signal_reject)
436// for all AliPHOSCradles.
437
438 for( int i=0; i<fCradles->GetEntries(); i++ )
439 GetCradle(i).Reconstruction(signal_step,min_signal_reject);
440}
441
442//______________________________________________________________________________
443
444void AliPHOS::ResetDigits(void)
445{
446 AliDetector::ResetDigits();
447
448 for( int i=0; i<fCradles->GetEntries(); i++ )
449 ((AliPHOSCradle*)(*fCradles)[i]) -> Clear();
450}
451
452//______________________________________________________________________________
453
454void AliPHOS::FinishEvent(void)
455{
456// Called at the end of each 'galice' event.
457
458 if( NULL!=fTreePHOS )
459 fTreePHOS->Fill();
460}
461
462//______________________________________________________________________________
463
464void AliPHOS::FinishRun(void)
465{
466}
467
468//______________________________________________________________________________
469
470void AliPHOS::Print(Option_t *opt)
471{
472// Print PHOS information.
473// For each AliPHOSCradle the function AliPHOSCradle::Print(opt) is called.
474
475 AliPHOS &PHOS = *(AliPHOS *)this; // Removing 'const'...
476
477 for( int i=0; i<fCradles->GetEntries(); i++ )
478 {
479 printf("PHOS cradle %d from %d\n",i+1, fCradles->GetEntries());
480 PHOS.GetCradle(i).Print(opt);
481 printf( "---------------------------------------------------\n");
482 }
483}
484
485//______________________________________________________________________________
486void AliPHOS::SetFlags(Float_t p1,Float_t p2,Float_t p3,Float_t p4,
487 Float_t p5,Float_t p6,Float_t p7,Float_t p8,Float_t p9)
488{
489 PHOSflags[0]=p1;
490 PHOSflags[1]=p2;
491 PHOSflags[2]=p3;
492 PHOSflags[3]=p4;
493 PHOSflags[4]=p5;
494 PHOSflags[5]=p6;
495 PHOSflags[6]=p7;
496 PHOSflags[7]=p8;
497 PHOSflags[8]=p9;
498}
499
500//______________________________________________________________________________
501void AliPHOS::SetCell(Float_t p1,Float_t p2,Float_t p3,Float_t p4,
502 Float_t p5,Float_t p6,Float_t p7,Float_t p8,Float_t p9)
503{
504 PHOScell[0]=p1;
505 PHOScell[1]=p2;
506 PHOScell[2]=p3;
507 PHOScell[3]=p4;
508 PHOScell[4]=p5;
509 PHOScell[5]=p6;
510 PHOScell[6]=p7;
511 PHOScell[7]=p8;
512 PHOScell[8]=p9;
513}
514
515//______________________________________________________________________________
516void AliPHOS::SetRadius(Float_t radius)
517{
518 PHOSradius=radius;
519}
520
521//______________________________________________________________________________
522void AliPHOS::SetCradleSize(Int_t nz, Int_t nphi, Int_t ncradles)
523{
524 PHOSsize[0]=nz;
525 PHOSsize[1]=nphi;
526 PHOSsize[2]=ncradles;
527}
528
529//______________________________________________________________________________
530void AliPHOS::SetCradleA(Float_t angle)
531{
532 PHOScradlesA=angle;
533}
534
535//______________________________________________________________________________
536void AliPHOS::SetCPV(Float_t p1,Float_t p2,Float_t p3,Float_t p4,
537 Float_t p5,Float_t p6,Float_t p7,Float_t p8,Float_t p9)
538{
539 PHOSCPV[0] = p1;
540 PHOSCPV[1] = p2;
541 PHOSCPV[2] = p3;
542 PHOSCPV[3] = p4;
543 PHOSCPV[4] = p5;
544 PHOSCPV[5] = p6;
545 PHOSCPV[6] = p7;
546 PHOSCPV[7] = p8;
547 PHOSCPV[8] = p9;
548}
549
550//______________________________________________________________________________
551void AliPHOS::SetExtra(Float_t p1,Float_t p2,Float_t p3,Float_t p4,
552 Float_t p5,Float_t p6,Float_t p7,Float_t p8,Float_t p9)
553{
554 PHOSextra[0] = p1;
555 PHOSextra[1] = p2;
556 PHOSextra[2] = p3;
557 PHOSextra[3] = p4;
558 PHOSextra[4] = p5;
559 PHOSextra[5] = p6;
560 PHOSextra[6] = p7;
561 PHOSextra[7] = p8;
562 PHOSextra[8] = p9;
563}
564
565//______________________________________________________________________________
566void AliPHOS::SetTextolitWall(Float_t dx, Float_t dy, Float_t dz)
567{
568 PHOSTXW[0] = dx;
569 PHOSTXW[1] = dy;
570 PHOSTXW[2] = dz;
571}
572
573//______________________________________________________________________________
574void AliPHOS::SetInnerAir(Float_t dx, Float_t dy, Float_t dz)
575{
576 PHOSAIR[0] = dx;
577 PHOSAIR[1] = dy;
578 PHOSAIR[2] = dz;
579}
580
581//______________________________________________________________________________
582void AliPHOS::SetFoam(Float_t dx, Float_t dy, Float_t dz, Float_t dr)
583{
584 PHOSFTI[0] = dx;
585 PHOSFTI[1] = dy;
586 PHOSFTI[2] = dz;
587 PHOSFTI[3] = dr;
588}
589
590ClassImp(AliPHOSCradle)
591
592//______________________________________________________________________________
593
594AliPHOSCradle::AliPHOSCradle(void) {}
595
596//______________________________________________________________________________
597
598AliPHOSCradle::AliPHOSCradle( int Geometry ,
599 float CrystalSideSize ,
600 float CrystalLength ,
601 float WrapThickness ,
602 float AirThickness ,
603 float PIN_SideSize ,
604 float PIN_Length ,
605 float Radius ,
606 float CPV_Thickness ,
607 float CPV_PHOS_Distance ,
608 int Nz ,
609 int Nphi ,
610 float Angle ) :
611 fGeometry (Geometry),
612// fCellEnergy (),
613// fChargedTracksInPIN (),
614// fCPV_hitsX (),
615// fCPV_hitsY (),
616 fCrystalSideSize (CrystalSideSize),
617 fCrystalLength (CrystalLength),
618 fWrapThickness (WrapThickness),
619 fAirThickness (AirThickness),
620 fPIN_SideSize (PIN_SideSize),
621 fPIN_Length (PIN_Length),
622 fRadius (Radius),
623 fCPV_PHOS_Distance (CPV_PHOS_Distance),
624 fCPV_Thickness (CPV_Thickness),
625 fNz (Nz),
626 fNphi (Nphi),
627 fPhi (Angle)
628{
629 fCellEnergy = TH2F("CellE","Energy deposition in a cells",fNz,0,fNz,fNphi,0,fNphi);
630 fCellEnergy .SetDirectory(0);
631 fChargedTracksInPIN = TH2S("PINCtracks","Amount of charged tracks in PIN",fNz,0,fNz,fNphi,0,fNphi);
632 fChargedTracksInPIN .SetDirectory(0);
633}
634
635//______________________________________________________________________________
636
bc9ab547 637AliPHOSCradle::~AliPHOSCradle(void) // 28.12.1998
638{
639 fGammasReconstructed.Delete();
640 fParticles .Delete();
641}
642
643//______________________________________________________________________________
644
fe4da5cc 645void AliPHOSCradle::Clear(Option_t *)
646{
647// Clear digit. information.
648
649 fCellEnergy .Reset();
650 fChargedTracksInPIN .Reset();
651 GetParticles() .Delete();
652 GetParticles() .Compress();
653 GetGammasReconstructed() .Delete();
654 GetGammasReconstructed() .Compress();
655
656 fCPV_hitsX.Set(0);
657 fCPV_hitsY.Set(0);
658}
659
660//______________________________________________________________________________
661
662void AliPHOSCradle::AddCPVHit(float x,float y)
663{
664// Add this hit to the hits list in CPV detector.
665
666 TArrayF a(fCPV_hitsX.GetSize()+1);
667
668 memcpy(a.GetArray(),fCPV_hitsX.GetArray(),sizeof(Float_t)*fCPV_hitsX.GetSize());
669 a[fCPV_hitsX.GetSize()] = x;
670 fCPV_hitsX = a;
671
672 // It must be: fCPV_hitsX.GetSize() == fCPV_hitsY.GetSize()
673
674 memcpy(a.GetArray(),fCPV_hitsY.GetArray(),sizeof(Float_t)*fCPV_hitsY.GetSize());
675 a[fCPV_hitsY.GetSize()] = y;
676 fCPV_hitsY = a;
677}
678
679//______________________________________________________________________________
680
bc9ab547 681void AliPHOSCradle::GetXY(const TVector3 &p,const TVector3 &v,float R,float &x,float &y,float &l) const
fe4da5cc 682{
683// This function calculates hit position (x,y) in the CRADLE cells plain from particle in
684// the direction given by 'p' (not required to be normalized) and start point
685// given by 3-vector 'v'. So the particle trajectory is t(l) = v + p*l
686// were 'l' is a number (distance from 'v' to CRADLE cells plain) and 't' is resulting
687// three-vector of trajectory point.
688//
689// After the call to this function user should test that l>=0 (the particle HITED the
690// plain) and (x,y) are in the region of CRADLE:
691//
692// Example:
693// AliPHOSCradle cradle(......);
bc9ab547 694// TVector3 p(....), v(....);
fe4da5cc 695// Float_t x,y,l;
696// cradle.GetXY(p,v,x,y,l);
bc9ab547 697// if( l<0 || TMath::Abs(x)>cradle.GetNz() *cradle.GetCellSideSize()/2
698// || TMath::Abs(y)>cradle.GetNphi()*cradle.GetCellSideSize()/2 )
fe4da5cc 699// cout << "Outside the CRADLE.\n";
700
701 // We have to create three vectors:
702 // s - central point on the PHOS surface
703 // n1 - first vector in CRADLE plain
704 // n2 - second vector in CRADLE plain
705 // This three vectors are orthonormalized.
706
bc9ab547 707 double phi = fPhi/180*TMath::Pi();
708 TVector3 n1( 0.0 , 0.0 , 1.0 ), // Z direction (X)
fe4da5cc 709 n2( -sin(phi) , cos(phi) , 0 ), // around beam (Y)
710 s ( R*cos(phi) , R*sin(phi) , 0 ); // central point
711
712 const double l1_min = 1e-2;
713 double l1,
bc9ab547 714 p_n1 = p*n1, // * - scalar product.
715 p_n2 = p*n2,
716 v_n1 = v*n1,
717 v_n2 = v*n2,
718 s_n1 = s*n1, // 0
719 s_n2 = s*n2; // 0
fe4da5cc 720
bc9ab547 721 if ( TMath::Abs(l1=p.X()-n1.X()*p_n1-n2.X()*p_n2)>l1_min )
722 { l = (-v.X()+s.X()+n1.X()*(v_n1-s_n1)+n2.X()*(v_n2-s_n2))/l1; }
723 else if ( TMath::Abs(l1=p.Y()-n1.Y()*p_n1-n2.Y()*p_n2)>l1_min )
724 { l = (-v.Y()+s.Y()+n1.Y()*(v_n1-s_n1)+n2.Y()*(v_n2-s_n2))/l1; }
725 else if ( TMath::Abs(l1=p.Z()-n1.Z()*p_n1-n2.Z()*p_n2)>l1_min )
726 { l = (-v.Z()+s.Z()+n1.Z()*(v_n1-s_n1)+n2.Z()*(v_n2-s_n2))/l1; }
727
728// double lx = (-v.X()+s.X()+n1.X()*(v.dot(n1)-s.dot(n1))+n2.X()*(v.dot(n2)-s.dot(n2)))/
729// (p.X()-n1.X()*p.dot(n1)-n2.X()*p.dot(n2)),
730// ly = (-v.Y()+s.Y()+n1.Y()*(v.dot(n1)-s.dot(n1))+n2.Y()*(v.dot(n2)-s.dot(n2)))/
731// (p.Y()-n1.Y()*p.dot(n1)-n2.Y()*p.dot(n2)),
732// lz = (-v.Z()+s.Z()+n1.Z()*(v.dot(n1)-s.dot(n1))+n2.Z()*(v.dot(n2)-s.dot(n2)))/
733// (p.Z()-n1.Z()*p.dot(n1)-n2.Z()*p.dot(n2));
734// cout.form("x: %g %g %g %g\n",lx,-v.X()+s.X()+n1.X()*(v.dot(n1)-s.dot(n1))+n2.X()*(v.dot(n2)-s.dot(n2)),p.X()-n1.X()*p.dot(n1)-n2.X()*p.dot(n2));
735// cout.form("y: %g %g %g %g\n",lx,-v.Y()+s.Y()+n1.Y()*(v.dot(n1)-s.dot(n1))+n2.Y()*(v.dot(n2)-s.dot(n2)),p.Y()-n1.Y()*p.dot(n1)-n2.Y()*p.dot(n2));
736// cout.form("z: %g %g %g %g\n",lx,-v.Z()+s.Z()+n1.Z()*(v.dot(n1)-s.dot(n1))+n2.Z()*(v.dot(n2)-s.dot(n2)),p.Z()-n1.Z()*p.dot(n1)-n2.Z()*p.dot(n2));
fe4da5cc 737// cout.form("lx,ly,lz = %g,%g,%g\n",lx,ly,lz);
738
739 x = p_n1*l + v_n1 - s_n1;
740 y = p_n2*l + v_n2 - s_n2;
741}
742
743//______________________________________________________________________________
744
745void AliPHOSCradle::Print(Option_t *opt)
746{
747// Print AliPHOSCradle information.
748//
749// options: 'd' - print energy deposition for EVERY cell
750// 'p' - print particles list that hit the cradle
751// 'r' - print list of reconstructed particles
752
753 AliPHOSCradle *cr = (AliPHOSCradle *)this; // Removing 'const'...
754
755 printf("AliPHOSCradle: Nz=%d Nphi=%d, fPhi=%f, E=%g, CPV hits amount = %d\n",fNz,fNphi,fPhi,
756 cr->fCellEnergy.GetSumOfWeights(),fCPV_hitsX.GetSize());
757
758 if( NULL!=strchr(opt,'d') )
759 {
760 printf("\n\nCells Energy (in MeV):\n\n |");
761 for( int x=0; x<fNz; x++ )
762 printf(" %4d|",x+1);
763 printf("\n");
764
765 for( int y=fNphi-1; y>=0; y-- )
766 {
767 printf("%3d|",y+1);
768 for( int x=0; x<fNz; x++ )
769 printf("%6d",(int)(cr->fCellEnergy.GetBinContent(cr->fCellEnergy.GetBin(x,y))*1000));
770 printf("\n");
771 }
772 printf("\n");
773 }
774
775 if( NULL!=strchr(opt,'p') )
776 {
777 printf("This cradle was hit by %d particles\n",
778 ((AliPHOSCradle*)this)->GetParticles().GetEntries());
779 TObjArray &p=((AliPHOSCradle*)this)->GetParticles();
780 for( int i=0; i<p.GetEntries(); i++ )
781 ((AliPHOSgamma*)(p[i]))->Print();
782 }
783
784 if( NULL!=strchr(opt,'p') )
785 {
786 printf("Amount of reconstructed gammas is %d\n",
787 ((AliPHOSCradle*)this)->GetGammasReconstructed().GetEntries());
788
789 TObjArray &p=((AliPHOSCradle*)this)->GetGammasReconstructed();
790 for( int i=0; i<p.GetEntries(); i++ )
791 ((AliPHOSgamma*)(p[i]))->Print();
792 }
793}
794
795//______________________________________________________________________________
796
797void AliPHOSCradle::Distortion(const TH2F *Noise, const TH2F *Stochastic, const TH2F *Calibration)
798{
799// This function changes histogram of cell energies fCellEnergy on the base of input
800// histograms Noise, Stochastic, Calibration. The histograms must have
801// size Nz x Nphi.
802
803 //////////////////////////////////
804 // Testing the histograms size. //
805 //////////////////////////////////
806
807 if( fNz!=fCellEnergy.GetNbinsX() || fNphi!=fCellEnergy.GetNbinsY() )
808 {
809 printf ("Bad size of CellEnergy! Must be: Nz x Nphi = %d x %d\n"
810 "but size of CellEnergy is: %d x %d\n",
811 fNz,fNphi,fCellEnergy.GetNbinsX(),fCellEnergy.GetNbinsY());
812 exit(1);
813 }
814
815 if( fNz!=fChargedTracksInPIN.GetNbinsX() || fNphi!=fChargedTracksInPIN.GetNbinsY() )
816 {
817 printf ("Bad size of ChargedTracksInPIN! Must be: Nz x Nphi = %d x %d\n"
818 "but size of ChargedTracksInPIN is: %d x %d\n",
819 fNz,fNphi,fChargedTracksInPIN.GetNbinsX(),fChargedTracksInPIN.GetNbinsY());
820 exit(1);
821 }
822
823 if( NULL!=Noise && (fNz!=Noise->GetNbinsX() || fNphi!=Noise->GetNbinsX()) )
824 {
825 printf ("Bad size of Noise! Must be: Nz x Nphi = %d x %d\n"
826 "but size of Noise is: %d x %d\n",
827 fNz,fNphi,fChargedTracksInPIN.GetNbinsX(),fChargedTracksInPIN.GetNbinsY());
828 exit(1);
829 }
830
831 if( NULL!=Stochastic && (fNz!=Stochastic->GetNbinsX() || fNphi!=Stochastic->GetNbinsX()) )
832 {
833 printf ("Bad size of Stochastic! Must be: Nz x Nphi = %d x %d\n"
834 "but size of Stochastic is: %d x %d\n",
835 fNz,fNphi,fChargedTracksInPIN.GetNbinsX(),fChargedTracksInPIN.GetNbinsY());
836 exit(1);
837 }
838
839 if( NULL!=Calibration && (fNz!=Calibration->GetNbinsX() || fNphi!=Calibration->GetNbinsX()) )
840 {
841 printf ("Bad size of Calibration! Must be: Nz x Nphi = %d x %d\n"
842 "but size of Calibration is: %d x %d\n",
843 fNz,fNphi,fChargedTracksInPIN.GetNbinsX(),fChargedTracksInPIN.GetNbinsY());
844 exit(1);
845 }
846
847 ////////////////////
848 // Do distortion! //
849 ////////////////////
850
851 for( int y=0; y<fNphi; y++ )
852 for( int x=0; x<fNz; x++ )
853 {
854 const int n = fCellEnergy.GetBin(x,y); // Bin number
855 static TRandom r;
856
857 Float_t E_old=fCellEnergy.GetBinContent(n), E_new=E_old;
858
859 if( NULL!=Stochastic )
860 E_new = r.Gaus(E_old,sqrt(E_old)*GetDistortedValue(Stochastic,n));
861
862 if( NULL!=Calibration )
863 E_new *= GetDistortedValue(Calibration,n);
864
865 if( NULL!=Noise )
866 E_new += GetDistortedValue(Noise,n);
867
868 fCellEnergy.SetBinContent(n,E_new);
869 }
870}
871
872////////////////////////////////////////////////////////////////////////////////
873
874TH2F* AliPHOSCradle::CreateHistForDistortion(const char *name, const char *title,
875 Int_t Nx, Int_t Ny,
876 Float_t MU_mu, Float_t MU_sigma,
877 Float_t SIGMA_mu, Float_t SIGMA_sigma)
878{
879// Create (new TH2F(...)) histogram with information (for every bin) that will
880// be used for VALUE creation.
881// Two values will be created for each bin:
882// MU = TRandom::Gaus(MU_mu,MU_sigma)
883// and
884// SIGMA = TRandom::Gaus(SIGMA_mu,SIGMA_sigma)
885// The VALUE in a particluar bin will be equal
886// VALUE = TRandom::Gaus(MU,SIGMA)
887//
888// Do not forget to delete the histogram at the end of the work.
889
890 TH2F *h = new TH2F( name,title, Nx,1,Nx, Ny,1,Ny );
891 if( h==NULL )
892 {
893 Error("CreateHistForDistortion","Can not create the histogram");
894 exit(1);
895 }
896 h->SetDirectory(0);
897
898 for( int y=0; y<Ny; y++ )
899 for( int x=0; x<Nx; x++ )
900 {
901 const int n = h->GetBin(x,y);
902 h->SetBinContent(n,r.Gaus( MU_mu, MU_sigma));
903 h->SetBinError (n,r.Gaus(SIGMA_mu,SIGMA_sigma));
904 }
905
906 return h;
907}
908
909////////////////////////////////////////////////////////////////////////////////
910
911Float_t AliPHOSCradle::GetDistortedValue(const TH2F *h, UInt_t n)
912{
913 return r.Gaus(((TH2F*)h)->GetBinContent(n),n);
914}
915
916////////////////////////////////////////////////////////////////////////////////
917//______________________________________________________________________________
918
919#ifdef WIN32
920 #define common_for_event_storing COMMON_FOR_EVENT_STORING
921#else
922 #define common_for_event_storing common_for_event_storing_
923#endif
924
925extern "C" struct
926{
927 enum { crystals_matrix_amount_max=4, crystals_in_matrix_amount_max=40000 };
928
929 // Event-independent information
930 UShort_t crystals_matrix_amount_PHOS,
931 crystal_matrix_type,
932 amount_of_crystals_on_Z,
933 amount_of_crystals_on_PHI;
934 Float_t radius,
935 crystal_size,
936 crystal_length,
937 matrix_coordinate_Z [crystals_matrix_amount_max],
938 matrix_coordinate_PHI [crystals_matrix_amount_max];
939 UInt_t event_number;
940 UShort_t crystals_amount_with_amplitudes [crystals_matrix_amount_max],
941 crystals_amplitudes_Iad [crystals_matrix_amount_max]
942 [crystals_in_matrix_amount_max][2];
943} common_for_event_storing;
944
945// integer*4 crystals_amount_max,crystals_in_matrix_amount_max,
946// + crystals_matrix_amount_max
947// parameter (crystals_matrix_amount_max=4)
948// parameter (crystals_in_matrix_amount_max=40000)
949// parameter (crystals_amount_max =crystals_matrix_amount_max*
950// + crystals_in_matrix_amount_max)
951//
952// * All units are in GeV, cm, radian
953// real crystal_amplitudes_unit, radius_unit,
954// + crystal_size_unit, crystal_length_unit,
955// + matrix_coordinate_Z_unit, matrix_coordinate_PHI_unit
956// integer crystal_amplitudes_in_units_min
957// parameter (crystal_amplitudes_in_units_min = 1)
958// parameter (crystal_amplitudes_unit = 0.001 ) ! 1.0 MeV
959// parameter (radius_unit = 0.1 ) ! 0.1 cm
960// parameter (crystal_size_unit = 0.01 ) ! 0.01 cm
961// parameter (crystal_length_unit = 0.01 ) ! 0.01 cm
962// parameter (matrix_coordinate_Z_unit = 0.1 ) ! 0.1 cm
963// parameter (matrix_coordinate_PHI_unit = 1e-4 ) ! 1e-4 radian
964//
965// integer*2 crystals_matrix_amount_PHOS, crystal_matrix_type,
966// + amount_of_crystals_on_Z, amount_of_crystals_on_PHI,
967// + crystals_amount_with_amplitudes, crystals_amplitudes_Iad
968// integer*4 event_number
969//
970// real radius, crystal_size, crystal_length,
971// + matrix_coordinate_Z, matrix_coordinate_PHI
972//
973// real crystals_amplitudes, crystals_energy_total
974// integer event_file_unit_number
975//
976// common /common_for_event_storing/
977// + ! Event-independent information
978// + crystals_matrix_amount_PHOS,
979// + crystal_matrix_type,
980// + amount_of_crystals_on_Z,
981// + amount_of_crystals_on_PHI,
982// + radius,
983// + crystal_size,
984// + crystal_length,
985// + matrix_coordinate_Z (crystals_matrix_amount_max),
986// + matrix_coordinate_PHI (crystals_matrix_amount_max),
987// +
988// + ! Event-dependent information
989// + event_number,
990// + crystals_amount_with_amplitudes
991// + (crystals_matrix_amount_max),
992// + crystals_amplitudes_Iad (2,crystals_in_matrix_amount_max,
993// + crystals_matrix_amount_max),
994// +
995// + ! These information don't store in data file
996// + crystals_amplitudes (crystals_amount_max),
997// + crystals_energy_total,
998// + event_file_unit_number
999
1000
1001// parameter (NGp=1000,nsps=10,nvertmax=1000)
1002// COMMON /GAMMA/KG,MW(ngp),ID(ngp),JD(ngp),E(ngp),E4(ngp),
1003// , XW(ngp),YW(ngp),ES(nsps,ngp),ET(nsps,ngp),ISsd(ngp),
1004// , IGDEV(ngp),ZGDEV(ngp),sigexy(3,ngp),Emimx(2,nsps,ngp),
1005// , kgfix,igfix(ngp),cgfix(3,ngp),sgfix(3,ngp),hiw(ngp),
1006// , wsw(nsps,ngp),h1w(ngp),h0w(ngp),raxay(5,ngp),
1007// , sigmaes0(nsps,ngp),dispeces(nsps,ngp),
1008// , igamvert(ngp)
1009
1010
1011#ifdef WIN32
1012#define rcgamma RCGAMMA
1013#else
1014#define rcgamma rcgamma_
1015#endif
1016
1017extern "C" struct
1018{
1019 enum {NGP=1000, nsps=10, nvertmax=1000};
1020 int recons_gammas_amount, mw[NGP],ID[NGP],JD[NGP];
1021 float E[NGP], E4[NGP], XW[NGP], YW[NGP], ES[NGP][nsps],ET[NGP][nsps],ISsd[NGP],
1022 igdev[NGP],Zgdev[NGP];
1023// sigexy(3,ngp),Emimx(2,nsps,ngp),
1024// , kgfix,igfix(ngp),cgfix(3,ngp),sgfix(3,ngp),hiw(ngp),
1025// , wsw(nsps,ngp),h1w(ngp),h0w(ngp),raxay(5,ngp),
1026// , sigmaes0(nsps,ngp),dispeces(nsps,ngp),
1027// , igamvert(ngp)
1028} rcgamma;
1029
1030#ifdef WIN32
1031#define reconsfirst RECONSFIRST
1032#define type_of_call _stdcall
1033#else
1034#define reconsfirst reconsfirst_
1035#define type_of_call
1036#endif
1037
1038extern "C" void type_of_call reconsfirst(const float &,const float &);
1039
1040void AliPHOSCradle::Reconstruction(Float_t signal_step, UInt_t min_signal_reject)
1041{
1042// Call of PHOS reconstruction program.
1043// signal_step=0.001 GeV (1MeV)
1044// min_signal_reject = 15 or 30 MeV
1045
bc9ab547 1046
fe4da5cc 1047 common_for_event_storing.event_number = 0; // We do not know event number?
1048 common_for_event_storing.crystals_matrix_amount_PHOS = 1;
1049 common_for_event_storing.crystal_matrix_type = 1; // 1 - rectangular
1050 common_for_event_storing.amount_of_crystals_on_Z = fNz;
1051 common_for_event_storing.amount_of_crystals_on_PHI = fNphi;
1052
1053 common_for_event_storing.radius = fRadius;
1054 common_for_event_storing.crystal_size = GetCellSideSize();
1055 common_for_event_storing.crystal_length = fCrystalLength;
1056
1057 common_for_event_storing.matrix_coordinate_Z [0] = 0;
1058 common_for_event_storing.matrix_coordinate_PHI [0] = fPhi;
1059
1060 #define k common_for_event_storing.crystals_amount_with_amplitudes[0]
1061 k=0;
1062
1063 for( int y=0; y<fNphi; y++ )
1064 for( int x=0; x<fNz; x++ )
1065 {
1066 UInt_t n = fCellEnergy.GetBin(x,y);
1067 UInt_t signal = (int) (fCellEnergy.GetBinContent(n)/signal_step);
1068 if( signal>=min_signal_reject )
1069 {
1070 common_for_event_storing.crystals_amplitudes_Iad[0][k][0] = signal;
1071 common_for_event_storing.crystals_amplitudes_Iad[0][k][1] = x + y*fNz;
1072 k++;
1073 }
1074 }
1075 #undef k
1076
1077 GetGammasReconstructed().Delete();
1078 GetGammasReconstructed().Compress();
1079
1080 const float stochastic_term = 0.03, // per cents over sqrt(E); E in GeV
1081 electronic_noise = 0.01; // GeV
1082 reconsfirst(stochastic_term,electronic_noise); // Call of reconstruction program.
fe4da5cc 1083
1084 for( int i=0; i<rcgamma.recons_gammas_amount; i++ )
1085 {
1086// new (GetGammasReconstructed().UncheckedAt(i) ) AliPHOSgamma;
1087// AliPHOSgamma &g = *(AliPHOSgamma*)(GetGammasReconstructed().UncheckedAt(i));
1088
1089 AliPHOSgamma *gggg = new AliPHOSgamma;
1090 if( NULL==gggg )
1091 {
1092 Error("Reconstruction","Can not create AliPHOSgamma");
1093 exit(1);
1094 }
1095
1096 GetGammasReconstructed().Add(gggg);
1097 AliPHOSgamma &g=*gggg;
1098
1099 Float_t thetta, alpha, betta, R=fRadius+rcgamma.Zgdev[i]/10;
1100
1101 g.fX = rcgamma.YW[i]/10;
fe4da5cc 1102 g.fY = rcgamma.XW[i]/10;
fe4da5cc 1103 g.fE = rcgamma.E [i];
fe4da5cc 1104
1105 thetta = atan(g.fX/R);
1106
1107 alpha = atan(g.fY/R);
bc9ab547 1108 betta = fPhi/180*TMath::Pi() + alpha;
fe4da5cc 1109
1110 g.fPx = g.fE * cos(thetta) * cos(betta);
1111 g.fPy = g.fE * cos(thetta) * sin(betta);
1112 g.fPz = g.fE * sin(thetta);
1113 }
1114}
1115
1116//______________________________________________________________________________
1117//______________________________________________________________________________
1118//______________________________________________________________________________
1119//______________________________________________________________________________
1120//______________________________________________________________________________
1121
1122ClassImp(AliPHOSgamma)
1123
1124//______________________________________________________________________________
1125
1126void AliPHOSgamma::Print(Option_t *)
1127{
1128 float mass = fE*fE - fPx*fPx - fPy*fPy - fPz*fPz;
1129
1130 if( mass>=0 )
1131 mass = sqrt( mass);
1132 else
1133 mass = -sqrt(-mass);
1134
bc9ab547 1135 printf("XY=(%+7.2f,%+7.2f) (%+7.2f,%+7.2f,%+7.2f;%7.2f) mass=%8.4f Ipart=%2d\n",
1136 fX,fY,fPx,fPy,fPz,fE,mass,fIpart);
fe4da5cc 1137}
1138
1139//______________________________________________________________________________
1140
1141AliPHOSgamma &AliPHOSgamma::operator=(const AliPHOSgamma &g)
1142{
1143 fX = g.fX;
fe4da5cc 1144 fY = g.fY;
fe4da5cc 1145 fE = g.fE;
fe4da5cc 1146 fPx = g.fPx;
1147 fPy = g.fPy;
1148 fPz = g.fPz;
bc9ab547 1149 fIpart = g.fIpart;
fe4da5cc 1150
1151 return *this;
1152}
1153
1154//______________________________________________________________________________
1155//______________________________________________________________________________
1156//______________________________________________________________________________
1157//______________________________________________________________________________
1158//______________________________________________________________________________
1159
1160ClassImp(AliPHOShit)
1161
1162//______________________________________________________________________________
1163
1164AliPHOShit::AliPHOShit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits):
1165AliHit(shunt, track)
1166{
1167 Int_t i;
1168 for (i=0;i<5;i++) fVolume[i] = vol[i];
1169 fX = hits[0];
1170 fY = hits[1];
1171 fZ = hits[2];
1172 fELOS = hits[3];
1173}
1174
1175//______________________________________________________________________________