X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=EMCAL%2FAliEMCAL.cxx;h=1e95f60c6b56adc636861ad056b9a2ba569398cd;hb=d49fe99ae1b943655ce70f572a557d9434e695fa;hp=aad20680a75e43a4efb399906f8dcdf3c662cbcd;hpb=1d7b3dd64ef1e53ab5b784e7e1606dc84907b3c3;p=u%2Fmrichter%2FAliRoot.git diff --git a/EMCAL/AliEMCAL.cxx b/EMCAL/AliEMCAL.cxx index aad20680a75..1e95f60c6b5 100644 --- a/EMCAL/AliEMCAL.cxx +++ b/EMCAL/AliEMCAL.cxx @@ -27,39 +27,42 @@ // --- ROOT system --- class TFile; -#include "TBranch.h" -#include "TClonesArray.h" -#include "TTree.h" +#include +#include +#include +#include +#include +#include // --- Standard library --- -#include // --- AliRoot header files --- +#include "AliMagF.h" #include "AliEMCAL.h" +#include "AliEMCALGetter.h" #include "AliRun.h" -#include "AliMagF.h" -#include "AliEMCALGeometry.h" -//#include "AliEMCALQAChecker.h" +#include "AliEMCALSDigitizer.h" +#include "AliEMCALDigitizer.h" +#include "AliAltroBuffer.h" ClassImp(AliEMCAL) //____________________________________________________________________________ AliEMCAL::AliEMCAL():AliDetector() { // Default ctor - fName="EMCAL"; - //fQATask = 0; - fTreeQA = 0; - fGeom = 0 ; + fName = "EMCAL" ; } //____________________________________________________________________________ AliEMCAL::AliEMCAL(const char* name, const char* title): AliDetector(name,title) { // ctor : title is used to identify the layout - - fGeom = AliEMCALGeometry::GetInstance(GetTitle(),"") ; - //fQATask = 0; - fTreeQA = 0; + + fTimeMax = 1.28E-5 ; + fTimePeak = 2.0E-6 ; + fTimeRes = 1.5E-6 ; + fHighGainFactor = 40; + fHighGainOffset = 0x200 ; } //____________________________________________________________________________ @@ -68,15 +71,30 @@ AliEMCAL::~AliEMCAL() } +//____________________________________________________________________________ +void AliEMCAL::Copy(AliEMCAL & emcal) +{ + TObject::Copy(emcal) ; +} + +//____________________________________________________________________________ +AliDigitizer* AliEMCAL::CreateDigitizer(AliRunDigitizer* manager) const +{ + return new AliEMCALDigitizer(manager); +} + //____________________________________________________________________________ void AliEMCAL::CreateMaterials() { // Definitions of materials to build EMCAL and associated tracking media. // media number in idtmed are 1599 to 1698. - // --- Air --- - AliMaterial(0, "Air$", 14.61, 7.3, 0.001205, 30420., 67500., 0, 0) ; - + // --- Air --- + Float_t aAir[4]={12.0107,14.0067,15.9994,39.948}; + Float_t zAir[4]={6.,7.,8.,18.}; + Float_t wAir[4]={0.000124,0.755267,0.231781,0.012827}; + Float_t dAir = 1.20479E-3; + AliMixture(0, "Air$", aAir, zAir, dAir, 4, wAir) ; // --- Lead --- AliMaterial(1, "Pb$", 207.2, 82, 11.35, 0.56, 0., 0, 0) ; @@ -94,38 +112,31 @@ void AliEMCAL::CreateMaterials() AliMaterial(3, "Al$", 26.98, 13., 2.7, 8.9, 999., 0, 0) ; // --- Absorption length is ignored ^ - - // DEFINITION OF THE TRACKING MEDIA // for EMCAL: idtmed[1599->1698] equivalent to fIdtmed[0->100] Int_t * idtmed = fIdtmed->GetArray() - 1599 ; Int_t isxfld = gAlice->Field()->Integ() ; Float_t sxmgmx = gAlice->Field()->Max() ; - - - // Air -> idtmed[1599] - AliMedium(0, "Air $", 0, 0, + // Air -> idtmed[1599] + AliMedium(0, "Air $", 0, 0, isxfld, sxmgmx, 10.0, 1.0, 0.1, 0.1, 10.0, 0, 0) ; - // The Lead -> idtmed[1600] + // The Lead -> idtmed[1600] AliMedium(1, "Lead $", 1, 0, isxfld, sxmgmx, 10.0, 0.1, 0.1, 0.1, 0.1, 0, 0) ; - // The scintillator of the CPV made of Polystyrene scintillator -> idtmed[1601] AliMedium(2, "CPV scint. $", 2, 1, isxfld, sxmgmx, 10.0, 0.001, 0.1, 0.001, 0.001, 0, 0) ; - // Various Aluminium parts made of Al -> idtmed[1602] + // Various Aluminium parts made of Al -> idtmed[1602] AliMedium(3, "Al parts $", 3, 0, isxfld, sxmgmx, 10.0, 0.1, 0.1, 0.001, 0.001, 0, 0) ; - - // --- Set decent energy thresholds for gamma and electron tracking // Tracking threshold for photons and electrons in Lead @@ -139,66 +150,213 @@ void AliEMCAL::CreateMaterials() gMC->Gstpar(idtmed[1600], "DCUTE",0.00001) ; gMC->Gstpar(idtmed[1600], "DCUTM",0.00001) ; -// --- and in aluminium parts --- +// --- in aluminium parts --- gMC->Gstpar(idtmed[1602], "LOSS",3.) ; gMC->Gstpar(idtmed[1602], "DRAY",1.) ; gMC->Gstpar(idtmed[1602], "DCUTE",0.00001) ; gMC->Gstpar(idtmed[1602], "DCUTM",0.00001) ; - // --- and finally thresholds for photons and electrons in the scintillator --- gMC->Gstpar(idtmed[1601],"CUTGAM",0.00008) ; gMC->Gstpar(idtmed[1601],"CUTELE",0.001) ; gMC->Gstpar(idtmed[1601],"BCUTE",0.0001) ; + //set constants for Birk's Law implentation + fBirkC0 = 1; + fBirkC1 = 0.013/dP; + fBirkC2 = 9.6e-6/(dP * dP); + } + +//____________________________________________________________________________ +void AliEMCAL::Digits2Raw() +{ +// convert digits of the current event to raw data + + // get the digits + AliEMCALGetter * gime = AliEMCALGetter::Instance(AliRunLoader::GetGAliceName()) ; + if (!gime) { + Error("Digits2Raw", "EMCAL Getter not instantiated") ; + return ; + } + gime->Event(gime->EventNumber(), "D") ; + TClonesArray* digits = gime->Digits() ; + + if (!digits) { + Error("Digits2Raw", "no digits found !"); + return; + } + + // get the geometry + AliEMCALGeometry* geom = gime->EMCALGeometry(); + if (!geom) { + Error("Digits2Raw", "no geometry found !"); + return; + } + + // some digitization constants + const Int_t kDDLOffset = 0x800; + const Int_t kThreshold = 3; + const Int_t kChannelsperDDL = 897 ; + + AliAltroBuffer* buffer = NULL; + Int_t prevDDL = -1; + Int_t adcValuesLow[fkTimeBins]; + Int_t adcValuesHigh[fkTimeBins]; + + // loop over digits (assume ordered digits) + for (Int_t iDigit = 0; iDigit < digits->GetEntries(); iDigit++) { + AliEMCALDigit* digit = gime->Digit(iDigit); + if (digit->GetAmp() < kThreshold) + continue; + Int_t iDDL = digit->GetId() / kChannelsperDDL ; + // for each DDL id is numbered from 1 to kChannelsperDDL -1 + Int_t idDDL = digit->GetId() - iDDL * ( kChannelsperDDL - 1 ) ; + // new DDL + if (iDDL != prevDDL) { + // write real header and close previous file + if (buffer) { + buffer->Flush(); + buffer->WriteDataHeader(kFALSE, kFALSE); + delete buffer; + } + + // open new file and write dummy header + TString fileName("EMCAL_") ; + fileName += (iDDL + kDDLOffset) ; + fileName += ".ddl" ; + buffer = new AliAltroBuffer(fileName.Data(), 1); + buffer->WriteDataHeader(kTRUE, kFALSE); //Dummy; + + prevDDL = iDDL; + } + + // out of time range signal (?) + if (digit->GetTimeR() > fTimeMax) { + buffer->FillBuffer(digit->GetAmp()); + buffer->FillBuffer(fkTimeBins); // time bin + buffer->FillBuffer(3); // bunch length + buffer->WriteTrailer(3, idDDL, 0, 0); // trailer + + // simulate linear rise and gaussian decay of signal + } else { + Bool_t highGain = kFALSE; + + // write low and eventually high gain channel + buffer->WriteChannel(iDDL, 0, 0, + fkTimeBins, adcValuesLow, kThreshold); + if (highGain) { + buffer->WriteChannel(iDDL, 0, fHighGainOffset, + fkTimeBins, adcValuesHigh, 1); + } + } + } + // write real header and close last file + if (buffer) { + buffer->Flush(); + buffer->WriteDataHeader(kFALSE, kFALSE); + delete buffer; + } + gime->EmcalLoader()->UnloadDigits(); +} + +//____________________________________________________________________________ +void AliEMCAL::Hits2SDigits() +{ +// create summable digits + + AliEMCALSDigitizer* emcalDigitizer = + new AliEMCALSDigitizer(fLoader->GetRunLoader()->GetFileName().Data()) ; + emcalDigitizer->SetEventRange(0, -1) ; // do all the events + emcalDigitizer->ExecuteTask() ; +} //____________________________________________________________________________ -AliEMCALGeometry * AliEMCAL::GetGeometry() const -{ - // gets the pointer to the AliEMCALGeometry unique instance - - if (fGeom) - return fGeom ; - else - return AliEMCALGeometry::GetInstance(GetTitle(),"") ; +AliLoader* AliEMCAL::MakeLoader(const char* topfoldername) +{ +//different behaviour than standard (singleton getter) +// --> to be discussed and made eventually coherent + fLoader = new AliEMCALLoader(GetName(),topfoldername); + return fLoader; +} + +//__________________________________________________________________ +Double_t AliEMCAL::RawResponseFunction(Double_t *x, Double_t *par) +{ + // Shape of the electronics raw reponse: + // 1. the signal rises linearly from par[4] to par[1] to reach the maximu par[3] + // 2. the signal decays with a gaussian shape for par[4]+par[1] with a sigma of par[2] + + Float_t xx = x[0] ; + + Double_t signal = 0. ; + + if (xx < par[4] + par[1]) // signal is rising + signal = (gRandom->Rndm() + par[3]) * (xx - par[4]) / (par[1] - par[4]) ; + else // signal is decaying + signal = (gRandom->Rndm() + par[3]) * TMath::Gaus(xx, par[4] + par[1], par[2]) ; + + return signal < 0. ? 0. : signal ; +} + +//__________________________________________________________________ +Bool_t AliEMCAL::RawSampledResponse(const Float_t dtime, const Int_t damp, Int_t * adcH, Int_t * adcL) const +{ + // for a start time dtime and an amplitude damp given by digit, + // calculates the raw sampled response + + const Int_t kRawSignalOverflow = 0x3FF ; + Bool_t highGain = kFALSE ; + + TF1 f1("signal", RawResponseFunction, 0, fkTimeBins, 5); + + f1.SetParNames("Time Max", "Peaking time", "Decay width", "Amp max", "Start time") ; + f1.SetParameter(0, fTimeMax) ; + f1.SetParameter(1, fTimePeak) ; + f1.SetParameter(2, fTimeRes) ; + f1.SetParameter(3, damp) ; + f1.SetParameter(4, dtime) ; + + for (Int_t iTime = 0; iTime < fkTimeBins; iTime++) { + Double_t time = iTime * fTimeMax/fkTimeBins; + Double_t signal = f1.Eval(time) ; + adcL[iTime] = static_cast(signal + 0.5) ; + if ( adcL[iTime] > kRawSignalOverflow) // larger than 10 bits + adcL[iTime] = kRawSignalOverflow ; + adcH[iTime] = static_cast(0.5 + (signal / fHighGainFactor)) ; + if (adcH[iTime] > 0) + highGain = kTRUE; + } + return highGain ; } //____________________________________________________________________________ void AliEMCAL::SetTreeAddress() { - + // Linking Hits in Tree to Hits array TBranch *branch; char branchname[20]; sprintf(branchname,"%s",GetName()); // Branch address for hit tree - TTree *treeH = gAlice->TreeH(); - if (treeH && fHits) { + TTree *treeH = TreeH(); + if (treeH) { branch = treeH->GetBranch(branchname); - if (branch) branch->SetAddress(&fHits); + if (branch) + { + if (fHits == 0x0) + fHits= new TClonesArray("AliEMCALHit",1000); + branch->SetAddress(&fHits); + } + else + { + Warning("SetTreeAddress","<%s> Failed",GetName()); + } } } -//____________________________________________________________________________ -void AliEMCAL::WriteQA() -{ - // Make TreeQA in the output file. - - if(fTreeQA == 0) - fTreeQA = new TTree("TreeQA", "QA Alarms") ; - // Create Alarms branches -// Int_t bufferSize = 32000 ; -// Int_t splitlevel = 0 ; -// TFolder * alarmsF = (TFolder*)gROOT->FindObjectAny("Folders/Run/Conditions/QA/PHOS") ; -// TString branchName(alarmsF->GetName()); -// TBranch * alarmsBranch = fTreeQA->Branch(branchName,"TFolder", &alarmsF, bufferSize, splitlevel); -// TString branchTitle = branchName + " QA alarms" ; -// alarmsBranch->SetTitle(branchTitle); -// alarmsBranch->Fill() ; - - //fTreeQA->Fill() ; -} + +