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