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