Brand-new implementation of Raw->SDigits. One first reconstructs the signal from...
[u/mrichter/AliRoot.git] / VZERO / AliVZERO.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 //                                                                       //
20 //                          V-Zero   Detector                            //
21 //  This class contains the base procedures for the VZERO  detector      //
22 //  Default geometry of November 2003 :   V0R box is 4.4 cm thick        //
23 //                                  scintillators are 2 cm thick         //
24 //  All comments should be sent to Brigitte CHEYNIS :                    //
25 //                                 b.cheynis@ipnl.in2p3.fr               //
26 //                                                                       //
27 //                                                                       //
28 ///////////////////////////////////////////////////////////////////////////
29
30
31 // --- Standard libraries ---
32 #include <Riostream.h>
33 #include <stdlib.h>
34
35 // --- ROOT libraries ---
36 #include <TNamed.h>
37 #include "TROOT.h"
38 #include "TFile.h"
39 #include "TNetFile.h"
40 #include "TRandom.h"
41 #include "TTree.h"
42 #include "TBranch.h"
43 #include "TClonesArray.h"
44 #include "TStopwatch.h"
45 #include "TParameter.h"
46 #include "TF1.h"
47
48 // --- AliRoot header files ---
49 #include "AliRun.h"
50 #include "AliMC.h"
51 #include "AliVZERO.h"
52 #include "AliVZEROLoader.h"
53 #include "AliVZERODigitizer.h"
54 #include "AliVZEROBuffer.h"
55 #include "AliRunDigitizer.h"
56 #include "AliVZEROdigit.h"
57 #include "AliVZEROSDigit.h"
58 #include "AliDAQ.h"
59 #include "AliRawReader.h"
60 #include "AliCDBManager.h"
61 #include "AliCDBEntry.h"
62 #include "AliVZERORawStream.h"
63 #include "AliVZEROCalibData.h"
64 #include "AliVZERORecoParam.h"
65 #include "AliVZEROReconstructor.h"
66
67 ClassImp(AliVZERO)
68  //__________________________________________________________________
69 AliVZERO::AliVZERO(): AliDetector(),
70           fIdSens1(0),
71           fThickness(0.),
72           fThickness1(0.),
73           fMaxStepQua(0.),
74           fMaxStepAlu(0.),
75           fMaxDestepQua(0.),
76           fMaxDestepAlu(0.),
77           fCalibData(NULL),
78           fTimeSlewing(NULL),
79           fSignalShape(NULL),
80           fRecoParam(NULL)
81 {
82 /// Default Constructor
83     
84     AliDebug(1,Form("default (empty) ctor this=%p",this));
85     fIshunt          = 0;
86
87     for(Int_t i = 0 ; i < 64; ++i) {
88       fNBins[i] = 0;
89       fBinSize[i] = 0;
90     }
91 }
92
93 //_____________________________________________________________________________
94 AliVZERO::AliVZERO(const char *name, const char *title)
95        : AliDetector(name,title),
96          fIdSens1(0),
97          fThickness(4.4),
98          fThickness1(2.0),
99          fMaxStepQua(0.05),
100          fMaxStepAlu(0.01),
101          fMaxDestepQua(-1.0),
102          fMaxDestepAlu(-1.0),
103          fCalibData(NULL),
104          fTimeSlewing(NULL),
105          fSignalShape(NULL),
106          fRecoParam(NULL)
107 {
108   
109   // Standard constructor for VZERO Detector
110   
111   AliDebug(1,Form("ctor this=%p",this));
112   
113   //  fIshunt       =  1;  // All hits are associated with primary particles  
114    
115   fHits         =  new TClonesArray("AliVZEROhit", 400);
116   fDigits       =  new TClonesArray("AliVZEROdigit",400); 
117    
118   gAlice->GetMCApp()->AddHitList(fHits);
119
120 //   fThickness    =  4.4;   // total thickness of the V0R box in cm
121 //   fThickness1   =  2.0;   // thickness of scintillating cells in cm
122 //   
123 //   fMaxStepQua   =  0.05; 
124 //   fMaxStepAlu   =  0.01; 
125 //   
126 //   fMaxDestepQua =  -1.0;
127 //   fMaxDestepAlu =  -1.0;
128   
129   for(Int_t i = 0 ; i < 64; ++i) {
130     fNBins[i] = 0;
131     fBinSize[i] = 0;
132   }
133 }
134
135 //_____________________________________________________________________________
136 AliVZERO::~AliVZERO()
137 {
138   //
139   // Default destructor for VZERO Detector
140   //
141   
142     if (fHits) {
143         fHits->Delete();
144         delete fHits;
145         fHits=0; }
146     
147     if (fDigits) {
148         fDigits->Delete();
149         delete fDigits;
150         fDigits=0; }
151     if (fSignalShape) {
152       delete fSignalShape;
153       fSignalShape = NULL;
154     }
155     if (fRecoParam) {
156       delete fRecoParam;
157       fRecoParam = NULL;
158     }
159 }
160
161 //_____________________________________________________________________________
162 void AliVZERO::CreateGeometry()
163 {
164   //
165   // Builds simple Geant3 geometry 
166   //
167 }
168 //_____________________________________________________________________________
169 void AliVZERO::CreateMaterials()
170 {
171   //
172   // Creates materials used for Geant3 geometry 
173   //
174 }
175
176 //_____________________________________________________________________________
177 void AliVZERO::Init()
178 {
179   //
180   // Initialises the VZERO  class after it has been built
181   //
182 }
183
184
185 //_____________________________________________________________________________
186 void AliVZERO::SetMaxStepQua(Float_t p1)
187 {
188   //
189   // Possible parametrisation of steps in active materials
190   //
191      fMaxStepQua = p1;
192 }
193
194 //_____________________________________________________________________________
195 void AliVZERO::SetMaxStepAlu(Float_t p1)
196 {
197   //
198   // Possible parametrisation of steps in Aluminum foils (not used in 
199   // version v2)
200   //
201     fMaxStepAlu = p1;
202 }
203
204 //_____________________________________________________________________________
205 void AliVZERO::SetMaxDestepQua(Float_t p1)
206 {
207   //
208   // Possible parametrisation of steps in active materials (quartz)
209   //
210     fMaxDestepQua = p1;
211 }
212
213 //_____________________________________________________________________________
214 void AliVZERO::SetMaxDestepAlu(Float_t p1)
215 {
216   //
217   // Possible parametrisation of steps in Aluminum (not used in 
218   // version v2)
219   //
220     fMaxDestepAlu = p1;
221 }
222
223 //_____________________________________________________________________________
224 AliLoader* AliVZERO::MakeLoader(const char* topfoldername)
225
226   //
227   // Builds VZEROgetter (AliLoader type)
228   // if detector wants to use customized getter, it must overload this method
229   //
230 //  Info("MakeLoader","Creating AliVZEROLoader. Top folder is %s.",topfoldername); 
231  
232   AliDebug(1,Form("Creating AliVZEROLoader, Top folder is %s ",topfoldername));
233   fLoader = new AliVZEROLoader(GetName(),topfoldername);
234   return fLoader;
235 }
236
237 //_____________________________________________________________________________
238 void AliVZERO::SetTreeAddress()
239 {
240   //
241   // Sets tree address for hits.
242   //
243   if (fLoader->TreeH() && (fHits == 0x0))
244     fHits = new  TClonesArray("AliVZEROhit", 400);
245
246   AliDetector::SetTreeAddress();
247 }
248
249 //_____________________________________________________________________________
250 AliDigitizer* AliVZERO::CreateDigitizer(AliRunDigitizer* manager) const
251 {
252   //
253   // Creates a digitizer for VZERO
254   //
255   return new AliVZERODigitizer(manager);
256 }
257
258 //_____________________________________________________________________________
259 void AliVZERO::Hits2Digits(){
260   //
261   // Converts hits to digits
262   //
263   // Creates the VZERO digitizer 
264   AliVZERODigitizer* dig = new AliVZERODigitizer(this,AliVZERODigitizer::kHits2Digits);
265
266   // Creates the digits
267   dig->Exec("");
268
269 }
270
271 //_____________________________________________________________________________
272 void AliVZERO::Hits2SDigits(){
273   //
274   // Converts hits to summable digits
275   //
276   // Creates the VZERO digitizer 
277   AliVZERODigitizer* dig = new AliVZERODigitizer(this,AliVZERODigitizer::kHits2SDigits);
278
279   // Creates the sdigits
280   dig->Exec("");
281
282 }
283
284 //_____________________________________________________________________________
285 void AliVZERO::Digits2Raw()
286 {
287    //
288    //  Converts digits of the current event to raw data
289    //
290   
291    AliVZERO *fVZERO = (AliVZERO*)gAlice->GetDetector("VZERO");
292    fLoader->LoadDigits();
293    TTree* digits = fLoader->TreeD();
294    if (!digits) {
295       Error("Digits2Raw", "no digits tree");
296       return;
297    }
298    TClonesArray * VZEROdigits = new TClonesArray("AliVZEROdigit",1000);
299    fVZERO->SetTreeAddress();            
300    digits->GetBranch("VZERODigit")->SetAddress(&VZEROdigits); 
301   
302    const char *fileName    = AliDAQ::DdlFileName("VZERO",0);
303    AliVZEROBuffer* buffer  = new AliVZEROBuffer(fileName);
304   
305    // Get Trigger information first
306    // Read trigger inputs from trigger-detector object
307    AliDataLoader * dataLoader = fLoader->GetDigitsDataLoader();
308    if( !dataLoader->IsFileOpen() ) 
309         dataLoader->OpenFile( "READ" );
310    AliTriggerDetector* trgdet = (AliTriggerDetector*)dataLoader->GetDirectory()->Get( "Trigger" );
311    UInt_t triggerInfo = 0;
312    if(trgdet) {
313       triggerInfo = trgdet->GetMask() & 0xffff;
314    }
315    else {
316       AliError(Form("There is no trigger object for %s",fLoader->GetName()));
317    }
318
319    buffer->WriteTriggerInfo((UInt_t)triggerInfo); 
320    buffer->WriteTriggerScalers(); 
321    buffer->WriteBunchNumbers(); 
322   
323    Int_t aBBflagsV0A = 0;
324    Int_t aBBflagsV0C = 0;
325    Int_t aBGflagsV0A = 0;
326    Int_t aBGflagsV0C = 0;
327
328    if (digits->GetUserInfo()->FindObject("BBflagsV0A")) {
329      aBBflagsV0A = ((TParameter<int>*)digits->GetUserInfo()->FindObject("BBflagsV0A"))->GetVal();
330    }
331    else
332      AliWarning("V0A beam-beam flags not found in digits tree UserInfo! The flags will not be written to the raw-data stream!");
333
334    if (digits->GetUserInfo()->FindObject("BBflagsV0C")) {
335      aBBflagsV0C = ((TParameter<int>*)digits->GetUserInfo()->FindObject("BBflagsV0C"))->GetVal();
336    }
337    else
338      AliWarning("V0C beam-beam flags not found in digits tree UserInfo! The flags will not be written to the raw-data stream!");
339
340    if (digits->GetUserInfo()->FindObject("BGflagsV0A")) {
341      aBGflagsV0A = ((TParameter<int>*)digits->GetUserInfo()->FindObject("BGflagsV0A"))->GetVal();
342    }
343    else
344      AliWarning("V0A beam-gas flags not found in digits tree UserInfo! The flags will not be written to the raw-data stream!");
345
346    if (digits->GetUserInfo()->FindObject("BGflagsV0C")) {
347      aBGflagsV0C = ((TParameter<int>*)digits->GetUserInfo()->FindObject("BGflagsV0C"))->GetVal();
348    }
349    else
350      AliWarning("V0C beam-gas flags not found in digits tree UserInfo! The flags will not be written to the raw-data stream!");
351
352    // Now retrieve the channel information: charge smaples + time and 
353    // dump it into ADC and Time arrays
354    Int_t nEntries = Int_t(digits->GetEntries());
355    Short_t aADC[64][AliVZEROdigit::kNClocks];
356    Float_t aTime[64];
357    Float_t aWidth[64];
358    Bool_t  aIntegrator[64];
359    Bool_t  aBBflag[64];
360    Bool_t  aBGflag[64];
361   
362    for (Int_t i = 0; i < nEntries; i++) {
363      fVZERO->ResetDigits();
364      digits->GetEvent(i);
365      Int_t ndig = VZEROdigits->GetEntriesFast(); 
366    
367      if(ndig == 0) continue;
368      for(Int_t k=0; k<ndig; k++){
369          AliVZEROdigit* fVZERODigit = (AliVZEROdigit*) VZEROdigits->At(k);
370          // Convert aliroot channel k into FEE channel iChannel before writing data
371          Int_t iChannel       = AliVZEROCalibData::GetBoardNumber(fVZERODigit->PMNumber()) * 8 +
372            AliVZEROCalibData::GetFEEChannelNumber(fVZERODigit->PMNumber());
373          for(Int_t iClock = 0; iClock < AliVZEROdigit::kNClocks; ++iClock) aADC[iChannel][iClock] = fVZERODigit->ChargeADC(iClock);
374          aTime[iChannel]      = fVZERODigit->Time();
375          aWidth[iChannel]     = fVZERODigit->Width();
376          aIntegrator[iChannel]= fVZERODigit->Integrator();
377          if(fVZERODigit->PMNumber() < 32) {
378            aBBflag[iChannel] = (aBBflagsV0C >> fVZERODigit->PMNumber()) & 0x1;
379            aBGflag[iChannel] = (aBGflagsV0C >> fVZERODigit->PMNumber()) & 0x1;
380          }
381          else {
382            aBBflag[iChannel] = (aBBflagsV0A >> (fVZERODigit->PMNumber()-32)) & 0x1;
383            aBGflag[iChannel] = (aBGflagsV0A >> (fVZERODigit->PMNumber()-32)) & 0x1;
384          }
385          AliDebug(1,Form("DDL: %s\tdigit number: %d\tPM number: %d\tADC: %f\tTime: %f",
386                          fileName,k,fVZERODigit->PMNumber(),aADC[k],aTime[k])); 
387      }        
388    }
389
390    // Now fill raw data
391           
392    for (Int_t  iCIU = 0; iCIU < 8; iCIU++) { 
393  
394    // decoding of one Channel Interface Unit numbered iCIU - there are 8 channels per CIU (and 8 CIUs) :
395   
396       for(Int_t iChannel_Offset = iCIU*8; iChannel_Offset < (iCIU*8)+8; iChannel_Offset=iChannel_Offset+4) { 
397          for(Int_t iChannel = iChannel_Offset; iChannel < iChannel_Offset+4; iChannel++) {
398              buffer->WriteChannel(iChannel, aADC[iChannel], aIntegrator[iChannel]);       
399          }
400          buffer->WriteBeamFlags(&aBBflag[iChannel_Offset],&aBGflag[iChannel_Offset]); 
401          buffer->WriteMBInfo(); 
402          buffer->WriteMBFlags();   
403          buffer->WriteBeamScalers(); 
404       } 
405
406       for(Int_t iChannel = iCIU*8 + 7; iChannel >= iCIU*8; iChannel--) {
407           buffer->WriteTiming(aTime[iChannel], aWidth[iChannel]); 
408       }
409
410     // End of decoding of one CIU card
411     
412   } // end of decoding the eight CIUs
413      
414   delete buffer;
415   fLoader->UnloadDigits();  
416 }
417
418 //_____________________________________________________________________________
419 Bool_t AliVZERO::Raw2SDigits(AliRawReader* rawReader){
420   // Converts the VZERO raw data into digits
421   // The method is used for merging simulated and
422   // real data events
423   TStopwatch timer;
424   timer.Start();
425
426   if(!fLoader) {
427     AliError("no VZERO loader found");
428     return kFALSE;
429   }
430   fLoader->LoadSDigits("UPDATE");
431
432   if (!fLoader->TreeS()) fLoader->MakeTree("S");
433   fLoader->MakeSDigitsContainer();
434   TTree* treeS  = fLoader->TreeS();
435
436   TClonesArray *sdigits = new TClonesArray("AliVZEROSDigit", 64);
437   treeS->Branch("VZEROSDigit", &sdigits); 
438
439   {
440     rawReader->Reset();
441     AliVZERORawStream rawStream(rawReader);    
442      
443     if (!rawStream.Next()) return kFALSE; // No VZERO data found
444
445     GetCalibData();
446
447     Int_t nSDigits = 0;
448     Float_t *charges = NULL;
449     Int_t nbins = 0;
450     for(Int_t iChannel=0; iChannel < 64; ++iChannel) {
451       Int_t offlineCh = rawStream.GetOfflineChannel(iChannel);
452       Short_t chargeADC[AliVZEROdigit::kNClocks];
453       for(Int_t iClock=0; iClock < AliVZEROdigit::kNClocks; ++iClock) {
454         chargeADC[iClock] = rawStream.GetPedestal(iChannel,iClock);
455       }
456       // Integrator flag
457       Bool_t integrator = rawStream.GetIntegratorFlag(iChannel,AliVZEROdigit::kNClocks/2);
458       // HPTDC data (leading time and width)
459       Int_t board = AliVZEROCalibData::GetBoardNumber(offlineCh);
460       Float_t time = rawStream.GetTime(iChannel)*fCalibData->GetTimeResolution(board);
461       //      Float_t width = rawStream.GetWidth(iChannel)*fCalibData->GetWidthResolution(board);
462       Float_t adc = 0;
463
464       // Pedestal retrieval and suppression
465       Float_t maxadc = 0;
466       Int_t imax = -1;
467       Float_t adcPedSub[AliVZEROdigit::kNClocks];
468       Float_t integral = fSignalShape->Integral(0,200);
469       for(Int_t iClock=0; iClock < AliVZEROdigit::kNClocks; ++iClock) {
470         Bool_t iIntegrator = (iClock%2 == 0) ? integrator : !integrator;
471         Int_t k = offlineCh + 64*iIntegrator;
472         adcPedSub[iClock] = (Float_t)chargeADC[iClock] - fCalibData->GetPedestal(k);
473         if(adcPedSub[iClock] <= fRecoParam->GetNSigmaPed()*fCalibData->GetSigma(k)) {
474           adcPedSub[iClock] = 0;
475           continue;
476         }
477         if(iClock < fRecoParam->GetStartClock() || iClock > fRecoParam->GetEndClock()) continue;
478         if(adcPedSub[iClock] > maxadc) {
479           maxadc = adcPedSub[iClock];
480           imax   = iClock;
481         }
482       }
483       if (imax != -1) {
484         Int_t start = imax - fRecoParam->GetNPreClocks();
485         if (start < 0) start = 0;
486         Int_t end = imax + fRecoParam->GetNPostClocks();
487         if (end > 20) end = 20;
488         for(Int_t iClock = start; iClock <= end; iClock++) {
489           adc += adcPedSub[iClock];
490         }
491       }
492       Float_t correctedTime = CorrectLeadingTime(offlineCh,time,adc);
493
494       if (!charges) {
495         nbins = fNBins[offlineCh];
496         charges = new Float_t[nbins];
497       }
498       else if (nbins != fNBins[offlineCh]) {
499         delete [] charges;
500         nbins = fNBins[offlineCh];
501         charges = new Float_t[nbins];
502       }
503       memset(charges,0,nbins*sizeof(Float_t));
504
505       // Now lets produce SDigit
506       if ((correctedTime > (AliVZEROReconstructor::kInvalidTime + 1e-6)) &&
507           (adc > 1e-6)) {
508         for(Int_t iBin = 0; iBin < nbins; ++iBin) {
509           Float_t t = fBinSize[offlineCh]*Float_t(iBin);
510           if ((t < correctedTime) ||
511               (t > (correctedTime+200.))) continue;
512           charges[iBin] = kChargePerADC*adc*(fSignalShape->Eval(t-correctedTime)*fBinSize[offlineCh]/integral);
513         }
514       }
515
516       TClonesArray &sdigitsref = *sdigits;  
517       new (sdigitsref[nSDigits++]) AliVZEROSDigit(offlineCh,fNBins[offlineCh],charges);
518     }
519     if (charges) delete [] charges;
520   }
521   
522   treeS->Fill();
523   fLoader->WriteSDigits("OVERWRITE");
524   fLoader->UnloadSDigits();     
525         
526   timer.Stop();
527   timer.Print();
528   return kTRUE;
529 }
530
531 //_____________________________________________________________________________
532 void AliVZERO::GetCalibData()
533 {
534   // Gets calibration object for VZERO set
535   // Do nothing in case it is already loaded
536   if (fCalibData) return;
537
538   AliCDBEntry *entry = AliCDBManager::Instance()->Get("VZERO/Calib/Data");
539   if (entry) fCalibData = (AliVZEROCalibData*) entry->GetObject();
540   if (!fCalibData)  AliFatal("No calibration data from calibration database !");
541
542   AliCDBEntry *entry2 = AliCDBManager::Instance()->Get("VZERO/Calib/TimeSlewing");
543   if (!entry2) AliFatal("VZERO time slewing function is not found in OCDB !");
544   fTimeSlewing = (TF1*)entry2->GetObject();
545
546   for(Int_t i = 0 ; i < 64; ++i) {
547     Int_t board = AliVZEROCalibData::GetBoardNumber(i);
548     fNBins[i] = TMath::Nint(((Float_t)(fCalibData->GetMatchWindow(board)+1)*25.0+
549                              (Float_t)kMaxTDCWidth*fCalibData->GetWidthResolution(board))/
550                             fCalibData->GetTimeResolution(board));
551     fBinSize[i] = fCalibData->GetTimeResolution(board);
552   }
553
554   fSignalShape = new TF1("VZEROSDigitSignalShape",this,&AliVZERO::SignalShape,0,200,6,"AliVZERO","SignalShape");
555   fSignalShape->SetParameters(0,1.57345e1,-4.25603e-1,
556                               2.9,6.40982,3.69339e-01);
557
558   fRecoParam = new AliVZERORecoParam;
559
560   return;
561 }
562
563 Float_t AliVZERO::CorrectLeadingTime(Int_t i, Float_t time, Float_t adc) const
564 {
565   // Correct the leading time
566   // for slewing effect and
567   // misalignment of the channels
568   if (time < 1e-6) return -1024;
569
570   // In case of pathological signals
571   if (adc < 1e-6) return time;
572
573   // Slewing correction
574   Float_t thr = fCalibData->GetDiscriThr(i);
575   time -= fTimeSlewing->Eval(adc/thr);
576
577   return time;
578 }
579
580 double AliVZERO::SignalShape(double *x, double *par)
581 {
582   // this function simulates the signal
583   // shape used in Raw->SDigits method
584
585   Double_t xx = x[0];
586   if (xx <= par[0]) return 0;
587   Double_t a = 1./TMath::Power((xx-par[0])/par[1],1./par[2]);
588   if (xx <= par[3]) return a;
589   Double_t b = 1./TMath::Power((xx-par[3])/par[4],1./par[5]);
590   Double_t f = a*b/(a+b);
591   AliDebug(100,Form("x=%f func=%f",xx,f));
592   return f;
593 }