]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TRD/AliTRDv1.cxx
Some unuseful print commented
[u/mrichter/AliRoot.git] / TRD / AliTRDv1.cxx
index 0c238f9c51b1b11165d97fbed4d692f133ca0b11..77208541b8e3b437bfb2809393782e35121a93d4 100644 (file)
 #include <TVector.h>
 #include <TVirtualMC.h>
 #include <TGeoManager.h>
+#include <TGeoMatrix.h>
+#include <TGeoPhysicalNode.h>
 
 #include "AliConst.h"
 #include "AliLog.h"
+#include "AliTrackReference.h"
 #include "AliMC.h"
 #include "AliRun.h"
+#include "AliGeomManager.h"
 
 #include "AliTRDgeometry.h"
+#include "AliTRDCommonParam.h"
 #include "AliTRDhit.h"
-#include "AliTRDsim.h"
+#include "AliTRDsimTR.h"
 #include "AliTRDv1.h"
 
 ClassImp(AliTRDv1)
@@ -50,6 +55,7 @@ AliTRDv1::AliTRDv1()
   ,fTR(NULL)
   ,fTypeOfStepManager(0)
   ,fStepSize(0)
+  ,fWion(0)
   ,fDeltaE(NULL)
   ,fDeltaG(NULL)
   ,fTrackLength0(0)
@@ -68,6 +74,7 @@ AliTRDv1::AliTRDv1(const char *name, const char *title)
   ,fTR(NULL)
   ,fTypeOfStepManager(2)
   ,fStepSize(0.1)
+  ,fWion(0)
   ,fDeltaE(NULL)
   ,fDeltaG(NULL)
   ,fTrackLength0(0)
@@ -79,6 +86,17 @@ AliTRDv1::AliTRDv1(const char *name, const char *title)
 
   SetBufferSize(128000);
 
+  if      (AliTRDCommonParam::Instance()->IsXenon()) {
+    fWion = 23.53; // Ionization energy XeCO2 (85/15)
+  }
+  else if (AliTRDCommonParam::Instance()->IsArgon()) {
+    fWion = 27.21; // Ionization energy ArCO2 (82/18)
+  }
+  else {
+    AliFatal("Wrong gas mixture");
+    exit(1);
+  }
+
 }
 
 //_____________________________________________________________________________
@@ -117,14 +135,16 @@ void AliTRDv1::AddAlignableVolumes() const
   TString volPath;
   TString symName;
 
-  TString vpStr  = "ALIC_1/B077_1/BSEGMO";
-  TString vpApp1 = "_1/BTRD";
-  TString vpApp2 = "_1";
-  TString vpApp3 = "/UTR1_1/UTS1_1/UTI1_1/UT";
+  TString vpStr   = "ALIC_1/B077_1/BSEGMO";
+  TString vpApp1  = "_1/BTRD";
+  TString vpApp2  = "_1";
+  TString vpApp3a = "/UTR1_1/UTS1_1/UTI1_1/UT";
+  TString vpApp3b = "/UTR2_1/UTS2_1/UTI2_1/UT";
+  TString vpApp3c = "/UTR3_1/UTS3_1/UTI3_1/UT";
 
-  TString snStr  = "TRD/sm";
-  TString snApp1 = "/st";
-  TString snApp2 = "/pl";
+  TString snStr   = "TRD/sm";
+  TString snApp1  = "/st";
+  TString snApp2  = "/pl";
 
   //
   // The super modules
@@ -132,16 +152,16 @@ void AliTRDv1::AddAlignableVolumes() const
   //                           ...
   //                         TRD/sm17
   //
-  for (Int_t isect = 0; isect < AliTRDgeometry::Nsect(); isect++) {
+  for (Int_t isector = 0; isector < AliTRDgeometry::Nsector(); isector++) {
 
     volPath  = vpStr;
-    volPath += isect;
+    volPath += isector;
     volPath += vpApp1;
-    volPath += isect;
+    volPath += isector;
     volPath += vpApp2;
 
     symName  = snStr;
-    symName += Form("%02d",isect);
+    symName += Form("%02d",isector);
 
     gGeoManager->SetAlignableEntry(symName.Data(),volPath.Data());
 
@@ -153,29 +173,67 @@ void AliTRDv1::AddAlignableVolumes() const
   //                           ...
   //                         TRD/sm17/st4/pl5
   //
-  for (Int_t isect = 0; isect < AliTRDgeometry::Nsect(); isect++) {
-    for (Int_t icham = 0; icham < AliTRDgeometry::Ncham(); icham++) {
-      for (Int_t iplan = 0; iplan < AliTRDgeometry::Nplan(); iplan++) {
+  AliGeomManager::ELayerID idTRD1 = AliGeomManager::kTRD1;
+  Int_t layer, modUID;
+  
+  for (Int_t isector = 0; isector < AliTRDgeometry::Nsector(); isector++) {
 
-        Int_t idet = AliTRDgeometry::GetDetectorSec(iplan,icham);
+    if (fGeometry->GetSMstatus(isector) == 0) continue;
+
+    for (Int_t istack = 0; istack < AliTRDgeometry::Nstack(); istack++) {
+      for (Int_t ilayer = 0; ilayer < AliTRDgeometry::Nlayer(); ilayer++) {
+
+       layer = idTRD1 + ilayer;
+       modUID = AliGeomManager::LayerToVolUIDSafe(layer,isector*5+istack);
+
+        Int_t idet = AliTRDgeometry::GetDetectorSec(ilayer,istack);
 
         volPath  = vpStr;
-        volPath += isect;
+        volPath += isector;
         volPath += vpApp1;
-        volPath += isect;
+        volPath += isector;
         volPath += vpApp2;
-        volPath += vpApp3;
+        switch (isector) {
+        case 13:
+        case 14:
+        case 15:
+          if (istack == 2) {
+            continue;
+         }
+          volPath += vpApp3c;
+          break;
+        case 11:
+        case 12:
+          volPath += vpApp3b;
+          break;
+        default:
+          volPath += vpApp3a;
+       };
         volPath += Form("%02d",idet);
         volPath += vpApp2;
 
         symName  = snStr;
-        symName += Form("%02d",isect);
+        symName += Form("%02d",isector);
         symName += snApp1;
-        symName += icham;
+        symName += istack;
         symName += snApp2;
-        symName += iplan;
-
-        gGeoManager->SetAlignableEntry(symName.Data(),volPath.Data());
+        symName += ilayer;
+
+        TGeoPNEntry *alignableEntry = 
+         gGeoManager->SetAlignableEntry(symName.Data(),volPath.Data(),modUID);
+
+       // Add the tracking to local matrix following the TPC example
+       if (alignableEntry) {
+         TGeoHMatrix *globMatrix = alignableEntry->GetGlobalOrig();
+         Double_t sectorAngle = 20.0 * (isector % 18) + 10.0;
+         TGeoHMatrix *t2lMatrix  = new TGeoHMatrix();
+         t2lMatrix->RotateZ(sectorAngle);
+         t2lMatrix->MultiplyLeft(&(globMatrix->Inverse()));
+         alignableEntry->SetMatrix(t2lMatrix);
+       }
+       else {
+         AliError(Form("Alignable entry %s is not valid!",symName.Data()));
+       }
 
       }
     }
@@ -225,9 +283,6 @@ void AliTRDv1::CreateTRhit(Int_t det)
   // volume.
   //
 
-  // Ionization energy
-  const Float_t kWion        = 23.53;
-
   // Maximum number of TR photons per track
   const Int_t   kNTR         = 50;
 
@@ -268,9 +323,17 @@ void AliTRDv1::CreateTRhit(Int_t det)
 
     // The absorbtion cross sections in the drift gas
     // Gas-mixture (Xe/CO2)
-    Double_t muXe = fTR->GetMuXe(energyMeV);
+    Double_t muNo = 0.0;
+    if      (AliTRDCommonParam::Instance()->IsXenon()) {
+      muNo = fTR->GetMuXe(energyMeV);
+    }
+    else if (AliTRDCommonParam::Instance()->IsArgon()) {
+      muNo = fTR->GetMuAr(energyMeV);
+    }
     Double_t muCO = fTR->GetMuCO(energyMeV);
-    sigma = (0.85 * muXe + 0.15 * muCO) * fGasDensity * fTR->GetTemp();
+    sigma = (fGasNobleFraction * muNo + (1.0 - fGasNobleFraction) * muCO) 
+          * fGasDensity 
+          * fTR->GetTemp();
 
     // The distance after which the energy of the TR photon
     // is deposited.
@@ -293,16 +356,15 @@ void AliTRDv1::CreateTRhit(Int_t det)
     posHit[2] = pos[2] + mom[2] / pTot * absLength;
 
     // Create the charge 
-    Int_t q = ((Int_t) (energyeV / kWion));
+    Int_t q = ((Int_t) (energyeV / fWion));
 
     // Add the hit to the array. TR photon hits are marked 
     // by negative charge
-    // The hit time is needed for pile-up events
     AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
           ,det
           ,posHit
           ,-q
-         ,gMC->TrackTime()*1.0e06
+          ,gMC->TrackTime()*1.0e06
           ,kTRUE);
 
   }
@@ -325,10 +387,10 @@ void AliTRDv1::Init()
     AliInfo("TR simulation off");
   }
   else {
-    fTR = new AliTRDsim();
+    fTR = new AliTRDsimTR();
   }
 
-  // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
+  // First ionization potential (eV) for the gas mixture
   const Float_t kPoti = 12.1;
   // Maximum energy (50 keV);
   const Float_t kEend = 50000.0;
@@ -394,13 +456,15 @@ void AliTRDv1::StepManagerGeant()
   // to Bethe-Bloch. The energy distribution of the delta electrons follows
   // a spectrum taken from Geant3.
   //
+  // Works only for Xe/CO2!!
+  //
   // Version by A. Bercuci
   //
 
-  Int_t    pla = 0;
-  Int_t    cha = 0;
-  Int_t    sec = 0;
-  Int_t    det = 0;
+  Int_t    layer  = 0;
+  Int_t    stack  = 0;
+  Int_t    sector = 0;
+  Int_t    det    = 0;
   Int_t    iPdg;
   Int_t    qTot;
 
@@ -432,12 +496,11 @@ void AliTRDv1::StepManagerGeant()
 
   TArrayI        processes;
 
-  const Int_t    kNplan       = AliTRDgeometry::Nplan();
-  const Int_t    kNcham       = AliTRDgeometry::Ncham();
-  const Int_t    kNdetsec     = kNplan * kNcham;
+  const Int_t    kNlayer      = AliTRDgeometry::Nlayer();
+  const Int_t    kNstack      = AliTRDgeometry::Nstack();
+  const Int_t    kNdetsec     = kNlayer * kNstack;
 
   const Double_t kBig         = 1.0e+12; // Infinitely big
-  const Float_t  kWion        = 23.53;   // Ionization energy
   const Float_t  kPTotMaxEl   = 0.002;   // Maximum momentum for e+ e- g
 
   // Minimum energy for the step size adjustment
@@ -492,17 +555,17 @@ void AliTRDv1::StepManagerGeant()
       cIdPath      = gGeoManager->GetPath();
       cIdSector[0] = cIdPath[21];
       cIdSector[1] = cIdPath[22];
-      sec = atoi(cIdSector);
+      sector = atoi(cIdSector);
 
-      // The plane and chamber number
+      // The layer and stack number
       cIdChamber[0]   = cIdCurrent[2];
       cIdChamber[1]   = cIdCurrent[3];
       Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
-      cha = ((Int_t) idChamber / kNplan);
-      pla = ((Int_t) idChamber % kNplan);
+      stack = ((Int_t) idChamber / kNlayer);
+      layer = ((Int_t) idChamber % kNlayer);
 
       // The detector number
-      det = fGeometry->GetDetector(pla,cha,sec);
+      det = fGeometry->GetDetector(layer,stack,sector);
 
       // Special hits only in the drift region
       if      ((drRegion) &&
@@ -511,13 +574,13 @@ void AliTRDv1::StepManagerGeant()
         // Create a track reference at the entrance of each
         // chamber that contains the momentum components of the particle
         gMC->TrackMomentum(mom);
-        AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
+        AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber(), AliTrackReference::kTRD);
 
         // Create the hits from TR photons if electron/positron is
         // entering the drift volume
         if ((fTR) && 
             (TMath::Abs(gMC->TrackPid()) == kPdgElectron)) {
-          CreateTRhit(det);    
+          CreateTRhit(det);
         }
 
       }
@@ -527,7 +590,7 @@ void AliTRDv1::StepManagerGeant()
         // Create a track reference at the exit of each
         // chamber that contains the momentum components of the particle
         gMC->TrackMomentum(mom);
-        AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
+        AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber(), AliTrackReference::kTRD);
 
       }
 
@@ -546,7 +609,7 @@ void AliTRDv1::StepManagerGeant()
         pid = 0;
       }
       else {
-        pid =  processes[nofprocesses-1];              
+        pid = processes[nofprocesses-1];               
       }                
                
       // Generate Edep according to GEANT parametrisation
@@ -592,15 +655,14 @@ void AliTRDv1::StepManagerGeant()
       // Generate the electron cluster size
       if (eDelta > 0.0) {
 
-        qTot = ((Int_t) (eDelta / kWion) + 1);
+        qTot = ((Int_t) (eDelta / fWion) + 1);
 
         // Create a new dEdx hit
-        // The hit time is needed for pile-up events
         AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
               ,det
               ,hits
               ,qTot
-             ,gMC->TrackTime()*1.0e06
+              ,gMC->TrackTime()*1.0e06
               ,drRegion);
 
       }
@@ -654,11 +716,13 @@ void AliTRDv1::StepManagerErmilova()
   // to Bethe-Bloch. The energy distribution of the delta electrons follows
   // a spectrum taken from Ermilova et al.
   //
+  // Works only for Xe/CO2!!
+  //
 
-  Int_t    pla = 0;
-  Int_t    cha = 0;
-  Int_t    sec = 0;
-  Int_t    det = 0;
+  Int_t    layer  = 0;
+  Int_t    stack  = 0;
+  Int_t    sector = 0;
+  Int_t    det    = 0;
   Int_t    iPdg;
   Int_t    qTot;
 
@@ -689,12 +753,11 @@ void AliTRDv1::StepManagerErmilova()
   TLorentzVector pos;
   TLorentzVector mom;
 
-  const Int_t    kNplan       = AliTRDgeometry::Nplan();
-  const Int_t    kNcham       = AliTRDgeometry::Ncham();
-  const Int_t    kNdetsec     = kNplan * kNcham;
+  const Int_t    kNlayer      = AliTRDgeometry::Nlayer();
+  const Int_t    kNstack      = AliTRDgeometry::Nstack();
+  const Int_t    kNdetsec     = kNlayer * kNstack;
 
   const Double_t kBig         = 1.0e+12; // Infinitely big
-  const Float_t  kWion        = 23.53;   // Ionization energy
   const Float_t  kPTotMaxEl   = 0.002;   // Maximum momentum for e+ e- g 
 
   // Minimum energy for the step size adjustment
@@ -740,17 +803,17 @@ void AliTRDv1::StepManagerErmilova()
       cIdPath      = gGeoManager->GetPath();
       cIdSector[0] = cIdPath[21];
       cIdSector[1] = cIdPath[22];
-      sec = atoi(cIdSector);
+      sector = atoi(cIdSector);
 
       // The plane and chamber number
       cIdChamber[0] = cIdCurrent[2];
       cIdChamber[1] = cIdCurrent[3];
       Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
-      cha = ((Int_t) idChamber / kNplan);
-      pla = ((Int_t) idChamber % kNplan);
+      stack = ((Int_t) idChamber / kNlayer);
+      layer = ((Int_t) idChamber % kNlayer);
 
       // The detector number
-      det = fGeometry->GetDetector(pla,cha,sec);
+      det = fGeometry->GetDetector(layer,stack,sector);
 
       // Special hits only in the drift region
       if      ((drRegion) &&
@@ -759,13 +822,13 @@ void AliTRDv1::StepManagerErmilova()
         // Create a track reference at the entrance of each
         // chamber that contains the momentum components of the particle
         gMC->TrackMomentum(mom);
-        AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
+        AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber(), AliTrackReference::kTRD);
 
         // Create the hits from TR photons if electron/positron is
         // entering the drift volume
         if ((fTR) && 
             (TMath::Abs(gMC->TrackPid()) == kPdgElectron)) {
-          CreateTRhit(det);    
+          CreateTRhit(det);
         }
 
       }
@@ -775,7 +838,7 @@ void AliTRDv1::StepManagerErmilova()
         // Create a track reference at the exit of each
         // chamber that contains the momentum components of the particle
         gMC->TrackMomentum(mom);
-        AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
+        AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber(), AliTrackReference::kTRD);
 
       }
 
@@ -786,16 +849,15 @@ void AliTRDv1::StepManagerErmilova()
       // Generate the electron cluster size
       if (eDelta > 0.0) {
 
-        qTot = ((Int_t) (eDelta / kWion) + 1);
+        qTot = ((Int_t) (eDelta / fWion) + 1);
 
        // Create a new dEdx hit
-        // The hit time is needed for pile-up events
         if (drRegion) {
           AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
                 ,det
                 ,hits
                 ,qTot
-               ,gMC->TrackTime()*1.0e06
+                ,gMC->TrackTime()*1.0e06
                 ,kTRUE);
        }
         else {
@@ -803,7 +865,7 @@ void AliTRDv1::StepManagerErmilova()
                 ,det
                 ,hits
                 ,qTot
-               ,gMC->TrackTime()*1.0e06
+                ,gMC->TrackTime()*1.0e06
                 ,kFALSE);
        }
 
@@ -860,14 +922,16 @@ void AliTRDv1::StepManagerFixedStep()
   // along its path across the drift volume. The step size is fixed in
   // this version of the step manager.
   //
+  // Works for Xe/CO2 as well as Ar/CO2
+  //
 
   // PDG code electron
   const Int_t   kPdgElectron = 11;
 
-  Int_t    pla = 0;
-  Int_t    cha = 0;
-  Int_t    sec = 0;
-  Int_t    det = 0;
+  Int_t    layer  = 0;
+  Int_t    stack  = 0;
+  Int_t    sector = 0;
+  Int_t    det    = 0;
   Int_t    qTot;
 
   Float_t  hits[3];
@@ -889,13 +953,11 @@ void AliTRDv1::StepManagerFixedStep()
   TLorentzVector pos;
   TLorentzVector mom;
 
-  const Int_t    kNplan       = AliTRDgeometry::Nplan();
-  const Int_t    kNcham       = AliTRDgeometry::Ncham();
-  const Int_t    kNdetsec     = kNplan * kNcham;
+  const Int_t    kNlayer      = AliTRDgeometry::Nlayer();
+  const Int_t    kNstack      = AliTRDgeometry::Nstack();
+  const Int_t    kNdetsec     = kNlayer * kNstack;
 
   const Double_t kBig         = 1.0e+12;
-
-  const Float_t  kWion        = 23.53;   // Ionization energy
   const Float_t  kEkinMinStep = 1.0e-5;  // Minimum energy for the step size adjustment
 
   // Set the maximum step size to a very large number for all 
@@ -933,19 +995,19 @@ void AliTRDv1::StepManagerFixedStep()
   cIdPath      = gGeoManager->GetPath();
   cIdSector[0] = cIdPath[21];
   cIdSector[1] = cIdPath[22];
-  sec = atoi(cIdSector);
+  sector = atoi(cIdSector);
 
   // The plane and chamber number
   cIdChamber[0]   = cIdCurrent[2];
   cIdChamber[1]   = cIdCurrent[3];
   Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
-  cha = ((Int_t) idChamber / kNplan);
-  pla = ((Int_t) idChamber % kNplan);
+  stack = ((Int_t) idChamber / kNlayer);
+  layer = ((Int_t) idChamber % kNlayer);
 
   // The detector number
-  det = fGeometry->GetDetector(pla,cha,sec);
+  det = fGeometry->GetDetector(layer,stack,sector);
 
-  // 0:InFlight 1:Entering 2:Exiting
+  // 0: InFlight 1:Entering 2:Exiting
   Int_t trkStat = 0;
 
   // Special hits only in the drift region
@@ -955,14 +1017,14 @@ void AliTRDv1::StepManagerFixedStep()
     // Create a track reference at the entrance of each
     // chamber that contains the momentum components of the particle
     gMC->TrackMomentum(mom);
-    AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
+    AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber(), AliTrackReference::kTRD);
     trkStat = 1;
 
     // Create the hits from TR photons if electron/positron is
     // entering the drift volume
     if ((fTR) && 
         (TMath::Abs(gMC->TrackPid()) == kPdgElectron)) {
-      CreateTRhit(det);    
+      CreateTRhit(det);
     }
 
   }
@@ -972,23 +1034,22 @@ void AliTRDv1::StepManagerFixedStep()
     // Create a track reference at the exit of each
     // chamber that contains the momentum components of the particle
     gMC->TrackMomentum(mom);
-    AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
+    AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber(), AliTrackReference::kTRD);
     trkStat = 2;
 
   }
   
   // Calculate the charge according to GEANT Edep
   // Create a new dEdx hit
-  // The hit time is needed for pile-up events
   eDep = TMath::Max(gMC->Edep(),0.0) * 1.0e+09;
-  qTot = (Int_t) (eDep / kWion);
+  qTot = (Int_t) (eDep / fWion);
   if ((qTot) ||
       (trkStat)) {
     AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
           ,det
           ,hits
           ,qTot
-         ,gMC->TrackTime()*1.0e06
+          ,gMC->TrackTime()*1.0e06
           ,drRegion);
   }
 
@@ -1189,6 +1250,10 @@ Double_t IntSpecGeant(Double_t *x, Double_t *)
   Int_t    i;
   Double_t energy = x[0];
 
+  if (energy >= arre[npts-1]) {
+    return 0.0;
+  }
+
   for (i = 0; i < npts; i++) {
     if (energy < arre[i]) {
       break;