#include <TClonesArray.h>
#include <TTree.h>
#include <TFile.h>
#include <TSystem.h>
#include <TRandom.h>
#include <TParticle.h>
#include "AliDetector.h"
#include "AliRawDataHeaderSim.h"
#include "AliRawReader.h"
#include "AliLoader.h"
#include "AliRun.h"
#include "AliMC.h"
#include "AliLog.h"
#include "AliDAQ.h"
#include "AliZDC.h"
#include "AliZDCHit.h"
#include "AliZDCSDigit.h"
#include "AliZDCDigit.h"
#include "AliZDCDigitizer.h"
#include "AliZDCRawStream.h"
#include "AliZDCPedestals.h"
#include "AliZDCEnCalib.h"
#include "AliZDCTowerCalib.h"
#include "AliFstream.h"
ClassImp(AliZDC)
AliZDC::AliZDC() :
AliDetector(),
fNoShower(0),
fPedCalib(0),
fEnCalibData(0),
fTowCalibData(0),
fZDCCalibFName(""),
fSpectatorTracked(1),
fBeamEnergy(0.),
fIspASystem(kFALSE),
fIsRELDISgen(kFALSE),
fOnlyZEM(kFALSE),
fFindMother(kFALSE)
{
fIshunt = 1;
fNhits = 0;
fHits = 0;
fDigits = 0;
fNdigits = 0;
}
AliZDC::AliZDC(const char *name, const char *title) :
AliDetector(name,title),
fNoShower(0),
fPedCalib(0),
fEnCalibData(0),
fTowCalibData(0),
fZDCCalibFName(""),
fSpectatorTracked(1),
fBeamEnergy(0.),
fIspASystem(kFALSE),
fIsRELDISgen(kFALSE),
fOnlyZEM(kFALSE),
fFindMother(kFALSE)
{
fIshunt = 1;
fNhits = 0;
fDigits = 0;
fNdigits = 0;
fHits = new TClonesArray("AliZDCHit",1000);
gAlice->GetMCApp()->AddHitList(fHits);
SetName("ZDC"); SetTitle("ZDC");
}
AliZDC::~AliZDC()
{
fIshunt = 0;
if(fPedCalib) delete fPedCalib;
if(fEnCalibData) delete fEnCalibData;
if(fEnCalibData) delete fEnCalibData;
}
AliZDC::AliZDC(const AliZDC& ZDC) :
AliDetector("ZDC","ZDC"),
fNoShower(ZDC.fNoShower),
fPedCalib(ZDC.fPedCalib),
fEnCalibData(ZDC.fEnCalibData),
fTowCalibData(ZDC.fTowCalibData),
fZDCCalibFName(ZDC.fZDCCalibFName),
fSpectatorTracked(ZDC.fSpectatorTracked),
fBeamEnergy(ZDC.fBeamEnergy),
fIspASystem(ZDC.fIspASystem),
fIsRELDISgen(ZDC.fIsRELDISgen),
fOnlyZEM(ZDC.fOnlyZEM),
fFindMother(ZDC.fFindMother)
{
}
AliZDC& AliZDC::operator=(const AliZDC& ZDC)
{
if(this!=&ZDC){
fNoShower = ZDC.fNoShower;
fPedCalib = ZDC.fPedCalib;
fEnCalibData = ZDC.fEnCalibData;
fTowCalibData = ZDC.fTowCalibData;
fZDCCalibFName = ZDC.fZDCCalibFName;
fBeamEnergy = ZDC.fBeamEnergy;
fIspASystem = ZDC.fIspASystem;
fIsRELDISgen = ZDC.fIsRELDISgen;
fOnlyZEM = ZDC.fOnlyZEM;
fFindMother = ZDC.fFindMother;
} return *this;
}
void AliZDC::AddHit(Int_t track, Int_t *vol, Float_t *hits)
{
static Float_t trackTime=0., trackEta=0., primKinEn=0., xImpact=0., yImpact=0., sFlag=0.;
static Int_t pcPDGcode=0, motPDGcode=0;
AliZDCHit *newquad, *curprimquad;
newquad = new AliZDCHit(fIshunt, track, vol, hits);
TClonesArray &lhits = *fHits;
if(fNhits==0){
sFlag = newquad->GetSFlag();
primKinEn = newquad->GetPrimKinEn();
xImpact = newquad->GetXImpact();
yImpact = newquad->GetYImpact();
pcPDGcode = newquad->GetPDGCode();
motPDGcode = newquad->GetMotherPDGCode();
trackTime = newquad->GetTrackTOF();
trackEta = newquad->GetTrackEta();
}
else{
newquad->SetPrimKinEn(primKinEn);
newquad->SetXImpact(xImpact);
newquad->SetYImpact(yImpact);
newquad->SetSFlag(sFlag);
newquad->SetPDGCode(pcPDGcode);
newquad->SetMotherPDGCode(motPDGcode);
newquad->SetTrackTOF(trackTime);
newquad->SetTrackEta(trackEta);
}
Int_t j;
for(j=0; j<fNhits; j++){
curprimquad = (AliZDCHit*) lhits[j];
if(*curprimquad == *newquad){
*curprimquad = *curprimquad+*newquad;
delete newquad;
return;
}
}
new(lhits[fNhits]) AliZDCHit(*newquad);
fNhits++;
delete newquad;
}
Float_t AliZDC::ZMin(void) const
{
return -11600.;
}
Float_t AliZDC::ZMax(void) const
{
return 11750.;
}
void AliZDC::MakeBranch(Option_t *opt)
{
char branchname[10];
snprintf(branchname, 10, "%s", GetName());
const char *cH = strstr(opt,"H");
if(cH && fLoader->TreeH()) {
if (fHits) {
fHits->Clear();
fNhits = 0;
}
else {
fHits = new TClonesArray("AliZDCHit",1000);
if (gAlice && gAlice->GetMCApp())
gAlice->GetMCApp()->AddHitList(fHits);
}
}
AliDetector::MakeBranch(opt);
}
void AliZDC::Hits2SDigits()
{
AliDebug(1,"\n AliZDC::Hits2SDigits() ");
fLoader->LoadHits("read");
fLoader->LoadSDigits("recreate");
AliRunLoader* runLoader = fLoader->GetRunLoader();
AliZDCSDigit sdigit;
AliZDCSDigit* psdigit = &sdigit;
for(Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
Float_t pmZNC[5], pmZPC[5], pmZNA[5], pmZPA[5], pmZEM1=0., pmZEM2=0.;
for(Int_t i=0; i<5; i++) pmZNC[i] = pmZPC[i] = pmZNA[i] = pmZPA[i] = 0;
runLoader->GetEvent(iEvent);
TTree* treeH = fLoader->TreeH();
Int_t ntracks = (Int_t) treeH->GetEntries();
ResetHits();
Int_t sector[2]; Float_t trackTime = 0.;
for(Int_t itrack = 0; itrack < ntracks; itrack++) {
treeH->GetEntry(itrack);
for(AliZDCHit* zdcHit = (AliZDCHit*)FirstHit(-1); zdcHit;
zdcHit = (AliZDCHit*)NextHit()) {
sector[0] = zdcHit->GetVolume(0);
sector[1] = zdcHit->GetVolume(1);
if((sector[1] < 1) || (sector[1]>5)) {
Error("Hits2SDigits", "sector[0] = %d, sector[1] = %d", sector[0], sector[1]);
continue;
}
Float_t lightQ = zdcHit->GetLightPMQ();
Float_t lightC = zdcHit->GetLightPMC();
trackTime = zdcHit->GetTrackTOF();
if(sector[0] == 3) trackTime += 320;
if(sector[0] == 1) {
pmZNC[0] += lightC;
pmZNC[sector[1]] += lightQ;
}
else if(sector[0] == 2) {
pmZPC[0] += lightC;
pmZPC[sector[1]] += lightQ;
}
else if(sector[0] == 3) {
if(sector[1] == 1) pmZEM1 += lightC;
else pmZEM2 += lightQ;
}
if(sector[0] == 4) {
pmZNA[0] += lightC;
pmZNA[sector[1]] += lightQ;
}
else if(sector[0] == 5) {
pmZPA[0] += lightC;
pmZPA[sector[1]] += lightQ;
}
}
}
fLoader->MakeTree("S");
TTree* treeS = fLoader->TreeS();
const Int_t kBufferSize = 4000;
treeS->Branch(GetName(), "AliZDCSDigit", &psdigit, kBufferSize);
sector[0] = 1;
for(Int_t j = 0; j < 5; j++) {
sector[1] = j;
if(pmZNC[j]>0){
new(psdigit) AliZDCSDigit(sector, pmZNC[j], trackTime);
treeS->Fill();
}
}
sector[0] = 2;
for(Int_t j = 0; j < 5; j++) {
sector[1] = j;
if(pmZPC[j]>0){
new(psdigit) AliZDCSDigit(sector, pmZPC[j], trackTime);
treeS->Fill();
}
}
sector[0] = 3;
sector[1] = 1;
if(pmZEM1>0){
new(psdigit) AliZDCSDigit(sector, pmZEM1, trackTime);
treeS->Fill();
}
sector[1] = 2;
if(pmZEM2>0){
new(psdigit) AliZDCSDigit(sector, pmZEM2, trackTime);
treeS->Fill();
}
sector[0] = 4;
for(Int_t j = 0; j < 5; j++) {
sector[1] = j;
if(pmZNA[j]>0){
new(psdigit) AliZDCSDigit(sector, pmZNA[j], trackTime);
treeS->Fill();
}
}
sector[0] = 5;
sector[1] = 0;
for(Int_t j = 0; j < 5; j++) {
sector[1] = j;
if(pmZPA[j]>0){
new(psdigit) AliZDCSDigit(sector, pmZPA[j], trackTime);
treeS->Fill();
}
}
fLoader->WriteSDigits("OVERWRITE");
}
fLoader->UnloadHits();
fLoader->UnloadSDigits();
}
AliDigitizer* AliZDC::CreateDigitizer(AliDigitizationInput* digInput) const{
AliZDCDigitizer *zdcDigitizer = new AliZDCDigitizer(digInput);
if(fSpectatorTracked==0) zdcDigitizer->SetSpectators2Track();
if(fBeamEnergy>0.01) zdcDigitizer->SetBeamEnergy(fBeamEnergy);
if(fIspASystem==kTRUE) zdcDigitizer->SetpAsystem();
if(fIsRELDISgen==kTRUE) zdcDigitizer->SetRELDISGenerator();
return zdcDigitizer;
}
void AliZDC::Digits2Raw()
{
const int knADCData1=12, knADCData2=12;
const int knADCData3=12, knADCData4=12;
UInt_t lADCHeader1;
UInt_t lADCHeader2;
UInt_t lADCHeader3;
UInt_t lADCHeader4;
UInt_t lADCData1[2*knADCData1];
UInt_t lADCData2[2*knADCData2];
UInt_t lADCData3[2*knADCData3];
UInt_t lADCData4[2*knADCData4];
UInt_t lADCEndBlock;
fLoader->LoadDigits("read");
AliZDCDigit digit;
AliZDCDigit* pdigit = &digit;
TTree* treeD = fLoader->TreeD();
if(!treeD) return;
treeD->SetBranchAddress("ZDC", &pdigit);
AliZDCChMap * chMap = GetChMap();
const int nCh = knADCData1+knADCData2+knADCData3+knADCData4;
Int_t mapADC[nCh][4];
for(Int_t i=0; i<nCh; i++){
mapADC[i][0] = chMap->GetADCModule(i);
mapADC[i][1] = chMap->GetADCChannel(i);
mapADC[i][2] = chMap->GetDetector(i);
mapADC[i][3] = chMap->GetSector(i);
}
UInt_t lADCHeaderGEO1 = 0;
UInt_t lADCHeaderGEO2 = 1;
UInt_t lADCHeaderGEO3 = 2;
UInt_t lADCHeaderGEO4 = 3;
UInt_t lADCHeaderCRATE = 0;
UInt_t lADCHeaderCNT1 = knADCData1;
UInt_t lADCHeaderCNT2 = knADCData2;
UInt_t lADCHeaderCNT3 = knADCData3;
UInt_t lADCHeaderCNT4 = knADCData4;
lADCHeader1 = lADCHeaderGEO1 << 27 | 0x1 << 25 | lADCHeaderCRATE << 16 |
lADCHeaderCNT1 << 8 ;
lADCHeader2 = lADCHeaderGEO2 << 27 | 0x1 << 25 | lADCHeaderCRATE << 16 |
lADCHeaderCNT2 << 8 ;
lADCHeader3 = lADCHeaderGEO3 << 27 | 0x1 << 25 | lADCHeaderCRATE << 16 |
lADCHeaderCNT3 << 8 ;
lADCHeader4 = lADCHeaderGEO4 << 27 | 0x1 << 25 | lADCHeaderCRATE << 16 |
lADCHeaderCNT4 << 8 ;
UInt_t lADCDataGEO = 0;
UInt_t lADCDataValue1[2*knADCData1];
UInt_t lADCDataValue2[2*knADCData2];
UInt_t lADCDataValue3[2*knADCData3];
UInt_t lADCDataValue4[2*knADCData4];
UInt_t lADCDataOvFlwHG = 0;
UInt_t lADCDataOvFlwLG = 0;
for(Int_t i=0; i<2*knADCData1 ; i++) lADCDataValue1[i] = 0;
for(Int_t i=0; i<2*knADCData2 ; i++) lADCDataValue2[i] = 0;
for(Int_t i=0; i<2*knADCData3 ; i++) lADCDataValue3[i] = 0;
for(Int_t i=0; i<2*knADCData4 ; i++) lADCDataValue4[i] = 0;
UInt_t lADCDataChannel = 0;
Int_t indADC0=0, indADC1=0, indADC2=0, indADC3=0;
for(Int_t iDigit=0; iDigit<(Int_t) (treeD->GetEntries()); iDigit++){
treeD->GetEntry(iDigit);
if(!pdigit) continue;
for(Int_t k=0; k<nCh; k++){
if(iDigit<knADCData1+knADCData2){
if(digit.GetSector(0)==mapADC[k][2] && digit.GetSector(1)==mapADC[k][3]){
lADCDataGEO = (UInt_t) mapADC[k][0];
lADCDataChannel = (UInt_t) mapADC[k][1];
break;
}
}
else{
if(digit.GetSector(0)==mapADC[k][2] && digit.GetSector(1)==mapADC[k][3]){
lADCDataGEO = (UInt_t) mapADC[k][0];
lADCDataChannel = (UInt_t) mapADC[k][1];
if(k>knADCData1+knADCData2) break;
}
}
}
if(lADCDataGEO==0){
if(indADC0>=knADCData1){
AliWarning(Form(" Problem with digit index %d for ADC0\n", indADC0));
return;
}
Int_t indLG = indADC0+knADCData1;
if(digit.GetADCValue(0) > 2047) lADCDataOvFlwHG = 1;
lADCDataValue1[indADC0] = digit.GetADCValue(0);
lADCData1[indADC0] = lADCDataGEO << 27 | lADCDataChannel << 17 |
lADCDataOvFlwHG << 12 | (lADCDataValue1[indADC0] & 0xfff);
if(digit.GetADCValue(1) > 2047) lADCDataOvFlwLG = 1;
lADCDataValue1[indLG] = digit.GetADCValue(1);
lADCData1[indLG] = lADCDataGEO << 27 | lADCDataChannel << 17 | 0x1 << 16 |
lADCDataOvFlwLG << 12 | (lADCDataValue1[indLG] & 0xfff);
indADC0++;
}
else if(lADCDataGEO==1){
if(indADC1>=knADCData2){
AliWarning(Form(" Problem with digit index %d for ADC1\n", indADC1));
return;
}
Int_t indLG = indADC1+knADCData2;
if(digit.GetADCValue(0) > 2047) lADCDataOvFlwHG = 1;
lADCDataValue2[indADC1] = digit.GetADCValue(0);
lADCData2[indADC1] = lADCDataGEO << 27 | lADCDataChannel << 17 |
lADCDataOvFlwHG << 12 | (lADCDataValue2[indADC1] & 0xfff);
if(digit.GetADCValue(1) > 2047) lADCDataOvFlwLG = 1;
lADCDataValue2[indLG] = digit.GetADCValue(1);
lADCData2[indLG] = lADCDataGEO << 27 | lADCDataChannel << 17 | 0x1 << 16 |
lADCDataOvFlwLG << 12 | (lADCDataValue2[indLG] & 0xfff);
indADC1++;
}
else if(lADCDataGEO==2){
if(indADC2>=knADCData3){
AliWarning(Form(" Problem with digit index %d for ADC2\n", indADC2));
return;
}
Int_t indLG = indADC2+knADCData3;
if(digit.GetADCValue(0) > 2047) lADCDataOvFlwHG = 1;
lADCDataValue3[indADC1] = digit.GetADCValue(0);
lADCData3[indADC1] = lADCDataGEO << 27 | lADCDataChannel << 17 |
lADCDataOvFlwHG << 12 | (lADCDataValue3[indADC2] & 0xfff);
if(digit.GetADCValue(1) > 2047) lADCDataOvFlwLG = 1;
lADCDataValue3[indLG] = digit.GetADCValue(1);
lADCData3[indLG] = lADCDataGEO << 27 | lADCDataChannel << 17 | 0x1 << 16 |
lADCDataOvFlwLG << 12 | (lADCDataValue3[indLG] & 0xfff);
indADC2++;
}
else if(lADCDataGEO==3){
if(indADC3>=knADCData4){
AliWarning(Form(" Problem with digit index %d for ADC2\n", indADC3));
return;
}
Int_t indLG = indADC3+knADCData4;
if(digit.GetADCValue(0) > 2047) lADCDataOvFlwHG = 1;
lADCDataValue4[indADC3] = digit.GetADCValue(0);
lADCData4[indADC3] = lADCDataGEO << 27 | lADCDataChannel << 17 |
lADCDataOvFlwHG << 12 | (lADCDataValue4[indADC3] & 0xfff);
if(digit.GetADCValue(1) > 2047) lADCDataOvFlwLG = 1;
lADCDataValue4[indLG] = digit.GetADCValue(1);
lADCData4[indLG] = lADCDataGEO << 27 | lADCDataChannel << 17 | 0x1 << 16 |
lADCDataOvFlwLG << 12 | (lADCDataValue4[indLG] & 0xfff);
indADC3++;
}
}
UInt_t lADCEndBlockGEO = 0;
AliRunLoader* runLoader = fLoader->GetRunLoader();
UInt_t lADCEndBlockEvCount = runLoader->GetEventNumber();
lADCEndBlock = lADCEndBlockGEO << 27 | 0x1 << 26 | lADCEndBlockEvCount;
TString fileName;
fileName.Form("%s",AliDAQ::DdlFileName("ZDC",0));
AliFstream* file = new AliFstream(fileName.Data());
AliRawDataHeaderSim header;
header.fSize = sizeof(header) +
sizeof(lADCHeader1) + sizeof(lADCData1) + sizeof(lADCEndBlock) +
sizeof(lADCHeader2) + sizeof(lADCData2) + sizeof(lADCEndBlock) +
sizeof(lADCHeader3) + sizeof(lADCData3) + sizeof(lADCEndBlock) +
sizeof(lADCHeader4) + sizeof(lADCData4) + sizeof(lADCEndBlock);
header.SetAttribute(0);
file->WriteBuffer((char*)(&header), sizeof(header));
file->WriteBuffer((char*) &lADCHeader1, sizeof (lADCHeader1));
file->WriteBuffer((char*) &lADCData1, sizeof(lADCData1));
file->WriteBuffer((char*) &lADCEndBlock, sizeof(lADCEndBlock));
file->WriteBuffer((char*) &lADCHeader2, sizeof (lADCHeader2));
file->WriteBuffer((char*) (lADCData2), sizeof(lADCData2));
file->WriteBuffer((char*) &lADCEndBlock, sizeof(lADCEndBlock));
file->WriteBuffer((char*) &lADCHeader3, sizeof (lADCHeader3));
file->WriteBuffer((char*) (lADCData3), sizeof(lADCData3));
file->WriteBuffer((char*) &lADCEndBlock, sizeof(lADCEndBlock));
file->WriteBuffer((char*) &lADCHeader4, sizeof (lADCHeader4));
file->WriteBuffer((char*) (lADCData4), sizeof(lADCData4));
file->WriteBuffer((char*) &lADCEndBlock, sizeof(lADCEndBlock));
delete file;
fLoader->UnloadDigits();
}
Bool_t AliZDC::Raw2SDigits(AliRawReader* rawReader)
{
const int kNch = 48;
AliLoader* loader = (AliRunLoader::Instance())->GetLoader("ZDCLoader");
if(!loader) {
AliError("no ZDC loader found");
return kFALSE;
}
TTree* treeS = loader->TreeS();
if(!treeS){
loader->MakeTree("S");
treeS = loader->TreeS();
}
AliZDCSDigit sdigit;
AliZDCSDigit* psdigit = &sdigit;
const Int_t kBufferSize = 4000;
treeS->Branch("ZDC", "AliZDCSDigit", &psdigit, kBufferSize);
AliZDCRawStream rawStream(rawReader);
Int_t sector[2], resADC, rawADC, corrADC, nPheVal;
Int_t jcount = 0;
while(rawStream.Next()){
if(rawStream.IsADCDataWord()){
if(jcount < kNch){
for(Int_t j=0; j<2; j++) sector[j] = rawStream.GetSector(j);
rawADC = rawStream.GetADCValue();
resADC = rawStream.GetADCGain();
corrADC = rawADC - Pedestal(sector[0], sector[1], resADC);
if(corrADC<0) corrADC=0;
nPheVal = ADCch2Phe(sector[0], sector[1], corrADC, resADC);
new(psdigit) AliZDCSDigit(sector, (Float_t) nPheVal, 0.);
treeS->Fill();
jcount++;
}
}
}
fLoader->WriteSDigits("OVERWRITE");
fLoader->UnloadSDigits();
return kTRUE;
}
Int_t AliZDC::Pedestal(Int_t Det, Int_t Quad, Int_t Res) const
{
Float_t pedValue = 0.;
AliCDBManager *man = AliCDBManager::Instance();
AliCDBEntry *entry = man->Get("ZDC/Calib/Pedestals");
AliZDCPedestals *calibPed = 0x0;
if(!entry) AliFatal("No calibration data loaded!");
else{
calibPed = (AliZDCPedestals*) entry->GetObject();
if(!calibPed){
printf("\t No calibration object found for ZDC!");
return -1;
}
Int_t index=0, kNch=24;
if(Quad!=5){
if(Det==1) index = Quad+kNch*Res;
else if(Det==2) index = Quad+5+kNch*Res;
else if(Det==3) index = Quad+9+kNch*Res;
else if(Det==4) index = Quad+12+kNch*Res;
else if(Det==5) index = Quad+17+kNch*Res;
}
else index = (Det-1)/3+22+kNch*Res;
Float_t meanPed = calibPed->GetMeanPed(index);
Float_t pedWidth = calibPed->GetMeanPedWidth(index);
pedValue = gRandom->Gaus(meanPed,pedWidth);
}
return (Int_t) pedValue;
}
Int_t AliZDC::ADCch2Phe(Int_t Det, Int_t Quad, Int_t ADCVal, Int_t Res) const
{
Float_t pmGain[6][5];
Float_t resADC[2];
for(Int_t j = 0; j < 5; j++){
pmGain[0][j] = 50000.;
pmGain[1][j] = 100000.;
pmGain[2][j] = 100000.;
pmGain[3][j] = 50000.;
pmGain[4][j] = 100000.;
pmGain[5][j] = 100000.;
}
resADC[0] = 0.0000008;
resADC[1] = 0.0000064;
Int_t nPhe = (Int_t) (ADCVal / (pmGain[Det-1][Quad] * resADC[Res]));
return nPhe;
}
void AliZDC::SetTreeAddress(){
if(fLoader->TreeH() && (fHits == 0x0))
fHits = new TClonesArray("AliZDCHit",1000);
AliDetector::SetTreeAddress();
}
AliZDCChMap* AliZDC::GetChMap() const
{
AliZDCChMap *calibdata = 0x0;
AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/ChMap");
if(!entry) AliFatal("No calibration data loaded!");
else{
calibdata = dynamic_cast<AliZDCChMap*> (entry->GetObject());
if(!calibdata) AliFatal("Wrong calibration object in calibration file!");
}
return calibdata;
}