X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=EMCAL%2FAliEMCAL.cxx;h=58fcab1c4b0a5eec31addc1c106b8ba27d8ef53c;hb=defd951885ee71623935defb75437af7c947bced;hp=2822ed5ff77a883852ee6b38c8b22f1589c22267;hpb=4d33c79705e5d70cfcb46afd0a19194b584996c0;p=u%2Fmrichter%2FAliRoot.git diff --git a/EMCAL/AliEMCAL.cxx b/EMCAL/AliEMCAL.cxx index 2822ed5ff77..58fcab1c4b0 100644 --- a/EMCAL/AliEMCAL.cxx +++ b/EMCAL/AliEMCAL.cxx @@ -14,7 +14,6 @@ **************************************************************************/ /* $Id$ */ - //_________________________________________________________________________ // Base Class for EMCAL description: // This class contains material definitions @@ -22,74 +21,180 @@ //*-- Author: Yves Schutz (SUBATECH) // //*-- Additional Contributions: Sahal Yacoob (LBNL/UCT) +// : Alexei Pavlinov (WSU) // ////////////////////////////////////////////////////////////////////////////// // --- ROOT system --- class TFile; #include -#include -#include +#include +#include #include -#include #include +#include +#include // --- Standard library --- // --- AliRoot header files --- #include "AliMagF.h" +#include "AliLog.h" #include "AliEMCAL.h" -#include "AliEMCALGetter.h" #include "AliRun.h" +#include "AliRunLoader.h" +#include "AliCDBManager.h" +#include "AliEMCALLoader.h" #include "AliEMCALSDigitizer.h" #include "AliEMCALDigitizer.h" -#include "AliAltroBuffer.h" +#include "AliEMCALDigit.h" +#include "AliEMCALRawUtils.h" +#include "AliCDBManager.h" +#include "AliCDBEntry.h" +#include "AliEMCALRawUtils.h" +#include "AliRawReader.h" +#include "AliEMCALTriggerData.h" +#include "AliEMCALRecParam.h" +#include "AliRawEventHeaderBase.h" ClassImp(AliEMCAL) -Double_t AliEMCAL::fgCapa = 1.; // 1pF -Int_t AliEMCAL::fgOrder = 2 ; -Double_t AliEMCAL::fgTimeMax = 2.56E-5 ; // each sample is over 100 ns fTimeMax/fTimeBins -Double_t AliEMCAL::fgTimePeak = 4.1E-6 ; // 4 micro seconds -Double_t AliEMCAL::fgTimeTrigger = 100E-9 ; // 100ns, just for a reference - + +//for embedding +AliEMCALRawUtils* AliEMCAL::fgRawUtils = 0; // EMCAL raw utilities class + //____________________________________________________________________________ -AliEMCAL::AliEMCAL():AliDetector() +AliEMCAL::AliEMCAL() + : AliDetector(), + fBirkC0(0), + fBirkC1(0.), + fBirkC2(0.), + fGeometry(0), + fCheckRunNumberAndGeoVersion(kTRUE), + fTriggerData(0x0) { // Default ctor fName = "EMCAL" ; + InitConstants(); + + // Should call AliEMCALGeometry::GetInstance(EMCAL->GetTitle(),"") for getting EMCAL geometry } //____________________________________________________________________________ -AliEMCAL::AliEMCAL(const char* name, const char* title): AliDetector(name,title) +AliEMCAL::AliEMCAL(const char* name, const char* title, + const Bool_t checkGeoAndRun) + : AliDetector(name,title), + fBirkC0(0), + fBirkC1(0.), + fBirkC2(0.), + fGeometry(0), + fCheckRunNumberAndGeoVersion(checkGeoAndRun), + fTriggerData(0x0) { // ctor : title is used to identify the layout - - fHighCharge = 8.2 ; // adjusted for a high gain range of 5.12 GeV (10 bits) - fHighGain = 6.64 ; - fHighLowGainFactor = 16. ; // adjusted for a low gain range of 82 GeV (10 bits) - fLowGainOffset = 1 ; // offset added to the module id to distinguish high and low gain data + InitConstants(); + } //____________________________________________________________________________ AliEMCAL::~AliEMCAL() { - + //dtor + delete fgRawUtils; + delete fTriggerData; } //____________________________________________________________________________ -void AliEMCAL::Copy(AliEMCAL & emcal) +void AliEMCAL::InitConstants() { - TObject::Copy(emcal) ; - emcal.fHighCharge = fHighCharge ; - emcal.fHighGain = fHighGain ; - emcal.fHighLowGainFactor = fHighLowGainFactor ; - emcal.fLowGainOffset = fLowGainOffset; + //initialize EMCAL values + fBirkC0 = 1; + fBirkC1 = 0.013/1.032; + fBirkC2 = 9.6e-6/(1.032 * 1.032); } +//Not needed, modify $ALICE_ROOT/data/galice.cuts instead. +//Load the modified one in the configuration file with SetTransPar +// //____________________________________________________________________________ +// void AliEMCAL::DefineMediumParameters() +// { +// // +// // EMCAL cuts (Geant3) +// // +// Int_t * idtmed = fIdtmed->GetArray() - 1599 ; +// // --- Set decent energy thresholds for gamma and electron tracking + +// // Tracking threshold for photons and electrons in Lead +// Float_t cutgam=10.e-5; // 100 kev; +// Float_t cutele=10.e-5; // 100 kev; +// TString ntmp(GetTitle()); +// ntmp.ToUpper(); +// if(ntmp.Contains("10KEV")) { +// cutele = cutgam = 1.e-5; +// } else if(ntmp.Contains("50KEV")) { +// cutele = cutgam = 5.e-5; +// } else if(ntmp.Contains("100KEV")) { +// cutele = cutgam = 1.e-4; +// } else if(ntmp.Contains("200KEV")) { +// cutele = cutgam = 2.e-4; +// } else if(ntmp.Contains("500KEV")) { +// cutele = cutgam = 5.e-4; +// } + +// TVirtualMC::GetMC()->Gstpar(idtmed[1600],"CUTGAM", cutgam); +// TVirtualMC::GetMC()->Gstpar(idtmed[1600],"CUTELE", cutele); // 1MEV -> 0.1MEV; 15-aug-05 +// TVirtualMC::GetMC()->Gstpar(idtmed[1600],"BCUTE", cutgam); // BCUTE and BCUTM start from GUTGUM +// TVirtualMC::GetMC()->Gstpar(idtmed[1600],"BCUTM", cutgam); // BCUTE and BCUTM start from GUTGUM +// // --- Generate explicitly delta rays in Lead --- +// TVirtualMC::GetMC()->Gstpar(idtmed[1600], "LOSS", 3) ; +// TVirtualMC::GetMC()->Gstpar(idtmed[1600], "DRAY", 1) ; +// TVirtualMC::GetMC()->Gstpar(idtmed[1600], "DCUTE", cutele) ; +// TVirtualMC::GetMC()->Gstpar(idtmed[1600], "DCUTM", cutele) ; + +// // --- in aluminium parts --- +// TVirtualMC::GetMC()->Gstpar(idtmed[1602],"CUTGAM", cutgam) ; +// TVirtualMC::GetMC()->Gstpar(idtmed[1602],"CUTELE", cutele) ; +// TVirtualMC::GetMC()->Gstpar(idtmed[1602],"BCUTE", cutgam); // BCUTE and BCUTM start from GUTGUM +// TVirtualMC::GetMC()->Gstpar(idtmed[1602],"BCUTM", cutgam); // BCUTE and BCUTM start from GUTGUM +// TVirtualMC::GetMC()->Gstpar(idtmed[1602], "LOSS",3.) ; +// TVirtualMC::GetMC()->Gstpar(idtmed[1602], "DRAY",1.) ; +// TVirtualMC::GetMC()->Gstpar(idtmed[1602], "DCUTE", cutele) ; +// TVirtualMC::GetMC()->Gstpar(idtmed[1602], "DCUTM", cutele) ; + +// // --- and finally thresholds for photons and electrons in the scintillator --- +// TVirtualMC::GetMC()->Gstpar(idtmed[1601],"CUTGAM", cutgam) ; +// TVirtualMC::GetMC()->Gstpar(idtmed[1601],"CUTELE", cutele) ;// 1MEV -> 0.1MEV; 15-aug-05 +// TVirtualMC::GetMC()->Gstpar(idtmed[1601],"BCUTE", cutgam); // BCUTE and BCUTM start from GUTGUM +// TVirtualMC::GetMC()->Gstpar(idtmed[1601],"BCUTM", cutgam); // BCUTE and BCUTM start from GUTGUM +// TVirtualMC::GetMC()->Gstpar(idtmed[1601], "LOSS",3) ; // generate delta rays +// TVirtualMC::GetMC()->Gstpar(idtmed[1601], "DRAY",1) ; +// TVirtualMC::GetMC()->Gstpar(idtmed[1601], "DCUTE", cutele) ; +// TVirtualMC::GetMC()->Gstpar(idtmed[1601], "DCUTM", cutele) ; + +// // S steel - +// TVirtualMC::GetMC()->Gstpar(idtmed[1603],"CUTGAM", cutgam); +// TVirtualMC::GetMC()->Gstpar(idtmed[1603],"CUTELE", cutele); +// TVirtualMC::GetMC()->Gstpar(idtmed[1603],"BCUTE", cutgam); // BCUTE and BCUTM start from GUTGUM +// TVirtualMC::GetMC()->Gstpar(idtmed[1603],"BCUTM", cutgam); // BCUTE and BCUTM start from GUTGUM +// // --- Generate explicitly delta rays +// TVirtualMC::GetMC()->Gstpar(idtmed[1603], "LOSS",3); +// TVirtualMC::GetMC()->Gstpar(idtmed[1603], "DRAY",1); +// TVirtualMC::GetMC()->Gstpar(idtmed[1603], "DCUTE", cutele) ; +// TVirtualMC::GetMC()->Gstpar(idtmed[1603], "DCUTM", cutele) ; + +// AliEMCALGeometry* geom = GetGeometry(); +// if(geom->GetILOSS()>=0) { +// for(int i=1600; i<=1603; i++) TVirtualMC::GetMC()->Gstpar(idtmed[i], "LOSS", geom->GetILOSS()) ; +// } +// if(geom->GetIHADR()>=0) { +// for(int i=1600; i<=1603; i++) TVirtualMC::GetMC()->Gstpar(idtmed[i], "HADR", geom->GetIHADR()) ; +// } +// } + //____________________________________________________________________________ -AliDigitizer* AliEMCAL::CreateDigitizer(AliRunDigitizer* manager) const +AliDigitizer* AliEMCAL::CreateDigitizer(AliDigitizationInput* digInput) const { - return new AliEMCALDigitizer(manager); + //create and return the digitizer + return new AliEMCALDigitizer(digInput); } //____________________________________________________________________________ @@ -97,7 +202,6 @@ void AliEMCAL::CreateMaterials() { // Definitions of materials to build EMCAL and associated tracking media. // media number in idtmed are 1599 to 1698. - // --- Air --- Float_t aAir[4]={12.0107,14.0067,15.9994,39.948}; Float_t zAir[4]={6.,7.,8.,18.}; @@ -121,170 +225,203 @@ void AliEMCAL::CreateMaterials() AliMaterial(3, "Al$", 26.98, 13., 2.7, 8.9, 999., 0, 0) ; // --- Absorption length is ignored ^ + // 25-aug-04 by PAI - see PMD/AliPMDv0.cxx for STEEL definition + Float_t asteel[4] = { 55.847,51.9961,58.6934,28.0855 }; + Float_t zsteel[4] = { 26.,24.,28.,14. }; + Float_t wsteel[4] = { .715,.18,.1,.005 }; + AliMixture(4, "STAINLESS STEEL$", asteel, zsteel, 7.88, 4, wsteel); + + // Oct 26,2010 : Multipurpose Copy Paper UNV-21200), weiht 75 g/m**2. + // *Cellulose C6H10O5 + // Component C A=12.01 Z=6. W=6./21. + // Component H A=1. Z=1. W=10./21. + // Component O A=16. Z=8. W=5./21. + Float_t apaper[3] = { 12.01, 1.0, 16.0}; + Float_t zpaper[3] = { 6.0, 1.0, 8.0}; + Float_t wpaper[3] = {6./21., 10./21., 5./21.}; + AliMixture(5, "BondPaper$", apaper, zpaper, 0.75, 3, wpaper); + // DEFINITION OF THE TRACKING MEDIA + // Look to the $ALICE_ROOT/data/galice.cuts for particular values + // of cuts. + // Don't forget to add a new tracking medium with non-default cuts // 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() ; + Int_t isxfld = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->Integ() ; + Float_t sxmgmx = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->Max() ; // Air -> idtmed[1599] - AliMedium(0, "Air $", 0, 0, + AliMedium(0, "Air$", 0, 0, isxfld, sxmgmx, 10.0, 1.0, 0.1, 0.1, 10.0, 0, 0) ; // The Lead -> idtmed[1600] - AliMedium(1, "Lead $", 1, 0, + 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) ; + float deemax = 0.1; // maximum fractional energy loss in one step (0 < DEEMAX < deemax ) + AliMedium(2, "Scintillator$", 2, 1, + isxfld, sxmgmx, 10.0, 0.001, deemax, 0.001, 0.001, 0, 0) ; // Various Aluminium parts made of Al -> idtmed[1602] - AliMedium(3, "Al parts $", 3, 0, + AliMedium(3, "Al$", 3, 0, isxfld, sxmgmx, 10.0, 0.1, 0.1, 0.001, 0.001, 0, 0) ; + // 25-aug-04 by PAI : see PMD/AliPMDv0.cxx for STEEL definition -> idtmed[1603] + AliMedium(4, "S steel$", 4, 0, + isxfld, sxmgmx, 10.0, 0.01, 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 - gMC->Gstpar(idtmed[1600],"CUTGAM",0.00008) ; - gMC->Gstpar(idtmed[1600],"CUTELE",0.001) ; - gMC->Gstpar(idtmed[1600],"BCUTE",0.0001) ; - - // --- Generate explicitly delta rays in Lead --- - gMC->Gstpar(idtmed[1600], "LOSS",3.) ; - gMC->Gstpar(idtmed[1600], "DRAY",1.) ; - gMC->Gstpar(idtmed[1600], "DCUTE",0.00001) ; - gMC->Gstpar(idtmed[1600], "DCUTM",0.00001) ; + // Oct 26,2010; Nov 24,2010 -> idtmed[1604] + deemax = 0.01; + AliMedium(5, "Paper$", 5, 0, + isxfld, sxmgmx, 10.0, deemax, 0.1, 0.001, 0.001, 0, 0) ; -// --- 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::Init() +{ + // Init + //Not needed, modify $ALICE_ROOT/data/galice.cuts instead. + //Load the modified one in the configuration file with SetTransPar + //DefineMediumParameters(); +} + +//____________________________________________________________________________ +void AliEMCAL::Digits2Raw() { + + static AliEMCALRawUtils rawUtils; + rawUtils.Digits2Raw(); } - //____________________________________________________________________________ -void AliEMCAL::Digits2Raw() -{ - // convert digits of the current event to raw data - AliEMCALLoader * loader = dynamic_cast(fLoader) ; +void AliEMCAL::Hits2SDigits() +{ +// create summable digits - // get the digits - loader->LoadDigits(); - TClonesArray* digits = loader->Digits() ; + GetGeometry(); + AliEMCALSDigitizer emcalDigitizer(fLoader->GetRunLoader()->GetFileName().Data()) ; + emcalDigitizer.SetEventRange(0, -1) ; // do all the events + emcalDigitizer.Digitize() ; +} - if (!digits) { - Error("Digits2Raw", "no digits found !"); - return; +//______________________________________________________________________ +Bool_t AliEMCAL::Raw2SDigits(AliRawReader* rawReader){ + + // Conversion from raw data to EMCAL sdigits. + // Does the same as AliEMCALReconstructor::ConvertDigits() + // Needed to embed real data and simulation + // Works on a single-event basis + + rawReader->Reset() ; + + //Get/create the sdigits tree and array + AliRunLoader *rl = AliRunLoader::Instance(); + AliEMCALLoader *emcalLoader = dynamic_cast(rl->GetDetectorLoader("EMCAL")); + + if(!emcalLoader){ + AliFatal("NULL loader"); + return kFALSE; } + + emcalLoader->GetEvent(); + emcalLoader->LoadSDigits("UPDATE"); - // get the digitizer - loader->LoadDigitizer(); - AliEMCALDigitizer * digitizer = dynamic_cast(loader->Digitizer()) ; + TTree * treeS = emcalLoader->TreeS(); + if ( !treeS ) { + emcalLoader->MakeSDigitsContainer(); + treeS = emcalLoader->TreeS(); + } - // get the geometry - AliEMCALGeometry* geom = GetGeometry(); - if (!geom) { - Error("Digits2Raw", "no geometry found !"); - return; + if(!emcalLoader->SDigits()) { + AliFatal("No sdigits array available\n"); + return kFALSE; } - - // some digitization constants - const Int_t kDDLOffset = 0x800; - const Int_t kThreshold = 1; - 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 = dynamic_cast(digits->At(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; + TClonesArray * sdigits = emcalLoader->SDigits(); + sdigits->Clear("C"); + + //Trigger sdigits + if(!fTriggerData)fTriggerData = new AliEMCALTriggerData(); + fTriggerData->SetMode(1); + const Int_t nTRU = GetGeometry()->GetNTotalTRU(); + TClonesArray *digitsTrg = new TClonesArray("AliEMCALTriggerRawDigit", nTRU * 96); + Int_t bufsize = 32000; + treeS->Branch("EMTRG", &digitsTrg, bufsize); + + + //Only physics events + if (rawReader->GetType()== AliRawEventHeaderBase::kPhysicsEvent) { + + if(!fgRawUtils) fgRawUtils = new AliEMCALRawUtils; + //must be done here because, in constructor, option is not yet known + fgRawUtils->SetOption(GetOption()); + + // Set parameters from OCDB to raw utils + AliEMCALRecParam* recpar = emcalLoader->ReconstructionParameters(0); + // fgRawUtils->SetRawFormatHighLowGainFactor(recpar->GetHighLowGainFactor()); + // fgRawUtils->SetRawFormatOrder(recpar->GetOrderParameter()); + // fgRawUtils->SetRawFormatTau(recpar->GetTau()); + fgRawUtils->SetNoiseThreshold(recpar->GetNoiseThreshold()); + fgRawUtils->SetNPedSamples(recpar->GetNPedSamples()); + fgRawUtils->SetRemoveBadChannels(recpar->GetRemoveBadChannels()); + fgRawUtils->SetFittingAlgorithm(recpar->GetFittingAlgorithm()); + fgRawUtils->SetFALTROUsage(recpar->UseFALTRO()); + // fgRawUtils->SetTimeMin(recpar->GetTimeMin()); + // fgRawUtils->SetTimeMax(recpar->GetTimeMax()); + + //Fit + fgRawUtils->Raw2Digits(rawReader,sdigits,emcalLoader->PedestalData(),digitsTrg,fTriggerData); + + }//skip calibration event + else{ + AliDebug(1," Calibration Event, skip!"); + } + + //Final arrangements of the array, set all sdigits as embedded + sdigits->Sort() ; + for (Int_t iSDigit = 0 ; iSDigit < sdigits->GetEntriesFast() ; iSDigit++) { + AliEMCALDigit * sdigit = dynamic_cast(sdigits->At(iSDigit)) ; + if(sdigit){ + sdigit->SetIndexInList(iSDigit) ; + sdigit->SetType(AliEMCALDigit::kEmbedded); } - - // out of time range signal (?) - if (digit->GetTimeR() > GetRawFormatTimeMax() ) { - buffer->FillBuffer(digit->GetAmp()); - buffer->FillBuffer(GetRawFormatTimeBins() ); // time bin - buffer->FillBuffer(3); // bunch length - buffer->WriteTrailer(3, idDDL, 0, 0); // trailer - - // calculate the time response function - } else { - Double_t energy = 0 ; - energy = digit->GetAmp() * digitizer->GetECAchannel() + digitizer->GetECApedestal() ; - - Bool_t lowgain = RawSampledResponse(digit->GetTimeR(), energy, adcValuesHigh, adcValuesLow) ; - - if (lowgain) - buffer->WriteChannel(iDDL, 0, fLowGainOffset, - GetRawFormatTimeBins(), adcValuesLow, kThreshold); - else - buffer->WriteChannel(iDDL, 0, 0, - GetRawFormatTimeBins(), adcValuesHigh, kThreshold); - + else { + AliFatal("sdigit is NULL!"); } - } + } + + AliDebug(1,Form("Embedded sdigits entries %d \n",sdigits->GetEntriesFast())); + + //Write array, clean arrays, unload .. + + Int_t bufferSize = 32000 ; + TBranch * sdigitsBranch = treeS->GetBranch("EMCAL"); + if (sdigitsBranch) + sdigitsBranch->SetAddress(&sdigits); + else + treeS->Branch("EMCAL",&sdigits,bufferSize); + + treeS->Fill(); + emcalLoader->WriteSDigits("OVERWRITE"); + emcalLoader->UnloadSDigits(); + + digitsTrg->Delete(); + delete digitsTrg; + + return kTRUE; - // write real header and close last file - if (buffer) { - buffer->Flush(); - buffer->WriteDataHeader(kFALSE, kFALSE); - delete buffer; - } - - loader->UnloadDigits(); } //____________________________________________________________________________ -void AliEMCAL::Hits2SDigits() -{ -// create summable digits - - AliEMCALSDigitizer emcalDigitizer(fLoader->GetRunLoader()->GetFileName().Data()) ; - emcalDigitizer.SetEventRange(0, -1) ; // do all the events - emcalDigitizer.ExecuteTask() ; -} -//____________________________________________________________________________ AliLoader* AliEMCAL::MakeLoader(const char* topfoldername) { //different behaviour than standard (singleton getter) @@ -293,99 +430,115 @@ AliLoader* AliEMCAL::MakeLoader(const char* topfoldername) return fLoader; } -//__________________________________________________________________ -Double_t AliEMCAL::RawResponseFunction(Double_t *x, Double_t *par) +//____________________________________________________________________________ + +AliEMCALGeometry* AliEMCAL::GetGeometry() const { - // Shape of the electronics raw reponse: - // It is a semi-gaussian, 2nd order Gamma function of the general form - // v(t) = n**n * Q * A**n / C *(t/tp)**n * exp(-n * t/tp) with - // tp : peaking time par[0] - // n : order of the function - // C : integrating capacitor in the preamplifier - // A : open loop gain of the preamplifier - // Q : the total APD charge to be measured Q = C * energy + //Initializes and returns geometry - Double_t signal ; - Double_t xx = x[0] - ( fgTimeTrigger + par[3] ) ; - - if (xx < 0 || xx > fgTimeMax) - signal = 0. ; - else { - Double_t fac = par[0] * TMath::Power(fgOrder, fgOrder) * TMath::Power(par[1], fgOrder) / fgCapa ; - signal = fac * par[2] * TMath::Power(xx / fgTimePeak, fgOrder) * TMath::Exp(-fgOrder * (xx / fgTimePeak)) ; + // Pass the transpor model name (Geant3, Geant4, Fluka) and title to the geometry + TString mcname = ""; + TString mctitle = ""; + if(TVirtualMC::GetMC()){ + mcname = TVirtualMC::GetMC()->GetName() ; + mctitle = TVirtualMC::GetMC()->GetTitle() ; } - return signal ; -} - -//__________________________________________________________________ -Double_t AliEMCAL::RawResponseFunctionMax(Double_t charge, Double_t gain) -{ - return ( charge * TMath::Power(fgOrder, fgOrder) * TMath::Power(gain, fgOrder) - / ( fgCapa * TMath::Exp(fgOrder) ) ); - -} -//__________________________________________________________________ -Bool_t AliEMCAL::RawSampledResponse( -const Double_t dtime, const Double_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 AliEMCAL::RawResponseFunction - - const Int_t kRawSignalOverflow = 0x3FF ; - Bool_t lowGain = kFALSE ; - - TF1 signalF("signal", RawResponseFunction, 0, GetRawFormatTimeMax(), 4); - - for (Int_t iTime = 0; iTime < GetRawFormatTimeBins(); iTime++) { - signalF.SetParameter(0, GetRawFormatHighCharge() ) ; - signalF.SetParameter(1, GetRawFormatHighGain() ) ; - signalF.SetParameter(2, damp) ; - signalF.SetParameter(3, dtime) ; - Double_t time = iTime * GetRawFormatTimeMax() / GetRawFormatTimeBins() ; - Double_t signal = signalF.Eval(time) ; - if ( static_cast(signal+0.5) > kRawSignalOverflow ){ // larger than 10 bits - signal = kRawSignalOverflow ; - lowGain = kTRUE ; + + TString geoName(GetTitle()); + + //Check if run number and requested geometry correspond to the same geometry as + //in real data taking. To prevent errors in official simulation productions + if(!(AliEMCALGeometry::GetInstance())) + { + // Check the transport model name and option, set sampling fraction depending on it + if(!fCheckRunNumberAndGeoVersion){// Set geometry with the name used in the configuration file + AliInfo(Form("Geometry name in use <<%s>>, requested via Config file", geoName.Data())); + return AliEMCALGeometry::GetInstance(GetTitle(),"EMCAL",mcname,mctitle) ; } - adcH[iTime] = static_cast(signal + 0.5) ; - - signalF.SetParameter(0, GetRawFormatLowCharge() ) ; - signalF.SetParameter(1, GetRawFormatLowGain() ) ; - signal = signalF.Eval(time) ; - if ( static_cast(signal+0.5) > kRawSignalOverflow) // larger than 10 bits - signal = kRawSignalOverflow ; - adcL[iTime] = static_cast(0.5 + signal ) ; - - } - return lowGain ; -} + else + {//Check run number and version and set the corresponding one. + //Get run number + //AliRunLoader *rl = AliRunLoader::Instance(); + //Int_t runNumber = rl->GetRunNumber(); + + AliCDBManager* man = AliCDBManager::Instance(); + Int_t runNumber = man->GetRun(); + + //Instanciate geometry depending on the run number + if(runNumber >= 104064 && runNumber < 140000 ){//2009-2010 runs + //First year geometry, 4 SM. + + if(!geoName.Contains("FIRSTYEARV1")) + { + AliInfo(Form("*** ATTENTION *** \n \t Specified geometry name <<%s>> for run %d is not considered! \n \t In use <>, check run number and year \n ", + geoName.Data(),runNumber)); + } + else + { + AliDebug(1,"Initialized geometry with name <>"); + } + + return AliEMCALGeometry::GetInstance("EMCAL_FIRSTYEARV1","EMCAL",mcname,mctitle) ;// Set geometry with the name used in the configuration file + } + else if(runNumber >= 140000 && runNumber <= 170593) + { + //Almost complete EMCAL geometry, 10 SM. Year 2011 configuration + + if(!geoName.Contains("COMPLETEV1")) + { + AliInfo(Form("*** ATTENTION *** \n \t Specified geometry name <<%s>> for run %d is not considered! \n \t In use <>, check run number and year \n ", + geoName.Data(),runNumber)); + } + else + { + AliDebug(1,"Initialized geometry with name <>"); + } + + return AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1","EMCAL",mcname,mctitle) ;// Set geometry with the name used in the configuration file + } + else if(runNumber > 176000 && runNumber <= 197692) + { + //Complete EMCAL geometry, 12 SM. Year 2012 and on + if(!geoName.Contains("COMPLETE12SMV1")) + { + AliInfo(Form("*** ATTENTION *** \n \t Specified geometry name <<%s>> for run %d is not considered! \n \t In use <>, check run number and year \n ", geoName.Data(),runNumber)); + } else { + AliDebug(1,"Initialized geometry with name <>"); + } + + return AliEMCALGeometry::GetInstance("EMCAL_COMPLETE12SMV1","EMCAL",mcname,mctitle) ;// Set geometry with the name used in the configuration file -//____________________________________________________________________________ -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 = TreeH(); - if (treeH) { - branch = treeH->GetBranch(branchname); - if (branch) - { - if (fHits == 0x0) - fHits= new TClonesArray("AliEMCALHit",1000); - branch->SetAddress(&fHits); } - else + else { - Warning("SetTreeAddress","<%s> Failed",GetName()); + //EMCAL + DCAL geometry, 20 SM. Year 2014 and on + + if(geoName.Contains("DCAL_DEV")) + { + AliInfo("*** ATTENTION *** \n \t Set geometry name <> \n "); + return AliEMCALGeometry::GetInstance("EMCAL_COMPLETE12SMV1_DCAL_DEV","EMCAL",mcname,mctitle) ; + } + else if(geoName.Contains("DCAL_8SM")) + { + AliInfo("*** ATTENTION *** \n \t Set geometry name <> \n "); + return AliEMCALGeometry::GetInstance("EMCAL_COMPLETE12SMV1_DCAL_8SM","EMCAL",mcname,mctitle) ; + } + + if(!geoName.Contains("COMPLETE12SMV1")) + { + AliInfo(Form("*** ATTENTION *** \n \t Specified geometry name <<%s>> for run %d is not considered! \n \t In use <>, check run number and year \n ", geoName.Data(),runNumber)); + } + else + { + AliDebug(1,"Initialized geometry with name <>"); + } + + return AliEMCALGeometry::GetInstance("EMCAL_COMPLETE12SMV1_DCAL","EMCAL",mcname,mctitle) ;// Set geometry with the name used in the configuration file } - } + } + }// Init geometry for the first time + + + return AliEMCALGeometry::GetInstance(); + } - - - - -