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