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