X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=VZERO%2FAliVZERODigitizer.cxx;h=193206e76354c3bd8bd91015bd7c78edf61d2924;hb=a39f5b706fdf439da800dfe272e80fe7061ff504;hp=abd88987f4405f1b6d8262073744ae261bc6d948;hpb=df791951d1717244a50f752e90695ef60825624f;p=u%2Fmrichter%2FAliRoot.git diff --git a/VZERO/AliVZERODigitizer.cxx b/VZERO/AliVZERODigitizer.cxx index abd88987f44..193206e7635 100644 --- a/VZERO/AliVZERODigitizer.cxx +++ b/VZERO/AliVZERODigitizer.cxx @@ -26,70 +26,91 @@ // --- ROOT system --- #include #include +#include +#include +#include +#include #include +#include +#include // --- AliRoot header files --- -#include "AliVZEROConst.h" #include "AliRun.h" #include "AliVZERO.h" #include "AliVZEROhit.h" #include "AliRunLoader.h" #include "AliLoader.h" +#include "AliGRPObject.h" #include "AliRunDigitizer.h" #include "AliCDBManager.h" #include "AliCDBStorage.h" #include "AliCDBEntry.h" #include "AliVZEROCalibData.h" - +#include "AliCTPTimeParams.h" +#include "AliLHCClockPhase.h" #include "AliVZEROdigit.h" #include "AliVZERODigitizer.h" +#include "AliVZEROSDigit.h" ClassImp(AliVZERODigitizer) AliVZERODigitizer::AliVZERODigitizer() - :AliDigitizer(), - fCalibData(GetCalibData()), - fPhotoCathodeEfficiency(0.18), - fPMVoltage(768.0), - fPMGain(TMath::Power((fPMVoltage / 112.5) ,7.04277)), - fNdigits(0), - fDigits(0) - + :AliDigitizer(), + fCalibData(GetCalibData()), + fPhotoCathodeEfficiency(0.18), + fNdigits(0), + fDigits(0), + fSignalShape(NULL), + fPMResponse(NULL), + fSinglePhESpectrum(NULL), + fEvenOrOdd(kFALSE), + fTask(kHits2Digits), + fVZERO(NULL) { // default constructor + // Initialize OCDB and containers used in the digitization -// fNdigits = 0; -// fDigits = 0; -// -// fPhotoCathodeEfficiency = 0.18; -// fPMVoltage = 768.0; -// fPMGain = TMath::Power((fPMVoltage / 112.5) ,7.04277); - -// fCalibData = GetCalibData(); + Init(); } +//____________________________________________________________________________ + AliVZERODigitizer::AliVZERODigitizer(AliVZERO *vzero, DigiTask_t task) + :AliDigitizer(), + fCalibData(GetCalibData()), + fPhotoCathodeEfficiency(0.18), + fNdigits(0), + fDigits(0), + fSignalShape(NULL), + fPMResponse(NULL), + fSinglePhESpectrum(NULL), + fEvenOrOdd(kFALSE), + fTask(task), + fVZERO(vzero) +{ + // constructor + // Initialize OCDB and containers used in the digitization + + Init(); +} + //____________________________________________________________________________ AliVZERODigitizer::AliVZERODigitizer(AliRunDigitizer* manager) :AliDigitizer(manager), fCalibData(GetCalibData()), fPhotoCathodeEfficiency(0.18), - fPMVoltage(768.0), - fPMGain(TMath::Power((fPMVoltage / 112.5) ,7.04277)), fNdigits(0), - fDigits(0) - + fDigits(0), + fSignalShape(NULL), + fPMResponse(NULL), + fSinglePhESpectrum(NULL), + fEvenOrOdd(kFALSE), + fTask(kHits2Digits), + fVZERO(NULL) { // constructor - -// fNdigits = 0; -// fDigits = 0; -// -// fPhotoCathodeEfficiency = 0.18; -// fPMVoltage = 768.0; -// fPMGain = TMath::Power( (fPMVoltage / 112.5) ,7.04277 ); - -// fCalibData = GetCalibData(); - + // Initialize OCDB and containers used in the digitization + + Init(); } //____________________________________________________________________________ @@ -102,15 +123,97 @@ ClassImp(AliVZERODigitizer) delete fDigits; fDigits=0; } + + if (fSignalShape) { + delete fSignalShape; + fSignalShape = NULL; + } + if (fPMResponse) { + delete fPMResponse; + fPMResponse = NULL; + } + if (fSinglePhESpectrum) { + delete fSinglePhESpectrum; + fSinglePhESpectrum = NULL; + } + + for(Int_t i = 0 ; i < 64; ++i) { + if (fTime[i]) delete [] fTime[i]; + } } //_____________________________________________________________________________ Bool_t AliVZERODigitizer::Init() { // Initialises the digitizer - - // Initialises the Digit array - fDigits = new TClonesArray ("AliVZEROdigit", 1000); + // Initialize OCDB and containers used in the digitization + + // check if the digitizer was already initialized + if (fSignalShape) return kTRUE; + + fSignalShape = new TF1("VZEROSignalShape",this,&AliVZERODigitizer::SignalShape,0,200,6,"AliVZERODigitizer","SignalShape"); + // fSignalShape->SetParameters(0,1.57345e1,-4.25603e-1,2.9,6.40982,3.69339e-01); + // fSignalShape->SetParameters(1.34330e+00,1.13007e+02,-4.95705e-01, + // 3.68911e+00,1.01040e+00, 3.94675e-01); + fSignalShape->SetParameters(-1.07335e+00,2.16002e+01,-1.26133e-01, + 1.41619e+00,5.50334e-01,3.86111e-01); + fPMResponse = new TF1("VZEROPMResponse",this,&AliVZERODigitizer::PMResponse,-kPMRespTime,2.*kPMRespTime,0,"AliVZERODigitizer","PMResponse"); + fSinglePhESpectrum = new TF1("VZEROSinglePhESpectrum",this,&AliVZERODigitizer::SinglePhESpectrum,0,20,0,"AliVZERODigitizer","SinglePhESpectrum"); + + // Now get the CTP L0->L1 delay + AliCDBEntry *entry = AliCDBManager::Instance()->Get("GRP/CTP/CTPtiming"); + if (!entry) AliFatal("CTP timing parameters are not found in OCDB !"); + AliCTPTimeParams *ctpParams = (AliCTPTimeParams*)entry->GetObject(); + Float_t l1Delay = (Float_t)ctpParams->GetDelayL1L0()*25.0; + + AliCDBEntry *entry1 = AliCDBManager::Instance()->Get("GRP/CTP/TimeAlign"); + if (!entry1) AliFatal("CTP time-alignment is not found in OCDB !"); + AliCTPTimeParams *ctpTimeAlign = (AliCTPTimeParams*)entry1->GetObject(); + l1Delay += ((Float_t)ctpTimeAlign->GetDelayL1L0()*25.0); + + AliCDBEntry *entry2 = AliCDBManager::Instance()->Get("VZERO/Calib/TimeDelays"); + if (!entry2) AliFatal("VZERO time delays are not found in OCDB !"); + TH1F *delays = (TH1F*)entry2->GetObject(); + + AliCDBEntry *entry3 = AliCDBManager::Instance()->Get("GRP/Calib/LHCClockPhase"); + if (!entry3) AliFatal("LHC clock-phase shift is not found in OCDB !"); + AliLHCClockPhase *phase = (AliLHCClockPhase*)entry3->GetObject(); + + for(Int_t i = 0 ; i < 64; ++i) { + + for(Int_t j = 0; j < kNClocks; ++j) fAdc[i][j] = 0; + fLeadingTime[i] = fTimeWidth[i] = 0; + + fPmGain[i] = fCalibData->GetGain(i); + + fAdcPedestal[i][0] = fCalibData->GetPedestal(i); + fAdcSigma[i][0] = fCalibData->GetSigma(i); + fAdcPedestal[i][1] = fCalibData->GetPedestal(i+64); + fAdcSigma[i][1] = fCalibData->GetSigma(i+64); + + Int_t board = AliVZEROCalibData::GetBoardNumber(i); + fNBins[i] = TMath::Nint(((Float_t)(fCalibData->GetMatchWindow(board)+1)*25.0+ + (Float_t)kMaxTDCWidth*fCalibData->GetWidthResolution(board))/ + fCalibData->GetTimeResolution(board)); + fNBinsLT[i] = TMath::Nint(((Float_t)(fCalibData->GetMatchWindow(board)+1)*25.0)/ + fCalibData->GetTimeResolution(board)); + fBinSize[i] = fCalibData->GetTimeResolution(board); + fHptdcOffset[i] = (((Float_t)fCalibData->GetRollOver(board)- + (Float_t)fCalibData->GetTriggerCountOffset(board))*25.0+ + fCalibData->GetTimeOffset(i)- + l1Delay- + phase->GetMeanPhase()+ + delays->GetBinContent(i+1)+ + kV0Offset); + fClockOffset[i] = (((Float_t)fCalibData->GetRollOver(board)- + (Float_t)fCalibData->GetTriggerCountOffset(board))*25.0+ + fCalibData->GetTimeOffset(i)- + l1Delay+ + kV0Offset); + + fTime[i] = new Float_t[fNBins[i]]; + memset(fTime[i],0,fNBins[i]*sizeof(Float_t)); + } return kTRUE; } @@ -119,165 +222,74 @@ Bool_t AliVZERODigitizer::Init() void AliVZERODigitizer::Exec(Option_t* /*option*/) { // Creates digits from hits - - Int_t map[80]; // 48 values on V0C + 32 on V0A - Int_t adc[64]; // 32 PMs on V0C + 32 PMs on V0A - Float_t time[80], time_ref[80], time2[64]; - Float_t adc_gain[80]; - fNdigits = 0; - Float_t cPM = fPhotoCathodeEfficiency * fPMGain; - - // Retrieval of ADC gain values from local CDB - // I use only the first 64th values of the calibration array in CDB - // as I have no beam burst structure - odd or even beam burst number - - // Reminder : We have 16 scintillating cells mounted on 8 PMs - // on Ring 3 and Ring 4 in V0C - added to produce ADC outputs - // on these rings... - - for(Int_t i=0; i<16; i++) { adc_gain[i] = fCalibData->GetGain(i) ; }; - - for(Int_t j=16; j<48; j=j+2) { - Int_t i=(j+17)/2; - adc_gain[j] = fCalibData->GetGain(i) ; - adc_gain[j+1] = fCalibData->GetGain(i) ; } - for(Int_t i=48; i<80; i++) { adc_gain[i] = fCalibData->GetGain(i-16) ; } - -// for(Int_t i=0; i<80; i++) { printf(" i = %d gain = %f\n\n", i, adc_gain[i] );} - - AliRunLoader* outRunLoader = - AliRunLoader::GetRunLoader(fManager->GetOutputFolderName()); - if (!outRunLoader) { - Error("Exec", "Can not get output Run Loader"); - return;} - - AliLoader* outLoader = outRunLoader->GetLoader("VZEROLoader"); - if (!outLoader) { - Error("Exec", "Can not get output VZERO Loader"); - return;} - - outLoader->LoadDigits("update"); - if (!outLoader->TreeD()) outLoader->MakeTree("D"); - outLoader->MakeDigitsContainer(); - TTree* treeD = outLoader->TreeD(); - Int_t bufsize = 16000; - treeD->Branch("VZERODigit", &fDigits, bufsize); - - for (Int_t iInput = 0; iInput < fManager->GetNinputs(); iInput++) { - AliRunLoader* runLoader = - AliRunLoader::GetRunLoader(fManager->GetInputFolderName(iInput)); - AliLoader* loader = runLoader->GetLoader("VZEROLoader"); - if (!loader) { - Error("Exec", "Can not get VZERO Loader for input %d", iInput); - continue;} - - if (!runLoader->GetAliRun()) runLoader->LoadgAlice(); - - AliVZERO* vzero = (AliVZERO*) runLoader->GetAliRun()->GetDetector("VZERO"); - if (!vzero) { - Error("Exec", "No VZERO detector for input %d", iInput); - continue;} - - loader->LoadHits(); - TTree* treeH = loader->TreeH(); - if (!treeH) { - Error("Exec", "Cannot get TreeH for input %d", iInput); - continue; } - - for(Int_t i=0; i<80; i++) {map[i] = 0; time[i] = 0.0;} - - TClonesArray* hits = vzero->Hits(); - -// Now makes Digits from hits - - Int_t nTracks = (Int_t) treeH->GetEntries(); - for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) { - for (Int_t i=0; i<80; i++) {time_ref[i] = 999999.0;} - vzero->ResetHits(); - treeH->GetEvent(iTrack); - Int_t nHits = hits->GetEntriesFast(); - for (Int_t iHit = 0; iHit < nHits; iHit++) { - AliVZEROhit* hit = (AliVZEROhit *)hits->UncheckedAt(iHit); - Int_t nPhot = hit->Nphot(); - Int_t cell = hit->Cell(); - map[cell] += nPhot; - Float_t dt_scintillator = gRandom->Gaus(0,0.7); - Float_t t = dt_scintillator + 1e9*hit->Tof(); - if (t > 0.0) { - if(t < time_ref[cell]) time_ref[cell] = t; - time[cell] = TMath::Min(t,time_ref[cell]); - } - } // hit loop - } // track loop + fNdigits = 0; - loader->UnloadHits(); - - } // input loop - -// Now builds the scintillator cell response (80 cells i.e. 80 responses) - - for (Int_t i=0; i<80; i++) { - Float_t q1 = Float_t ( map[i] )* cPM * kQe; - Float_t noise = gRandom->Gaus(10.5,3.22); - Float_t pmResponse = q1/kC*TMath::Power(ktheta/kthau,1/(1-ktheta/kthau)) - + noise*1e-3; - map[i] = Int_t( pmResponse * adc_gain[i]); - } - -// Now transforms 80 cell responses into 64 photomultiplier responses - - for (Int_t j=0; j<32; j++){ - adc[j] = map [j]; - time2[j]= time[j];} - - for (Int_t j=48; j<80; j++){ - adc[j-16] = map [j]; - time2[j-16]= time[j];} - - for (Int_t j=0; j<16; j++){ - adc[16+j] = map [16+2*j]+ map [16+2*j+1]; - Float_t min_time = TMath::Min(time [16+2*j],time [16+2*j+1]); - time2[16+j] = min_time; - if(min_time==0.0) { - time2[16+j] = TMath::Max(time[16+2*j],time[16+2*j+1]); - } - } - - -// Now add digits to the digit Tree - - for (Int_t i=0; i<64; i++) { - if(adc[i] > 0) { -// printf(" Event, cell, adc, tof = %d %d %d %f\n", -// outRunLoader->GetEventNumber(),i, map[i], time2[i]*10.0); -// multiply by 10 to have 100 ps per channel : - AddDigit(i, adc[i], Int_t(time2[i]*10.0) ); - } - } - - treeD->Fill(); - outLoader->WriteDigits("OVERWRITE"); - outLoader->UnloadDigits(); - ResetDigit(); + if (fVZERO && !fManager) { + AliLoader *loader = fVZERO->GetLoader(); + if (!loader) { + AliError("Can not get VZERO Loader via AliVZERO object!"); + return; + } + AliRunLoader* runLoader = AliRunLoader::Instance(); + for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); ++iEvent) { + runLoader->GetEvent(iEvent); + if (fTask == kHits2Digits) { + DigitizeHits(); + DigitizeSDigits(); + WriteDigits(loader); + } + else { + DigitizeHits(); + WriteSDigits(loader); + } + } + } + else if (fManager) { + ReadSDigits(); + DigitizeSDigits(); + AliRunLoader *currentLoader = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName()); + AliLoader *loader = currentLoader->GetLoader("VZEROLoader"); + if (!loader) { + AliError("Cannot get VZERO Loader via RunDigitizer!"); + return; + } + WriteDigits(loader); + } + else { + AliFatal("Invalid digitization task! Exiting!"); + } } //____________________________________________________________________________ -void AliVZERODigitizer::AddDigit(Int_t PMnumber, Int_t adc, Int_t time) +void AliVZERODigitizer::AddDigit(Int_t pmnumber, Float_t time, Float_t width, Bool_t integrator, Short_t *chargeADC, Int_t *labels) { // Adds Digit TClonesArray &ldigits = *fDigits; - new(ldigits[fNdigits++]) AliVZEROdigit(PMnumber,adc,time); + + new(ldigits[fNdigits++]) AliVZEROdigit(pmnumber,time,width,integrator,chargeADC,labels); + +} +//____________________________________________________________________________ +void AliVZERODigitizer::AddSDigit(Int_t pmnumber, Int_t nbins, Float_t *charges, Int_t *labels) + { + +// Adds SDigit + + TClonesArray &ldigits = *fDigits; + + new(ldigits[fNdigits++]) AliVZEROSDigit(pmnumber,nbins,charges,labels); + } //____________________________________________________________________________ -void AliVZERODigitizer::ResetDigit() +void AliVZERODigitizer::ResetDigits() { -// + // Clears Digits -// + fNdigits = 0; - if (fDigits) fDigits->Delete(); + if (fDigits) fDigits->Clear(); } //____________________________________________________________________________ @@ -294,7 +306,7 @@ AliVZEROCalibData* AliVZERODigitizer::GetCalibData() const // AliWarning("Load of calibration data from default storage failed!"); // AliWarning("Calibration data will be loaded from local storage ($ALICE_ROOT)"); // Int_t runNumber = man->GetRun(); -// entry = man->GetStorage("local://$ALICE_ROOT") +// entry = man->GetStorage("local://$ALICE_ROOT/OCDB") // ->Get("VZERO/Calib/Data",runNumber); // // } @@ -311,3 +323,330 @@ AliVZEROCalibData* AliVZERODigitizer::GetCalibData() const } +double AliVZERODigitizer::SignalShape(double *x, double *par) +{ + // this function simulates the time + // of arrival of the photons at the + // photocathode + Double_t xx = x[0]; + if (xx <= par[0]) return 0; + Double_t a = 1./TMath::Power((xx-par[0])/par[1],1./par[2]); + if (xx <= par[3]) return a; + Double_t b = 1./TMath::Power((xx-par[3])/par[4],1./par[5]); + Double_t f = a*b/(a+b); + AliDebug(100,Form("x=%f func=%f",xx,f)); + return f; +} + +double AliVZERODigitizer::PMResponse(double *x, double * /* par */) +{ + // this function describes the + // PM time response to a single + // photoelectron + Double_t xx = x[0]+kPMRespTime; + return xx*xx*TMath::Exp(-xx*xx/(kPMRespTime*kPMRespTime)); +} + +double AliVZERODigitizer::SinglePhESpectrum(double *x, double * /* par */) +{ + // this function describes the + // PM amplitude response to a single + // photoelectron + Double_t xx = x[0]; + if (xx < 0) return 0; + return (TMath::Poisson(xx,kPMNbOfSecElec)+kPMTransparency*TMath::Poisson(xx,1.0)); +} + +Int_t AliVZERODigitizer::Cell2Pmt(Int_t cell) const +{ + // The method maps the scintillator + // indexes to the PM ones + if (cell < 0 || cell >= 80) { + AliError(Form("Wrong VZERO cell index %d",cell)); + return -1; + } + if (cell < 16) return cell; + if (cell < 48) return 8 + cell/2; + return cell - 16; +} + +void AliVZERODigitizer::DigitizeHits() +{ + // Digitize the hits to the level of + // SDigits (fTime arrays) + + for(Int_t i = 0 ; i < 64; ++i) { + memset(fTime[i],0,fNBins[i]*sizeof(Float_t)); + fLabels[i][0] = fLabels[i][1] = fLabels[i][2] = -1; + } + Float_t integral = fPMResponse->Integral(-kPMRespTime,2.*kPMRespTime); + Float_t meansPhE = fSinglePhESpectrum->Mean(0,20); + + AliLoader* loader = fVZERO->GetLoader(); + if (!loader) { + AliError("Can not get VZERO Loader!"); + return; + } + loader->LoadHits(); + TTree* treeH = loader->TreeH(); + if (!treeH) { + AliError("Cannot get TreeH!"); + return; + } + TClonesArray* hits = fVZERO->Hits(); + +// Now makes Digits from hits + Int_t nTracks = (Int_t) treeH->GetEntries(); + for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) { + fVZERO->ResetHits(); + treeH->GetEvent(iTrack); + Int_t nHits = hits->GetEntriesFast(); + for (Int_t iHit = 0; iHit < nHits; iHit++) { + AliVZEROhit* hit = (AliVZEROhit *)hits->UncheckedAt(iHit); + Int_t nPhot = hit->Nphot(); + Int_t cell = hit->Cell(); + Int_t pmt = Cell2Pmt(cell); + if (pmt < 0) continue; + Int_t trackLabel = hit->GetTrack(); + for(Int_t l = 0; l < 3; ++l) { + if (fLabels[pmt][l] < 0) { + fLabels[pmt][l] = trackLabel; + break; + } + } + Float_t dt_scintillator = gRandom->Gaus(0,kIntTimeRes); + Float_t t = dt_scintillator + 1e9*hit->Tof(); + if (pmt < 32) t += kV0CDelayCables; + t += fHptdcOffset[pmt]; + Int_t nPhE; + Float_t prob = fCalibData->GetLightYields(pmt)*fPhotoCathodeEfficiency; // Optical losses included! + if (nPhot > 100) + nPhE = (Int_t)gRandom->Gaus(prob*Float_t(nPhot)+0.5, + sqrt(Float_t(nPhot)*prob*(1.-prob))); + else + nPhE = gRandom->Binomial(nPhot,prob); + Float_t charge = TMath::Qe()*fPmGain[pmt]*fBinSize[pmt]/integral; + for (Int_t iPhE = 0; iPhE < nPhE; ++iPhE) { + Float_t tPhE = t + fSignalShape->GetRandom(0,fBinSize[pmt]*Float_t(fNBins[pmt])); + Float_t gainVar = fSinglePhESpectrum->GetRandom(0,20)/meansPhE; + Int_t firstBin = TMath::Max(0,(Int_t)((tPhE-kPMRespTime)/fBinSize[pmt])); + Int_t lastBin = TMath::Min(fNBins[pmt]-1,(Int_t)((tPhE+2.*kPMRespTime)/fBinSize[pmt])); + for(Int_t iBin = firstBin; iBin <= lastBin; ++iBin) { + Float_t tempT = fBinSize[pmt]*(0.5+iBin)-tPhE; + fTime[pmt][iBin] += gainVar*charge*fPMResponse->Eval(tempT); + } + } // ph.e. loop + } // hit loop + } // track loop + loader->UnloadHits(); +} + + +void AliVZERODigitizer::DigitizeSDigits() +{ + // Digitize the fTime arrays (SDigits) to the level of + // Digits (fAdc arrays) + for(Int_t i = 0 ; i < 64; ++i) { + for(Int_t j = 0; j < kNClocks; ++j) fAdc[i][j] = 0; + fLeadingTime[i] = fTimeWidth[i] = 0; + } + + Float_t maximum = 0.9*fSignalShape->GetMaximum(0,200); // Not exact, one needs to do this on the convoluted + Float_t integral2 = fSignalShape->Integral(0,200); // function. Anyway the effect is small <10% on the 2.5 ADC thr + for (Int_t ipmt = 0; ipmt < 64; ++ipmt) { + Float_t thr = fCalibData->GetCalibDiscriThr(ipmt,kFALSE)*kChargePerADC*maximum*fBinSize[ipmt]/integral2; + Bool_t ltFound = kFALSE, ttFound = kFALSE; + for (Int_t iBin = 0; iBin < fNBins[ipmt]; ++iBin) { + Float_t t = fBinSize[ipmt]*Float_t(iBin); + if (fTime[ipmt][iBin] > thr) { + if (!ltFound && (iBin < fNBinsLT[ipmt])) { + ltFound = kTRUE; + fLeadingTime[ipmt] = t; + } + } + else { + if (ltFound) { + if (!ttFound) { + ttFound = kTRUE; + fTimeWidth[ipmt] = t - fLeadingTime[ipmt]; + } + } + } + Float_t tadc = t - fClockOffset[ipmt]; + Int_t clock = kNClocks/2 - Int_t(tadc/25.0); + if (clock >= 0 && clock < kNClocks) + fAdc[ipmt][clock] += fTime[ipmt][iBin]/kChargePerADC; + } + AliDebug(1,Form("Channel %d Offset %f Time %f",ipmt,fClockOffset[ipmt],fLeadingTime[ipmt])); + Int_t board = AliVZEROCalibData::GetBoardNumber(ipmt); + if (ltFound && ttFound) { + fTimeWidth[ipmt] = fCalibData->GetWidthResolution(board)* + Float_t(Int_t(fTimeWidth[ipmt]/fCalibData->GetWidthResolution(board))); + if (fTimeWidth[ipmt] < Float_t(kMinTDCWidth)*fCalibData->GetWidthResolution(board)) + fTimeWidth[ipmt] = Float_t(kMinTDCWidth)*fCalibData->GetWidthResolution(board); + if (fTimeWidth[ipmt] > Float_t(kMaxTDCWidth)*fCalibData->GetWidthResolution(board)) + fTimeWidth[ipmt] = Float_t(kMaxTDCWidth)*fCalibData->GetWidthResolution(board); + } + } + + fEvenOrOdd = gRandom->Integer(2); + for (Int_t j=0; j<64; ++j){ + for (Int_t iClock = 0; iClock < kNClocks; ++iClock) { + Int_t integrator = (iClock + fEvenOrOdd) % 2; + AliDebug(1,Form("ADC %d %d %f",j,iClock,fAdc[j][iClock])); + fAdc[j][iClock] += gRandom->Gaus(fAdcPedestal[j][integrator], fAdcSigma[j][integrator]); + } + } + +} + +void AliVZERODigitizer::WriteDigits(AliLoader *loader) +{ + // Take fAdc arrays filled by the previous + // method and produce and add digits to the digit Tree + + loader->LoadDigits("UPDATE"); + + if (!loader->TreeD()) loader->MakeTree("D"); + loader->MakeDigitsContainer(); + TTree* treeD = loader->TreeD(); + DigitsArray(); + treeD->Branch("VZERODigit", &fDigits); + + Short_t *chargeADC = new Short_t[kNClocks]; + for (Int_t i=0; i<64; i++) { + for (Int_t j = 0; j < kNClocks; ++j) { + Int_t tempadc = Int_t(fAdc[i][j]); + if (tempadc > 1023) tempadc = 1023; + chargeADC[j] = tempadc; + } + AddDigit(i, fLeadingTime[i], fTimeWidth[i], Bool_t((10+fEvenOrOdd)%2), chargeADC, fLabels[i]); + } + delete [] chargeADC; + + treeD->Fill(); + loader->WriteDigits("OVERWRITE"); + loader->UnloadDigits(); + ResetDigits(); +} + +void AliVZERODigitizer::WriteSDigits(AliLoader *loader) +{ + // Take fTime arrays filled by the previous + // method and produce and add sdigits to the sdigit Tree + + loader->LoadSDigits("UPDATE"); + + if (!loader->TreeS()) loader->MakeTree("S"); + loader->MakeSDigitsContainer(); + TTree* treeS = loader->TreeS(); + SDigitsArray(); + treeS->Branch("VZEROSDigit", &fDigits); + + for (Int_t ipmt = 0; ipmt < 64; ++ipmt) { + AddSDigit(ipmt,fNBins[ipmt],fTime[ipmt],fLabels[ipmt]); + } + + treeS->Fill(); + loader->WriteSDigits("OVERWRITE"); + loader->UnloadSDigits(); + ResetDigits(); +} + +void AliVZERODigitizer::ReadSDigits() +{ + // Read SDigits which are then to precessed + // in the following method + for(Int_t i = 0 ; i < 64; ++i) { + memset(fTime[i],0,fNBins[i]*sizeof(Float_t)); + fLabels[i][0] = fLabels[i][1] = fLabels[i][2] = -1; + } + + // Loop over input files + Int_t nFiles= fManager->GetNinputs(); + for (Int_t inputFile = 0; inputFile < nFiles; inputFile++) { + // Get the current loader + AliRunLoader* currentLoader = + AliRunLoader::GetRunLoader(fManager->GetInputFolderName(inputFile)); + + AliLoader *loader = currentLoader->GetLoader("VZEROLoader"); + loader->LoadSDigits("READ"); + + // Get the tree of summable digits + TTree* sdigitsTree = loader->TreeS(); + if (!sdigitsTree) { + AliError("No sdigit tree from manager"); + continue; + } + + // Get the branch + TBranch* sdigitsBranch = sdigitsTree->GetBranch("VZEROSDigit"); + if (!sdigitsBranch) { + AliError("Failed to get sdigit branch"); + return; + } + + // Set the branch address + TClonesArray *sdigitsArray = NULL; + sdigitsBranch->SetAddress(&sdigitsArray); + + // Sum contributions from the sdigits + // Get number of entries in the tree + Int_t nentries = Int_t(sdigitsBranch->GetEntries()); + for (Int_t entry = 0; entry < nentries; ++entry) { + sdigitsBranch->GetEntry(entry); + // Get the number of sdigits + Int_t nsdigits = sdigitsArray->GetEntries(); + for (Int_t sdigit = 0; sdigit < nsdigits; sdigit++) { + AliVZEROSDigit* sDigit = static_cast(sdigitsArray->UncheckedAt(sdigit)); + Int_t pmNumber = sDigit->PMNumber(); + Int_t nbins = sDigit->GetNBins(); + if (nbins != fNBins[pmNumber]) { + AliError(Form("Incompatible number of bins between digitizer (%d) and sdigit (%d) for PM %d! Skipping sdigit!", + fNBins[pmNumber],nbins,pmNumber)); + continue; + } + // Sum the charges + Float_t *charges = sDigit->GetCharges(); + for(Int_t iBin = 0; iBin < nbins; ++iBin) fTime[pmNumber][iBin] += charges[iBin]; + // and the labels + Int_t *labels = sDigit->GetTracks(); + Int_t j = 0; + for(Int_t i = 0; i < 3; ++i) { + if (fLabels[pmNumber][i] < 0) { + if (labels[j] < 0) break; + fLabels[pmNumber][i] = labels[j]; + j++; + } + } + } + } + loader->UnloadSDigits(); + } +} + +//____________________________________________________________________ +TClonesArray* +AliVZERODigitizer::DigitsArray() +{ + // Initialize digit array if not already and + // return pointer to it. + if (!fDigits) { + fDigits = new TClonesArray("AliVZEROdigit", 64); + fNdigits = 0; + } + return fDigits; +} + +//____________________________________________________________________ +TClonesArray* +AliVZERODigitizer::SDigitsArray() +{ + // Initialize sdigit array if not already and + // return pointer to it. + if (!fDigits) { + fDigits = new TClonesArray("AliVZEROSDigit", 64); + fNdigits = 0; + } + return fDigits; +}