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