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