MUON/AliMUONDigitizerV3.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
  16 // $Id$
  17
  18
  19 #include "AliMUONDigitizerV3.h"
  20
  21 #include "AliMUON.h"
  22 #include "AliMUONCalibrationData.h"
  23 #include "AliMUONConstants.h"
  24 #include "AliMUONDigit.h"
  25 #include "AliMUONLogger.h"
  26 #include "AliMUONTriggerElectronics.h"
  27 #include "AliMUONTriggerStoreV1.h"
  28 #include "AliMUONVCalibParam.h"
  29 #include "AliMUONVDigitStore.h"
  30 #include "AliMUONGeometryTransformer.h" //ADDED for trigger noise
  31 #include "AliMUONRecoParam.h"
  32 #include "AliMUONTriggerChamberEfficiency.h"
  33 #include "AliMUONTriggerUtilities.h"
  34
  35 #include "AliMpCDB.h"
  36 #include "AliMpSegmentation.h"
  37 #include "AliMpCathodType.h"
  38 #include "AliMpConstants.h"
  39 #include "AliMpDEIterator.h"
  40 #include "AliMpDEManager.h"
  41 #include "AliMpPad.h"
  42 #include "AliMpStationType.h"
  43 #include "AliMpVSegmentation.h"
  44 #include "AliMpDDLStore.h"
  45
  46 #include "AliCDBManager.h"
  47 #include "AliCodeTimer.h"
  48 #include "AliLog.h"
  49 #include "AliRun.h"
  50 #include "AliDigitizationInput.h"
  51 #include "AliLoader.h"
  52 #include "AliRunLoader.h"
  53
  54 #include <Riostream.h>
  55 #include <TF1.h>
  56 #include <TFile.h>
  57 #include <TMath.h>
  58 #include <TRandom.h>
  59 #include <TString.h>
  60 #include <TSystem.h>
  61 #include <TTree.h>
  62
  63 //-----------------------------------------------------------------------------
  64 /// \class AliMUONDigitizerV3
  65 ///
  66 /// The digitizer is performing the transformation to go from SDigits (digits
  67 /// w/o any electronic noise) to Digits (w/ electronic noise, and decalibration)
  68 ///
  69 /// The decalibration is performed by doing the reverse operation of the
  70 /// calibration, that is we do (Signal+pedestal)/gain -> ADC
  71 ///
  72 /// Note also that the digitizer takes care of merging sdigits that belongs
  73 /// to the same pad, either because we're merging several input sdigit files
  74 /// or with a single file because the sdigitizer does not merge sdigits itself
  75 /// (for performance reason mainly, and because anyway we know we have to do it
  76 /// here, at the digitization level).
  77 ///
  78 /// August 2011. In order to remove the need for specific MC OCDB storages,
  79 /// we're introducing a dependence of simulation on AliMUONRecoParam (stored
  80 /// in MUON/Calib/RecoParam in OCDB), which is normally (or conceptually, if
  81 /// you will) only a reconstruction object. That's not a pretty solution, but,
  82 /// well, we have to do it...
  83 /// This dependence comes from the fact that we must know how to decalibrate
  84 /// the digits, so that the decalibration (done here) - calibration (done during
  85 /// reco) process is (as much as possible) neutral.
  86 ///
  87 ///
  88 /// \author Laurent Aphecetche
  89 ///
  90 //-----------------------------------------------------------------------------
  91
  92 namespace
  93 {
  94   AliMUON* muon()
  95   {
  96     return static_cast<AliMUON*>(gAlice->GetModule("MUON"));
  97   }
  98
  99   //ADDED for trigger noise
 100   const AliMUONGeometryTransformer* GetTransformer()
 101   {
 102       return muon()->GetGeometryTransformer();
 103   }
 104 }
 105
 106 Double_t AliMUONDigitizerV3::fgNSigmas = 4.0;
 107
 108 /// \cond CLASSIMP
 109 ClassImp(AliMUONDigitizerV3)
 110 /// \endcond
 111
 112 //_____________________________________________________________________________
 113 AliMUONDigitizerV3::AliMUONDigitizerV3(AliDigitizationInput* digInput,
 114                                        Int_t generateNoisyDigits)
 115 : AliDigitizer(digInput),
 116 fIsInitialized(kFALSE),
 117 fCalibrationData(0x0),
 118 fTriggerProcessor(0x0),
 119 fNoiseFunctionTrig(0x0),
 120 fGenerateNoisyDigits(generateNoisyDigits),
 121 fLogger(new AliMUONLogger(1000)),
 122 fTriggerStore(new AliMUONTriggerStoreV1),
 123 fDigitStore(0x0),
 124 fOutputDigitStore(0x0),
 125 fInputDigitStores(0x0),
 126 fRecoParam(0x0),
 127 fTriggerEfficiency(0x0),
 128 fTriggerUtilities(0x0),
 129 fEfficiencyResponse(2*AliMUONConstants::NTriggerCh()*AliMUONConstants::NTriggerCircuit())
 130 {
 131   /// Ctor.
 132
 133   AliDebug(1,Form("AliDigitizationInput=%p",fDigInput));
 134
 135 }
 136
 137 //_____________________________________________________________________________
 138 AliMUONDigitizerV3::~AliMUONDigitizerV3()
 139 {
 140   /// Dtor. Note we're the owner of some pointers.
 141
 142   AliDebug(1,"dtor");
 143
 144  // delete fCalibrationData;
 145   delete fTriggerProcessor;
 146   delete fNoiseFunctionTrig;
 147   delete fTriggerStore;
 148   delete fDigitStore;
 149   delete fOutputDigitStore;
 150   delete fInputDigitStores;
 151   delete fTriggerUtilities;
 152
 153   AliInfo("Summary of messages");
 154   fLogger->Print();
 155
 156   delete fLogger;
 157 }
 158
 159 //_____________________________________________________________________________
 160 void
 161 AliMUONDigitizerV3::ApplyResponseToTrackerDigit(AliMUONVDigit& digit, Bool_t addNoise)
 162 {
 163   /// For tracking digits, starting from an ideal digit's charge, we :
 164   ///
 165   /// - "divide" by a gain (thus decalibrating the digit)
 166   /// - add a pedestal (thus decalibrating the digit)
 167   /// - add some electronics noise (thus leading to a realistic adc), if requested to do so
 168   /// - sets the signal to zero if below 3*sigma of the noise
 169
 170   Float_t charge = digit.Charge();
 171
 172   if (!digit.IsChargeInFC())
 173   {
 174     charge *= AliMUONConstants::DefaultADC2MV()*AliMUONConstants::DefaultA0()*AliMUONConstants::DefaultCapa();
 175     fLogger->Log("CHECK ME ! WAS NOT SUPPOSED TO BE HERE !!! ARE YOU RECONSTRUCTING OLD SIMULATIONS ? ");
 176     AliError("CHECK ME ! WAS NOT SUPPOSED TO BE HERE !!! ARE YOU RECONSTRUCTING OLD SIMULATIONS ? ");
 177   }
 178
 179   // We set the charge to 0, as the only relevant piece of information
 180   // after Digitization is the ADC value.
 181   digit.SetCharge(0);
 182
 183   Int_t detElemId = digit.DetElemId();
 184   Int_t manuId = digit.ManuId();
 185
 186   AliMUONVCalibParam* pedestal = fCalibrationData->Pedestals(detElemId,manuId);
 187   if (!pedestal)
 188   {
 189     fLogger->Log(Form("%s:%d:Could not get pedestal for DE=%4d manuId=%4d. Disabling.",
 190                       __FILE__,__LINE__,
 191                       detElemId,manuId));
 192     digit.SetADC(0);
 193     return;
 194   }
 195
 196   Int_t manuChannel = digit.ManuChannel();
 197
 198   if ( pedestal->ValueAsFloat(manuChannel,0) == AliMUONVCalibParam::InvalidFloatValue() ||
 199       pedestal->ValueAsFloat(manuChannel,1) == AliMUONVCalibParam::InvalidFloatValue() )
 200   {
 201     // protection against invalid pedestal value
 202     digit.SetADC(0);
 203     return;
 204   }
 205
 206   TString calibrationMode(fRecoParam->GetCalibrationMode());
 207   calibrationMode.ToUpper();
 208
 209   AliMUONVCalibParam* gain = fCalibrationData->Gains(detElemId,manuId);
 210   if (!gain)
 211   {
 212
 213     if (!calibrationMode.Contains("NOGAIN") )
 214     {
 215       fLogger->Log(Form("%s:%d:Could not get gain for DE=%4d manuId=%4d. Disabling.",
 216                         __FILE__,__LINE__,
 217                         detElemId,manuId));
 218       digit.SetADC(0);
 219       return;
 220     }
 221   }
 222
 223   Int_t adc = DecalibrateTrackerDigit(*pedestal,gain,manuChannel,charge,addNoise,
 224                                       digit.IsNoiseOnly(),
 225                                       calibrationMode);
 226
 227   digit.SetADC(adc);
 228 }
 229
 230
 231 //_____________________________________________________________________________
 232 void
 233 AliMUONDigitizerV3::ApplyResponseToTriggerDigit(AliMUONVDigit& digit)
 234 {
 235   /// For trigger digits, starting from an ideal digit, we :
 236   ///
 237   /// - apply efficiency (on demand)
 238   /// - apply trigger masks
 239
 240   Int_t detElemId = digit.DetElemId();
 241   Int_t localCircuit = digit.ManuId();
 242   Int_t strip = digit.ManuChannel();
 243   Int_t cathode = digit.Cathode();
 244   Int_t trigCh = detElemId/100 - 11;
 245
 246   Int_t arrayIndex = GetArrayIndex(cathode, trigCh, localCircuit);
 247
 248   // Trigger chamber efficiency
 249   if ( fTriggerEfficiency ) {
 250     if ( fEfficiencyResponse[arrayIndex] < 0 ) {
 251       Bool_t isTrig[2] = {kTRUE, kTRUE};
 252       fTriggerEfficiency->IsTriggered(detElemId, localCircuit, isTrig[0], isTrig[1]);
 253       Int_t arrayIndexBend = GetArrayIndex(0, trigCh, localCircuit);
 254       Int_t arrayIndexNonBend = GetArrayIndex(1, trigCh, localCircuit);
 255       fEfficiencyResponse[arrayIndexBend] = isTrig[0];
 256       fEfficiencyResponse[arrayIndexNonBend] = isTrig[1];
 257     }
 258     AliDebug(1,Form("ch %i  cath %i  board %i  strip %i  efficiency %i\n", trigCh, cathode, localCircuit, strip, fEfficiencyResponse[arrayIndex]));
 259     if ( fEfficiencyResponse[arrayIndex] == 0 ) {
 260       digit.SetCharge(0);
 261       digit.SetADC(0);
 262       //AliDebug(1,Form("ch %i  cath %i  board %i  strip %i  NOT efficient\n", trigCh, cathode, localCircuit, strip));
 263       return;
 264     }
 265   }
 266
 267   // Masked channels
 268   Bool_t isMasked = fTriggerUtilities->IsMasked(digit);
 269   AliDebug(1,Form("ch %i  cath %i  board %i  strip %i  mask %i\n", trigCh, cathode, localCircuit, strip, !isMasked));
 270   if ( isMasked ) {
 271     digit.SetCharge(0);
 272     digit.SetADC(0);
 273     //AliDebug(1,Form("ch %i  cath %i  board %i  strip %i  masked\n", trigCh, cathode, localCircuit, strip));
 274     return;
 275   }
 276 }
 277
 278
 279
 280 //_____________________________________________________________________________
 281 void
 282 AliMUONDigitizerV3::ApplyResponse(const AliMUONVDigitStore& store,
 283                                   AliMUONVDigitStore& filteredStore)
 284 {
 285   /// Loop over all chamber digits, and apply the response to them
 286   /// Note that this method may remove digits.
 287
 288   filteredStore.Clear();
 289
 290   const Bool_t kAddNoise = kTRUE;
 291
 292   TIter next(store.CreateIterator());
 293   AliMUONVDigit* digit;
 294
 295   if ( fTriggerEfficiency ) fEfficiencyResponse.Reset(-1);
 296
 297   while ( ( digit = static_cast<AliMUONVDigit*>(next()) ) )
 298   {
 299     AliMp::StationType stationType = AliMpDEManager::GetStationType(digit->DetElemId());
 300
 301     if ( stationType != AliMp::kStationTrigger )
 302     {
 303       Bool_t addNoise = kAddNoise;
 304       if (digit->IsConverted()) addNoise = kFALSE; // No need to add extra noise to a converted real digit
 305       ApplyResponseToTrackerDigit(*digit,addNoise);
 306     }
 307     else {
 308       ApplyResponseToTriggerDigit(*digit);
 309     }
 310
 311     if ( digit->ADC() > 0  || digit->Charge() > 0 )
 312     {
 313       filteredStore.Add(*digit,AliMUONVDigitStore::kIgnore);
 314     }
 315   }
 316 }
 317
 318 //_____________________________________________________________________________
 319 Int_t
 320 AliMUONDigitizerV3::DecalibrateTrackerDigit(const AliMUONVCalibParam& pedestals,
 321                                             const AliMUONVCalibParam* gains,
 322                                             Int_t channel,
 323                                             Float_t charge,
 324                                             Bool_t addNoise,
 325                                             Bool_t noiseOnly,
 326                                             const TString& calibrationMode)
 327 {
 328   /// Decalibrate (i.e. go from charge to adc) a tracker digit, given its
 329   /// pedestal and gain parameters.
 330   /// Must insure before calling that channel is valid (i.e. between 0 and
 331   /// pedestals or gains->GetSize()-1, but also corresponding to a valid channel
 332   /// otherwise results are not predictible...)
 333   ///
 334   /// This method is completely tied to what happens in its sister method :
 335   /// AliMUONDigitCalibrator::CalibrateDigit, which is doing the reverse work...
 336   ///
 337
 338   static const Int_t kMaxADC = (1<<12)-1; // We code the charge on a 12 bits ADC.
 339
 340   Bool_t nogain = calibrationMode.Contains("NOGAIN");
 341
 342   Float_t a1(0.0);
 343   Int_t thres(4095);
 344   Int_t qual(0xF);
 345   Float_t capa(AliMUONConstants::DefaultCapa()); // capa = 0.2 and a0 = 1.25
 346   Float_t a0(AliMUONConstants::DefaultA0());  // is equivalent to gain = 4 mV/fC
 347   Float_t adc2mv(AliMUONConstants::DefaultADC2MV()); // 1 ADC channel = 0.61 mV
 348
 349   if ( ! nogain )
 350   {
 351     if  (!gains)
 352     {
 353       AliFatalClass("Cannot make gain decalibration without gain values !");
 354     }
 355     a0 = gains->ValueAsFloat(channel,0);
 356     a1 = gains->ValueAsFloat(channel,1);
 357     thres = gains->ValueAsInt(channel,2);
 358     qual = gains->ValueAsInt(channel,3);
 359   }
 360
 361   Float_t pedestalMean = pedestals.ValueAsFloat(channel,0);
 362   Float_t pedestalSigma = pedestals.ValueAsFloat(channel,1);
 363
 364   AliDebugClass(1,Form("DE %04d MANU %04d CH %02d PEDMEAN %7.2f PEDSIGMA %7.2f",
 365                        pedestals.ID0(),pedestals.ID1(),channel,pedestalMean,pedestalSigma));
 366
 367   if ( qual <= 0 ) return 0;
 368
 369   Float_t chargeThres = a0*thres;
 370
 371   Float_t padc(0); // (adc - ped) value
 372
 373   if ( nogain || charge <= chargeThres || TMath::Abs(a1) < 1E-12 )
 374   {
 375     // linear part only
 376
 377     if ( TMath::Abs(a0) > 1E-12 )
 378     {
 379       padc = charge/a0;
 380     }
 381   }
 382   else
 383   {
 384     // FIXME: when we'll use capacitances and real gains, revise this part
 385     // to take capa properly into account...
 386
 387     AliWarningClass("YOU PROBABLY NEED TO REVISE THIS PART OF CODE !!!");
 388
 389     // linear + parabolic part
 390     Double_t qt = chargeThres - charge;
 391     Double_t delta = a0*a0-4*a1*qt;
 392     if ( delta < 0 )
 393     {
 394       AliErrorClass(Form("delta=%e DE %d Manu %d Channel %d "
 395                     " charge %e a0 %e a1 %e thres %d ped %e pedsig %e",
 396                     delta,pedestals.ID0(),pedestals.ID1(),
 397                     channel, charge, a0, a1, thres, pedestalMean,
 398                     pedestalSigma));
 399     }
 400     else
 401     {
 402       delta = TMath::Sqrt(delta);
 403
 404       padc = ( ( -a0 + delta ) > 0 ? ( -a0 + delta ) : ( -a0 - delta ) );
 405
 406       padc /= 2*a1;
 407
 408       if ( padc < 0 )
 409       {
 410         if ( TMath::Abs(padc) > 1E-3)
 411         {
 412           // this is more than a precision problem : let's signal it !
 413           AliErrorClass(Form("padc=%e DE %d Manu %d Channel %d "
 414                              " charge %e a0 %e a1 %e thres %d ped %e pedsig %e delta %e",
 415                              padc,pedestals.ID0(),pedestals.ID1(),
 416                              channel, charge, a0, a1, thres, pedestalMean,
 417                              pedestalSigma,delta));
 418         }
 419
 420         // ok. consider we're just at thres, let it be zero.
 421         padc = 0;
 422       }
 423
 424       padc += thres;
 425
 426     }
 427   }
 428
 429   padc /= capa*adc2mv;
 430
 431   Int_t adc(0);
 432
 433   Float_t adcNoise = 0.0;
 434
 435   if ( addNoise )
 436   {
 437     if ( noiseOnly )
 438     {
 439       adcNoise = NoiseFunction()->GetRandom()*pedestalSigma;
 440     }
 441     else
 442     {
 443       adcNoise = gRandom->Gaus(0.0,pedestalSigma);
 444     }
 445   }
 446
 447   adc = TMath::Nint(padc + pedestalMean + adcNoise + 0.5);
 448
 449   if ( adc < TMath::Nint(pedestalMean + fgNSigmas*pedestalSigma + 0.5) )
 450   {
 451     // this is an error only in specific cases
 452     if ( !addNoise || (addNoise && noiseOnly) )
 453     {
 454       AliDebugClass(1,Form(" DE %04d Manu %04d Channel %02d "
 455                                                                                                          " a0 %7.2f a1 %7.2f thres %04d ped %7.2f pedsig %7.2f adcNoise %7.2f "
 456                                                                                                          " charge=%7.2f padc=%7.2f adc=%04d ZS=%04d fgNSigmas=%e addNoise %d noiseOnly %d ",
 457                                                                                                          pedestals.ID0(),pedestals.ID1(),channel,
 458                                                                                                          a0, a1, thres, pedestalMean, pedestalSigma, adcNoise,
 459                                                                                                          charge, padc, adc,
 460                                                                                                          TMath::Nint(pedestalMean + fgNSigmas*pedestalSigma + 0.5),
 461                                                                                                          fgNSigmas,addNoise,noiseOnly));
 462     }
 463
 464     adc = 0;
 465   }
 466
 467   // be sure we stick to 12 bits.
 468   if ( adc > kMaxADC )
 469   {
 470     adc = kMaxADC;
 471   }
 472
 473   return adc;
 474 }
 475
 476 //_____________________________________________________________________________
 477 void
 478 AliMUONDigitizerV3::CreateInputDigitStores()
 479 {
 480   /// Create input digit stores
 481   ///
 482
 483   if (fInputDigitStores)
 484   {
 485     AliFatal("Should be called only once !");
 486   }
 487
 488   fInputDigitStores = new TObjArray;
 489
 490   fInputDigitStores->SetOwner(kTRUE);
 491
 492   for ( Int_t iFile = 0; iFile < fDigInput->GetNinputs(); ++iFile )
 493   {
 494     AliLoader* inputLoader = GetLoader(fDigInput->GetInputFolderName(iFile));
 495
 496     inputLoader->LoadSDigits("READ");
 497
 498     TTree* iTreeS = inputLoader->TreeS();
 499     if (!iTreeS)
 500     {
 501       AliFatal(Form("Could not get access to input file #%d",iFile));
 502     }
 503
 504     fInputDigitStores->AddAt(AliMUONVDigitStore::Create(*iTreeS),iFile);
 505   }
 506 }
 507
 508 //_____________________________________________________________________________
 509 void
 510 AliMUONDigitizerV3::Digitize(Option_t*)
 511 {
 512   /// Main method.
 513   /// We first loop over input files, and merge the sdigits we found there.
 514   /// Second, we digitize all the resulting sdigits
 515   /// Then we generate noise-only digits (for tracker only)
 516   /// And we finally generate the trigger outputs.
 517
 518   AliCodeTimerAuto("",0)
 519
 520   if ( fDigInput->GetNinputs() == 0 )
 521   {
 522     AliWarning("No input set. Nothing to do.");
 523     return;
 524   }
 525
 526   if ( !fIsInitialized )
 527   {
 528     AliError("Not initialized. Cannot perform the work. Sorry");
 529     return;
 530   }
 531
 532   Int_t nInputFiles = fDigInput->GetNinputs();
 533
 534   AliLoader* outputLoader = GetLoader(fDigInput->GetOutputFolderName());
 535
 536   outputLoader->MakeDigitsContainer();
 537
 538   TTree* oTreeD = outputLoader->TreeD();
 539
 540   if (!oTreeD)
 541   {
 542     AliFatal("Cannot create output TreeD");
 543   }
 544
 545   // Loop over all the input files, and merge the sdigits found in those
 546   // files.
 547
 548   for ( Int_t iFile = 0; iFile < nInputFiles; ++iFile )
 549   {
 550     AliLoader* inputLoader = GetLoader(fDigInput->GetInputFolderName(iFile));
 551
 552     inputLoader->LoadSDigits("READ");
 553
 554     TTree* iTreeS = inputLoader->TreeS();
 555     if (!iTreeS)
 556     {
 557       AliFatal(Form("Could not get access to input file #%d",iFile));
 558     }
 559
 560     if (!fInputDigitStores)
 561     {
 562       CreateInputDigitStores();
 563     }
 564
 565     AliMUONVDigitStore* dstore = static_cast<AliMUONVDigitStore*>(fInputDigitStores->At(iFile));
 566
 567     dstore->Connect(*iTreeS);
 568
 569     iTreeS->GetEvent(0);
 570
 571     MergeWithSDigits(fDigitStore,*dstore,fDigInput->GetMask(iFile));
 572
 573     inputLoader->UnloadSDigits();
 574
 575     dstore->Clear();
 576   }
 577
 578
 579   // At this point, we do have digit arrays (one per chamber) which contains
 580   // the merging of all the sdigits of the input file(s).
 581   // We now massage them to apply the detector response, i.e. this
 582   // is here that we do the "digitization" work.
 583
 584   if (!fOutputDigitStore)
 585   {
 586     fOutputDigitStore = fDigitStore->Create();
 587   }
 588
 589   if ( fGenerateNoisyDigits>=2 )
 590   {
 591     // Generate noise-only digits for trigger.
 592     GenerateNoisyDigitsForTrigger(*fDigitStore);
 593   }
 594   ApplyResponse(*fDigitStore,*fOutputDigitStore);
 595
 596   if ( fGenerateNoisyDigits )
 597   {
 598     // Generate noise-only digits for tracker.
 599     GenerateNoisyDigits(*fOutputDigitStore);
 600   }
 601
 602   // We generate the global and local trigger decisions.
 603   fTriggerProcessor->Digits2Trigger(*fOutputDigitStore,*fTriggerStore);
 604
 605   // Prepare output tree
 606   Bool_t okD = fOutputDigitStore->Connect(*oTreeD,kFALSE);
 607   Bool_t okT = fTriggerStore->Connect(*oTreeD,kFALSE);
 608   if (!okD || !okT)
 609   {
 610     AliError(Form("Could not make branch : Digit %d Trigger %d",okD,okT));
 611     return;
 612   }
 613
 614   // Fill the output treeD
 615   oTreeD->Fill();
 616
 617   // Write to the output tree(D).
 618   // Please note that as GlobalTrigger, LocalTrigger and Digits are in the same
 619   // tree (=TreeD) in different branches, this WriteDigits in fact writes all of
 620   // the 3 branches.
 621   outputLoader->WriteDigits("OVERWRITE");
 622
 623   outputLoader->UnloadDigits();
 624
 625   // Finally, we clean up after ourselves.
 626   fTriggerStore->Clear();
 627   fDigitStore->Clear();
 628   fOutputDigitStore->Clear();
 629 }
 630
 631
 632 //_____________________________________________________________________________
 633 void
 634 AliMUONDigitizerV3::GenerateNoisyDigits(AliMUONVDigitStore& digitStore)
 635 {
 636   /// According to a given probability, generate digits that
 637   /// have a signal above the noise cut (ped+n*sigma_ped), i.e. digits
 638   /// that are "only noise".
 639
 640   AliCodeTimerAuto("",0)
 641
 642   for ( Int_t i = 0; i < AliMUONConstants::NTrackingCh(); ++i )
 643   {
 644     AliMpDEIterator it;
 645
 646     it.First(i);
 647
 648     while ( !it.IsDone() )
 649     {
 650       for ( Int_t cathode = 0; cathode < 2; ++cathode )
 651       {
 652         GenerateNoisyDigitsForOneCathode(digitStore,it.CurrentDEId(),cathode);
 653       }
 654       it.Next();
 655     }
 656   }
 657 }
 658
 659 //_____________________________________________________________________________
 660 void
 661 AliMUONDigitizerV3::GenerateNoisyDigitsForOneCathode(AliMUONVDigitStore& digitStore,
 662                                                      Int_t detElemId, Int_t cathode)
 663 {
 664   /// Generate noise-only digits for one cathode of one detection element.
 665   /// Called by GenerateNoisyDigits()
 666
 667   const AliMpVSegmentation* seg
 668     = AliMpSegmentation::Instance()->GetMpSegmentation(detElemId,AliMp::GetCathodType(cathode));
 669   Int_t nofPads = seg->NofPads();
 670
 671   Int_t maxIx = seg->MaxPadIndexX();
 672   Int_t maxIy = seg->MaxPadIndexY();
 673
 674   static const Double_t kProbToBeOutsideNsigmas = TMath::Erfc(fgNSigmas/TMath::Sqrt(2.0)) / 2. ;
 675
 676   Int_t nofNoisyPads = TMath::Nint(kProbToBeOutsideNsigmas*nofPads);
 677   if ( !nofNoisyPads ) return;
 678
 679   nofNoisyPads =
 680     TMath::Nint(gRandom->Gaus(nofNoisyPads,
 681                               nofNoisyPads/TMath::Sqrt(nofNoisyPads)));
 682
 683   AliDebug(3,Form("DE %d cath %d nofNoisyPads %d",detElemId,cathode,nofNoisyPads));
 684
 685   for ( Int_t i = 0; i < nofNoisyPads; ++i )
 686   {
 687     Int_t ix(-1);
 688     Int_t iy(-1);
 689     AliMpPad pad;
 690
 691     do {
 692       ix = gRandom->Integer(maxIx+1);
 693       iy = gRandom->Integer(maxIy+1);
 694       pad = seg->PadByIndices(ix,iy,kFALSE);
 695     } while ( !pad.IsValid() );
 696
 697     Int_t manuId = pad.GetManuId();
 698     Int_t manuChannel = pad.GetManuChannel();
 699
 700     AliMUONVCalibParam* pedestals = fCalibrationData->Pedestals(detElemId,manuId);
 701
 702     if (!pedestals)
 703     {
 704       // no pedestal available for this channel, simply give up
 705       continue;
 706     }
 707
 708     AliMUONVDigit* d = digitStore.CreateDigit(detElemId,manuId,manuChannel,cathode);
 709
 710     d->SetPadXY(ix,iy);
 711
 712     d->SetCharge(0.0); // charge is zero, the ApplyResponseToTrackerDigit will add the noise
 713     d->NoiseOnly(kTRUE);
 714     ApplyResponseToTrackerDigit(*d,kTRUE);
 715     if ( d->ADC() > 0 )
 716     {
 717       Bool_t ok = digitStore.Add(*d,AliMUONVDigitStore::kDeny);
 718       // this can happen (that we randomly chose a digit that is
 719       // already there). We simply ignore this, but log the occurence
 720       // to cross-check that it's not too frequent.
 721       if (!ok)
 722       {
 723         fLogger->Log("Collision while adding noiseOnly digit");
 724       }
 725       else
 726       {
 727         fLogger->Log("Added noiseOnly digit");
 728       }
 729     }
 730     delete d;
 731   }
 732 }
 733
 734
 735 //_____________________________________________________________________________
 736 void
 737 AliMUONDigitizerV3::GenerateNoisyDigitsForTrigger(AliMUONVDigitStore& digitStore)
 738 {
 739   /// Generate noise-only digits for one cathode of one detection element.
 740   /// Called by GenerateNoisyDigits()
 741
 742   if ( !fNoiseFunctionTrig )
 743   {
 744     fNoiseFunctionTrig = new TF1("AliMUONDigitizerV3::fNoiseFunctionTrig","landau",
 745                                  50.,270.);
 746
 747     fNoiseFunctionTrig->SetParameters(3.91070e+02, 9.85026, 9.35881e-02);
 748   }
 749
 750   AliMpPad pad[2];
 751   AliMUONVDigit *d[2]={0x0};
 752
 753   for ( Int_t chamberId = AliMUONConstants::NTrackingCh(); chamberId < AliMUONConstants::NCh(); ++chamberId )
 754   {
 755
 756     Int_t nofNoisyPads = 50;
 757
 758     Float_t r=-1, fi = 0., gx, gy, x, y, z, xg01, yg01, zg, xg02, yg02;
 759     AliMpDEIterator it;
 760
 761     AliDebug(3,Form("Chamber %d nofNoisyPads %d",chamberId,nofNoisyPads));
 762
 763     for ( Int_t i = 0; i < nofNoisyPads; ++i )
 764     {
 765       //printf("Generating noise %i\n",i);
 766         Int_t ix(-1);
 767         Int_t iy(-1);
 768         Bool_t isOk = kFALSE;
 769         Int_t detElemId = -1;
 770         do {
 771           //r = gRandom->Landau(9.85026, 9.35881e-02);
 772             r = fNoiseFunctionTrig->GetRandom();
 773             fi = 2. * TMath::Pi() * gRandom->Rndm();
 774             //printf("r = %f\tfi = %f\n", r, fi);
 775             gx = r * TMath::Cos(fi);
 776             gy = r * TMath::Sin(fi);
 777
 778             for ( it.First(chamberId); ! it.IsDone(); it.Next() ){
 779                 Int_t currDetElemId = it.CurrentDEId();
 780                 const AliMpVSegmentation* seg
 781                     = AliMpSegmentation::Instance()->GetMpSegmentation(currDetElemId,AliMp::GetCathodType(0));
 782                 if (!seg) continue;
 783                 Float_t deltax = seg->GetDimensionX();
 784                 Float_t deltay = seg->GetDimensionY();
 785                 GetTransformer()->Local2Global(currDetElemId, -deltax, -deltay, 0, xg01, yg01, zg);
 786                 GetTransformer()->Local2Global(currDetElemId,  deltax,  deltay, 0, xg02, yg02, zg);
 787                 Float_t xg1 = xg01, xg2 = xg02, yg1 = yg01, yg2 = yg02;
 788                 if(xg01>xg02){
 789                     xg1 = xg02;
 790                     xg2 = xg01;
 791                 }
 792                 if(yg01>yg02){
 793                     yg1 = yg02;
 794                     yg2 = yg01;
 795                 }
 796                 if(gx>=xg1 && gx<=xg2 && gy>=yg1 && gy<=yg2){
 797                     detElemId = currDetElemId;
 798                     GetTransformer()->Global2Local(detElemId, gx, gy, 0, x, y, z);
 799                     pad[0] = seg->PadByPosition(x,y,kFALSE);
 800                     if(!pad[0].IsValid()) continue;
 801                     isOk = kTRUE;
 802                     break;
 803                 }
 804             } // loop on slats
 805         } while ( !isOk );
 806
 807         const AliMpVSegmentation* seg1
 808             = AliMpSegmentation::Instance()->GetMpSegmentation(detElemId,AliMp::GetCathodType(1));
 809         pad[1] = seg1->PadByPosition(x,y,kFALSE);
 810
 811         for ( Int_t cathode = 0; cathode < 2; ++cathode ){
 812           Int_t manuId = pad[cathode].GetLocalBoardId(0);
 813           Int_t manuChannel = pad[cathode].GetLocalBoardChannel(0);
 814           d[cathode] = digitStore.CreateDigit(detElemId,manuId,manuChannel,cathode);
 815           ix = pad[cathode].GetIx();
 816           iy = pad[cathode].GetIy();
 817           d[cathode]->SetPadXY(ix,iy);
 818           //d[cathode].SetSignal(1);
 819           //d[cathode].SetPhysicsSignal(0);
 820           d[cathode]->SetCharge(1);
 821           d[cathode]->NoiseOnly(kTRUE);
 822           AliDebug(3,Form("Adding a pure noise digit :"));
 823
 824           Bool_t ok = digitStore.Add(*d[cathode],AliMUONVDigitStore::kDeny);
 825           if (!ok)
 826           {
 827               fLogger->Log("Collision while adding TriggerNoise digit");
 828           }
 829           else
 830           {
 831               fLogger->Log("Added triggerNoise digit");
 832           }
 833         } //loop on cathodes
 834     } // loop on noisy pads
 835   } // loop on chambers
 836 }
 837
 838
 839 //_____________________________________________________________________________
 840 AliLoader*
 841 AliMUONDigitizerV3::GetLoader(const TString& folderName)
 842 {
 843   /// Get a MUON loader
 844
 845   AliDebug(2,Form("Getting access to folder %s",folderName.Data()));
 846   AliLoader* loader = AliRunLoader::GetDetectorLoader("MUON",folderName.Data());
 847   if (!loader)
 848   {
 849     AliError(Form("Could not get MuonLoader from folder %s",folderName.Data()));
 850     return 0x0;
 851   }
 852   return loader;
 853 }
 854
 855 //_____________________________________________________________________________
 856 Bool_t
 857 AliMUONDigitizerV3::Init()
 858 {
 859   /// Initialization of the digitization :
 860   /// a) create the calibrationData, according to run number
 861   /// b) create the trigger processing task
 862
 863   AliDebug(2,"");
 864
 865   if ( fIsInitialized )
 866   {
 867     AliError("Object already initialized.");
 868     return kFALSE;
 869   }
 870
 871   if (!fDigInput)
 872   {
 873     AliError("fDigInput is null !");
 874     return kFALSE;
 875   }
 876
 877   // Load mapping
 878   if ( ! AliMpCDB::LoadDDLStore() ) {
 879     AliFatal("Could not access mapping from OCDB !");
 880   }
 881
 882   if (!fCalibrationData)
 883       AliFatal("Calibration data object not defined");
 884
 885   if ( !fCalibrationData->Pedestals() )
 886   {
 887     AliFatal("Could not access pedestals from OCDB !");
 888   }
 889   if ( !fCalibrationData->Gains() )
 890   {
 891     AliFatal("Could not access gains from OCDB !");
 892   }
 893
 894
 895    AliInfo("Using trigger configuration from CDB");
 896
 897   fTriggerProcessor = new AliMUONTriggerElectronics(fCalibrationData);
 898
 899   AliDebug(1, Form("Will %s generate noise-only digits for tracker",
 900                      (fGenerateNoisyDigits ? "":"NOT")));
 901
 902   fTriggerUtilities = new AliMUONTriggerUtilities(fCalibrationData);
 903
 904   if ( muon()->GetTriggerEffCells() ) {
 905     // Apply trigger efficiency
 906     AliDebug(1, "Will apply trigger efficiency");
 907     fTriggerEfficiency = new AliMUONTriggerChamberEfficiency(fCalibrationData->TriggerEfficiency());
 908   }
 909
 910   fIsInitialized = kTRUE;
 911   return kTRUE;
 912 }
 913
 914
 915 //_____________________________________________________________________________
 916 Int_t AliMUONDigitizerV3::GetArrayIndex(Int_t cathode, Int_t trigCh, Int_t localCircuit)
 917 {
 918   /// Get index of array with trigger status map or efficiency
 919   return
 920     AliMUONConstants::NTriggerCircuit() * AliMUONConstants::NTriggerCh() * cathode +
 921     AliMUONConstants::NTriggerCircuit() * trigCh + localCircuit-1;
 922 }
 923
 924
 925 //_____________________________________________________________________________
 926 void
 927 AliMUONDigitizerV3::MergeWithSDigits(AliMUONVDigitStore*& outputStore,
 928                                      const AliMUONVDigitStore& input,
 929                                      Int_t mask)
 930 {
 931   /// Merge the sdigits in inputData with the digits already present in outputData
 932
 933   if ( !outputStore ) outputStore = input.Create();
 934
 935   TIter next(input.CreateIterator());
 936   AliMUONVDigit* sdigit;
 937
 938   while ( ( sdigit = static_cast<AliMUONVDigit*>(next()) ) )
 939   {
 940     // Update the track references using the mask.
 941     // FIXME: this is dirty, for backward compatibility only.
 942     // Should re-design all this way of keeping track of MC information...
 943     if ( mask ) sdigit->PatchTracks(mask);
 944     // Then add or update the digit to the output.
 945     AliMUONVDigit* added = outputStore->Add(*sdigit,AliMUONVDigitStore::kMerge);
 946     if (!added)
 947     {
 948       AliError("Could not add digit in merge mode");
 949     }
 950   }
 951 }
 952
 953 //_____________________________________________________________________________
 954 TF1*
 955 AliMUONDigitizerV3::NoiseFunction()
 956 {
 957   /// Return noise function
 958   static TF1* f = 0x0;
 959   if (!f)
 960   {
 961     f = new TF1("AliMUONDigitizerV3::NoiseFunction","gaus",fgNSigmas,fgNSigmas*10);
 962     f->SetParameters(1,0,1);
 963   }
 964   return f;
 965 }
 966
 967 //_____________________________________________________________________________
 968 void AliMUONDigitizerV3::SetCalibrationData(AliMUONCalibrationData* calibrationData,
 969                                             AliMUONRecoParam* recoParam)
 970 {
 971   fCalibrationData = calibrationData;
 972   fRecoParam = recoParam;
 973   if (!fRecoParam)
 974   {
 975     AliError("Cannot work (e.g. decalibrate) without recoparams !");
 976   }
 977 }
 978