ITS/UPGRADE/AliITSupgrade.cxx

   1 // **************************************************************************
   2 // * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3 // *                                                                        *
   4 // * Author: The ALICE Off-line Project.                                    *
   5 // * Contributors are mentioned in the code where appropriate.              *
   6 // *                                                                        *
   7 // * Permission to use, copy, modify and distribute this software and its   *
   8 // * documentation strictly for non-commercial purposes is hereby granted   *
   9 // * without fee, provided that the above copyright notice appears in all   *
  10 // * copies and that both the copyright notice and this permission notice   *
  11 // * appear in the supporting documentation. The authors make no claims     *
  12 // * about the suitability of this software for any purpose. It is          *
  13 // * provided "as is" without express or implied warranty.                  *
  14 // **************************************************************************
  15
  16 /* $Id$ */
  17
  18 #include <TArrayD.h>            //new constructor
  19 #include <TFile.h>
  20 #include <TGeoManager.h>
  21 #include <TGeoVolume.h>         //CreateGeometry()
  22 #include <TGeoMatrix.h>
  23 #include <TVirtualMC.h>         //->gMC in StepManager
  24 #include <TPDGCode.h>           //StepHistory
  25 #include <TClonesArray.h>
  26 #include <TGeoGlobalMagField.h>
  27 #include "AliRun.h"             //CreateMaterials()
  28 #include "AliMC.h"             //CreateMaterials()
  29 #include "AliMagF.h"            //CreateMaterials()
  30 #include "AliCDBManager.h"      //CreateMaterials()
  31 #include "AliCDBEntry.h"        //CreateMaterials()
  32 #include "AliITSupgrade.h"      //class header
  33 #include "AliITShit.h"           // StepManager()
  34 #include "AliITSDigitUpgrade.h"
  35 #include "AliTrackReference.h"   // StepManager()
  36 #include "AliITSDetTypeSim.h"
  37 ClassImp(AliITSupgrade)
  38 //__________________________________________________________________________________________________
  39   AliITSupgrade::AliITSupgrade():
  40     AliITS(),
  41     fWidths(0),
  42     fRadii(0),
  43     fNSectors(20),
  44     fRadiiCu(0),
  45     fWidthsCu(0),
  46     fCopper(0),
  47     fBeampipe(0),
  48     fRadiusBP(0),
  49     fWidthBP(0),
  50     fHalfLengthBP(0),
  51     fNlayers(0),
  52     fHalfLength(0),
  53     fSdigits(0),
  54     fDigits(0),
  55     fSegmentation(0x0)
  56 {
  57   //
  58   // default ctor
  59   //
  60 }
  61
  62 //__________________________________________________________________________________________________
  63 AliITSupgrade::AliITSupgrade(const char *name,const char *title, Bool_t isBeamPipe):
  64   AliITS(name,title),
  65   fWidths(0),
  66   fRadii(0),
  67   fNSectors(20),
  68   fRadiiCu(0),
  69   fWidthsCu(0),
  70   fCopper(0),
  71   fBeampipe(isBeamPipe),
  72   fRadiusBP(0),
  73   fWidthBP(0),
  74   fHalfLengthBP(0),
  75   fNlayers(0),
  76   fHalfLength(0),
  77   fSdigits(0),
  78   fDigits(0),
  79   fSegmentation(0x0)
  80 {
  81   fNlayers = 6;
  82   fWidths.Set(fNlayers);
  83   fRadii.Set(fNlayers);
  84   fRadiiCu.Set(fNlayers);
  85   fWidthsCu.Set(fNlayers);
  86   fCopper.Set(fNlayers);
  87   fHalfLength.Set(fNlayers);
  88
  89   // Default values are used in order to simulate the standard ITS material budget (see The ALICE Collaboration et al 2008 JINST 3 S08002 - Figure 3.2).
  90   // The cell segmentation is chosen to achieve the tracking resolution as described in Table 3.2 of The ALICE Collaboration et al 2008 JINST 3 S08002,
  91   // apart from SPD0 where the 9 um value is considered : resolution*sqrt(12)
  92   // Cilinder lenghts are chosen according to Table 1 of the Alignment paper : http://iopscience.iop.org/1748-0221/5/03/P03003
  93
  94   Double_t xsz[6]={31.18*1e-04,41.6*1e-04,131.6*1e-04,131.6*1e-04,69.3*1e-04,69.3*1e-04};
  95   Double_t zsz[6]={416*1e-04,416*1e-04,97*1e-04,97*1e-04,2875*1e-04,2875*1e-04};
  96   Double_t Lhalf[6]={14.1,14.1,22.2,29.7,43.1,48.9};
  97   Double_t r[6]={4.,7.6,14.9,23.8,39.1,43.6};
  98   Double_t thick[6]={75.*1e-04,150.*1e-04,150.*1e-04,150.*1e-04,150.*1e-04,150.*1e-04};
  99
 100   Int_t  nlayers =6;
 101   TArrayD xsizeSi(nlayers);
 102   TArrayD zsizeSi(nlayers);
 103
 104   // adding beam pipe upgrade
 105   if(isBeamPipe){
 106     fRadiusBP = 2.9;
 107     fWidthBP = 0.08;
 108     fHalfLengthBP = Lhalf[0];
 109   }
 110
 111   // setting the geonetry parameters and the segmentation
 112   Int_t npixHalf[6], npixR[6];
 113   Double_t HalfL[6], xszInt[6];
 114   Int_t c[6]={1,1,1,1,1,1};
 115
 116   for(Int_t i=0; i<nlayers; i++){
 117     // recalc length
 118     npixHalf[i]=(Int_t)(Lhalf[i]/zsz[i]);
 119     HalfL[i]=npixHalf[i]*zsz[i];
 120     // recalc segmentation
 121     npixR[i] = (Int_t)(2*TMath::Pi()*r[i]/xsz[i]);
 122     xszInt[i]= 2*TMath::Pi()*r[i]/npixR[i];
 123     xsizeSi.AddAt(xszInt[i],i);
 124     zsizeSi.AddAt(zsz[i],i);
 125
 126     fHalfLength.AddAt(HalfL[i],i);
 127
 128     fRadii.AddAt(r[i],i);
 129     fWidths.AddAt(thick[i],i);
 130
 131     fRadiiCu.AddAt(r[i]+thick[i],i);
 132     fWidthsCu.AddAt(0.015,i);//micron
 133     fCopper.AddAt(c[i],i);
 134   }
 135
 136   SetFullSegmentation(xsizeSi,zsizeSi);
 137 }
 138 //__________________________________________________________________________________________________
 139 AliITSupgrade::AliITSupgrade(const char *name,const char *title, TArrayD widths, TArrayD radii,TArrayD halfLengths, TArrayD radiiCu, TArrayD widthsCu, TArrayS copper,Bool_t bp,Double_t radiusBP, Double_t widthBP, Double_t halfLengthBP):
 140   AliITS(name,title),
 141   fWidths(0),
 142   fRadii(0),
 143   fNSectors(20),
 144   fRadiiCu(radiiCu),
 145   fWidthsCu(widthsCu),
 146   fCopper(copper),
 147   fBeampipe(bp),
 148   fRadiusBP(0),
 149   fWidthBP(0),
 150   fHalfLengthBP(0),
 151   fNlayers(widths.GetSize()),
 152   fHalfLength(halfLengths),
 153   fSdigits(0),
 154   fDigits(0),
 155   fSegmentation(0x0)
 156 {
 157
 158   //
 159   // constructor
 160   for(Int_t i=0;i<fNlayers;i++){
 161     fWidths.Set(i+1);
 162     fWidths.AddAt(widths.At(i),i);
 163     fRadii.Set(i+1);
 164     fRadii.AddAt(radii.At(i),i);
 165     AliDebug(1,"Creating Volume");
 166   }
 167
 168   if(bp){
 169     fRadiusBP=radiusBP;
 170     fWidthBP=widthBP;
 171     fHalfLengthBP=halfLengthBP;
 172   }
 173   Init();
 174
 175 }
 176
 177 AliITSupgrade::~AliITSupgrade(){
 178
 179   if(fSdigits) {fSdigits->Delete(); delete fSdigits;}
 180   if(fDigits)  {fDigits->Delete(); delete fDigits;}
 181 }
 182
 183
 184
 185 //_________________________________________________________________________________________________
 186 void AliITSupgrade::AddAlignableVolumes()const
 187 {
 188   //not needed
 189   return;
 190 }
 191
 192 //__________________________________________________________________________________________________
 193
 194 void AliITSupgrade::CreateMaterials()
 195 {
 196   //
 197   // Definition of ITS materials
 198   //
 199
 200   AliInfo("Start ITS materials");
 201   //data from PDG booklet 2002     density [gr/cm^3] rad len [cm] abs len [cm]
 202   Float_t   aAir[4]={12,14,16,36} ,   zAir[4]={6,7,8,18} ,   wAir[4]={0.000124,0.755267,0.231781,0.012827} , dAir=0.00120479; Int_t nAir=4;//mixture
 203   Float_t       aBe = 9.012 ,          zBe   = 4 ,            dBe   =  1.848    ,   radBe   =  65.19/dBe , absBe   = 77.8/dBe  ; // UPGRADE -> beryllium beampipe
 204   Float_t       aSi = 28.085 ,         zSi   = 14 ,           dSi   =  2.329    ,   radSi   =  21.82/dSi , absSi   = 108.4/dSi  ; // UPGRADE -> Si tube
 205   Float_t aCu = 63.54, zCu = 29 , dCu= 8.96 , radCu = 12.86/dCu, absCu= 137.3/dCu;//Upgrade -> Copper Tube
 206
 207   Int_t   matId=0;                           //tmp material id number
 208   Int_t   unsens = 0, sens=1;                //sensitive or unsensitive medium
 209   Int_t   itgfld = 3;                        //type of field intergration 0 no field -1 user in guswim 1 Runge Kutta 2 helix 3 const field along z
 210   Float_t maxfld = 5.;                       //max field value
 211   Float_t tmaxfd = -10.0;                    //max deflection angle due to magnetic field in one step
 212   Float_t deemax = - 0.2;                    //max fractional energy loss in one step
 213   Float_t stemax = - 0.1;                    //max step allowed [cm]
 214   Float_t epsil  =   0.001;                  //abs tracking precision [cm]
 215   Float_t stmin  = - 0.001;                  //min step size [cm] in continius process transport, negative value: choose it automatically
 216   Float_t tmaxfdSi = 0.1; // .10000E+01; // Degree
 217   Float_t stemaxSi = 0.0075; //  .10000E+01; // cm
 218   Float_t deemaxSi = 0.1; // 0.30000E-02; // Fraction of particle's energy 0<deemax<=1
 219   Float_t epsilSi  = 1.0E-4;// .10000E+01;
 220   Float_t stminSi  = 0.0; // cm "Default value used"
 221
 222   Float_t epsilBe  = .001;    // Tracking precision,
 223   Float_t stemaxBe = -0.01;   // Maximum displacement for multiple scat
 224   Float_t tmaxfdBe = -20.;    // Maximum angle due to field deflection
 225   Float_t deemaxBe = -.3;     // Maximum fractional energy loss, DLS
 226   Float_t stminBe  = -.8;
 227   Int_t   isxfld = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->Integ();//from CreateMaterials in STRUCT/AliPIPEv3.cxx
 228   Float_t sxmgmx = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->Max();//from CreateMaterials in STRUCT/AliPIPEv3.cxx
 229
 230
 231   AliMixture(++matId,"UpgradeAir"  ,aAir  ,zAir  ,dAir  ,nAir  ,wAir  );
 232   AliMedium(kAir  ,"UpgradeAir"  ,matId, unsens, itgfld, maxfld, tmaxfd, stemax, deemax, epsil, stmin);
 233
 234   AliMaterial(++matId,"UpgradeSi"  ,aSi  ,zSi  ,dSi  ,radSi  ,absSi  );
 235   AliMedium(kSi  ,"UpgradeSi"  , matId, sens, itgfld, maxfld, tmaxfdSi, stemaxSi, deemaxSi, epsilSi, stminSi);
 236
 237   AliMaterial(++matId,"UpgradeBe"  ,aBe  ,zBe  ,dBe  ,radBe  ,absBe  );
 238   AliMedium(kBe  ,"UpgradeBe"  , matId, unsens, isxfld, sxmgmx, tmaxfdBe, stemaxBe, deemaxBe, epsilBe, stminBe);
 239
 240   AliMaterial(++matId, "UpgradeCu", aCu, zCu, dCu, radCu, absCu);
 241   AliMedium(kCu, "UpgradeCu", matId, unsens, itgfld, maxfld, tmaxfd, stemax, deemax, epsil, stmin);
 242
 243   AliInfo("End ITS materials");
 244
 245 }//void AliITS::CreateMaterials()
 246 //__________________________________________________________________________________________________
 247 void AliITSupgrade::CreateGeometry()
 248 {
 249   //
 250   //Creates detailed geometry simulation (currently GEANT volumes tree)
 251   //
 252   AliInfo("Start ITS upgrade preliminary version building");
 253   if(!gMC->IsRootGeometrySupported()) return;
 254   TGeoVolumeAssembly *vol= CreateVol();
 255   gGeoManager->GetVolume("ALIC")->AddNode(vol,0);
 256   AliInfo("Stop ITS upgrade preliminary version building");
 257   PrintSummary();
 258 }
 259 //__________________________________________________________________________________________________
 260 void AliITSupgrade::StepManager()
 261 {
 262   // Full Step Manager.
 263   // Arguments: none
 264   //   Returns: none
 265   //  StepHistory(); return; //uncomment to print tracks history
 266   //  StepCount(); return;   //uncomment to count photons
 267
 268   if(!fSegmentation) AliFatal("No segmentation available");
 269
 270   if(!(this->IsActive())) return;
 271   if(!(gMC->TrackCharge())) return;
 272   TString volumeName=gMC->CurrentVolName();
 273   if(volumeName.Contains("Be")) return;
 274   if(volumeName.Contains("Cu")) return;
 275   if(gMC->IsTrackExiting()) {
 276     AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber(), AliTrackReference::kITS);
 277   } // if Outer ITS mother Volume
 278
 279   static TLorentzVector position, momentum; // Saves on calls to construtors
 280   static AliITShit hit;// Saves on calls to constructors
 281
 282
 283   Int_t  status = 0;
 284
 285   //
 286   // Track status
 287   if(gMC->IsTrackInside())      status +=  1;
 288   if(gMC->IsTrackEntering())    status +=  2;
 289   if(gMC->IsTrackExiting())     status +=  4;
 290   if(gMC->IsTrackOut())         status +=  8;
 291   if(gMC->IsTrackDisappeared()) status += 16;
 292   if(gMC->IsTrackStop())        status += 32;
 293   if(gMC->IsTrackAlive())       status += 64;
 294
 295   //
 296   // Fill hit structure.
 297   //
 298   Int_t copy=-1;
 299   gMC->CurrentVolID(copy);
 300
 301   volumeName.Remove(0,12);          // remove letters to get the layer number
 302   hit.SetModule(fSegmentation->GetIdIndex(volumeName.Atoi(),copy)); // layer and sector information are together in the IdIndex (if copy=0 the idIndex is the layer));
 303   hit.SetTrack(gAlice->GetMCApp()->GetCurrentTrackNumber());
 304
 305   gMC->TrackPosition(position);
 306   gMC->TrackMomentum(momentum);
 307   hit.SetPosition(position);
 308
 309   hit.SetTime(gMC->TrackTime());
 310   hit.SetMomentum(momentum);
 311   hit.SetStatus(status);
 312   hit.SetEdep(gMC->Edep());
 313   hit.SetShunt(GetIshunt());
 314   if(gMC->IsTrackEntering()){
 315     hit.SetStartPosition(position);
 316     hit.SetStartTime(gMC->TrackTime());
 317     hit.SetStartStatus(status);
 318     return; // don't save entering hit.
 319   }
 320   // Fill hit structure with this new hit.
 321
 322   new((*fHits)[fNhits++]) AliITShit(hit); // Use Copy Construtor.
 323   // Save old position... for next hit.
 324   hit.SetStartPosition(position);
 325   hit.SetStartTime(gMC->TrackTime());
 326   hit.SetStartStatus(status);
 327   return;
 328
 329 }
 330 //__________________________________________________________________________________________________
 331 TGeoVolumeAssembly * AliITSupgrade::CreateVol()
 332 {
 333
 334    TGeoVolumeAssembly *vol = new TGeoVolumeAssembly("ITSupgrade");
 335   TGeoMedium *si   =gGeoManager->GetMedium("ITS_UpgradeSi");
 336   TGeoMedium *cu   =gGeoManager->GetMedium("ITS_UpgradeCu");
 337   TGeoMedium *be   =gGeoManager->GetMedium("ITS_UpgradeBe");
 338   for(Int_t ivol=0;ivol<fNlayers;ivol++){
 339
 340     if(fNSectors<1){
 341       TGeoVolume *layer=gGeoManager->MakeTube(Form("LayerSilicon%i",ivol),si   ,    fRadii.At(ivol)   ,   fRadii.At(ivol)+fWidths.At(ivol)  ,   fHalfLength.At(ivol)); //upgraded situation
 342
 343       TGeoVolume *layerCu=gGeoManager->MakeTube(Form("LayerCu%i",ivol),cu   ,    fRadiiCu.At(ivol)   ,   fRadiiCu.At(ivol)+fWidthsCu.At(ivol) ,  fHalfLength.At(ivol)  ); //upgraded situation
 344
 345       vol ->AddNode(layer,ivol);
 346       if(fCopper.At(ivol)){
 347         vol->AddNode(layerCu,ivol);
 348       }
 349
 350     }else{
 351
 352       TGeoVolume *layer = gGeoManager->MakeTubs(Form("LayerSilicon%i",ivol),si,  fRadii.At(ivol),   fRadii.At(ivol)+fWidths.At(ivol)  ,fHalfLength.At(ivol),0,(360./fNSectors));
 353       TGeoVolume *layerCu = gGeoManager->MakeTubs(Form("LayerCu%i",ivol),cu   ,    fRadiiCu.At(ivol)   ,   fRadiiCu.At(ivol)+fWidthsCu.At(ivol) ,  fHalfLength.At(ivol) ,0,(360./fNSectors));
 354
 355
 356       for(Int_t i=0;i<fNSectors;i++){
 357         TGeoRotation *rot1 = new TGeoRotation(" ",0.0,0.0,360.*i/fNSectors);//sector rotation
 358         TGeoRotation *rot2 = new TGeoRotation(" ",0.0,0.0,360.*i/fNSectors);//
 359         vol->AddNode(layer,i,rot1);
 360         if(fCopper.At(ivol)){
 361           vol->AddNode(layerCu,i,rot2);
 362         }
 363       }
 364     }
 365   }
 366
 367
 368   if(fBeampipe) {
 369     TGeoVolume *beampipe=gGeoManager->MakeTube("BeamPipe", be   ,    fRadiusBP   ,  fRadiusBP+ fWidthBP ,  fHalfLengthBP  ); //upgraded situation
 370      vol->AddNode(beampipe,0);
 371    }
 372   return vol;
 373
 374
 375 }
 376 //_________________________________________________________________________________________________
 377 void AliITSupgrade::SetFullSegmentation(TArrayD xsize,TArrayD zsize){
 378
 379   Bool_t Check=kFALSE;
 380   for(Int_t lay = 0; lay< xsize.GetSize(); lay++){
 381     Double_t arch = fRadii.At(lay)*(TMath::Pi()*2/fNSectors);
 382     Int_t nPixRPhi = (Int_t)(arch/xsize.At(lay));
 383     Int_t nPixZed = (Int_t)((2*fHalfLength.At(lay))/zsize.At(lay));
 384     if(nPixRPhi>9999)Check=kTRUE;
 385     if(nPixZed>99999)Check=kTRUE;
 386   }
 387   if(Check) AliFatal(" Segmentation is too small!! ");
 388   if(fSegmentation) fSegmentation->SetNSectors(fNSectors);
 389   TArrayD nSect(1);
 390   nSect.AddAt(fNSectors,0);
 391   TFile *file= TFile::Open("Segmentation.root","recreate");
 392   file->WriteObjectAny(&xsize,"TArrayD","CellSizeX");
 393   file->WriteObjectAny(&zsize,"TArrayD","CellSizeZ");
 394   file->WriteObjectAny(&nSect,"TArrayD","nSectors");
 395   file->Close();
 396 }
 397 //_________________________________________________________________________________________________
 398 void AliITSupgrade::StepHistory()
 399 {
 400   // This methode is invoked from StepManager() in order to print out
 401   static Int_t iStepN;
 402   const char *sParticle;
 403   switch(gMC->TrackPid()){
 404   case kProton:      sParticle="PROTON"    ;break;
 405   case kNeutron:     sParticle="neutron"   ;break;
 406   case kGamma:       sParticle="gamma"     ;break;
 407   case kPi0:         sParticle="Pi0"       ;break;
 408   case kPiPlus:      sParticle="Pi+"       ;break;
 409   case kPiMinus:     sParticle="Pi-"       ;break;
 410   case kElectron:    sParticle="electron"  ;break;
 411   default:           sParticle="not known" ;break;
 412   }
 413
 414   TString flag="fanny combination";
 415   if(gMC->IsTrackAlive())
 416     if(gMC->IsTrackEntering())      flag="enters to";
 417     else if(gMC->IsTrackExiting())  flag="exits from";
 418     else if(gMC->IsTrackInside())   flag="inside";
 419     else if(gMC->IsTrackStop())     flag="stopped in";
 420
 421   Int_t vid=0,copy=0;
 422   TString path=gMC->CurrentVolName(); path.Prepend("-");path.Prepend(gMC->CurrentVolOffName(1));//current volume and his mother are always there
 423   vid=gMC->CurrentVolOffID(2,copy);  if(vid) {path.Prepend("-");path.Prepend(gMC->VolName(vid));}
 424   vid=gMC->CurrentVolOffID(3,copy);  if(vid) {path.Prepend("-");path.Prepend(gMC->VolName(vid));}
 425
 426   AliDebug(10, Form("Step %i: %s (%i) %s %s m=%.6f GeV q=%.1f dEdX=%.4f Etot=%.4f",iStepN,sParticle,gMC->TrackPid(),flag.Data(),path.Data(),gMC->TrackMass(),gMC->TrackCharge(),gMC->Edep()*1e9,gMC->Etot()));
 427
 428   Double_t gMcTrackPos[3]; gMC->TrackPosition(gMcTrackPos[0],gMcTrackPos[1],gMcTrackPos[2]);
 429   Double_t  gMcTrackPosLoc[3]; gMC->Gmtod(gMcTrackPos,gMcTrackPosLoc,1);
 430   AliDebug(10,Form("gMC Track Position (MARS) x: %5.3lf, y: %5.3lf, z: %5.3lf (r: %5.3lf) ---> (LOC) x: %5.3f, y: %5.3f, z: %5.3f",gMcTrackPos[0],gMcTrackPos[1],gMcTrackPos[2],TMath::Sqrt(gMcTrackPos[0]*gMcTrackPos[0]+gMcTrackPos[1]*gMcTrackPos[1]+gMcTrackPos[2]*gMcTrackPos[2]),gMcTrackPosLoc[0],gMcTrackPosLoc[1],gMcTrackPosLoc[2]));
 431
 432
 433   AliDebug(10,Form("Step %i: tid=%i flags alive=%i disap=%i enter=%i exit=%i inside=%i out=%i stop=%i new=%i",
 434                    iStepN, gAlice->GetMCApp()->GetCurrentTrackNumber(),
 435                    gMC->IsTrackAlive(), gMC->IsTrackDisappeared(),gMC->IsTrackEntering(), gMC->IsTrackExiting(),
 436                    gMC->IsTrackInside(),gMC->IsTrackOut(),        gMC->IsTrackStop(),     gMC->IsNewTrack()));
 437
 438   Float_t a,z,den,rad,abs; a=z=den=rad=abs=-1;
 439   Int_t mid=gMC->CurrentMaterial(a,z,den,rad,abs);
 440   AliDebug(10, Form("Step %i: mid=%i a=%7.2f z=%7.2f den=%9.4f rad=%9.2f abs=%9.2f\n\n",iStepN,mid,a,z,den,rad,abs));
 441
 442   TArrayI proc;  gMC->StepProcesses(proc);
 443   AliDebug(1,"Processes in this step:");
 444   for ( int i = 0 ; i < proc.GetSize(); i++)
 445     {
 446       AliDebug(10, Form("%s",TMCProcessName[proc.At(i)]));
 447     }
 448   AliDebug(1,"End process list");
 449   iStepN++;
 450 }//StepHistory()
 451 //______________________________________________________
 452 void AliITSupgrade::Hits2SDigits(){
 453
 454   // Interface method ivoked from AliSimulation to create a list of sdigits corresponding to list of hits. Every hit generates one or more sdigits.
 455   // Arguments: none
 456   //   Returns: none
 457   AliDebug(1,"Start Hits2SDigits.");
 458
 459   if(!fSegmentation)fSegmentation=new AliITSsegmentationUpgrade();
 460
 461   if(!GetLoader()->TreeH()) {GetLoader()->LoadHits();}
 462   if(!GetLoader()->TreeS())
 463
 464     for(Int_t iEvt=0;iEvt < GetLoader()->GetRunLoader()->GetNumberOfEvents();iEvt++){//events loop
 465       GetLoader()->GetRunLoader()->GetEvent(iEvt);                          //get next event
 466       {
 467         GetLoader()->MakeTree("S");
 468         MakeBranch("S");
 469       }
 470       SetTreeAddress();
 471       Int_t nSdigit[10]; for(Int_t i=0;i<10;i++)  nSdigit[i] =0;
 472       for(Int_t iEnt=0;iEnt<GetLoader()->TreeH()->GetEntries();iEnt++){//prims loop
 473         GetLoader()->TreeH()->GetEntry(iEnt);
 474         Hit2SumDig(Hits(),SDigitsList(),nSdigit);//convert this hit to list of sdigits
 475       }//prims loop
 476       GetLoader()->TreeS()->Fill();
 477       GetLoader()->WriteSDigits("OVERWRITE");
 478       SDigitsReset();
 479     }//events loop
 480   GetLoader()->UnloadHits();
 481   GetLoader()->UnloadSDigits();
 482   AliDebug(1,"Stop Hits2SDigits.");
 483
 484 }
 485 //____________________________
 486 void AliITSupgrade::Hit2SumDig(TClonesArray *hits,TObjArray *pSDig, Int_t *nSdigit)
 487 {
 488   // Adds  sdigits of this hit to the list
 489   //   Returns: none
 490
 491   AliDebug(1,"Start Hit2SumDig");
 492
 493
 494   if(!fSegmentation){
 495     AliDebug(1,"Segmentation Not inizialized!!");
 496     return ;
 497   }
 498
 499   TClonesArray *pSdigList[10];
 500
 501   for(Int_t i=0;i<fNlayers;i++){
 502     pSdigList[i]=(TClonesArray*)(*pSDig)[i];
 503     if(pSdigList[i]->GetEntries()!=0) AliErrorClass("Some of sdigits lists is not empty");         //in principle those lists should be empty
 504   }
 505
 506   for(Int_t iHit=0;iHit<hits->GetEntries();iHit++){         //hits loop
 507     AliITShit *hit = (AliITShit*)hits->At(iHit);
 508     Double_t xz[2];
 509
 510     Int_t module;
 511     if(!fSegmentation->GlobalToDet(fSegmentation->GetLayerFromIdIndex(hit->GetModule()),hit->GetXG(),hit->GetYG(),hit->GetZG(),xz[0],xz[1],module)) continue;
 512     AliITSDigitUpgrade digit;
 513     digit.SetSignal(hit->GetIonization());
 514     digit.SetNelectrons(hit->GetIonization()/(3.62*1e-09));
 515     digit.SetLayer(fSegmentation->GetLayerFromIdIndex(hit->GetModule()));
 516     digit.SetModule(fSegmentation->GetModuleFromIdIndex(hit->GetModule()));//set the module (=sector) of ITSupgrade
 517     digit.SetTrackID(hit->GetTrack());
 518
 519     Int_t xpix = 999;
 520     Int_t zpix = 999; // shift at the next line to have zpix always positive. Such numbers are used to build the Pixel Id in the layer (> 0!)
 521     fSegmentation->DetToPixID(xz[0], xz[1],fSegmentation->GetLayerFromIdIndex(hit->GetModule()), xpix, zpix);
 522     digit.SetPixId(xpix,zpix);
 523     new((*pSdigList[fSegmentation->GetLayerFromIdIndex(hit->GetModule())])[nSdigit[fSegmentation->GetLayerFromIdIndex(hit->GetModule())]++]) AliITSDigitUpgrade(digit);
 524
 525   }
 526
 527   AliDebug(1,"Stop Hit2SumDig.");
 528 }
 529 //_______________________________________________________________________________________________
 530 void AliITSupgrade::MakeBranch(Option_t *option){
 531   //Create Tree branches
 532   AliDebug(1,Form("Start with option= %s.",option));
 533
 534   const Int_t kBufSize = 4000;
 535
 536   const char *cH = strstr(option,"H");
 537   const char *cD = strstr(option,"D");
 538   //const char *cR = strstr(option,"R");
 539   const char *cS = strstr(option,"S");
 540
 541   if(cH&&fLoader->TreeH()){
 542
 543     HitCreate();
 544     MakeBranchInTree(fLoader->TreeH(),"ITSupgrade",&fHits,kBufSize,0);
 545   }
 546
 547
 548   if(cS&&fLoader->TreeS()){
 549     SDigitsCreate();
 550     for(Int_t i=0;i<fNlayers;i++) MakeBranchInTree(fLoader->TreeS(),Form("Layer%d",i),&((*fSdigits)[i]),kBufSize,0);
 551   }
 552
 553
 554   if(cD&&fLoader->TreeD()){
 555     DigitsCreate();
 556     for(Int_t i=0;i<fNlayers;i++) MakeBranchInTree(fLoader->TreeD(),Form("Layer%d",i),&((*fDigits)[i]),kBufSize,0);
 557   }
 558   AliDebug(1,"Stop.");
 559 }
 560 //____________________________________________________________________________________________________
 561 void AliITSupgrade::SetTreeAddress()
 562 {
 563   //Set branch address for the Hits and Digits Tree.
 564   AliDebug(1,"Start.");
 565   if(fLoader->TreeH() && fLoader->TreeH()->GetBranch("ITSupgrade" )){
 566     HitCreate();
 567     fLoader->TreeH()->SetBranchAddress("ITSupgrade",&fHits);
 568
 569   }
 570
 571   if(fLoader->TreeS() && fLoader->TreeS()->GetBranch("Layer0" )){
 572     SDigitsCreate();
 573     for(int i=0;i<fNlayers;i++){
 574       fLoader->TreeS()->SetBranchAddress(Form("Layer%d",i),&((*fSdigits)[i]));
 575     }
 576   }
 577
 578   if(fLoader->TreeD() && fLoader->TreeD()->GetBranch("Layer0")){
 579     DigitsCreate();
 580     for(int i=0;i<fNlayers;i++){
 581       fLoader->TreeD()->SetBranchAddress(Form("Layer%d",i),&((*fDigits)[i]));
 582     }
 583   }
 584   AliDebug(1,"Stop.");
 585 }
 586 //____________________________________________________________________________________________________
 587 void AliITSupgrade::PrintSummary()
 588 {
 589
 590   // pdg booklet or http://pdg.lbl.gov/2010/reviews/rpp2010-rev-atomic-nuclear-prop.pdf
 591
 592   // X0 as X0[g/cm^2]/density
 593   Double_t beX0= 65.19/1.848;
 594   Double_t  cuX0= 12.86/8.96;
 595   Double_t siX0= 21.82/2.329;
 596
 597   if(!fSegmentation) fSegmentation= new AliITSsegmentationUpgrade();
 598   printf(" %10s %10s %10s %10s %10s %10s %10s %10s\n","Name","R [cm]","x/X0 [%]","Phi res[um]","Z res[um]","X segm[um]","Y segm[um]","widths [um]");
 599
 600    if(fBeampipe) printf(" Beampipe   %10.3f %10f %10s %10s %10s %10s %10.1f\n", fRadiusBP, 100.*fWidthBP/beX0,"-","-","-","-",fWidthBP*1.e+4);
 601    else printf(" * No Beam Pipe Upgrade  \n");
 602   for(Int_t i=0; i< fNlayers; i++){
 603     printf(" Si Layer %1i %10.3f %10f %10.2f %10.2f %10.2f %10.2f %10.1f\n",i, fRadii.At(i), 100.*fWidths.At(i)/siX0, (fSegmentation->GetCellSizeX(i) *1.e+4)/TMath::Sqrt(12.), (fSegmentation->GetCellSizeZ(i)*1.e+4)/TMath::Sqrt(12.),fSegmentation->GetCellSizeX(i) *1.e+4, fSegmentation->GetCellSizeZ(i)*1.e+4,fWidths.At(i)*1.e+4);
 604     printf(" Cu Layer   %10.3f %10f %10s %10s %10s %10s %10.1f\n", fRadiiCu.At(i), 100.*fWidthsCu.At(i)/cuX0,"-","-","-","-",fWidthsCu.At(i)*1.e+4);
 605   }
 606 }
 607 //____________________________________________________________________________________________________
 608 void AliITSupgrade::SetSegmentationX(Double_t x, Int_t lay){
 609 TFile *f = TFile::Open("Segmentation.root","UPDATE");
 610 if(!f) {
 611   AliError("\n\n\n Segmentation.root file does not exist. The segmentation is 0! \n\n\n");
 612   return;
 613  }
 614 else {
 615   TArrayD *a = (TArrayD*)f->Get("CellSizeX");
 616   if(!a)  {
 617    AliError("CessSizeX array does not exist!");
 618    return;
 619    }
 620   a->AddAt(x,lay);
 621   f->WriteObjectAny(a,"TArrayD","CellSizeX");
 622   f->Close();
 623   delete f;
 624  }
 625 }
 626 //____________________________________________________________________________________________________
 627 void AliITSupgrade::SetSegmentationZ(Double_t z, Int_t lay){
 628
 629 TFile *f = TFile::Open("Segmentation.root","UPDATE");
 630 if(!f) {
 631   AliError("\n\n\n Segmentation.root file does not exist. The segmentation is 0! \n\n\n");
 632   return;
 633  }
 634 else {
 635   TArrayD *a = (TArrayD*)f->Get("CellSizeZ");
 636   if(!a)  {
 637    AliError("CellSizeZ array does not exist!");
 638    return;
 639    }
 640   a->AddAt(z,lay);
 641   f->WriteObjectAny(a,"TArrayD","CellSizeZ");
 642   f->Close();
 643   delete f;
 644  }
 645
 646 }