Completing removal of old PHOS code
[u/mrichter/AliRoot.git] / PHOS / AliPHOS.cxx
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$
18 Revision 1.10  1999/12/09 14:57:51  fca
19 Removal of obsolete PHOS code
20
21 Revision 1.9  1999/11/08 07:12:31  fca
22 Minor corrections thanks to I.Hrivnacova
23
24 Revision 1.8  1999/09/29 09:24:23  fca
25 Introduction of the Copyright and cvs Log
26
27 */
28
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"
40 #include "TMath.h"
41
42 // --- Standard library ---
43 #include <stdio.h>
44 #include <string.h>
45 #include <stdlib.h>
46 #include <iostream.h>
47
48 // --- galice header files ---
49 #include "AliPHOS.h"
50 #include "AliRun.h"
51
52 //______________________________________________________________________________
53
54
55 ClassImp(AliPHOS)
56
57 //______________________________________________________________________________
58
59 AliPHOS::~AliPHOS(void)
60 {
61   delete fHits;                 // 28.12.1998
62   delete fTreePHOS;             // 28.12.1998
63   fCradles->Delete();
64   delete fCradles;
65 }
66
67 //______________________________________________________________________________
68
69 AliPHOS::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
87 AliPHOS::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 /*
96 <img src="picts/aliphos.gif">
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
117 void 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.;
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
164 void 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 //___________________________________________
171 void 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 //___________________________________________
207 void AliPHOS::CreateMaterials()
208 {
209 // *** DEFINITION OF AVAILABLE PHOS MATERIALS *** 
210
211 // CALLED BY : PHOS_MEDIA 
212 // ORIGIN    : NICK VAN EIJNDHOVEN 
213
214
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
255     Int_t *idtmed = fIdtmed->GetArray()-699;
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
269     AliMedium(0, "PHOS Xtal    $", 0, 1, ISXFLD, SXMGMX, 10., .1, .1, .1, .1);
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);
278
279 // --- Generate explicitly delta rays in the steel cover --- 
280     gMC->Gstpar(idtmed[704], "LOSS", 3.);
281     gMC->Gstpar(idtmed[704], "DRAY", 1.);
282 // --- and in aluminium parts --- 
283     gMC->Gstpar(idtmed[701], "LOSS", 3.);
284     gMC->Gstpar(idtmed[701], "DRAY", 1.);
285 }
286  
287 //______________________________________________________________________________
288
289 void 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             (),
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
317 Int_t AliPHOS::DistancetoPrimitive(Int_t , Int_t )
318 {
319    return 9999;
320 }
321  
322 //___________________________________________
323 void 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
340 void 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
391 void 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   {
407     Error("SetTreeAddress","Can not find Tree \"%s\"\n",fTreeName.Data());
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
423 AliPHOSCradle *AliPHOS::GetCradleOfTheParticle(const TVector3 &p,const TVector3 &v) const
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;
434     const float d = cradle->GetRadius();
435     cradle->GetXY(p,v,d,x,y,l);
436
437     if( l>0 && TMath::Abs(x)<cradle->GetNz  ()*cradle->GetCellSideSize()/2 
438             && TMath::Abs(y)<cradle->GetNphi()*cradle->GetCellSideSize()/2 )
439       return cradle;
440   }
441
442   return NULL;
443 }
444
445 //______________________________________________________________________________
446
447 void 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
458 void 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
468 void 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
478 void AliPHOS::FinishRun(void)
479 {
480 }
481
482 //______________________________________________________________________________
483
484 void 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 //______________________________________________________________________________
500 void 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 //______________________________________________________________________________
515 void 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 //______________________________________________________________________________
530 void AliPHOS::SetRadius(Float_t radius)
531 {
532    PHOSradius=radius;
533 }
534
535 //______________________________________________________________________________
536 void 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 //______________________________________________________________________________
544 void AliPHOS::SetCradleA(Float_t angle)
545 {
546    PHOScradlesA=angle;
547 }
548
549 //______________________________________________________________________________
550 void 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 //______________________________________________________________________________
565 void 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 //______________________________________________________________________________
573 void 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 //______________________________________________________________________________
581 void 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
589 ClassImp(AliPHOSCradle)
590
591 //______________________________________________________________________________
592
593 AliPHOSCradle::AliPHOSCradle(void) {}
594
595 //______________________________________________________________________________
596
597 AliPHOSCradle::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             ,
605                               int   Nz                 ,
606                               int   Nphi               ,
607                               float Angle              ) :
608     fGeometry                   (Geometry),
609 //  fCellEnergy                 (),
610 //  fChargedTracksInPIN         (),
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),
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
630 AliPHOSCradle::~AliPHOSCradle(void)        // 28.12.1998
631 {
632   fGammasReconstructed.Delete();
633   fParticles          .Delete();
634 }
635
636 //______________________________________________________________________________
637
638 void 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
649 }
650
651 //______________________________________________________________________________
652
653 void AliPHOSCradle::GetXY(const TVector3 &p,const TVector3 &v,float R,float &x,float &y,float &l) const
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(......);
666 //   TVector3 p(....), v(....);
667 //   Float_t x,y,l;
668 //   cradle.GetXY(p,v,x,y,l);
669 //   if( l<0 || TMath::Abs(x)>cradle.GetNz()  *cradle.GetCellSideSize()/2
670 //           || TMath::Abs(y)>cradle.GetNphi()*cradle.GetCellSideSize()/2 )
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
679   double phi = fPhi/180*TMath::Pi();
680   TVector3        n1(   0.0      ,   0.0      , 1.0 ),   // Z direction (X)
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,
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
692   
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));
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
717 void 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
727   printf("AliPHOSCradle:  Nz=%d  Nphi=%d, fPhi=%f, E=%g\n",fNz,fNphi,fPhi,
728        cr->fCellEnergy.GetSumOfWeights());
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
769 void 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
846 TH2F* 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
883 Float_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
897 /* extern "C" */ struct
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
989 /* extern "C" */ struct
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
1002 /*
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
1011 extern "C" void type_of_call reconsfirst(const float &,const float &);
1012 */
1013
1014 void 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
1020
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
1056 //  reconsfirst(stochastic_term,electronic_noise); // Call of reconstruction program.
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;
1076     g.fY      = rcgamma.XW[i]/10;
1077     g.fE      = rcgamma.E [i];
1078
1079     thetta      = atan(g.fX/R);
1080
1081     alpha = atan(g.fY/R);
1082     betta = fPhi/180*TMath::Pi() + alpha;
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
1096 ClassImp(AliPHOSgamma)
1097
1098 //______________________________________________________________________________
1099
1100 void 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
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);
1111 }
1112
1113 //______________________________________________________________________________
1114
1115 AliPHOSgamma &AliPHOSgamma::operator=(const AliPHOSgamma &g)
1116 {
1117   fX           = g.fX;
1118   fY           = g.fY;
1119   fE           = g.fE;
1120   fPx          = g.fPx;
1121   fPy          = g.fPy;
1122   fPz          = g.fPz;
1123   fIpart       = g.fIpart;
1124
1125   return *this;
1126 }
1127
1128 //______________________________________________________________________________
1129 //______________________________________________________________________________
1130 //______________________________________________________________________________
1131 //______________________________________________________________________________
1132 //______________________________________________________________________________
1133
1134 ClassImp(AliPHOShit)
1135
1136 //______________________________________________________________________________
1137
1138 AliPHOShit::AliPHOShit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits):
1139 AliHit(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 //______________________________________________________________________________