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