]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PHOS/AliPHOSv1.cxx
Writing histograms to file is temporarily disabled.
[u/mrichter/AliRoot.git] / PHOS / AliPHOSv1.cxx
index 0c268721b703b29845aa8a505864ba99ac55fe41..b99e522710377f6031b0aae1b98de7515809c64a 100644 (file)
 /* History of cvs commits:
  *
  * $Log$
+ * Revision 1.111  2007/07/24 09:41:19  morsch
+ * AliStack included for kKeepBit.
+ *
+ * Revision 1.110  2007/03/10 08:58:52  kharlov
+ * Protection for noCPV geometry
+ *
+ * Revision 1.109  2007/03/01 11:37:37  kharlov
+ * Strip units changed from 8x1 to 8x2 (T.Pocheptsov)
+ *
+ * Revision 1.108  2007/02/02 09:40:50  alibrary
+ * Includes required by ROOT head
+ *
+ * Revision 1.107  2007/02/01 10:34:47  hristov
+ * Removing warnings on Solaris x86
+ *
  * Revision 1.106  2006/11/14 17:11:15  hristov
  * Removing inheritances from TAttLine, TAttMarker and AliRndm in AliModule. The copy constructor and assignment operators are moved to the private part of the class and not implemented. The corresponding changes are propagated to the detectors
  *
@@ -41,6 +56,7 @@
 
 
 // --- ROOT system ---
+#include <TClonesArray.h>
 #include <TParticle.h>
 #include <TVirtualMC.h>
 
 #include "AliPHOSv1.h"
 #include "AliRun.h"
 #include "AliMC.h"
+#include "AliStack.h"
+#include "AliPHOSSimParam.h"
 
 ClassImp(AliPHOSv1)
 
 //____________________________________________________________________________
-AliPHOSv1::AliPHOSv1():
-  fLightYieldMean(0.),
-  fIntrinsicPINEfficiency(0.),
-  fLightYieldAttenuation(0.),
-  fRecalibrationFactor(0.),
-  fElectronsPerGeV(0.),
-  fAPDGain(0.),
-  fLightFactor(0.),
-  fAPDFactor(0.)
+AliPHOSv1::AliPHOSv1() : fCPVDigits("AliPHOSCPVDigit",20)
 {
   //Def ctor.
 }
 
 //____________________________________________________________________________
 AliPHOSv1::AliPHOSv1(const char *name, const char *title):
-  AliPHOSv0(name,title),
-  fLightYieldMean(0.),
-  fIntrinsicPINEfficiency(0.),
-  fLightYieldAttenuation(0.),
-  fRecalibrationFactor(0.),
-  fElectronsPerGeV(0.),
-  fAPDGain(0.),
-  fLightFactor(0.),
-  fAPDFactor(0.)
+  AliPHOSv0(name,title), fCPVDigits("AliPHOSCPVDigit",20)
 {
   //
   // We store hits :
@@ -98,31 +100,12 @@ AliPHOSv1::AliPHOSv1(const char *name, const char *title):
   // and the TreeD at the end of the event (branch is set in FinishEvent() ). 
   
   fHits= new TClonesArray("AliPHOSHit",1000) ;
+//   fCPVDigits("AliPHOSCPVDigit",20);
   gAlice->GetMCApp()->AddHitList(fHits) ; 
 
   fNhits = 0 ;
 
   fIshunt     =  2 ; // All hits are associated with primary particles
-
-  //Photoelectron statistics:
-  // The light yield is a poissonian distribution of the number of
-  // photons created in the PbWo4 crystal, calculated using following formula
-  // NumberOfPhotons = EnergyLost * LightYieldMean* APDEfficiency *
-  //              exp (-LightYieldAttenuation * DistanceToPINdiodeFromTheHit);
-  // LightYieldMean is parameter calculated to be over 47000 photons per GeV
-  // APDEfficiency is 0.02655
-  // k_0 is 0.0045 from Valery Antonenko
-  // The number of electrons created in the APD is
-  // NumberOfElectrons = APDGain * LightYield
-  // The APD Gain is 300
-  fLightYieldMean = 47000;
-  fIntrinsicPINEfficiency = 0.02655 ; //APD= 0.1875/0.1271 * 0.018 (PIN)
-  fLightYieldAttenuation  = 0.0045 ; 
-  fRecalibrationFactor    = 13.418/ fLightYieldMean ;
-  fElectronsPerGeV        = 2.77e+8 ;
-  fAPDGain                = 300. ;
-  fLightFactor            = fLightYieldMean * fIntrinsicPINEfficiency ; 
-  fAPDFactor              = (fRecalibrationFactor/100.) * fAPDGain ;   
 }
 
 //____________________________________________________________________________
@@ -201,22 +184,24 @@ void AliPHOSv1::StepManager(void)
   TLorentzVector pos      ;           // Lorentz vector of the track current position
   Int_t          copy     ;
 
-  TString name      =  GetGeometry()->GetName() ; 
-
   Int_t moduleNumber ;
   
-  static Int_t idPCPQ = gMC->VolId("PCPQ");
+  static Int_t idPCPQ = -1;
+  if (strstr(fTitle.Data(),"noCPV") == 0) 
+    idPCPQ = gMC->VolId("PCPQ");
+
   if( gMC->CurrentVolID(copy) == idPCPQ &&
-      (gMC->IsTrackEntering() ) &&
+     (gMC->IsTrackEntering() ) &&
       gMC->TrackCharge() != 0) {      
     
     gMC -> TrackPosition(pos);
-    
+
     Float_t xyzm[3], xyzd[3] ;
     Int_t i;
     for (i=0; i<3; i++) xyzm[i] = pos[i];
-    gMC -> Gmtod (xyzm, xyzd, 1);    // transform coordinate from master to daughter system
-    
+    gMC -> Gmtod (xyzm, xyzd, 1);    // transform coordinate from master to daughter system    
+
+
     Float_t        xyd[3]={0,0,0}   ;   //local position of the entering
     xyd[0]  = xyzd[0];
     xyd[1]  =-xyzd[2];
@@ -240,8 +225,8 @@ void AliPHOSv1::StepManager(void)
     gMC->CurrentVolOffID(3,moduleNumber);
     moduleNumber--;
     
-    TClonesArray *cpvDigits = new TClonesArray("AliPHOSCPVDigit",0);   // array of digits for current hit
-    CPVDigitize(pmom,xyd,cpvDigits);
+//     TClonesArray *cpvDigits = new TClonesArray("AliPHOSCPVDigit",0);   // array of digits for current hit
+    CPVDigitize(pmom,xyd,&fCPVDigits);
       
     Float_t xmean = 0;
     Float_t zmean = 0;
@@ -250,29 +235,29 @@ void AliPHOSv1::StepManager(void)
     
     // 2. go through the current digit list and sum digits in pads
     
-    ndigits = cpvDigits->GetEntriesFast();
+    ndigits = fCPVDigits.GetEntriesFast();
     for (idigit=0; idigit<ndigits-1; idigit++) {
-      AliPHOSCPVDigit  *cpvDigit1 = dynamic_cast<AliPHOSCPVDigit*>(cpvDigits->UncheckedAt(idigit));
+      AliPHOSCPVDigit  *cpvDigit1 = dynamic_cast<AliPHOSCPVDigit*>(fCPVDigits.UncheckedAt(idigit));
       Float_t x1 = cpvDigit1->GetXpad() ;
       Float_t z1 = cpvDigit1->GetYpad() ;
       for (Int_t jdigit=idigit+1; jdigit<ndigits; jdigit++) {
-       AliPHOSCPVDigit  *cpvDigit2 = dynamic_cast<AliPHOSCPVDigit*>(cpvDigits->UncheckedAt(jdigit));
+       AliPHOSCPVDigit  *cpvDigit2 = dynamic_cast<AliPHOSCPVDigit*>(fCPVDigits.UncheckedAt(jdigit));
        Float_t x2 = cpvDigit2->GetXpad() ;
        Float_t z2 = cpvDigit2->GetYpad() ;
        if (x1==x2 && z1==z2) {
          Float_t qsumpad = cpvDigit1->GetQpad() + cpvDigit2->GetQpad() ;
          cpvDigit2->SetQpad(qsumpad) ;
-         cpvDigits->RemoveAt(idigit) ;
+         fCPVDigits.RemoveAt(idigit) ;
        }
       }
     }
-    cpvDigits->Compress() ;
+    fCPVDigits.Compress() ;
     
     // 3. add digits to temporary hit list fTmpHits
     
-    ndigits = cpvDigits->GetEntriesFast();
+    ndigits = fCPVDigits.GetEntriesFast();
     for (idigit=0; idigit<ndigits; idigit++) {
-      AliPHOSCPVDigit  *cpvDigit = dynamic_cast<AliPHOSCPVDigit*>(cpvDigits->UncheckedAt(idigit));
+      AliPHOSCPVDigit  *cpvDigit = dynamic_cast<AliPHOSCPVDigit*>(fCPVDigits.UncheckedAt(idigit));
       relid[0] = moduleNumber + 1 ;                             // CPV (or PHOS) module number
       relid[1] =-1 ;                                            // means CPV
       relid[2] = cpvDigit->GetXpad() ;                          // column number of a pad
@@ -295,11 +280,7 @@ void AliPHOSv1::StepManager(void)
        qsum  += cpvDigit->GetQpad();
       }
     }
-    if (cpvDigits) {
-      cpvDigits->Delete();
-      delete cpvDigits;
-      cpvDigits=0;
-    }
+    fCPVDigits.Clear();
   }
 
  
@@ -311,10 +292,6 @@ void AliPHOSv1::StepManager(void)
     xyzte[1] = pos[1] ;
     xyzte[2] = pos[2] ;
 
-    Float_t global[3], local[3] ;
-    global[0] = pos[0] ;
-    global[1] = pos[1] ;
-    global[2] = pos[2] ;
     Float_t lostenergy = gMC->Edep(); 
     
     //Put in the TreeK particle entering PHOS and all its parents
@@ -349,24 +326,25 @@ void AliPHOSv1::StepManager(void)
       gMC->CurrentVolOffID(3, strip);
       Int_t cell ;
       gMC->CurrentVolOffID(2, cell);
-      
-      Int_t row = 1 + GetGeometry()->GetNZ() - strip % GetGeometry()->GetNZ() ;
-      Int_t col = (Int_t) TMath::Ceil((Double_t) strip/GetGeometry()->GetNZ()) -1 ;
-      
+
+      //Old formula for row is wrong. For example, I have strip 56 (28 for 2 x 8), row must be 1.
+      //But row == 1 + 56 - 56 % 56 == 57 (row == 1 + 28 - 28 % 28 == 29)
+      //Int_t row = 1 + GetGeometry()->GetEMCAGeometry()->GetNStripZ() - strip % (GetGeometry()->GetEMCAGeometry()->GetNStripZ()) ;
+      Int_t row = GetGeometry()->GetEMCAGeometry()->GetNStripZ() - (strip - 1) % (GetGeometry()->GetEMCAGeometry()->GetNStripZ()) ;
+      Int_t col = (Int_t) TMath::Ceil((Double_t) strip/(GetGeometry()->GetEMCAGeometry()->GetNStripZ())) -1 ;
+
+      // Absid for 8x2-strips. Looks nice :) 
       absid = (moduleNumber-1)*GetGeometry()->GetNCristalsInModule() + 
-       row + (col*GetGeometry()->GetEMCAGeometry()->GetNCellsInStrip() + cell-1)*GetGeometry()->GetNZ() ;
-      
-      gMC->Gmtod(global, local, 1) ;
+                   row * 2 + (col*GetGeometry()->GetEMCAGeometry()->GetNCellsXInStrip() + (cell - 1) / 2)*GetGeometry()->GetNZ() - (cell & 1 ? 1 : 0);
+
       
       //Calculates the light yield, the number of photons produced in the
       //crystal 
-      Float_t lightYield = gRandom->Poisson(fLightFactor * lostenergy *
-                                           exp(-fLightYieldAttenuation *
-                                               (local[1]+GetGeometry()->GetCrystalSize(1)/2.0 ))
-                                           ) ;
+      //There is no dependence of reponce on distance from energy deposition to APD
+      Float_t lightYield = gRandom->Poisson(AliPHOSSimParam::GetInstance()->GetLightFactor() * lostenergy) ;
 
       //Calculates de energy deposited in the crystal  
-      xyzte[4] = fAPDFactor * lightYield  ;
+      xyzte[4] = AliPHOSSimParam::GetInstance()->GetAPDFactor() * lightYield ;
       
       Int_t primary ;
       if(fIshunt == 2){
@@ -384,10 +362,9 @@ void AliPHOSv1::StepManager(void)
          part = gAlice->GetMCApp()->Particle(primary) ;
        }
       }
-      else
+      else{
        primary  =  gAlice->GetMCApp()->GetPrimary( gAlice->GetMCApp()->GetCurrentTrackNumber() ); 
-
-      
+      }
       
       // add current hit to the hit list
       // Info("StepManager","%d %d", primary, tracknumber) ; 
@@ -432,8 +409,6 @@ void AliPHOSv1::CPVDigitize (TLorentzVector p, Float_t *zxhit, TClonesArray *cpv
   Float_t pNorm = p.Py();
   Float_t eloss = kdEdx;
 
-//Info("CPVDigitize", "YVK : %f %f | %f %f %d", hitX, hitZ, pX, pZ, pNorm) ;
-
   Float_t dZY   = pZ/pNorm * GetGeometry()->GetCPVGasThickness();
   Float_t dXY   = pX/pNorm * GetGeometry()->GetCPVGasThickness();
   gRandom->Rannor(rnor1,rnor2);