DA output for test shuttle.
[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   // deletes the digitizer
270   delete dig;
271 }
272
273 //_____________________________________________________________________________
274 void AliVZERO::Hits2SDigits(){
275   //
276   // Converts hits to summable digits
277   //
278   // Creates the VZERO digitizer 
279   AliVZERODigitizer* dig = new AliVZERODigitizer(this,AliVZERODigitizer::kHits2SDigits);
280
281   // Creates the sdigits
282   dig->Exec("");
283
284   // deletes the digitizer
285   delete dig;
286 }
287
288 //_____________________________________________________________________________
289 void AliVZERO::Digits2Raw()
290 {
291    //
292    //  Converts digits of the current event to raw data
293    //
294   
295    AliVZERO *fVZERO = (AliVZERO*)gAlice->GetDetector("VZERO");
296    fLoader->LoadDigits();
297    TTree* digits = fLoader->TreeD();
298    if (!digits) {
299       Error("Digits2Raw", "no digits tree");
300       return;
301    }
302    TClonesArray * VZEROdigits = new TClonesArray("AliVZEROdigit",1000);
303    fVZERO->SetTreeAddress();            
304    digits->GetBranch("VZERODigit")->SetAddress(&VZEROdigits); 
305   
306    const char *fileName    = AliDAQ::DdlFileName("VZERO",0);
307    AliVZEROBuffer* buffer  = new AliVZEROBuffer(fileName);
308   
309    // Get Trigger information first
310    // Read trigger inputs from trigger-detector object
311    AliDataLoader * dataLoader = fLoader->GetDigitsDataLoader();
312    if( !dataLoader->IsFileOpen() ) 
313         dataLoader->OpenFile( "READ" );
314    AliTriggerDetector* trgdet = (AliTriggerDetector*)dataLoader->GetDirectory()->Get( "Trigger" );
315    UInt_t triggerInfo = 0;
316    if(trgdet) {
317       triggerInfo = trgdet->GetMask() & 0xffff;
318    }
319    else {
320       AliError(Form("There is no trigger object for %s",fLoader->GetName()));
321    }
322
323    buffer->WriteTriggerInfo((UInt_t)triggerInfo); 
324    buffer->WriteTriggerScalers(); 
325    buffer->WriteBunchNumbers(); 
326   
327    Int_t aBBflagsV0A = 0;
328    Int_t aBBflagsV0C = 0;
329    Int_t aBGflagsV0A = 0;
330    Int_t aBGflagsV0C = 0;
331
332    if (digits->GetUserInfo()->FindObject("BBflagsV0A")) {
333      aBBflagsV0A = ((TParameter<int>*)digits->GetUserInfo()->FindObject("BBflagsV0A"))->GetVal();
334    }
335    else
336      AliWarning("V0A beam-beam flags not found in digits tree UserInfo! The flags will not be written to the raw-data stream!");
337
338    if (digits->GetUserInfo()->FindObject("BBflagsV0C")) {
339      aBBflagsV0C = ((TParameter<int>*)digits->GetUserInfo()->FindObject("BBflagsV0C"))->GetVal();
340    }
341    else
342      AliWarning("V0C beam-beam flags not found in digits tree UserInfo! The flags will not be written to the raw-data stream!");
343
344    if (digits->GetUserInfo()->FindObject("BGflagsV0A")) {
345      aBGflagsV0A = ((TParameter<int>*)digits->GetUserInfo()->FindObject("BGflagsV0A"))->GetVal();
346    }
347    else
348      AliWarning("V0A beam-gas flags not found in digits tree UserInfo! The flags will not be written to the raw-data stream!");
349
350    if (digits->GetUserInfo()->FindObject("BGflagsV0C")) {
351      aBGflagsV0C = ((TParameter<int>*)digits->GetUserInfo()->FindObject("BGflagsV0C"))->GetVal();
352    }
353    else
354      AliWarning("V0C beam-gas flags not found in digits tree UserInfo! The flags will not be written to the raw-data stream!");
355
356    // Now retrieve the channel information: charge smaples + time and 
357    // dump it into ADC and Time arrays
358    Int_t nEntries = Int_t(digits->GetEntries());
359    Short_t aADC[64][AliVZEROdigit::kNClocks];
360    Float_t aTime[64];
361    Float_t aWidth[64];
362    Bool_t  aIntegrator[64];
363    Bool_t  aBBflag[64];
364    Bool_t  aBGflag[64];
365   
366    for (Int_t i = 0; i < nEntries; i++) {
367      fVZERO->ResetDigits();
368      digits->GetEvent(i);
369      Int_t ndig = VZEROdigits->GetEntriesFast(); 
370    
371      if(ndig == 0) continue;
372      for(Int_t k=0; k<ndig; k++){
373          AliVZEROdigit* fVZERODigit = (AliVZEROdigit*) VZEROdigits->At(k);
374          // Convert aliroot channel k into FEE channel iChannel before writing data
375          Int_t iChannel       = AliVZEROCalibData::GetBoardNumber(fVZERODigit->PMNumber()) * 8 +
376            AliVZEROCalibData::GetFEEChannelNumber(fVZERODigit->PMNumber());
377          for(Int_t iClock = 0; iClock < AliVZEROdigit::kNClocks; ++iClock) aADC[iChannel][iClock] = fVZERODigit->ChargeADC(iClock);
378          aTime[iChannel]      = fVZERODigit->Time();
379          aWidth[iChannel]     = fVZERODigit->Width();
380          aIntegrator[iChannel]= fVZERODigit->Integrator();
381          if(fVZERODigit->PMNumber() < 32) {
382            aBBflag[iChannel] = (aBBflagsV0C >> fVZERODigit->PMNumber()) & 0x1;
383            aBGflag[iChannel] = (aBGflagsV0C >> fVZERODigit->PMNumber()) & 0x1;
384          }
385          else {
386            aBBflag[iChannel] = (aBBflagsV0A >> (fVZERODigit->PMNumber()-32)) & 0x1;
387            aBGflag[iChannel] = (aBGflagsV0A >> (fVZERODigit->PMNumber()-32)) & 0x1;
388          }
389          AliDebug(1,Form("DDL: %s\tdigit number: %d\tPM number: %d\tADC: %d\tTime: %f",
390                          fileName,k,fVZERODigit->PMNumber(),aADC[iChannel][AliVZEROdigit::kNClocks/2],aTime[iChannel])); 
391      }        
392    }
393
394    // Now fill raw data
395           
396    for (Int_t  iCIU = 0; iCIU < 8; iCIU++) { 
397  
398    // decoding of one Channel Interface Unit numbered iCIU - there are 8 channels per CIU (and 8 CIUs) :
399   
400       for(Int_t iChannel_Offset = iCIU*8; iChannel_Offset < (iCIU*8)+8; iChannel_Offset=iChannel_Offset+4) { 
401          for(Int_t iChannel = iChannel_Offset; iChannel < iChannel_Offset+4; iChannel++) {
402              buffer->WriteChannel(iChannel, aADC[iChannel], aIntegrator[iChannel]);       
403          }
404          buffer->WriteBeamFlags(&aBBflag[iChannel_Offset],&aBGflag[iChannel_Offset]); 
405          buffer->WriteMBInfo(); 
406          buffer->WriteMBFlags();   
407          buffer->WriteBeamScalers(); 
408       } 
409
410       for(Int_t iChannel = iCIU*8 + 7; iChannel >= iCIU*8; iChannel--) {
411           buffer->WriteTiming(aTime[iChannel], aWidth[iChannel]); 
412       }
413
414     // End of decoding of one CIU card
415     
416   } // end of decoding the eight CIUs
417      
418   delete buffer;
419   fLoader->UnloadDigits();  
420 }
421
422 //_____________________________________________________________________________
423 Bool_t AliVZERO::Raw2SDigits(AliRawReader* rawReader){
424   // Converts the VZERO raw data into digits
425   // The method is used for merging simulated and
426   // real data events
427   TStopwatch timer;
428   timer.Start();
429
430   if(!fLoader) {
431     AliError("no VZERO loader found");
432     return kFALSE;
433   }
434   fLoader->LoadSDigits("UPDATE");
435
436   if (!fLoader->TreeS()) fLoader->MakeTree("S");
437   fLoader->MakeSDigitsContainer();
438   TTree* treeS  = fLoader->TreeS();
439
440   TClonesArray *sdigits = new TClonesArray("AliVZEROSDigit", 64);
441   treeS->Branch("VZEROSDigit", &sdigits); 
442
443   {
444     rawReader->Reset();
445     AliVZERORawStream rawStream(rawReader);    
446      
447     if (!rawStream.Next()) return kFALSE; // No VZERO data found
448
449     GetCalibData();
450
451     Int_t nSDigits = 0;
452     Float_t *charges = NULL;
453     Int_t nbins = 0;
454     for(Int_t iChannel=0; iChannel < 64; ++iChannel) {
455       Int_t offlineCh = rawStream.GetOfflineChannel(iChannel);
456       Short_t chargeADC[AliVZEROdigit::kNClocks];
457       for(Int_t iClock=0; iClock < AliVZEROdigit::kNClocks; ++iClock) {
458         chargeADC[iClock] = rawStream.GetPedestal(iChannel,iClock);
459       }
460       // Integrator flag
461       Bool_t integrator = rawStream.GetIntegratorFlag(iChannel,AliVZEROdigit::kNClocks/2);
462       // HPTDC data (leading time and width)
463       Int_t board = AliVZEROCalibData::GetBoardNumber(offlineCh);
464       Float_t time = rawStream.GetTime(iChannel)*fCalibData->GetTimeResolution(board);
465       //      Float_t width = rawStream.GetWidth(iChannel)*fCalibData->GetWidthResolution(board);
466       Float_t adc = 0;
467
468       // Pedestal retrieval and suppression
469       Float_t maxadc = 0;
470       Int_t imax = -1;
471       Float_t adcPedSub[AliVZEROdigit::kNClocks];
472       Float_t integral = fSignalShape->Integral(0,200);
473       for(Int_t iClock=0; iClock < AliVZEROdigit::kNClocks; ++iClock) {
474         Bool_t iIntegrator = (iClock%2 == 0) ? integrator : !integrator;
475         Int_t k = offlineCh + 64*iIntegrator;
476         adcPedSub[iClock] = (Float_t)chargeADC[iClock] - fCalibData->GetPedestal(k);
477         if(adcPedSub[iClock] <= fRecoParam->GetNSigmaPed()*fCalibData->GetSigma(k)) {
478           adcPedSub[iClock] = 0;
479           continue;
480         }
481         if(iClock < fRecoParam->GetStartClock() || iClock > fRecoParam->GetEndClock()) continue;
482         if(adcPedSub[iClock] > maxadc) {
483           maxadc = adcPedSub[iClock];
484           imax   = iClock;
485         }
486       }
487       if (imax != -1) {
488         Int_t start = imax - fRecoParam->GetNPreClocks();
489         if (start < 0) start = 0;
490         Int_t end = imax + fRecoParam->GetNPostClocks();
491         if (end > 20) end = 20;
492         for(Int_t iClock = start; iClock <= end; iClock++) {
493           adc += adcPedSub[iClock];
494         }
495       }
496       Float_t correctedTime = CorrectLeadingTime(offlineCh,time,adc);
497
498       if (!charges) {
499         nbins = fNBins[offlineCh];
500         charges = new Float_t[nbins];
501       }
502       else if (nbins != fNBins[offlineCh]) {
503         delete [] charges;
504         nbins = fNBins[offlineCh];
505         charges = new Float_t[nbins];
506       }
507       memset(charges,0,nbins*sizeof(Float_t));
508
509       // Now lets produce SDigit
510       if ((correctedTime > (AliVZEROReconstructor::kInvalidTime + 1e-6)) &&
511           (adc > 1e-6)) {
512         for(Int_t iBin = 0; iBin < nbins; ++iBin) {
513           Float_t t = fBinSize[offlineCh]*Float_t(iBin);
514           if ((t < correctedTime) ||
515               (t > (correctedTime+200.))) continue;
516           charges[iBin] = kChargePerADC*adc*(fSignalShape->Eval(t-correctedTime)*fBinSize[offlineCh]/integral);
517         }
518       }
519
520       TClonesArray &sdigitsref = *sdigits;  
521       new (sdigitsref[nSDigits++]) AliVZEROSDigit(offlineCh,fNBins[offlineCh],charges);
522     }
523     if (charges) delete [] charges;
524   }
525   
526   treeS->Fill();
527   fLoader->WriteSDigits("OVERWRITE");
528   fLoader->UnloadSDigits();     
529         
530   timer.Stop();
531   timer.Print();
532   return kTRUE;
533 }
534
535 //_____________________________________________________________________________
536 void AliVZERO::GetCalibData()
537 {
538   // Gets calibration object for VZERO set
539   // Do nothing in case it is already loaded
540   if (fCalibData) return;
541
542   AliCDBEntry *entry = AliCDBManager::Instance()->Get("VZERO/Calib/Data");
543   if (entry) fCalibData = (AliVZEROCalibData*) entry->GetObject();
544   if (!fCalibData)  AliFatal("No calibration data from calibration database !");
545
546   AliCDBEntry *entry2 = AliCDBManager::Instance()->Get("VZERO/Calib/TimeSlewing");
547   if (!entry2) AliFatal("VZERO time slewing function is not found in OCDB !");
548   fTimeSlewing = (TF1*)entry2->GetObject();
549
550   for(Int_t i = 0 ; i < 64; ++i) {
551     Int_t board = AliVZEROCalibData::GetBoardNumber(i);
552     fNBins[i] = TMath::Nint(((Float_t)(fCalibData->GetMatchWindow(board)+1)*25.0+
553                              (Float_t)kMaxTDCWidth*fCalibData->GetWidthResolution(board))/
554                             fCalibData->GetTimeResolution(board));
555     fBinSize[i] = fCalibData->GetTimeResolution(board);
556   }
557
558   fSignalShape = new TF1("VZEROSDigitSignalShape",this,&AliVZERO::SignalShape,0,200,6,"AliVZERO","SignalShape");
559   fSignalShape->SetParameters(0,1.57345e1,-4.25603e-1,
560                               2.9,6.40982,3.69339e-01);
561
562   fRecoParam = new AliVZERORecoParam;
563
564   return;
565 }
566
567 Float_t AliVZERO::CorrectLeadingTime(Int_t i, Float_t time, Float_t adc) const
568 {
569   // Correct the leading time
570   // for slewing effect and
571   // misalignment of the channels
572   if (time < 1e-6) return -1024;
573
574   // In case of pathological signals
575   if (adc < 1e-6) return time;
576
577   // Slewing correction
578   Float_t thr = fCalibData->GetCalibDiscriThr(i,kTRUE);
579   time -= fTimeSlewing->Eval(adc/thr);
580
581   return time;
582 }
583
584 double AliVZERO::SignalShape(double *x, double *par)
585 {
586   // this function simulates the signal
587   // shape used in Raw->SDigits method
588
589   Double_t xx = x[0];
590   if (xx <= par[0]) return 0;
591   Double_t a = 1./TMath::Power((xx-par[0])/par[1],1./par[2]);
592   if (xx <= par[3]) return a;
593   Double_t b = 1./TMath::Power((xx-par[3])/par[4],1./par[5]);
594   Double_t f = a*b/(a+b);
595   AliDebug(100,Form("x=%f func=%f",xx,f));
596   return f;
597 }