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