]> git.uio.no Git - u/mrichter/AliRoot.git/blob - MUON/AliMUONDigitizerV3.cxx
Removing, as not used anymore (Laurent)
[u/mrichter/AliRoot.git] / 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 "AliCDBManager.h"
22 #include "AliLog.h"
23 #include "AliMUON.h"
24 #include "AliMUONCalibrationData.h"
25 #include "AliMUONConstants.h"
26 #include "AliMUONDigit.h"
27 #include "AliMUONLogger.h"
28 #include "AliMUONTriggerEfficiencyCells.h"
29 #include "AliMUONTriggerElectronics.h"
30 #include "AliMUONTriggerStoreV1.h"
31 #include "AliMUONVCalibParam.h"
32 #include "AliMUONVDigitStore.h"
33 #include "AliMpCathodType.h"
34 #include "AliMpConstants.h"
35 #include "AliMpDEIterator.h"
36 #include "AliMpDEManager.h"
37 #include "AliMpDEManager.h"
38 #include "AliMpIntPair.h"
39 #include "AliMpPad.h"
40 #include "AliMpSegmentation.h"
41 #include "AliMpStationType.h"
42 #include "AliMpVSegmentation.h"
43 #include "AliRun.h"
44 #include "AliRunDigitizer.h"
45 #include "AliRunLoader.h"
46 #include <Riostream.h>
47 #include <TF1.h>
48 #include <TFile.h>
49 #include <TMath.h>
50 #include <TRandom.h>
51 #include <TString.h>
52 #include <TSystem.h>
53
54 ///
55 /// \class AliMUONDigitizerV3
56 /// The digitizer is performing the transformation to go from SDigits (digits
57 /// w/o any electronic noise) to Digits (w/ electronic noise, and decalibration)
58 /// 
59 /// The decalibration is performed by doing the reverse operation of the
60 /// calibration, that is we do (Signal+pedestal)/gain -> ADC
61 ///
62 /// Note also that the digitizer takes care of merging sdigits that belongs
63 /// to the same pad, either because we're merging several input sdigit files
64 /// or with a single file because the sdigitizer does not merge sdigits itself
65 /// (for performance reason mainly, and because anyway we know we have to do it
66 /// here, at the digitization level).
67 ///
68 /// \author Laurent Aphecetche
69
70 namespace
71 {
72   AliMUON* muon()
73   {
74     return static_cast<AliMUON*>(gAlice->GetModule("MUON"));
75   }
76 }
77
78 const Double_t AliMUONDigitizerV3::fgkNSigmas=3;
79
80 /// \cond CLASSIMP
81 ClassImp(AliMUONDigitizerV3)
82 /// \endcond
83
84 //_____________________________________________________________________________
85 AliMUONDigitizerV3::AliMUONDigitizerV3(AliRunDigitizer* manager, 
86                                        Bool_t generateNoisyDigits)
87 : AliDigitizer(manager),
88 fIsInitialized(kFALSE),
89 fCalibrationData(0x0),
90 fTriggerProcessor(0x0),
91 fTriggerEfficiency(0x0),
92 fGenerateNoisyDigitsTimer(),
93 fExecTimer(),
94 fNoiseFunction(0x0),
95   fGenerateNoisyDigits(generateNoisyDigits),
96   fLogger(new AliMUONLogger(1000)),
97 fTriggerStore(new AliMUONTriggerStoreV1),
98 fDigitStore(0x0),
99 fOutputDigitStore(0x0)
100 {
101   /// Ctor.
102
103   AliDebug(1,Form("AliRunDigitizer=%p",fManager));
104   fGenerateNoisyDigitsTimer.Start(kTRUE); fGenerateNoisyDigitsTimer.Stop();
105   fExecTimer.Start(kTRUE); fExecTimer.Stop();
106 }
107
108 //_____________________________________________________________________________
109 AliMUONDigitizerV3::~AliMUONDigitizerV3()
110 {
111   /// Dtor. Note we're the owner of some pointers.
112
113   AliDebug(1,"dtor");
114
115   delete fCalibrationData;
116   delete fTriggerProcessor;
117   delete fNoiseFunction;
118   delete fTriggerStore;
119   delete fDigitStore;
120   delete fOutputDigitStore;
121   
122   if ( fGenerateNoisyDigits )
123   {
124     AliDebug(1, Form("Execution time for GenerateNoisyDigits() : R:%.2fs C:%.2fs",
125                  fGenerateNoisyDigitsTimer.RealTime(),
126                  fGenerateNoisyDigitsTimer.CpuTime()));
127   }
128   AliDebug(1, Form("Execution time for Exec() : R:%.2fs C:%.2fs",
129                fExecTimer.RealTime(),fExecTimer.CpuTime()));
130  
131   AliInfo("Summary of messages");
132   fLogger->Print();
133   
134   delete fLogger;
135 }
136
137 //_____________________________________________________________________________
138 void 
139 AliMUONDigitizerV3::ApplyResponseToTrackerDigit(AliMUONVDigit& digit, Bool_t addNoise)
140 {
141   /// For tracking digits, starting from an ideal digit's charge, we :
142   ///
143   /// - add some noise (thus leading to a realistic charge), if requested to do so
144   /// - divide by a gain (thus decalibrating the digit)
145   /// - add a pedestal (thus decalibrating the digit)
146   /// - sets the signal to zero if below 3*sigma of the noise
147
148   static const Int_t kMaxADC = (1<<12)-1; // We code the charge on a 12 bits ADC.
149   
150   Float_t signal = digit.Charge();
151   
152   if ( !addNoise )
153   {
154     digit.SetADC(TMath::Min(kMaxADC,TMath::Nint(signal)));
155     return;
156   }
157   
158   Int_t detElemId = digit.DetElemId();
159   
160   Int_t manuId = digit.ManuId();
161   Int_t manuChannel = digit.ManuChannel();
162   
163   AliMUONVCalibParam* pedestal = fCalibrationData->Pedestals(detElemId,manuId);
164   if (!pedestal)
165   {
166     fLogger->Log(Form("%s:%d:Could not get pedestal for DE=%4d manuId=%4d. Disabling.",
167                       __FILE__,__LINE__,
168                       detElemId,manuId));
169     digit.SetCharge(0);
170     digit.SetADC(0);
171     return;    
172   }
173   Float_t pedestalMean = pedestal->ValueAsFloat(manuChannel,0);
174   Float_t pedestalSigma = pedestal->ValueAsFloat(manuChannel,1);
175   
176   AliMUONVCalibParam* gain = fCalibrationData->Gains(detElemId,manuId);
177   if (!gain)
178   {
179     fLogger->Log(Form("%s:%d:Could not get gain for DE=%4d manuId=%4d. Disabling.",
180                       __FILE__,__LINE__,
181                       detElemId,manuId));
182     digit.SetCharge(0);
183     digit.SetADC(0);
184     return;        
185   }    
186   Float_t gainMean = gain->ValueAsFloat(manuChannel,0);
187   
188   Float_t adcNoise = gRandom->Gaus(0.0,pedestalSigma);
189   
190   Int_t adc;
191   
192   if ( gainMean < 1E-6 )
193   {
194     AliError(Form("Got a too small gain %e for DE=%d manuId=%d manuChannel=%d. "
195                   "Setting signal to 0.",
196                   gainMean,detElemId,manuId,manuChannel));
197     adc = 0;
198   }
199   else
200   {
201     adc = TMath::Nint( signal / gainMean + pedestalMean + adcNoise);///
202     
203     if ( adc <= pedestalMean + fgkNSigmas*pedestalSigma ) 
204     {
205       adc = 0;
206     }
207   }
208   
209   // be sure we stick to 12 bits.
210   if ( adc > kMaxADC )
211   {
212     adc = kMaxADC;
213   }
214   
215   digit.SetCharge(adc);
216   digit.SetADC(adc);
217 }
218
219 //_____________________________________________________________________________
220 void 
221 AliMUONDigitizerV3::ApplyResponseToTriggerDigit(const AliMUONVDigitStore& digitStore,
222                                                 AliMUONVDigit& digit)
223 {
224   /// \todo add comment
225
226   if ( !fTriggerEfficiency ) return;
227
228   if (digit.IsEfficiencyApplied()) return;
229
230   AliMUONVDigit* correspondingDigit = FindCorrespondingDigit(digitStore,digit);
231
232   if (!correspondingDigit) return; //reject bad correspondences
233
234   Int_t detElemId = digit.DetElemId();
235
236   AliMpSegmentation* segmentation = AliMpSegmentation::Instance();
237   const AliMpVSegmentation* segment[2] = 
238   {
239     segmentation->GetMpSegmentation(detElemId,AliMp::GetCathodType(digit.Cathode())), 
240     segmentation->GetMpSegmentation(detElemId,AliMp::GetCathodType(correspondingDigit->Cathode()))
241   };
242
243   AliMpPad pad[2] = 
244   {
245     segment[0]->PadByIndices(AliMpIntPair(digit.PadX(),digit.PadY()),kTRUE), 
246     segment[1]->PadByIndices(AliMpIntPair(correspondingDigit->PadX(),correspondingDigit->PadY()),kTRUE)
247   };
248
249   Int_t p0(1);
250   if (digit.Cathode()==0) p0=0;
251
252   AliMpIntPair location = pad[p0].GetLocation(0);
253   Int_t nboard = location.GetFirst();
254
255   Bool_t isTrig[2];
256
257   fTriggerEfficiency->IsTriggered(detElemId, nboard-1, 
258                                   isTrig[0], isTrig[1]);
259   digit.EfficiencyApplied(kTRUE);
260   correspondingDigit->EfficiencyApplied(kTRUE);
261
262   if (!isTrig[digit.Cathode()])
263   {
264           digit.SetCharge(0);
265   }
266   
267   if ( &digit != correspondingDigit )
268   {
269           if (!isTrig[correspondingDigit->Cathode()])
270     {
271       correspondingDigit->SetCharge(0);
272           }
273   }
274 }
275
276 //_____________________________________________________________________________
277 void
278 AliMUONDigitizerV3::ApplyResponse(const AliMUONVDigitStore& store,
279                                   AliMUONVDigitStore& filteredStore)
280 {
281   /// Loop over all chamber digits, and apply the response to them
282   /// Note that this method may remove digits.
283
284   filteredStore.Clear();
285   
286   const Bool_t kAddNoise = kTRUE;
287   
288   TIter next(store.CreateIterator());
289   AliMUONVDigit* digit;
290   
291   while ( ( digit = static_cast<AliMUONVDigit*>(next()) ) )
292   {
293     AliMp::StationType stationType = AliMpDEManager::GetStationType(digit->DetElemId());
294     
295     if ( stationType != AliMp::kStationTrigger )
296     {
297       ApplyResponseToTrackerDigit(*digit,kAddNoise);
298     }
299     else
300     {
301       ApplyResponseToTriggerDigit(store,*digit);
302     }
303     if ( digit->Charge() > 0  )
304     {
305       filteredStore.Add(*digit,AliMUONVDigitStore::kIgnore);
306     }
307   }
308 }    
309
310 //_____________________________________________________________________________
311 void
312 AliMUONDigitizerV3::Exec(Option_t*)
313 {
314   /// Main method.
315   /// We first loop over input files, and merge the sdigits we found there.
316   /// Second, we digitize all the resulting sdigits
317   /// Then we generate noise-only digits (for tracker only)
318   /// And we finally generate the trigger outputs.
319     
320   AliDebug(1, "Running digitizer.");
321   
322   if ( fManager->GetNinputs() == 0 )
323   {
324     AliWarning("No input set. Nothing to do.");
325     return;
326   }
327   
328   if ( !fIsInitialized )
329   {
330     AliError("Not initialized. Cannot perform the work. Sorry");
331     return;
332   }
333   
334   fExecTimer.Start(kFALSE);
335
336   Int_t nInputFiles = fManager->GetNinputs();
337   
338   AliLoader* outputLoader = GetLoader(fManager->GetOutputFolderName());
339   
340   outputLoader->MakeDigitsContainer();
341   
342   TTree* oTreeD = outputLoader->TreeD();
343   
344   if (!oTreeD) 
345   {
346     AliFatal("Cannot create output TreeD");
347   }
348
349   // Loop over all the input files, and merge the sdigits found in those
350   // files.
351   
352   for ( Int_t iFile = 0; iFile < nInputFiles; ++iFile )
353   {    
354     AliLoader* inputLoader = GetLoader(fManager->GetInputFolderName(iFile));
355
356     inputLoader->LoadSDigits("READ");
357
358     TTree* iTreeS = inputLoader->TreeS();
359     if (!iTreeS)
360     {
361       AliFatal(Form("Could not get access to input file #%d",iFile));
362     }
363     
364     AliMUONVDigitStore* inputStore = AliMUONVDigitStore::Create(*iTreeS);
365     inputStore->Connect(*iTreeS);
366     
367     iTreeS->GetEvent(0);
368     
369     MergeWithSDigits(fDigitStore,*inputStore,fManager->GetMask(iFile));
370
371     inputLoader->UnloadSDigits();
372     
373     delete inputStore;
374   }
375   
376   // At this point, we do have digit arrays (one per chamber) which contains 
377   // the merging of all the sdigits of the input file(s).
378   // We now massage them to apply the detector response, i.e. this
379   // is here that we do the "digitization" work.
380   
381   if (!fOutputDigitStore)
382   {
383     fOutputDigitStore = fDigitStore->Create();
384   }
385   
386   ApplyResponse(*fDigitStore,*fOutputDigitStore);
387   
388   if ( fGenerateNoisyDigits )
389   {
390     // Generate noise-only digits for tracker.
391     GenerateNoisyDigits(*fOutputDigitStore);
392   }
393   
394   // We generate the global and local trigger decisions.
395   fTriggerProcessor->Digits2Trigger(*fOutputDigitStore,*fTriggerStore);
396
397   // Prepare output tree
398   Bool_t okD = fOutputDigitStore->Connect(*oTreeD,kFALSE);
399   Bool_t okT = fTriggerStore->Connect(*oTreeD,kFALSE);
400   if (!okD || !okT)
401   {
402     AliError(Form("Could not make branch : Digit %d Trigger %d",okD,okT));
403     return;
404   }
405   
406   // Fill the output treeD
407   oTreeD->Fill();
408   
409   // Write to the output tree(D).
410   // Please note that as GlobalTrigger, LocalTrigger and Digits are in the same
411   // tree (=TreeD) in different branches, this WriteDigits in fact writes all of 
412   // the 3 branches.
413   outputLoader->WriteDigits("OVERWRITE");  
414   
415   outputLoader->UnloadDigits();
416   
417   // Finally, we clean up after ourselves.
418   fTriggerStore->Clear();
419   fDigitStore->Clear();
420   fOutputDigitStore->Clear();
421   fExecTimer.Stop();
422 }
423
424 //_____________________________________________________________________________
425 AliMUONVDigit* 
426 AliMUONDigitizerV3::FindCorrespondingDigit(const AliMUONVDigitStore& digitStore,
427                                            AliMUONVDigit& digit) const
428 {                                                
429   /// Find, if it exists, the digit corresponding to digit.Hit(), in the 
430   /// other cathode
431
432   TIter next(digitStore.CreateIterator());
433   AliMUONVDigit* d;
434   
435   while ( ( d = static_cast<AliMUONVDigit*>(next()) ) )
436   {
437     if ( d->DetElemId() == digit.DetElemId() &&
438          d->Hit() == digit.Hit() &&
439          d->Cathode() != digit.Cathode() )
440     {
441       return d;
442     }      
443   }    
444   return 0x0;
445 }
446
447
448 //_____________________________________________________________________________
449 void
450 AliMUONDigitizerV3::GenerateNoisyDigits(AliMUONVDigitStore& digitStore)
451 {
452   /// According to a given probability, generate digits that
453   /// have a signal above the noise cut (ped+n*sigma_ped), i.e. digits
454   /// that are "only noise".
455   
456   if ( !fNoiseFunction )
457   {
458     fNoiseFunction = new TF1("AliMUONDigitizerV3::fNoiseFunction","gaus",
459                              fgkNSigmas,fgkNSigmas*10);
460     
461     fNoiseFunction->SetParameters(1,0,1);
462   }
463   
464   fGenerateNoisyDigitsTimer.Start(kFALSE);
465   
466   for ( Int_t i = 0; i < AliMUONConstants::NTrackingCh(); ++i )
467   {
468     AliMpDEIterator it;
469   
470     it.First(i);
471   
472     while ( !it.IsDone() )
473     {
474       for ( Int_t cathode = 0; cathode < 2; ++cathode )
475       {
476         GenerateNoisyDigitsForOneCathode(digitStore,it.CurrentDEId(),cathode);
477       }
478       it.Next();
479     }
480   }
481   
482   fGenerateNoisyDigitsTimer.Stop();
483 }
484  
485 //_____________________________________________________________________________
486 void
487 AliMUONDigitizerV3::GenerateNoisyDigitsForOneCathode(AliMUONVDigitStore& digitStore,
488                                                      Int_t detElemId, Int_t cathode)
489 {
490   /// Generate noise-only digits for one cathode of one detection element.
491   /// Called by GenerateNoisyDigits()
492   
493   const AliMpVSegmentation* seg 
494     = AliMpSegmentation::Instance()->GetMpSegmentation(detElemId,AliMp::GetCathodType(cathode));
495   Int_t nofPads = seg->NofPads();
496   
497   Int_t maxIx = seg->MaxPadIndexX();
498   Int_t maxIy = seg->MaxPadIndexY();
499   
500   static const Double_t kProbToBeOutsideNsigmas = TMath::Erfc(fgkNSigmas/TMath::Sqrt(2.0)) / 2. ;
501   
502   Int_t nofNoisyPads = TMath::Nint(kProbToBeOutsideNsigmas*nofPads);
503   if ( !nofNoisyPads ) return;
504   
505   nofNoisyPads = 
506     TMath::Nint(gRandom->Gaus(nofNoisyPads,
507                               nofNoisyPads/TMath::Sqrt(nofNoisyPads)));
508   
509   AliDebug(3,Form("DE %d cath %d nofNoisyPads %d",detElemId,cathode,nofNoisyPads));
510   
511   for ( Int_t i = 0; i < nofNoisyPads; ++i ) 
512   {
513     Int_t ix(-1);
514     Int_t iy(-1);
515     AliMpPad pad;
516     
517     do {
518       ix = gRandom->Integer(maxIx+1);
519       iy = gRandom->Integer(maxIy+1);
520       pad = seg->PadByIndices(AliMpIntPair(ix,iy),kFALSE);
521     } while ( !pad.IsValid() );
522
523     Int_t manuId = pad.GetLocation().GetFirst();
524     Int_t manuChannel = pad.GetLocation().GetSecond();    
525
526     AliMUONVCalibParam* pedestals = fCalibrationData->Pedestals(detElemId,manuId);
527     
528     if (!pedestals) 
529     {
530       // no pedestal available for this channel, simply give up
531       continue;
532     }
533     
534     AliMUONVDigit* d = digitStore.CreateDigit(detElemId,manuId,manuChannel,cathode);
535     
536     d->SetPadXY(ix,iy);
537     
538     Float_t pedestalMean = pedestals->ValueAsFloat(manuChannel,0);
539     Float_t pedestalSigma = pedestals->ValueAsFloat(manuChannel,1);
540     
541     Double_t ped = fNoiseFunction->GetRandom()*pedestalSigma;
542
543     d->SetCharge(TMath::Nint(ped+pedestalMean+0.5));
544     d->NoiseOnly(kTRUE);
545     ApplyResponseToTrackerDigit(*d,kFALSE);
546     if ( d->Charge() > 0 )
547     {
548       Bool_t ok = digitStore.Add(*d,AliMUONVDigitStore::kDeny);
549       // this can happen (that we randomly chose a digit that is
550       // already there). We simply ignore this, but log the occurence
551       // to cross-check that it's not too frequent.
552       if (!ok)
553       {
554         fLogger->Log("Collision while adding noiseOnly digit");
555       }
556       else
557       {
558         fLogger->Log("Added noiseOnly digit");
559       }
560     }
561     else
562     {
563       AliError("Pure noise below threshold. This should not happen. Not adding "
564                "this digit.");
565     }
566     delete d;
567   }
568 }
569
570 //_____________________________________________________________________________
571 AliLoader*
572 AliMUONDigitizerV3::GetLoader(const TString& folderName)
573 {
574   /// Get a MUON loader
575
576   AliDebug(2,Form("Getting access to folder %s",folderName.Data()));
577   AliLoader* loader = AliRunLoader::GetDetectorLoader("MUON",folderName.Data());
578   if (!loader)
579   {
580     AliError(Form("Could not get MuonLoader from folder %s",folderName.Data()));
581     return 0x0;
582   }
583   return loader;
584 }
585
586 //_____________________________________________________________________________
587 Bool_t
588 AliMUONDigitizerV3::Init()
589 {
590   /// Initialization of the TTask :
591   /// a) create the calibrationData, according to run number
592   /// b) create the trigger processing task
593
594   AliDebug(2,"");
595   
596   if ( fIsInitialized )
597   {
598     AliError("Object already initialized.");
599     return kFALSE;
600   }
601   
602   if (!fManager)
603   {
604     AliError("fManager is null !");
605     return kFALSE;
606   }
607   
608   Int_t runnumber = AliCDBManager::Instance()->GetRun();
609   
610   fCalibrationData = new AliMUONCalibrationData(runnumber);
611   if ( !fCalibrationData->Pedestals() )
612   {
613     AliFatal("Could not access pedestals from OCDB !");
614   }
615   if ( !fCalibrationData->Gains() )
616   {
617     AliFatal("Could not access gains from OCDB !");
618   }
619   fTriggerProcessor = new AliMUONTriggerElectronics(fCalibrationData);
620   
621   if ( muon()->GetTriggerEffCells() )
622   {
623     fTriggerEfficiency = fCalibrationData->TriggerEfficiency();
624     if ( fTriggerEfficiency )
625     {
626       AliDebug(1, "Will apply trigger efficiency");
627     }
628     else
629     {
630       AliFatal("I was requested to apply trigger efficiency, but I could "
631                "not get it !");
632     }
633   }
634   
635   AliDebug(1, Form("Will %s generate noise-only digits for tracker",
636                      (fGenerateNoisyDigits ? "":"NOT")));
637
638   fIsInitialized = kTRUE;
639   return kTRUE;
640 }
641
642 //_____________________________________________________________________________
643 void 
644 AliMUONDigitizerV3::MergeWithSDigits(AliMUONVDigitStore*& outputStore,
645                                      const AliMUONVDigitStore& input,
646                                      Int_t mask)
647 {
648   /// Merge the sdigits in inputData with the digits already present in outputData
649   
650   if ( !outputStore ) outputStore = input.Create();
651   
652   TIter next(input.CreateIterator());
653   AliMUONVDigit* sdigit;
654   
655   while ( ( sdigit = static_cast<AliMUONVDigit*>(next()) ) )
656   {
657     // Update the track references using the mask.
658     // FIXME: this is dirty, for backward compatibility only.
659     // Should re-design all this way of keeping track of MC information...
660     if ( mask ) sdigit->PatchTracks(mask);
661     // Then add or update the digit to the output.
662     AliMUONVDigit* added = outputStore->Add(*sdigit,AliMUONVDigitStore::kMerge);
663     if (!added)
664     {
665       AliError("Could not add digit in merge mode");
666     }
667   }
668 }