]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - MUON/AliMUONDigitizerV3.cxx
Adding fClusterMap and fSharedMap of TPC clusters for the tracks
[u/mrichter/AliRoot.git] / MUON / AliMUONDigitizerV3.cxx
index 0f576e51209dcd7b7132fba354f350f69de3b405..29edc6a0cd79644165bca9cd405c5a077a383aa8 100644 (file)
 #include "AliMUONDigitizerV3.h"
 
 #include "AliCDBManager.h"
+#include "AliCodeTimer.h"
 #include "AliLog.h"
 #include "AliMUON.h"
 #include "AliMUONCalibrationData.h"
 #include "AliMUONConstants.h"
 #include "AliMUONDigit.h"
 #include "AliMUONLogger.h"
-#include "AliMUONTriggerEfficiencyCells.h"
 #include "AliMUONTriggerElectronics.h"
 #include "AliMUONTriggerStoreV1.h"
 #include "AliMUONVCalibParam.h"
 #include "AliMUONVDigitStore.h"
+#include "AliMpCDB.h"
 #include "AliMpCathodType.h"
 #include "AliMpConstants.h"
 #include "AliMpDEIterator.h"
 #include "AliMpDEManager.h"
-#include "AliMpDEManager.h"
 #include "AliMpIntPair.h"
 #include "AliMpPad.h"
 #include "AliMpSegmentation.h"
@@ -51,7 +51,9 @@
 #include <TString.h>
 #include <TSystem.h>
 
-///
+#include "AliMUONGeometryTransformer.h" //ADDED for trigger noise
+
+//-----------------------------------------------------------------------------
 /// \class AliMUONDigitizerV3
 /// The digitizer is performing the transformation to go from SDigits (digits
 /// w/o any electronic noise) to Digits (w/ electronic noise, and decalibration)
@@ -66,6 +68,7 @@
 /// here, at the digitization level).
 ///
 /// \author Laurent Aphecetche
+//-----------------------------------------------------------------------------
 
 namespace
 {
@@ -73,6 +76,12 @@ namespace
   {
     return static_cast<AliMUON*>(gAlice->GetModule("MUON"));
   }
+
+  //ADDED for trigger noise
+  const AliMUONGeometryTransformer* GetTransformer()
+  {
+      return muon()->GetGeometryTransformer();
+  }
 }
 
 const Double_t AliMUONDigitizerV3::fgkNSigmas=3;
@@ -83,15 +92,13 @@ ClassImp(AliMUONDigitizerV3)
 
 //_____________________________________________________________________________
 AliMUONDigitizerV3::AliMUONDigitizerV3(AliRunDigitizer* manager, 
-                                       Bool_t generateNoisyDigits)
+                                       Int_t generateNoisyDigits)
 : AliDigitizer(manager),
 fIsInitialized(kFALSE),
 fCalibrationData(0x0),
 fTriggerProcessor(0x0),
-fTriggerEfficiency(0x0),
-fGenerateNoisyDigitsTimer(),
-fExecTimer(),
 fNoiseFunction(0x0),
+fNoiseFunctionTrig(0x0),
   fGenerateNoisyDigits(generateNoisyDigits),
   fLogger(new AliMUONLogger(1000)),
 fTriggerStore(new AliMUONTriggerStoreV1),
@@ -101,8 +108,6 @@ fOutputDigitStore(0x0)
   /// Ctor.
 
   AliDebug(1,Form("AliRunDigitizer=%p",fManager));
-  fGenerateNoisyDigitsTimer.Start(kTRUE); fGenerateNoisyDigitsTimer.Stop();
-  fExecTimer.Start(kTRUE); fExecTimer.Stop();
 }
 
 //_____________________________________________________________________________
@@ -115,19 +120,11 @@ AliMUONDigitizerV3::~AliMUONDigitizerV3()
   delete fCalibrationData;
   delete fTriggerProcessor;
   delete fNoiseFunction;
+  delete fNoiseFunctionTrig;
   delete fTriggerStore;
   delete fDigitStore;
   delete fOutputDigitStore;
-  
-  if ( fGenerateNoisyDigits )
-  {
-    AliDebug(1, Form("Execution time for GenerateNoisyDigits() : R:%.2fs C:%.2fs",
-                 fGenerateNoisyDigitsTimer.RealTime(),
-                 fGenerateNoisyDigitsTimer.CpuTime()));
-  }
-  AliDebug(1, Form("Execution time for Exec() : R:%.2fs C:%.2fs",
-               fExecTimer.RealTime(),fExecTimer.CpuTime()));
+   
   AliInfo("Summary of messages");
   fLogger->Print();
   
@@ -140,25 +137,19 @@ AliMUONDigitizerV3::ApplyResponseToTrackerDigit(AliMUONVDigit& digit, Bool_t add
 {
   /// For tracking digits, starting from an ideal digit's charge, we :
   ///
-  /// - add some noise (thus leading to a realistic charge), if requested to do so
-  /// - divide by a gain (thus decalibrating the digit)
+  /// - "divide" by a gain (thus decalibrating the digit)
   /// - add a pedestal (thus decalibrating the digit)
+  /// - add some electronics noise (thus leading to a realistic adc), if requested to do so
   /// - sets the signal to zero if below 3*sigma of the noise
 
-  static const Int_t kMaxADC = (1<<12)-1; // We code the charge on a 12 bits ADC.
-  
-  Float_t signal = digit.Charge();
-  
-  if ( !addNoise )
-  {
-    digit.SetADC(TMath::Min(kMaxADC,TMath::Nint(signal)));
-    return;
-  }
+  Float_t charge = digit.Charge();
   
+  // We set the charge to 0, as the only relevant piece of information
+  // after Digitization is the ADC value.  
+  digit.SetCharge(0);
+    
   Int_t detElemId = digit.DetElemId();
-  
   Int_t manuId = digit.ManuId();
-  Int_t manuChannel = digit.ManuChannel();
   
   AliMUONVCalibParam* pedestal = fCalibrationData->Pedestals(detElemId,manuId);
   if (!pedestal)
@@ -166,12 +157,9 @@ AliMUONDigitizerV3::ApplyResponseToTrackerDigit(AliMUONVDigit& digit, Bool_t add
     fLogger->Log(Form("%s:%d:Could not get pedestal for DE=%4d manuId=%4d. Disabling.",
                       __FILE__,__LINE__,
                       detElemId,manuId));
-    digit.SetCharge(0);
     digit.SetADC(0);
     return;    
   }
-  Float_t pedestalMean = pedestal->ValueAsFloat(manuChannel,0);
-  Float_t pedestalSigma = pedestal->ValueAsFloat(manuChannel,1);
   
   AliMUONVCalibParam* gain = fCalibrationData->Gains(detElemId,manuId);
   if (!gain)
@@ -179,100 +167,19 @@ AliMUONDigitizerV3::ApplyResponseToTrackerDigit(AliMUONVDigit& digit, Bool_t add
     fLogger->Log(Form("%s:%d:Could not get gain for DE=%4d manuId=%4d. Disabling.",
                       __FILE__,__LINE__,
                       detElemId,manuId));
-    digit.SetCharge(0);
     digit.SetADC(0);
     return;        
   }    
-  Float_t gainMean = gain->ValueAsFloat(manuChannel,0);
-  
-  Float_t adcNoise = gRandom->Gaus(0.0,pedestalSigma);
-  
-  Int_t adc;
+
+  Int_t manuChannel = digit.ManuChannel();
   
-  if ( gainMean < 1E-6 )
-  {
-    AliError(Form("Got a too small gain %e for DE=%d manuId=%d manuChannel=%d. "
-                  "Setting signal to 0.",
-                  gainMean,detElemId,manuId,manuChannel));
-    adc = 0;
-  }
-  else
-  {
-    adc = TMath::Nint( signal / gainMean + pedestalMean + adcNoise);///
-    
-    if ( adc <= pedestalMean + fgkNSigmas*pedestalSigma ) 
-    {
-      adc = 0;
-    }
-  }
+
   
-  // be sure we stick to 12 bits.
-  if ( adc > kMaxADC )
-  {
-    adc = kMaxADC;
-  }
+  Int_t adc = DecalibrateTrackerDigit(*pedestal,*gain,manuChannel,charge,addNoise);
   
-  digit.SetCharge(adc);
   digit.SetADC(adc);
 }
 
-//_____________________________________________________________________________
-void 
-AliMUONDigitizerV3::ApplyResponseToTriggerDigit(const AliMUONVDigitStore& digitStore,
-                                                AliMUONVDigit& digit)
-{
-  /// \todo add comment
-
-  if ( !fTriggerEfficiency ) return;
-
-  if (digit.IsEfficiencyApplied()) return;
-
-  AliMUONVDigit* correspondingDigit = FindCorrespondingDigit(digitStore,digit);
-
-  if (!correspondingDigit) return; //reject bad correspondences
-
-  Int_t detElemId = digit.DetElemId();
-
-  AliMpSegmentation* segmentation = AliMpSegmentation::Instance();
-  const AliMpVSegmentation* segment[2] = 
-  {
-    segmentation->GetMpSegmentation(detElemId,AliMp::GetCathodType(digit.Cathode())), 
-    segmentation->GetMpSegmentation(detElemId,AliMp::GetCathodType(correspondingDigit->Cathode()))
-  };
-
-  AliMpPad pad[2] = 
-  {
-    segment[0]->PadByIndices(AliMpIntPair(digit.PadX(),digit.PadY()),kTRUE), 
-    segment[1]->PadByIndices(AliMpIntPair(correspondingDigit->PadX(),correspondingDigit->PadY()),kTRUE)
-  };
-
-  Int_t p0(1);
-  if (digit.Cathode()==0) p0=0;
-
-  AliMpIntPair location = pad[p0].GetLocation(0);
-  Int_t nboard = location.GetFirst();
-
-  Bool_t isTrig[2];
-
-  fTriggerEfficiency->IsTriggered(detElemId, nboard-1, 
-                                  isTrig[0], isTrig[1]);
-  digit.EfficiencyApplied(kTRUE);
-  correspondingDigit->EfficiencyApplied(kTRUE);
-
-  if (!isTrig[digit.Cathode()])
-  {
-         digit.SetCharge(0);
-  }
-  
-  if ( &digit != correspondingDigit )
-  {
-         if (!isTrig[correspondingDigit->Cathode()])
-    {
-      correspondingDigit->SetCharge(0);
-         }
-  }
-}
-
 //_____________________________________________________________________________
 void
 AliMUONDigitizerV3::ApplyResponse(const AliMUONVDigitStore& store,
@@ -296,17 +203,114 @@ AliMUONDigitizerV3::ApplyResponse(const AliMUONVDigitStore& store,
     {
       ApplyResponseToTrackerDigit(*digit,kAddNoise);
     }
-    else
-    {
-      ApplyResponseToTriggerDigit(store,*digit);
-    }
-    if ( digit->Charge() > 0  )
+
+    if ( digit->ADC() > 0  || digit->Charge() > 0 )
     {
       filteredStore.Add(*digit,AliMUONVDigitStore::kIgnore);
     }
   }
 }    
 
+//_____________________________________________________________________________
+Int_t 
+AliMUONDigitizerV3::DecalibrateTrackerDigit(const AliMUONVCalibParam& pedestals,
+                                            const AliMUONVCalibParam& gains,
+                                            Int_t channel,
+                                            Float_t charge,
+                                            Bool_t addNoise)
+{
+  /// Decalibrate (i.e. go from charge to adc) a tracker digit, given its
+  /// pedestal and gain parameters.
+  /// Must insure before calling that channel is valid (i.e. between 0 and
+  /// pedestals or gains->GetSize()-1, but also corresponding to a valid channel
+  /// otherwise results are not predictible...)
+
+  static const Int_t kMaxADC = (1<<12)-1; // We code the charge on a 12 bits ADC.
+  
+  Float_t pedestalMean = pedestals.ValueAsFloat(channel,0);
+  Float_t pedestalSigma = pedestals.ValueAsFloat(channel,1);
+  
+  Float_t a0 = gains.ValueAsFloat(channel,0);
+  Float_t a1 = gains.ValueAsFloat(channel,1);
+  Int_t thres = gains.ValueAsInt(channel,2);
+  Int_t qual = gains.ValueAsInt(channel,3);
+  if ( qual <= 0 ) return 0;
+  
+  Float_t chargeThres = a0*thres;
+  
+  Float_t padc(0); // (adc - ped) value
+  
+  if ( charge <= chargeThres || TMath::Abs(a1) < 1E-12 ) 
+  {
+    // linear part only
+    
+    if ( TMath::Abs(a0) > 1E-12 ) 
+    {
+      padc = charge/a0;    
+    }
+  }
+  else 
+  {
+    // linear + parabolic part
+    Double_t qt = chargeThres - charge;
+    Double_t delta = a0*a0-4*a1*qt;
+    if ( delta < 0 ) 
+    {
+      AliErrorClass(Form("delta=%e DE %d Manu %d Channel %d "
+                    " charge %e a0 %e a1 %e thres %d ped %e pedsig %e",
+                    delta,pedestals.ID0(),pedestals.ID1(),
+                    channel, charge, a0, a1, thres, pedestalMean, 
+                    pedestalSigma));      
+    }      
+    else
+    {
+      delta = TMath::Sqrt(delta);
+      
+      padc = ( ( -a0 + delta ) > 0 ? ( -a0 + delta ) : ( -a0 - delta ) );
+      
+      padc /= 2*a1;
+    
+      if ( padc < 0 )
+      {
+        if ( TMath::Abs(padc) > 1E-3) 
+        {
+          // this is more than a precision problem : let's signal it !
+          AliErrorClass(Form("padc=%e DE %d Manu %d Channel %d "
+                             " charge %e a0 %e a1 %e thres %d ped %e pedsig %e delta %e",
+                             padc,pedestals.ID0(),pedestals.ID1(),
+                             channel, charge, a0, a1, thres, pedestalMean, 
+                             pedestalSigma,delta));
+        }
+
+        // ok. consider we're just at thres, let it be zero.
+        padc = 0;
+      }
+
+      padc += thres;
+
+    }
+  }
+  
+  Int_t adc(0);
+  
+  if ( padc > 0 ) 
+  {
+    Float_t adcNoise = 0.0;
+    
+    if ( addNoise ) adcNoise = gRandom->Gaus(0.0,pedestalSigma);
+    
+    adc = TMath::Nint(padc + pedestalMean + adcNoise);
+  }
+  
+  // be sure we stick to 12 bits.
+  if ( adc > kMaxADC )
+  {
+    adc = kMaxADC;
+  }
+  
+  return adc;
+}
+
 //_____________________________________________________________________________
 void
 AliMUONDigitizerV3::Exec(Option_t*)
@@ -317,6 +321,8 @@ AliMUONDigitizerV3::Exec(Option_t*)
   /// Then we generate noise-only digits (for tracker only)
   /// And we finally generate the trigger outputs.
     
+  AliCodeTimerAuto("")
+  
   AliDebug(1, "Running digitizer.");
   
   if ( fManager->GetNinputs() == 0 )
@@ -331,8 +337,6 @@ AliMUONDigitizerV3::Exec(Option_t*)
     return;
   }
   
-  fExecTimer.Start(kFALSE);
-
   Int_t nInputFiles = fManager->GetNinputs();
   
   AliLoader* outputLoader = GetLoader(fManager->GetOutputFolderName());
@@ -383,6 +387,12 @@ AliMUONDigitizerV3::Exec(Option_t*)
     fOutputDigitStore = fDigitStore->Create();
   }
   
+  if ( fGenerateNoisyDigits>=2 )
+  {
+    // Generate noise-only digits for trigger.
+    GenerateNoisyDigitsForTrigger(*fDigitStore);
+  }
+
   ApplyResponse(*fDigitStore,*fOutputDigitStore);
   
   if ( fGenerateNoisyDigits )
@@ -418,30 +428,6 @@ AliMUONDigitizerV3::Exec(Option_t*)
   fTriggerStore->Clear();
   fDigitStore->Clear();
   fOutputDigitStore->Clear();
-  fExecTimer.Stop();
-}
-
-//_____________________________________________________________________________
-AliMUONVDigit* 
-AliMUONDigitizerV3::FindCorrespondingDigit(const AliMUONVDigitStore& digitStore,
-                                           AliMUONVDigit& digit) const
-{                                                
-  /// Find, if it exists, the digit corresponding to digit.Hit(), in the 
-  /// other cathode
-
-  TIter next(digitStore.CreateIterator());
-  AliMUONVDigit* d;
-  
-  while ( ( d = static_cast<AliMUONVDigit*>(next()) ) )
-  {
-    if ( d->DetElemId() == digit.DetElemId() &&
-         d->Hit() == digit.Hit() &&
-         d->Cathode() != digit.Cathode() )
-    {
-      return d;
-    }      
-  }    
-  return 0x0;
 }
 
 
@@ -453,6 +439,8 @@ AliMUONDigitizerV3::GenerateNoisyDigits(AliMUONVDigitStore& digitStore)
   /// have a signal above the noise cut (ped+n*sigma_ped), i.e. digits
   /// that are "only noise".
   
+  AliCodeTimerAuto("")
+  
   if ( !fNoiseFunction )
   {
     fNoiseFunction = new TF1("AliMUONDigitizerV3::fNoiseFunction","gaus",
@@ -461,8 +449,6 @@ AliMUONDigitizerV3::GenerateNoisyDigits(AliMUONVDigitStore& digitStore)
     fNoiseFunction->SetParameters(1,0,1);
   }
   
-  fGenerateNoisyDigitsTimer.Start(kFALSE);
-  
   for ( Int_t i = 0; i < AliMUONConstants::NTrackingCh(); ++i )
   {
     AliMpDEIterator it;
@@ -478,8 +464,6 @@ AliMUONDigitizerV3::GenerateNoisyDigits(AliMUONVDigitStore& digitStore)
       it.Next();
     }
   }
-  
-  fGenerateNoisyDigitsTimer.Stop();
 }
  
 //_____________________________________________________________________________
@@ -543,7 +527,7 @@ AliMUONDigitizerV3::GenerateNoisyDigitsForOneCathode(AliMUONVDigitStore& digitSt
     d->SetCharge(TMath::Nint(ped+pedestalMean+0.5));
     d->NoiseOnly(kTRUE);
     ApplyResponseToTrackerDigit(*d,kFALSE);
-    if ( d->Charge() > 0 )
+    if ( d->ADC() > 0 )
     {
       Bool_t ok = digitStore.Add(*d,AliMUONVDigitStore::kDeny);
       // this can happen (that we randomly chose a digit that is
@@ -567,6 +551,111 @@ AliMUONDigitizerV3::GenerateNoisyDigitsForOneCathode(AliMUONVDigitStore& digitSt
   }
 }
 
+
+//_____________________________________________________________________________
+void
+AliMUONDigitizerV3::GenerateNoisyDigitsForTrigger(AliMUONVDigitStore& digitStore)
+{
+  /// Generate noise-only digits for one cathode of one detection element.
+  /// Called by GenerateNoisyDigits()
+
+  if ( !fNoiseFunctionTrig )
+  {
+    fNoiseFunctionTrig = new TF1("AliMUONDigitizerV3::fNoiseFunctionTrig","landau",
+                                50.,270.);
+    
+    fNoiseFunctionTrig->SetParameters(3.91070e+02, 9.85026, 9.35881e-02);
+  }
+
+  AliMpPad pad[2];
+  AliMUONVDigit *d[2]={0x0};
+
+  for ( Int_t chamberId = AliMUONConstants::NTrackingCh(); chamberId < AliMUONConstants::NCh(); ++chamberId )
+  {
+  
+    Int_t nofNoisyPads = 50;
+
+    Float_t r=-1, fi = 0., gx, gy, x, y, z, xg01, yg01, zg, xg02, yg02;
+    AliMpDEIterator it;
+  
+    AliDebug(3,Form("Chamber %d nofNoisyPads %d",chamberId,nofNoisyPads));
+
+    for ( Int_t i = 0; i < nofNoisyPads; ++i )
+    {
+      //printf("Generating noise %i\n",i);
+       Int_t ix(-1);
+       Int_t iy(-1);
+       Bool_t isOk = kFALSE;
+       Int_t detElemId = -1;
+       do {
+         //r = gRandom->Landau(9.85026, 9.35881e-02);
+           r = fNoiseFunctionTrig->GetRandom();
+           fi = 2. * TMath::Pi() * gRandom->Rndm();
+           //printf("r = %f\tfi = %f\n", r, fi);
+           gx = r * TMath::Cos(fi);
+           gy = r * TMath::Sin(fi);
+
+           for ( it.First(chamberId); ! it.IsDone(); it.Next() ){
+               Int_t currDetElemId = it.CurrentDEId();
+               const AliMpVSegmentation* seg
+                   = AliMpSegmentation::Instance()->GetMpSegmentation(currDetElemId,AliMp::GetCathodType(0));
+               if (!seg) continue;
+               Float_t deltax = seg->Dimensions().X();
+               Float_t deltay = seg->Dimensions().Y();
+               GetTransformer()->Local2Global(currDetElemId, -deltax, -deltay, 0, xg01, yg01, zg);
+               GetTransformer()->Local2Global(currDetElemId,  deltax,  deltay, 0, xg02, yg02, zg);
+               Float_t xg1 = xg01, xg2 = xg02, yg1 = yg01, yg2 = yg02;
+               if(xg01>xg02){
+                   xg1 = xg02;
+                   xg2 = xg01;
+               }
+               if(yg01>yg02){
+                   yg1 = yg02;
+                   yg2 = yg01;
+               }
+               if(gx>=xg1 && gx<=xg2 && gy>=yg1 && gy<=yg2){
+                   detElemId = currDetElemId;
+                   GetTransformer()->Global2Local(detElemId, gx, gy, 0, x, y, z);
+                   pad[0] = seg->PadByPosition(TVector2(x,y),kFALSE);
+                   if(!pad[0].IsValid()) continue;
+                   isOk = kTRUE;
+                   break;
+               }
+           } // loop on slats
+       } while ( !isOk );
+
+       const AliMpVSegmentation* seg1
+           = AliMpSegmentation::Instance()->GetMpSegmentation(detElemId,AliMp::GetCathodType(1));
+       pad[1] = seg1->PadByPosition(TVector2(x,y),kFALSE);
+
+       for ( Int_t cathode = 0; cathode < 2; ++cathode ){
+         Int_t manuId = pad[cathode].GetLocation(0).GetFirst();
+         Int_t manuChannel = pad[cathode].GetLocation(0).GetSecond();    
+         d[cathode] = digitStore.CreateDigit(detElemId,manuId,manuChannel,cathode);
+         ix = pad[cathode].GetIndices().GetFirst();
+         iy = pad[cathode].GetIndices().GetSecond();
+         d[cathode]->SetPadXY(ix,iy);
+         //d[cathode].SetSignal(1);
+         //d[cathode].SetPhysicsSignal(0);
+         d[cathode]->SetCharge(1);
+         d[cathode]->NoiseOnly(kTRUE);
+         AliDebug(3,Form("Adding a pure noise digit :"));
+
+         Bool_t ok = digitStore.Add(*d[cathode],AliMUONVDigitStore::kDeny);
+         if (!ok)
+         {
+             fLogger->Log("Collision while adding TriggerNoise digit");
+         }
+         else
+         {
+             fLogger->Log("Added triggerNoise digit");
+         }
+       } //loop on cathodes
+    } // loop on noisy pads
+  } // loop on chambers
+}
+
+
 //_____________________________________________________________________________
 AliLoader*
 AliMUONDigitizerV3::GetLoader(const TString& folderName)
@@ -607,6 +696,16 @@ AliMUONDigitizerV3::Init()
   
   Int_t runnumber = AliCDBManager::Instance()->GetRun();
   
+  if ( ! AliMpCDB::LoadMpSegmentation()  ) 
+  {
+    AliFatal("Could not access mapping from OCDB !");
+  }
+  
+  if ( ! AliMpCDB::LoadDDLStore() ) 
+  {
+    AliFatal("Could not access DDL Store from OCDB !");
+  }
+  
   fCalibrationData = new AliMUONCalibrationData(runnumber);
   if ( !fCalibrationData->Pedestals() )
   {
@@ -618,20 +717,6 @@ AliMUONDigitizerV3::Init()
   }
   fTriggerProcessor = new AliMUONTriggerElectronics(fCalibrationData);
   
-  if ( muon()->GetTriggerEffCells() )
-  {
-    fTriggerEfficiency = fCalibrationData->TriggerEfficiency();
-    if ( fTriggerEfficiency )
-    {
-      AliDebug(1, "Will apply trigger efficiency");
-    }
-    else
-    {
-      AliFatal("I was requested to apply trigger efficiency, but I could "
-               "not get it !");
-    }
-  }
-  
   AliDebug(1, Form("Will %s generate noise-only digits for tracker",
                      (fGenerateNoisyDigits ? "":"NOT")));