FMD/AliFMD.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 /* $Id$ */
  16 /** @file    AliFMD.cxx
  17     @author  Christian Holm Christensen <cholm@nbi.dk>
  18     @date    Sun Mar 26 17:59:18 2006
  19     @brief   Implementation of AliFMD base class
  20 */
  21 //____________________________________________________________________
  22 //
  23 // Forward Multiplicity Detector based on Silicon wafers. This class
  24 // is the driver for especially simulation.
  25 //
  26 // The Forward Multiplicity Detector consists of 3 sub-detectors FMD1,
  27 // FMD2, and FMD3, each of which has 1 or 2 rings of silicon sensors.
  28 //
  29 // This is the base class for all FMD manager classes.
  30 //
  31 // The actual code is done by various separate classes.   Below is
  32 // diagram showing the relationship between the various FMD classes
  33 // that handles the simulation
  34 //
  35 //
  36 //       +----------+   +----------+
  37 //       | AliFMDv1 |   | AliFMDv0 |
  38 //       +----------+   +----------+
  39 //            |              |                    +-----------------+
  40 //       +----+--------------+                 +--| AliFMDDigitizer |
  41 //       |                                     |  +-----------------+
  42 //       |           +---------------------+   |
  43 //       |        +--| AliFMDBaseDigitizer |<--+
  44 //       V     1  |  +---------------------+   |
  45 //  +--------+<>--+                            |  +------------------+
  46 //  | AliFMD |                                 +--| AliFMDSDigitizer |
  47 //  +--------+<>--+                               +------------------+
  48 //             1  |  +---------------------+
  49 //                +--| AliFMDReconstructor |
  50 //                   +---------------------+
  51 //
  52 // *  AliFMD
  53 //    This defines the interface for the various parts of AliROOT that
  54 //    uses the FMD, like AliFMDSimulator, AliFMDDigitizer,
  55 //    AliFMDReconstructor, and so on.
  56 //
  57 // *  AliFMDv0
  58 //    This is a concrete implementation of the AliFMD interface.
  59 //    It is the responsibility of this class to create the FMD
  60 //    geometry.
  61 //
  62 // *  AliFMDv1
  63 //    This is a concrete implementation of the AliFMD interface.
  64 //    It is the responsibility of this class to create the FMD
  65 //    geometry, process hits in the FMD, and serve hits and digits to
  66 //    the various clients.
  67 //
  68 // *  AliFMDSimulator
  69 //    This is the base class for the FMD simulation tasks.   The
  70 //    simulator tasks are responsible to implment the geoemtry, and
  71 //    process hits.
  72 //
  73 // *  AliFMDReconstructor
  74 //    This is a concrete implementation of the AliReconstructor that
  75 //    reconstructs pseudo-inclusive-multiplicities from digits (raw or
  76 //    from simulation)
  77 //
  78 // Calibration and geometry parameters are managed by separate
  79 // singleton managers.  These are AliFMDGeometry and
  80 // AliFMDParameters.  Please refer to these classes for more
  81 // information on these.
  82 //
  83
  84 // These files are not in the same directory, so there's no reason to
  85 // ask the preprocessor to search in the current directory for these
  86 // files by including them with `#include "..."'
  87 #include <TBrowser.h>           // ROOT_TBrowser
  88 #include <TClonesArray.h>       // ROOT_TClonesArray
  89 #include <TGeoGlobalMagField.h> // ROOT_TGeoGlobalMagField
  90 #include <TGeoManager.h>        // ROOT_TGeoManager
  91 #include <TRotMatrix.h>         // ROOT_TRotMatrix
  92 #include <TTree.h>              // ROOT_TTree
  93 #include <TVector2.h>           // ROOT_TVector2
  94 #include <TVirtualMC.h>         // ROOT_TVirtualMC
  95 #include <cmath>                // __CMATH__
  96
  97 #include <AliDigitizationInput.h>       // ALIRUNDIGITIZER_H
  98 #include <AliLoader.h>          // ALILOADER_H
  99 #include <AliRun.h>             // ALIRUN_H
 100 #include <AliMC.h>              // ALIMC_H
 101 #include <AliMagF.h>            // ALIMAGF_H
 102 // #include <AliLog.h>          // ALILOG_H
 103 #include "AliFMDDebug.h" // Better debug macros
 104 #include "AliFMD.h"             // ALIFMD_H
 105 #include "AliFMDDigit.h"        // ALIFMDDIGIT_H
 106 #include "AliFMDSDigit.h"       // ALIFMDSDIGIT_H
 107 #include "AliFMDHit.h"          // ALIFMDHIT_H
 108 #include "AliFMDGeometry.h"     // ALIFMDGEOMETRY_H
 109 #include "AliFMDDetector.h"     // ALIFMDDETECTOR_H
 110 #include "AliFMDRing.h"         // ALIFMDRING_H
 111 #include "AliFMDDigitizer.h"    // ALIFMDDIGITIZER_H
 112 #include "AliFMDHitDigitizer.h" // ALIFMDSDIGITIZER_H
 113 // #define USE_SSDIGITIZER
 114 //#ifdef USE_SSDIGITIZER
 115 //# include "AliFMDSSDigitizer.h"       // ALIFMDSDIGITIZER_H
 116 //#endif
 117 // #include "AliFMDGeometryBuilder.h"
 118 #include "AliFMDRawWriter.h"    // ALIFMDRAWWRITER_H
 119 #include "AliFMDRawReader.h"    // ALIFMDRAWREADER_H
 120 #include "AliTrackReference.h"
 121 #include "AliFMDStripIndex.h"
 122 #include "AliFMDEncodedEdx.h"
 123 #include "AliFMDParameters.h"
 124 #include "AliFMDReconstructor.h"
 125
 126 //____________________________________________________________________
 127 ClassImp(AliFMD)
 128 #if 0
 129   ; // This is to keep Emacs from indenting the next line
 130 #endif
 131
 132 //____________________________________________________________________
 133 AliFMD::AliFMD()
 134   : AliDetector(),
 135     fSDigits(0),
 136     fNsdigits(0),
 137     fDetailed(kTRUE),
 138     fUseOld(kFALSE),
 139     fUseAssembly(kTRUE),
 140     fBad(0)
 141 {
 142   //
 143   // Default constructor for class AliFMD
 144   //
 145   AliFMDDebug(10, ("\tDefault CTOR"));
 146   fHits        = 0;
 147   fDigits      = 0;
 148   fIshunt      = 0;
 149   // fBad         = new TClonesArray("AliFMDHit");
 150 }
 151
 152 //____________________________________________________________________
 153 AliFMD::AliFMD(const char *name, const char *title)
 154   : AliDetector (name, title),
 155     fSDigits(0),
 156     fNsdigits(0),
 157     fDetailed(kTRUE),
 158     fUseOld(kFALSE),
 159     fUseAssembly(kFALSE),
 160     fBad(0)
 161 {
 162   //
 163   // Standard constructor for Forward Multiplicity Detector
 164   //
 165   AliFMDDebug(10, ("\tStandard CTOR"));
 166   // fBad         = new TClonesArray("AliFMDHit");
 167
 168   // Initialise Hit array
 169   // HitsArray();
 170   // gAlice->GetMCApp()->AddHitList(fHits);
 171
 172   // (S)Digits for the detectors disk
 173   // DigitsArray();
 174   // SDigitsArray();
 175
 176   // CHC: What is this?
 177   fIshunt = 0;
 178   //PH  SetMarkerColor(kRed);
 179   //PH  SetLineColor(kYellow);
 180 }
 181
 182 //____________________________________________________________________
 183 AliFMD::~AliFMD ()
 184 {
 185   // Destructor for base class AliFMD
 186   if (fHits) {
 187     fHits->Delete();
 188     delete fHits;
 189     fHits = 0;
 190   }
 191   if (fDigits) {
 192     fDigits->Delete();
 193     delete fDigits;
 194     fDigits = 0;
 195   }
 196   if (fSDigits) {
 197     fSDigits->Delete();
 198     delete fSDigits;
 199     fSDigits = 0;
 200   }
 201   if (fBad) {
 202     fBad->Delete();
 203     delete fBad;
 204     fBad = 0;
 205   }
 206 }
 207
 208
 209 //====================================================================
 210 //
 211 // GEometry ANd Traking
 212 //
 213 //____________________________________________________________________
 214 void
 215 AliFMD::CreateGeometry()
 216 {
 217   //
 218   // Create the geometry of Forward Multiplicity Detector.  The actual
 219   // construction of the geometry is delegated to the class
 220   // AliFMDGeometryBuilder, invoked by the singleton manager
 221   // AliFMDGeometry.
 222   //
 223   AliFMDGeometry*  fmd = AliFMDGeometry::Instance();
 224   fmd->SetDetailed(fDetailed);
 225   fmd->UseAssembly(fUseAssembly);
 226   fmd->Build();
 227 }
 228
 229 //____________________________________________________________________
 230 void AliFMD::CreateMaterials()
 231 {
 232   // Define the materials and tracking mediums needed by the FMD
 233   // simulation.   These mediums are made by sending the messages
 234   // AliMaterial, AliMixture, and AliMedium to the passed AliModule
 235   // object module.   The defined mediums are
 236   //
 237   //    FMD Si$         Silicon (active medium in sensors)
 238   //    FMD C$          Carbon fibre (support cone for FMD3 and vacuum pipe)
 239   //    FMD Al$         Aluminium (honeycomb support plates)
 240   //    FMD PCB$        Printed Circuit Board (FEE board with VA1_3)
 241   //    FMD Chip$       Electronics chips (currently not used)
 242   //    FMD Air$        Air (Air in the FMD)
 243   //    FMD Plastic$    Plastic (Support legs for the hybrid cards)
 244   //
 245   // The geometry builder should really be the one that creates the
 246   // materials, but the architecture of AliROOT makes that design
 247   // akward.  What should happen, was that the AliFMDGeometryBuilder
 248   // made the mediums, and that this class retrives pointers from the
 249   // TGeoManager, and registers the mediums here.  Alas, it's not
 250   // really that easy.
 251   //
 252   AliFMDDebug(10, ("\tCreating materials"));
 253   // Get pointer to geometry singleton object.
 254   AliFMDGeometry* geometry = AliFMDGeometry::Instance();
 255   geometry->Init();
 256 #if 0
 257   if (gGeoManager && gGeoManager->GetMedium("FMD Si$")) {
 258     // We need to figure out the some stuff about the geometry
 259     fmd->ExtractGeomInfo();
 260     return;
 261   }
 262 #endif
 263   Int_t    id;
 264   Double_t a                = 0;
 265   Double_t z                = 0;
 266   Double_t density          = 0;
 267   Double_t radiationLength  = 0;
 268   Double_t absorbtionLength = 999;
 269   Int_t    fieldType        = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->Integ();     // Field type
 270   Double_t maxField         = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->Max();     // Field max.
 271   Double_t maxBending       = 0;     // Max Angle
 272   Double_t maxStepSize      = 0.001; // Max step size
 273   Double_t maxEnergyLoss    = 1;     // Max Delta E
 274   Double_t precision        = 0.001; // Precision
 275   Double_t minStepSize      = 0.001; // Minimum step size
 276
 277   // Silicon
 278   a                = 28.0855;
 279   z                = 14.;
 280   density          = geometry->GetSiDensity();
 281   radiationLength  = 9.36;
 282   maxBending       = 1;
 283   maxStepSize      = .001;
 284   precision        = .001;
 285   minStepSize      = .001;
 286   id               = kSiId;
 287   AliMaterial(id, "Si$", a, z, density, radiationLength, absorbtionLength);
 288   AliMedium(kSiId, "Si$", id,1,fieldType,maxField,maxBending,
 289             maxStepSize,maxEnergyLoss,precision,minStepSize);
 290
 291
 292   // Carbon
 293   a                = 12.011;
 294   z                = 6.;
 295   density          = 2.265;
 296   radiationLength  = 18.8;
 297   maxBending       = 10;
 298   maxStepSize      = .01;
 299   precision        = .003;
 300   minStepSize      = .003;
 301   id               = kCarbonId;
 302   AliMaterial(id, "Carbon$", a, z, density, radiationLength, absorbtionLength);
 303   AliMedium(kCarbonId, "Carbon$", id,0,fieldType,maxField,maxBending,
 304                     maxStepSize,maxEnergyLoss,precision,minStepSize);
 305
 306   // Aluminum
 307   a                = 26.981539;
 308   z                = 13.;
 309   density          = 2.7;
 310   radiationLength  = 8.9;
 311   id               = kAlId;
 312   AliMaterial(id, "Aluminum$",a,z, density, radiationLength, absorbtionLength);
 313   AliMedium(kAlId, "Aluminum$", id, 0, fieldType, maxField, maxBending,
 314             maxStepSize, maxEnergyLoss, precision, minStepSize);
 315
 316
 317   // Copper
 318   a                = 63.546;
 319   z                = 29;
 320   density          =  8.96;
 321   radiationLength  =  1.43;
 322   id               = kCopperId;
 323   AliMaterial(id, "Copper$",
 324                       a, z, density, radiationLength, absorbtionLength);
 325   AliMedium(kCopperId, "Copper$", id, 0, fieldType, maxField, maxBending,
 326             maxStepSize, maxEnergyLoss, precision, minStepSize);
 327
 328
 329   // Silicon chip
 330   {
 331     Float_t as[] = { 12.0107,      14.0067,      15.9994,
 332                       1.00794,     28.0855,     107.8682 };
 333     Float_t zs[] = {  6.,           7.,           8.,
 334                       1.,          14.,          47. };
 335     Float_t ws[] = {  0.039730642,  0.001396798,  0.01169634,
 336                       0.004367771,  0.844665,     0.09814344903 };
 337     density          = 2.36436;
 338     maxBending       = 10;
 339     maxStepSize      = .01;
 340     precision        = .003;
 341     minStepSize      = .003;
 342     id               = kSiChipId;
 343     AliMixture(id, "Si Chip$", as, zs, density, 6, ws);
 344     AliMedium(kSiChipId, "Si Chip$",  id, 0, fieldType, maxField, maxBending,
 345               maxStepSize, maxEnergyLoss, precision, minStepSize);
 346   }
 347
 348   // Kaption
 349   {
 350     Float_t as[] = { 1.00794,  12.0107,  14.010,   15.9994};
 351     Float_t zs[] = { 1.,        6.,       7.,       8.};
 352     Float_t ws[] = { 0.026362,  0.69113,  0.07327,  0.209235};
 353     density          = 1.42;
 354     maxBending       = 1;
 355     maxStepSize      = .001;
 356     precision        = .001;
 357     minStepSize      = .001;
 358     id               = kKaptonId;
 359     AliMixture(id, "Kaption$", as, zs, density, 4, ws);
 360     AliMedium(kKaptonId, "Kaption$", id,0,fieldType,maxField,maxBending,
 361               maxStepSize,maxEnergyLoss,precision,minStepSize);
 362   }
 363
 364   // Air
 365   {
 366     Float_t as[] = { 12.0107, 14.0067,   15.9994,  39.948 };
 367     Float_t zs[] = {  6.,      7.,       8.,       18. };
 368     Float_t ws[] = { 0.000124, 0.755267, 0.231781, 0.012827 };
 369     density      = .00120479;
 370     maxBending   = 1;
 371     maxStepSize  = .001;
 372     precision    = .001;
 373     minStepSize  = .001;
 374     id           = kAirId;
 375     AliMixture(id, "Air$", as, zs, density, 4, ws);
 376     AliMedium(kAirId, "Air$", id,0,fieldType,maxField,maxBending,
 377               maxStepSize,maxEnergyLoss,precision,minStepSize);
 378   }
 379
 380   // PCB
 381   {
 382     Float_t zs[] = { 14.,         20.,         13.,         12.,
 383                       5.,         22.,         11.,         19.,
 384                      26.,          9.,          8.,          6.,
 385                       7.,          1.};
 386     Float_t as[] = { 28.0855,     40.078,      26.981538,   24.305,
 387                      10.811,      47.867,      22.98977,    39.0983,
 388                      55.845,      18.9984,     15.9994,     12.0107,
 389                      14.0067,      1.00794};
 390     Float_t ws[] = {  0.15144894,  0.08147477,  0.04128158,  0.00904554,
 391                       0.01397570,  0.00287685,  0.00445114,  0.00498089,
 392                       0.00209828,  0.00420000,  0.36043788,  0.27529426,
 393                       0.01415852,  0.03427566};
 394     density      = 1.8;
 395     maxBending   = 1;
 396     maxStepSize  = .001;
 397     precision    = .001;
 398     minStepSize  = .001;
 399     id           = kPcbId;
 400     AliMixture(id, "PCB$", as, zs, density, 14, ws);
 401     AliMedium(kPcbId, "PCB$", id,0,fieldType,maxField,maxBending,
 402               maxStepSize,maxEnergyLoss,precision,minStepSize);
 403   }
 404
 405   // Stainless steel
 406   {
 407     Float_t as[] = { 55.847, 51.9961, 58.6934, 28.0855 };
 408     Float_t zs[] = { 26.,    24.,     28.,     14.     };
 409     Float_t ws[] = { .715,   .18,     .1,      .005    };
 410     density      = 7.88;
 411     id           = kSteelId;
 412     AliMixture(id, "Steel$", as, zs, density, 4, ws);
 413     AliMedium(kSteelId, "Steel$", id, 0, fieldType, maxField, maxBending,
 414               maxStepSize, maxEnergyLoss, precision, minStepSize);
 415   }
 416   // Plastic
 417   {
 418     Float_t as[] = { 1.01, 12.01 };
 419     Float_t zs[] = { 1.,   6.    };
 420     Float_t ws[] = { 1.,   1.    };
 421     density      = 1.03;
 422     maxBending   = 10;
 423     maxStepSize  = .01;
 424     precision    = .003;
 425     minStepSize  = .003;
 426     id           = kPlasticId;
 427     AliMixture(id, "Plastic$", as, zs, density, -2, ws);
 428     AliMedium(kPlasticId, "Plastic$", id,0,fieldType,maxField,maxBending,
 429               maxStepSize,maxEnergyLoss,precision,minStepSize);
 430   }
 431
 432 }
 433
 434 #if 0
 435 //____________________________________________________________________
 436 void
 437 AliFMD::SetTrackingParameters(Int_t imed,
 438                               Float_t gamma,
 439                               Float_t electron,
 440                               Float_t neutral_hadron,
 441                               Float_t charged_hadron,
 442                               Float_t muon,
 443                               Float_t electron_bremstrahlung,
 444                               Float_t muon__bremstrahlung,
 445                               Float_t electron_delta,
 446                               Float_t muon_delta,
 447                               Float_t muon_pair,
 448                               Int_t   annihilation,
 449                               Int_t   bremstrahlung,
 450                               Int_t   compton_scattering,
 451                               Int_t   decay,
 452                               Int_t   delta_ray,
 453                               Int_t   hadronic,
 454                               Int_t   energy_loss,
 455                               Int_t   multiple_scattering,
 456                               Int_t   pair_production,
 457                               Int_t   photon_production,
 458                               Int_t   rayleigh_scattering)
 459 {
 460   // Disabled by request of FCA, kept for reference only
 461   if (!gMC) return;
 462   TArrayI& idtmed = *(GetIdtmed());
 463   Int_t    iimed  = idtmed[imed];
 464   // gMC->Gstpar(iimed, "CUTGAM",       gamma);
 465   // gMC->Gstpar(iimed, "CUTELE",       electron);
 466   // gMC->Gstpar(iimed, "CUTNEU",       neutral_hadron);
 467   // gMC->Gstpar(iimed, "CUTHAD",       charged_hadron);
 468   // gMC->Gstpar(iimed, "CUTMUO",       muon);
 469   // gMC->Gstpar(iimed, "BCUTE",        electron_bremstrahlung);
 470   // gMC->Gstpar(iimed, "BCUTM",        muon__bremstrahlung);
 471   // gMC->Gstpar(iimed, "DCUTE",        electron_delta);
 472   // gMC->Gstpar(iimed, "DCUTM",        muon_delta);
 473   // gMC->Gstpar(iimed, "PPCUTM",       muon_pair);
 474   // gMC->Gstpar(iimed, "ANNI", Float_t(annihilation));
 475   // gMC->Gstpar(iimed, "BREM", Float_t(bremstrahlung));
 476   // gMC->Gstpar(iimed, "COMP", Float_t(compton_scattering));
 477   // gMC->Gstpar(iimed, "DCAY", Float_t(decay));
 478   // gMC->Gstpar(iimed, "DRAY", Float_t(delta_ray));
 479   // gMC->Gstpar(iimed, "HADR", Float_t(hadronic));
 480   // gMC->Gstpar(iimed, "LOSS", Float_t(energy_loss));
 481   // gMC->Gstpar(iimed, "MULS", Float_t(multiple_scattering));
 482   // gMC->Gstpar(iimed, "PAIR", Float_t(pair_production));
 483   // gMC->Gstpar(iimed, "PHOT", Float_t(photon_production));
 484   // gMC->Gstpar(iimed, "RAYL", Float_t(rayleigh_scattering));
 485 }
 486 #endif
 487
 488 //____________________________________________________________________
 489 void
 490 AliFMD::Init()
 491 {
 492   // Initialize the detector
 493   //
 494   AliFMDDebug(1, ("Initialising FMD detector object"));
 495   TVirtualMC*      mc     = TVirtualMC::GetMC();
 496   AliFMDGeometry*  fmd    = AliFMDGeometry::Instance();
 497   TArrayI          actGeo = fmd->ActiveIds();
 498   bool             valid  = true;
 499   if (actGeo.fN <= 0) valid = false;
 500   else {
 501     for (int i = 0; i < actGeo.fN; i++) {
 502       if (actGeo[i] < 0) {
 503         valid = false;
 504         break;
 505       }
 506     }
 507   }
 508   if (!valid) {
 509     AliFMDDebug(1, ("Extracting geometry info from loaded geometry"));
 510     fmd->ExtractGeomInfo();
 511     actGeo = fmd->ActiveIds();
 512   }
 513   TArrayI          actVmc(actGeo.fN);
 514   for (Int_t i = 0; i < actGeo.fN; i++) {
 515     if (actGeo[i] < 0) {
 516       AliError(Form("Invalid id: %d", actGeo[i]));
 517       continue;
 518     }
 519     TGeoVolume *sens = gGeoManager->GetVolume(actGeo[i]);
 520     if (!sens) {
 521       AliError(Form("No TGeo volume for sensitive volume ID=%d",actGeo[i]));
 522       continue;
 523     }
 524     actVmc[i] = mc->VolId(sens->GetName());
 525     AliFMDDebug(1, ("Active vol id # %d: %d changed to %d",
 526                     i, actGeo[i], actVmc[i]));
 527   }
 528   fmd->SetActive(actVmc.fArray, actVmc.fN);
 529   // fmd->InitTransformations();
 530 }
 531
 532 //____________________________________________________________________
 533 void
 534 AliFMD::FinishEvent()
 535 {
 536   // Called at the end of the an event in simulations.  If the debug
 537   // level is high enough, then the `bad' hits are printed.
 538   //
 539   if (AliLog::GetDebugLevel("FMD", "AliFMD") < 10) return;
 540   if (fBad && fBad->GetEntries() > 0) {
 541     AliWarning(Form("got %d 'bad' hits", fBad->GetEntries()));
 542     TIter next(fBad);
 543     AliFMDHit* hit;
 544     while ((hit = static_cast<AliFMDHit*>(next()))) hit->Print("D");
 545     fBad->Clear();
 546   }
 547 }
 548
 549
 550
 551 //====================================================================
 552 //
 553 // Hit and Digit managment
 554 //
 555 //____________________________________________________________________
 556 void
 557 AliFMD::MakeBranch(Option_t * option)
 558 {
 559   // Create Tree branches for the FMD.
 560   //
 561   // Options:
 562   //
 563   //    H          Make a branch of TClonesArray of AliFMDHit's
 564   //    D          Make a branch of TClonesArray of AliFMDDigit's
 565   //    S          Make a branch of TClonesArray of AliFMDSDigit's
 566   //
 567   const Int_t kBufferSize = 16000;
 568   TString branchname(GetName());
 569   TString opt(option);
 570
 571   if (opt.Contains("H", TString::kIgnoreCase)) {
 572     HitsArray();
 573     AliDetector::MakeBranch(option);
 574   }
 575   if (opt.Contains("D", TString::kIgnoreCase)) {
 576     DigitsArray();
 577     MakeBranchInTree(fLoader->TreeD(), branchname.Data(),
 578                      &fDigits, kBufferSize, 0);
 579   }
 580   if (opt.Contains("S", TString::kIgnoreCase)) {
 581     SDigitsArray();
 582     MakeBranchInTree(fLoader->TreeS(), branchname.Data(),
 583                      &fSDigits, kBufferSize, 0);
 584   }
 585 }
 586
 587 //____________________________________________________________________
 588 void
 589 AliFMD::SetTreeAddress()
 590 {
 591   // Set branch address for the Hits, Digits, and SDigits Tree.
 592   if (fLoader->TreeH()) HitsArray();
 593   AliDetector::SetTreeAddress();
 594
 595   TTree *treeD = fLoader->TreeD();
 596   if (treeD) {
 597     DigitsArray();
 598     TBranch* branch = treeD->GetBranch ("FMD");
 599     if (branch) branch->SetAddress(&fDigits);
 600   }
 601
 602   TTree *treeS = fLoader->TreeS();
 603   if (treeS) {
 604     SDigitsArray();
 605     TBranch* branch = treeS->GetBranch ("FMD");
 606     if (branch) branch->SetAddress(&fSDigits);
 607   }
 608 }
 609
 610 //____________________________________________________________________
 611 void
 612 AliFMD::SetHitsAddressBranch(TBranch *b)
 613 {
 614   // Set the TClonesArray to read hits into.
 615   b->SetAddress(&fHits);
 616 }
 617 //____________________________________________________________________
 618 void
 619 AliFMD::SetSDigitsAddressBranch(TBranch *b)
 620 {
 621   // Set the TClonesArray to read hits into.
 622   b->SetAddress(&fSDigits);
 623 }
 624
 625 //____________________________________________________________________
 626 void
 627 AliFMD::AddHit(Int_t track, Int_t *vol, Float_t *hits)
 628 {
 629   // Add a hit to the hits tree
 630   //
 631   // The information of the two arrays are decoded as
 632   //
 633   // Parameters
 634   //    track                Track #
 635   //    ivol[0]  [UShort_t ] Detector #
 636   //    ivol[1]  [Char_t   ] Ring ID
 637   //    ivol[2]  [UShort_t ] Sector #
 638   //    ivol[3]  [UShort_t ] Strip #
 639   //    hits[0]  [Float_t  ] Track's X-coordinate at hit
 640   //    hits[1]  [Float_t  ] Track's Y-coordinate at hit
 641   //    hits[3]  [Float_t  ] Track's Z-coordinate at hit
 642   //    hits[4]  [Float_t  ] X-component of track's momentum
 643   //    hits[5]  [Float_t  ] Y-component of track's momentum
 644   //    hits[6]  [Float_t  ] Z-component of track's momentum
 645   //    hits[7]  [Float_t  ] Energy deposited by track
 646   //    hits[8]  [Int_t    ] Track's particle Id #
 647   //    hits[9]  [Float_t  ] Time when the track hit
 648   //
 649   //
 650   AddHitByFields(track,
 651                  UShort_t(vol[0]),  // Detector #
 652                  Char_t(vol[1]),    // Ring ID
 653                  UShort_t(vol[2]),  // Sector #
 654                  UShort_t(vol[3]),  // Strip #
 655                  hits[0],           // X
 656                  hits[1],           // Y
 657                  hits[2],           // Z
 658                  hits[3],           // Px
 659                  hits[4],           // Py
 660                  hits[5],           // Pz
 661                  hits[6],           // Energy loss
 662                  Int_t(hits[7]),    // PDG
 663                  hits[8]);          // Time
 664 }
 665
 666 //____________________________________________________________________
 667 AliFMDHit*
 668 AliFMD::AddHitByFields(Int_t    track,
 669                        UShort_t detector,
 670                        Char_t   ring,
 671                        UShort_t sector,
 672                        UShort_t strip,
 673                        Float_t  x,
 674                        Float_t  y,
 675                        Float_t  z,
 676                        Float_t  px,
 677                        Float_t  py,
 678                        Float_t  pz,
 679                        Float_t  edep,
 680                        Int_t    pdg,
 681                        Float_t  t,
 682                        Float_t  l,
 683                        Bool_t   stop)
 684 {
 685   // Add a hit to the list
 686   //
 687   // Parameters:
 688   //
 689   //    track     Track #
 690   //    detector  Detector # (1, 2, or 3)
 691   //    ring      Ring ID ('I' or 'O')
 692   //    sector    Sector # (For inner/outer rings: 0-19/0-39)
 693   //    strip     Strip # (For inner/outer rings: 0-511/0-255)
 694   //    x         Track's X-coordinate at hit
 695   //    y         Track's Y-coordinate at hit
 696   //    z         Track's Z-coordinate at hit
 697   //    px        X-component of track's momentum
 698   //    py        Y-component of track's momentum
 699   //    pz        Z-component of track's momentum
 700   //    edep      Energy deposited by track
 701   //    pdg       Track's particle Id #
 702   //    t         Time when the track hit
 703   //    l         Track length through the material.
 704   //    stop      Whether track was stopped or disappeared
 705   //
 706   TClonesArray& a = *(HitsArray());
 707   // Search through the list of already registered hits, and see if we
 708   // find a hit with the same parameters.  If we do, then don't create
 709   // a new hit, but rather update the energy deposited in the hit.
 710   // This is done, so that a FLUKA based simulation will get the
 711   // number of hits right, not just the enerrgy deposition.
 712   AliFMDHit* hit = 0;
 713   for (Int_t i = 0; i < fNhits; i++) {
 714     if (!a.At(i)) continue;
 715     hit = static_cast<AliFMDHit*>(a.At(i));
 716     if (hit->Detector() == detector
 717         && hit->Ring() == ring
 718         && hit->Sector() == sector
 719         && hit->Strip() == strip
 720         && hit->Track() == track) {
 721       AliFMDDebug(1, ("already had a hit in FMD%d%c[%2d,%3d] for track # %d,"
 722                        " adding energy (%f) to that hit (%f) -> %f",
 723                        detector, ring, sector, strip, track, edep, hit->Edep(),
 724                        hit->Edep() + edep));
 725       hit->SetEdep(hit->Edep() + edep);
 726       return hit;
 727     }
 728   }
 729   // If hit wasn't already registered, do so know.
 730   hit = new (a[fNhits]) AliFMDHit(fIshunt, track, detector, ring, sector,
 731                                   strip, x, y, z, px, py, pz, edep, pdg, t,
 732                                   l, stop);
 733   // gMC->AddTrackReference(track, 12);
 734   fNhits++;
 735
 736   //Reference track
 737
 738   AliMC *mcApplication = (AliMC*)gAlice->GetMCApp();
 739
 740   AliTrackReference* trackRef =
 741     AddTrackReference(mcApplication->GetCurrentTrackNumber(),
 742                       AliTrackReference::kFMD);
 743   UInt_t stripId = AliFMDStripIndex::Pack(detector,ring,sector,strip);
 744   UInt_t dedx    = AliFMDEncodedEdx::Encode(edep, l);
 745
 746   trackRef->SetUserId((dedx << 19) | stripId);
 747
 748
 749   return hit;
 750 }
 751
 752 //____________________________________________________________________
 753 void
 754 AliFMD::AddDigit(Int_t* digits, Int_t*)
 755 {
 756   // Add a digit to the Digit tree
 757   //
 758   // Paramters
 759   //
 760   //    digits[0]  [UShort_t] Detector #
 761   //    digits[1]  [Char_t]   Ring ID
 762   //    digits[2]  [UShort_t] Sector #
 763   //    digits[3]  [UShort_t] Strip #
 764   //    digits[4]  [UShort_t] ADC Count
 765   //    digits[5]  [Short_t]  ADC Count, -1 if not used
 766   //    digits[6]  [Short_t]  ADC Count, -1 if not used
 767   //
 768   AddDigitByFields(UShort_t(digits[0]),  // Detector #
 769                    Char_t(digits[1]),    // Ring ID
 770                    UShort_t(digits[2]),  // Sector #
 771                    UShort_t(digits[3]),  // Strip #
 772                    UShort_t(digits[4]),  // ADC Count1
 773                    Short_t(digits[5]),   // ADC Count2
 774                    Short_t(digits[6]),   // ADC Count3
 775                    Short_t(digits[7]));
 776 }
 777
 778 //____________________________________________________________________
 779 void
 780 AliFMD::AddDigitByFields(UShort_t       detector,
 781                          Char_t         ring,
 782                          UShort_t       sector,
 783                          UShort_t       strip,
 784                          UShort_t       count1,
 785                          Short_t        count2,
 786                          Short_t        count3,
 787                          Short_t        count4,
 788                          UShort_t       nrefs,
 789                          Int_t*         refs)
 790 {
 791   // add a real digit - as coming from data
 792   //
 793   // Parameters
 794   //
 795   //    detector  Detector # (1, 2, or 3)
 796   //    ring      Ring ID ('I' or 'O')
 797   //    sector    Sector # (For inner/outer rings: 0-19/0-39)
 798   //    strip     Strip # (For inner/outer rings: 0-511/0-255)
 799   //    count1    ADC count (a 10-bit word)
 800   //    count2    ADC count (a 10-bit word), or -1 if not used
 801   //    count3    ADC count (a 10-bit word), or -1 if not used
 802   TClonesArray& a = *(DigitsArray());
 803
 804   AliFMDDebug(15, ("Adding digit # %5d/%5d for FMD%d%c[%2d,%3d]"
 805                    "=(%d,%d,%d,%d) with %d tracks",
 806                    fNdigits-1, a.GetEntriesFast(),
 807                    detector, ring, sector, strip,
 808                    count1, count2, count3, count4, nrefs));
 809   new (a[fNdigits++])
 810     AliFMDDigit(detector, ring, sector, strip,
 811                 count1, count2, count3, count4, nrefs, refs);
 812
 813 }
 814
 815 //____________________________________________________________________
 816 void
 817 AliFMD::AddSDigit(Int_t* digits)
 818 {
 819   // Add a digit to the SDigit tree
 820   //
 821   // Paramters
 822   //
 823   //    digits[0]  [UShort_t] Detector #
 824   //    digits[1]  [Char_t]   Ring ID
 825   //    digits[2]  [UShort_t] Sector #
 826   //    digits[3]  [UShort_t] Strip #
 827   //    digits[4]  [Float_t]  Total energy deposited
 828   //    digits[5]  [UShort_t] ADC Count
 829   //    digits[6]  [Short_t]  ADC Count, -1 if not used
 830   //    digits[7]  [Short_t]  ADC Count, -1 if not used
 831   //
 832   AddSDigitByFields(UShort_t(digits[0]),   // Detector #
 833                     Char_t(digits[1]),     // Ring ID
 834                     UShort_t(digits[2]),   // Sector #
 835                     UShort_t(digits[3]),   // Strip #
 836                     Float_t(digits[4]),    // Edep
 837                     UShort_t(digits[5]),   // ADC Count1
 838                     Short_t(digits[6]),    // ADC Count2
 839                     Short_t(digits[7]),    // ADC Count3
 840                     Short_t(digits[8]),    // ADC Count4
 841                     UShort_t(digits[9]),   // N particles
 842                     UShort_t(digits[10])); // N primaries
 843 }
 844
 845 //____________________________________________________________________
 846 void
 847 AliFMD::AddSDigitByFields(UShort_t       detector,
 848                           Char_t         ring,
 849                           UShort_t       sector,
 850                           UShort_t       strip,
 851                           Float_t        edep,
 852                           UShort_t       count1,
 853                           Short_t        count2,
 854                           Short_t        count3,
 855                           Short_t        count4,
 856                           UShort_t       ntot,
 857                           UShort_t       nprim,
 858                           Int_t*         refs)
 859 {
 860   // add a summable digit
 861   //
 862   // Parameters
 863   //
 864   //    detector  Detector # (1, 2, or 3)
 865   //    ring      Ring ID ('I' or 'O')
 866   //    sector    Sector # (For inner/outer rings: 0-19/0-39)
 867   //    strip     Strip # (For inner/outer rings: 0-511/0-255)
 868   //    edep      Total energy deposited
 869   //    count1    ADC count (a 10-bit word)
 870   //    count2    ADC count (a 10-bit word), or -1 if not used
 871   //    count3    ADC count (a 10-bit word), or -1 if not used
 872   //
 873   TClonesArray& a = *(SDigitsArray());
 874   // AliFMDDebug(0, ("Adding sdigit # %d", fNsdigits));
 875
 876   AliFMDDebug(15, ("Adding sdigit # %5d/%5d for FMD%d%c[%2d,%3d]"
 877                    "=(%d,%d,%d,%d) with %d tracks %d primaries (%p)",
 878                    fNsdigits-1, a.GetEntriesFast(),
 879                    detector, ring, sector, strip,
 880                    count1, count2, count3, count4, ntot, nprim, refs));
 881   new (a[fNsdigits++])
 882     AliFMDSDigit(detector, ring, sector, strip, edep,
 883                  count1, count2, count3, count4, ntot, nprim, refs);
 884 }
 885
 886 //____________________________________________________________________
 887 void
 888 AliFMD::ResetSDigits()
 889 {
 890   // Reset number of digits and the digits array for this detector.
 891   //
 892   fNsdigits   = 0;
 893   if (fSDigits) fSDigits->Clear();
 894 }
 895
 896
 897 //____________________________________________________________________
 898 TClonesArray*
 899 AliFMD::HitsArray()
 900 {
 901   // Initialize hit array if not already, and return pointer to it.
 902   if (!fHits) {
 903     fHits = new TClonesArray("AliFMDHit", 1000);
 904     fNhits = 0;
 905     if (gAlice && gAlice->GetMCApp() && gAlice->GetMCApp()->GetHitLists())
 906       gAlice->GetMCApp()->AddHitList(fHits);
 907   }
 908   return fHits;
 909 }
 910
 911 //____________________________________________________________________
 912 TClonesArray*
 913 AliFMD::DigitsArray()
 914 {
 915   // Initialize digit array if not already, and return pointer to it.
 916   if (!fDigits) {
 917     fDigits = new TClonesArray("AliFMDDigit", 1000);
 918     fNdigits = 0;
 919   }
 920   return fDigits;
 921 }
 922
 923 //____________________________________________________________________
 924 TClonesArray*
 925 AliFMD::SDigitsArray()
 926 {
 927   // Initialize digit array if not already, and return pointer to it.
 928   if (!fSDigits) {
 929     fSDigits = new TClonesArray("AliFMDSDigit", 1000);
 930     fNsdigits = 0;
 931   }
 932   return fSDigits;
 933 }
 934
 935 //====================================================================
 936 //
 937 // Digitization
 938 //
 939 //____________________________________________________________________
 940 void
 941 AliFMD::Hits2Digits()
 942 {
 943   // Create AliFMDDigit's from AliFMDHit's.  This is done by making a
 944   // AliFMDDigitizer, and executing that code.
 945   //
 946   AliFMDHitDigitizer digitizer(this, AliFMDHitDigitizer::kDigits);
 947   digitizer.Init();
 948   digitizer.Digitize("");
 949 }
 950
 951 //____________________________________________________________________
 952 void
 953 AliFMD::Hits2SDigits()
 954 {
 955   // Create AliFMDSDigit's from AliFMDHit's.  This is done by creating
 956   // an AliFMDSDigitizer object, and executing it.
 957   //
 958   AliFMDHitDigitizer digitizer(this, AliFMDHitDigitizer::kSDigits);
 959   digitizer.Init();
 960   digitizer.Digitize("");
 961 }
 962
 963
 964 //____________________________________________________________________
 965 AliDigitizer*
 966 AliFMD::CreateDigitizer(AliDigitizationInput* digInput) const
 967 {
 968   // Create a digitizer object
 969
 970   /* This is what we probably _should_ do */
 971   AliFMDBaseDigitizer* digitizer = 0;
 972
 973 #ifdef USE_SSDIGITIZER
 974   digitizer = new AliFMDSSDigitizer(digInput);
 975 #else
 976   /* This is what we actually do, and will work */
 977 #if 0
 978   AliInfo("SDigit->Digit conversion not really supported, "
 979           "doing Hit->Digit conversion instead");
 980 #endif
 981   digitizer = new AliFMDDigitizer(digInput);
 982 #endif
 983   return digitizer;
 984 }
 985
 986 //====================================================================
 987 //
 988 // Raw data simulation
 989 //
 990 //__________________________________________________________________
 991 void
 992 AliFMD::Digits2Raw()
 993 {
 994   // Turn digits into raw data.
 995   //
 996   // This uses the class AliFMDRawWriter to do the job.   Please refer
 997   // to that class for more information.
 998   AliFMDRawWriter writer(this);
 999   writer.Exec();
1000 }
1001
1002 //====================================================================
1003 //
1004 // Raw data reading
1005 //
1006 //__________________________________________________________________
1007 Bool_t
1008 AliFMD::Raw2SDigits(AliRawReader* reader)
1009 {
1010   // Turn digits into raw data.
1011   //
1012   // This uses the class AliFMDRawWriter to do the job.   Please refer
1013   // to that class for more information.
1014   AliFMDParameters::Instance()->Init();
1015   MakeTree("S");
1016   MakeBranch("S");
1017
1018   TClonesArray*       sdigits = SDigits();
1019   AliFMDReconstructor rec;
1020
1021   // The two boolean arguments
1022   //   Make sdigits instead of digits
1023   //   Subtract the pedestal off the signal
1024   rec.Digitize(reader, sdigits);
1025   //
1026   // Bool_t ret = fmdReader.ReadAdcs(sdigits, kTRUE, kTRUE);
1027   // sdigits->ls();
1028   UShort_t ns = sdigits->GetEntriesFast();
1029   if (AliLog::GetDebugLevel("FMD", 0) > 5) {
1030     for (UShort_t i = 0; i < ns; i++)
1031       sdigits->At(i)->Print("pl");
1032   }
1033   AliFMDDebug(1, ("Got a total of %d SDigits", ns));
1034
1035   fLoader->TreeS()->Fill();
1036   ResetSDigits();
1037   fLoader->WriteSDigits("OVERWRITE");
1038
1039   return kTRUE;
1040 }
1041
1042
1043 //====================================================================
1044 //
1045 // Utility
1046 //
1047 //__________________________________________________________________
1048 void
1049 AliFMD::Browse(TBrowser* b)
1050 {
1051   // Browse this object.
1052   //
1053   AliFMDDebug(30, ("\tBrowsing the FMD"));
1054   AliDetector::Browse(b);
1055   b->Add(AliFMDGeometry::Instance());
1056 }
1057
1058 //____________________________________________________________________
1059 void
1060 AliFMD::AddAlignableVolumes() const
1061 {
1062   //
1063   // Create entries for alignable volumes associating the symbolic volume
1064   // name with the corresponding volume path. Needs to be syncronized with
1065   // eventual changes in the geometry.
1066   //
1067   // This code was made by Raffaele Grosso <rgrosso@mail.cern.ch>.  I
1068   // (cholm) will probably want to change it.   For one, I think it
1069   // should be the job of the geometry manager to deal with this.
1070   AliInfo("Add FMD alignable volumes");
1071   AliFMDGeometry::Instance()->SetAlignableVolumes();
1072 #if 0
1073   for(size_t f = 1; f <= 3; f++){ // Detector 1,2,3
1074     for(size_t tb =  0; tb <2 ; tb++){ // Top/Bottom
1075       char     stb = tb == 0 ? 'T' : 'B';
1076       unsigned min = tb == 0 ? 0   : 5;
1077
1078       TString halfVol(Form("/ALIC_1/F%dM%c_%d", f, stb, f));
1079       TString halfSym(halfVol);
1080       if(!gGeoManager->SetAlignableEntry(halfSym.Data(),halfVol.Data()))
1081         AliFatal(Form("Alignable entry %s not created. "
1082                       "Volume path %s not valid",
1083                       halfSym.Data(),halfVol.Data()));
1084       for(size_t io = 0; io < 2; io++){ // inner, outer
1085         if (f==1 && io==1) continue; // Only one ring in FMD1
1086         if(tb == 1 && io==1) min=10;
1087         char     sio = (io == 0 ? 'I' : 'O');
1088         unsigned nio = (io == 0 ? 3   : 9);
1089         unsigned max = (io == 0 ? 5   : 10) + min;
1090
1091         for(size_t i = min; i < max; i++) { // Modules
1092           TString modVol(Form("%s/F%c%cV_7%d/F%cSE_%d", halfVol.Data(),
1093                               sio, stb, nio, sio, i));
1094           TString modSym(modVol);
1095           if(!gGeoManager->SetAlignableEntry(modSym.Data(),modVol.Data()))
1096             AliFatal(Form("Alignable entry %s not created. "
1097                           "Volume path %s not valid",
1098                           modSym.Data(), modVol.Data()));
1099         }
1100       }
1101     }
1102   }
1103 #endif
1104 }
1105 //___________________________________________________________________
1106 //
1107 // EOF
1108 //