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