First implementation of neural network PID

[u/mrichter/AliRoot.git] / TRD / AliTRDv1.cxx

diff --git a/TRD/AliTRDv1.cxx b/TRD/AliTRDv1.cxx
index a963c21c05986fb5169585e6d0d852a6f661dcb8..23c411dc7215798c19c3a87b950e2f8f9ade136e 100644 (file)
--- a/TRD/AliTRDv1.cxx
+++ b/TRD/AliTRDv1.cxx
@@ -30,15 +30,17 @@
 #include <TVector.h>
 #include <TVirtualMC.h>
 #include <TGeoManager.h>
+#include <TGeoPhysicalNode.h>
 
 #include "AliConst.h"
 #include "AliLog.h"
+#include "AliTrackReference.h"
 #include "AliMC.h"
 #include "AliRun.h"
 
 #include "AliTRDgeometry.h"
 #include "AliTRDhit.h"
-#include "AliTRDsim.h"
+#include "AliTRDsimTR.h"
 #include "AliTRDv1.h"
 
 ClassImp(AliTRDv1)
@@ -154,6 +156,9 @@ void AliTRDv1::AddAlignableVolumes() const
   //                         TRD/sm17/st4/pl5
   //
   for (Int_t isect = 0; isect < AliTRDgeometry::Nsect(); isect++) {
+
+    if (fGeometry->GetSMstatus(isect) == 0) continue;
+
     for (Int_t icham = 0; icham < AliTRDgeometry::Ncham(); icham++) {
       for (Int_t iplan = 0; iplan < AliTRDgeometry::Nplan(); iplan++) {
 
@@ -175,7 +180,25 @@ void AliTRDv1::AddAlignableVolumes() const
         symName += snApp2;
         symName += iplan;
 
-        gGeoManager->SetAlignableEntry(symName.Data(),volPath.Data());
+        TGeoPNEntry *alignableEntry = 
+         gGeoManager->SetAlignableEntry(symName.Data(),volPath.Data());
+
+       // Add the tracking to local matrix following the TPC example
+       if (alignableEntry) {
+         const char *path = alignableEntry->GetTitle();
+         if (!gGeoManager->cd(path)) {
+           AliFatal(Form("Volume path %s not valid!",path));
+         }
+         TGeoHMatrix *globMatrix = gGeoManager->GetCurrentMatrix();
+         Double_t sectorAngle = 20.0 * (isect % 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()));
+       }
 
       }
     }
@@ -301,7 +324,8 @@ void AliTRDv1::CreateTRhit(Int_t det)
           ,det
           ,posHit
           ,-q
-          ,kTRUE); 
+          ,gMC->TrackTime()*1.0e06
+          ,kTRUE);
 
   }
 
@@ -323,7 +347,7 @@ 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)
@@ -415,6 +439,10 @@ void AliTRDv1::StepManagerGeant()
   Bool_t   drRegion = kFALSE;
   Bool_t   amRegion = kFALSE;
 
+  TString  cIdPath;
+  Char_t   cIdSector[3];
+           cIdSector[2]  = 0;
+
   TString  cIdCurrent;
   TString  cIdSensDr = "J";
   TString  cIdSensAm = "K";
@@ -482,22 +510,17 @@ void AliTRDv1::StepManagerGeant()
       hits[1] = pos[1];
       hits[2] = pos[2];
 
-      // The sector number (0 - 17)
-      // The numbering goes clockwise and starts at y = 0
-      Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
-      if (phi < 90.0) {
-        phi = phi + 270.0;
-      }
-      else {
-        phi = phi -  90.0;
-      }
-      sec = ((Int_t) (phi / 20.0));
+      // The sector number (0 - 17), according to standard coordinate system
+      cIdPath      = gGeoManager->GetPath();
+      cIdSector[0] = cIdPath[21];
+      cIdSector[1] = cIdPath[22];
+      sec = atoi(cIdSector);
 
       // The plane and chamber number
       cIdChamber[0]   = cIdCurrent[2];
       cIdChamber[1]   = cIdCurrent[3];
       Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
-      cha = kNcham - ((Int_t) idChamber / kNplan) - 1;
+      cha = ((Int_t) idChamber / kNplan);
       pla = ((Int_t) idChamber % kNplan);
 
       // The detector number
@@ -510,7 +533,7 @@ 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
@@ -526,7 +549,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);
 
       }
 
@@ -598,6 +621,7 @@ void AliTRDv1::StepManagerGeant()
               ,det
               ,hits
               ,qTot
+              ,gMC->TrackTime()*1.0e06
               ,drRegion);
 
       }
@@ -673,6 +697,10 @@ void AliTRDv1::StepManagerErmilova()
   Bool_t   drRegion = kFALSE;
   Bool_t   amRegion = kFALSE;
 
+  TString  cIdPath;
+  Char_t   cIdSector[3];
+           cIdSector[2]  = 0;
+
   TString  cIdCurrent;
   TString  cIdSensDr = "J";
   TString  cIdSensAm = "K";
@@ -729,22 +757,17 @@ void AliTRDv1::StepManagerErmilova()
       hits[1] = pos[1];
       hits[2] = pos[2];
 
-      // The sector number (0 - 17)
-      // The numbering goes clockwise and starts at y = 0
-      Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
-      if (phi < 90.0) { 
-        phi = phi + 270.0;
-      }
-      else {
-        phi = phi -  90.0;
-      }
-      sec = ((Int_t) (phi / 20.0));
+      // The sector number (0 - 17), according to standard coordinate system
+      cIdPath      = gGeoManager->GetPath();
+      cIdSector[0] = cIdPath[21];
+      cIdSector[1] = cIdPath[22];
+      sec = atoi(cIdSector);
 
       // The plane and chamber number
       cIdChamber[0] = cIdCurrent[2];
       cIdChamber[1] = cIdCurrent[3];
       Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
-      cha = kNcham - ((Int_t) idChamber / kNplan) - 1;
+      cha = ((Int_t) idChamber / kNplan);
       pla = ((Int_t) idChamber % kNplan);
 
       // The detector number
@@ -757,7 +780,7 @@ 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
@@ -773,7 +796,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);
 
       }
 
@@ -792,6 +815,7 @@ void AliTRDv1::StepManagerErmilova()
                 ,det
                 ,hits
                 ,qTot
+                ,gMC->TrackTime()*1.0e06
                 ,kTRUE);
        }
         else {
@@ -799,6 +823,7 @@ void AliTRDv1::StepManagerErmilova()
                 ,det
                 ,hits
                 ,qTot
+                ,gMC->TrackTime()*1.0e06
                 ,kFALSE);
        }
 
@@ -871,11 +896,15 @@ void AliTRDv1::StepManagerFixedStep()
   Bool_t   drRegion = kFALSE;
   Bool_t   amRegion = kFALSE;
 
+  TString  cIdPath;
+  Char_t   cIdSector[3];
+           cIdSector[2]  = 0;
+
   TString  cIdCurrent;
   TString  cIdSensDr = "J";
   TString  cIdSensAm = "K";
   Char_t   cIdChamber[3];
-  cIdChamber[2] = 0;
+           cIdChamber[2] = 0;
 
   TLorentzVector pos;
   TLorentzVector mom;
@@ -920,22 +949,17 @@ void AliTRDv1::StepManagerFixedStep()
   hits[1] = pos[1];
   hits[2] = pos[2];
 
-  // The sector number (0 - 17)
-  // The numbering goes clockwise and starts at y = 0
-  Float_t phi = kRaddeg * TMath::ATan2(pos[0],pos[1]);
-  if (phi < 90.0) {
-    phi = phi + 270.0;
-  }
-  else {          
-    phi = phi -  90.0;
-  }
-  sec = ((Int_t) (phi / 20.0));
+  // The sector number (0 - 17), according to standard coordinate system
+  cIdPath      = gGeoManager->GetPath();
+  cIdSector[0] = cIdPath[21];
+  cIdSector[1] = cIdPath[22];
+  sec = atoi(cIdSector);
 
   // The plane and chamber number
   cIdChamber[0]   = cIdCurrent[2];
   cIdChamber[1]   = cIdCurrent[3];
   Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
-  cha = kNcham - ((Int_t) idChamber / kNplan) - 1;
+  cha = ((Int_t) idChamber / kNplan);
   pla = ((Int_t) idChamber % kNplan);
 
   // The detector number
@@ -951,7 +975,7 @@ 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
@@ -968,7 +992,7 @@ 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;
 
   }
@@ -983,6 +1007,7 @@ void AliTRDv1::StepManagerFixedStep()
           ,det
           ,hits
           ,qTot
+          ,gMC->TrackTime()*1.0e06
           ,drRegion);
   }
 
@@ -1183,6 +1208,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;