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