X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=PHOS%2FAliPHOS.cxx;h=d0d882a62d5dce981e3ed87366641f24e32d166f;hp=6ab9856d93ad87ca9004fcc23bece2c07440a2d7;hb=48b634ff934842fabcb6cfedda2ada1875d97b1f;hpb=eb0b105198969ebdd08334fdbd7c57e4b9e98038 diff --git a/PHOS/AliPHOS.cxx b/PHOS/AliPHOS.cxx index 6ab9856d93a..d0d882a62d5 100644 --- a/PHOS/AliPHOS.cxx +++ b/PHOS/AliPHOS.cxx @@ -16,15 +16,41 @@ /* History of cvs commits: * * $Log$ - * Revision 1.91 2005/07/27 15:08:53 kharlov - * Mixture ArCO2 is corrected + * Revision 1.116 2007/10/10 09:05:10 schutz + * Changing name QualAss to QA * - * Revision 1.90 2005/06/17 07:39:07 hristov - * Removing GetDebug and SetDebug from AliRun and AliModule. Using AliLog for the messages + * Revision 1.115 2007/08/22 09:20:50 hristov + * Updated QA classes (Yves) * - * Revision 1.89 2005/05/28 12:10:07 schutz - * Copy constructor is corrected (by T.P.) + * Revision 1.114 2007/08/07 14:12:03 kharlov + * Quality assurance added (Yves Schutz) * + * Revision 1.113 2007/07/18 16:29:54 policheh + * Raw Sdigits energy converted to GeV. + * + * Revision 1.112 2007/02/25 22:59:13 policheh + * Digits2Raw(): ALTRO buffer and mapping created per each DDL. + * + * Revision 1.111 2007/02/18 15:21:47 kharlov + * Corrections for memory leak in Digits2Raw due to AliAltroMapping + * + * Revision 1.110 2007/02/13 10:52:08 policheh + * Raw2SDigits() implemented + * + * Revision 1.109 2007/02/05 10:43:25 hristov + * Changes for correct initialization of Geant4 (Mihaela) + * + * Revision 1.108 2007/02/01 10:34:47 hristov + * Removing warnings on Solaris x86 + * + * Revision 1.107 2007/01/29 16:29:37 kharlov + * Digits2Raw(): special workaround for digits with time out of range + * + * Revision 1.106 2007/01/17 17:28:56 kharlov + * Extract ALTRO sample generation to a separate class AliPHOSPulseGenerator + * + * Revision 1.105 2007/01/12 21:44:29 kharlov + * Simulate and reconstruct two gains simulaneouslsy */ //_________________________________________________________________________ @@ -40,36 +66,46 @@ // --- ROOT system --- class TFile; +#include #include -#include -#include +#include #include -#include #include +#include +#include // --- Standard library --- // --- AliRoot header files --- #include "AliMagF.h" #include "AliPHOS.h" -#include "AliPHOSGetter.h" +#include "AliPHOSLoader.h" #include "AliRun.h" +#include "AliRawReader.h" #include "AliPHOSDigitizer.h" #include "AliPHOSSDigitizer.h" #include "AliPHOSDigit.h" #include "AliAltroBuffer.h" +#include "AliAltroMapping.h" +#include "AliCaloAltroMapping.h" #include "AliLog.h" +#include "AliCDBManager.h" +#include "AliCDBEntry.h" +#include "AliCDBStorage.h" +#include "AliPHOSCalibData.h" +#include "AliPHOSPulseGenerator.h" +#include "AliDAQ.h" +#include "AliPHOSRawFitterv0.h" +#include "AliPHOSCalibData.h" +#include "AliPHOSRawDigiProducer.h" +#include "AliPHOSQAChecker.h" +#include "AliPHOSRecoParam.h" +#include "AliPHOSSimParam.h" ClassImp(AliPHOS) -Double_t AliPHOS::fgCapa = 1.; // 1pF -Int_t AliPHOS::fgOrder = 2 ; -Double_t AliPHOS::fgTimeMax = 2.56E-5 ; // each sample is over 100 ns fTimeMax/fTimeBins -Double_t AliPHOS::fgTimePeak = 4.1E-6 ; // 4 micro seconds -Double_t AliPHOS::fgTimeTrigger = 100E-9 ; // 100ns, just for a reference - //____________________________________________________________________________ - AliPHOS:: AliPHOS() : AliDetector() + AliPHOS:: AliPHOS() : AliDetector(),fgCalibData(0) { // Default ctor fName = "PHOS" ; @@ -77,34 +113,16 @@ Double_t AliPHOS::fgTimeTrigger = 100E-9 ; // 100ns, just for a reference } //____________________________________________________________________________ -AliPHOS::AliPHOS(const char* name, const char* title): AliDetector(name, title) +AliPHOS::AliPHOS(const char* name, const char* title): AliDetector(name, title), +fgCalibData(0) { // 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 = GetGeometry()->GetNModules() + 1 ; - // offset added to the module id to distinguish high and low gain data } //____________________________________________________________________________ AliPHOS::~AliPHOS() { -} - -//____________________________________________________________________________ -void AliPHOS::Copy(TObject &obj)const -{ - // copy method to be used byy the cpy ctor - TObject::Copy(obj); - - AliPHOS &phos = static_cast(obj); - - phos.fHighCharge = fHighCharge ; - phos.fHighGain = fHighGain ; - phos.fHighLowGainFactor = fHighLowGainFactor ; - phos.fLowGainOffset = fLowGainOffset; + if(fgCalibData) delete fgCalibData ; } //____________________________________________________________________________ @@ -208,7 +226,8 @@ void AliPHOS::CreateMaterials() Float_t aG10[4] = { 12., 1., 16., 28.} ; Float_t zG10[4] = { 6., 1., 8., 14.} ; Float_t wG10[4] = { .259, .288, .248, .205} ; - Float_t dG10 = 1.7 ; + Float_t dG10 = 1.7 ; + AliMixture(12, "G10$", aG10, zG10, dG10, -4, wG10); @@ -237,7 +256,8 @@ void AliPHOS::CreateMaterials() wArCO[1] = (1-arContent)*1; wArCO[2] = (1-arContent)*2; Float_t dArCO = arContent*dAr + (1-arContent)*dCO ; - AliMixture(16, "ArCO2$", aArCO, zArCO, dArCO, -3, wArCO) ; + AliMixture(16, "ArCO2$", aArCO, zArCO, dArCO, -3, wArCO) ; + // --- Stainless steel (let it be pure iron) --- AliMaterial(17, "Steel$", 55.845, 26, 7.87, 1.76, 0., 0, 0) ; @@ -273,9 +293,8 @@ void AliPHOS::CreateMaterials() // DEFINITION OF THE TRACKING MEDIA // for PHOS: idtmed[699->798] equivalent to fIdtmed[0->100] - Int_t * idtmed = fIdtmed->GetArray() - 699 ; - 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() ; // The scintillator of the calorimeter made of PBW04 -> idtmed[699] AliMedium(0, "PHOS Xtal $", 0, 1, @@ -356,40 +375,11 @@ void AliPHOS::CreateMaterials() // Air -> idtmed[798] AliMedium(99, "Air $", 99, 0, isxfld, sxmgmx, 10.0, 1.0, 0.1, 0.1, 10.0, 0, 0) ; +} - // --- Set decent energy thresholds for gamma and electron tracking - - // Tracking threshold for photons and electrons in the scintillator crystal - gMC->Gstpar(idtmed[699], "CUTGAM",0.5E-4) ; - gMC->Gstpar(idtmed[699], "CUTELE",1.0E-4) ; - - // --- Generate explicitly delta rays in the titan cover --- - gMC->Gstpar(idtmed[704], "LOSS",3.) ; - gMC->Gstpar(idtmed[704], "DRAY",1.) ; - // --- and in aluminium parts --- - gMC->Gstpar(idtmed[701], "LOSS",3.) ; - gMC->Gstpar(idtmed[701], "DRAY",1.) ; - // --- and in PIN diode - gMC->Gstpar(idtmed[705], "LOSS",3) ; - gMC->Gstpar(idtmed[705], "DRAY",1) ; - // --- and in the passive convertor - gMC->Gstpar(idtmed[712], "LOSS",3) ; - gMC->Gstpar(idtmed[712], "DRAY",1) ; - // Tracking threshold for photons and electrons in the gas ArC02 - gMC->Gstpar(idtmed[715], "CUTGAM",1.E-5) ; - gMC->Gstpar(idtmed[715], "CUTELE",1.E-5) ; - gMC->Gstpar(idtmed[715], "CUTNEU",1.E-5) ; - gMC->Gstpar(idtmed[715], "CUTHAD",1.E-5) ; - gMC->Gstpar(idtmed[715], "CUTMUO",1.E-5) ; - gMC->Gstpar(idtmed[715], "BCUTE",1.E-5) ; - gMC->Gstpar(idtmed[715], "BCUTM",1.E-5) ; - gMC->Gstpar(idtmed[715], "DCUTE",1.E-5) ; - gMC->Gstpar(idtmed[715], "DCUTM",1.E-5) ; - gMC->Gstpar(idtmed[715], "PPCUTM",1.E-5) ; - gMC->Gstpar(idtmed[715], "LOSS",2.) ; - gMC->Gstpar(idtmed[715], "DRAY",0.) ; - gMC->Gstpar(idtmed[715], "STRA",2.) ; - +//_____________________________________________________________________________ +void AliPHOS::Init() +{ } //____________________________________________________________________________ @@ -397,7 +387,11 @@ void AliPHOS::Digits2Raw() { // convert digits of the current event to raw data - AliPHOSLoader * loader = dynamic_cast(fLoader) ; + if(AliPHOSSimParam::GetInstance()->IsEDigitizationOn()){ + AliError("Energy digitization should be OFF if use Digits2Raw") ; + } + + AliPHOSLoader * loader = static_cast(fLoader) ; // get the digits loader->LoadDigits(); @@ -408,10 +402,6 @@ void AliPHOS::Digits2Raw() return; } - // get the digitizer - loader->LoadDigitizer(); - AliPHOSDigitizer * digitizer = dynamic_cast(loader->Digitizer()) ; - // get the geometry AliPHOSGeometry* geom = GetGeometry(); if (!geom) { @@ -419,85 +409,147 @@ void AliPHOS::Digits2Raw() return; } + // get mapping from OCDB + const TObjArray* maps = AliPHOSRecoParam::GetMappings(); + if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!"); + // some digitization constants - const Int_t kDDLOffset = 0x600; // assigned to PHOS - const Int_t kThreshold = 1; // skip digits below this threshold + const Float_t kThreshold = 1.; // skip digits below 1 ADC channel + const Int_t kAdcThreshold = 1; // Lower ADC threshold to write to raw data - AliAltroBuffer* buffer = NULL; Int_t prevDDL = -1; - Int_t adcValuesLow[fkTimeBins]; - Int_t adcValuesHigh[fkTimeBins]; + + if(fgCalibData==0) + fgCalibData= new AliPHOSCalibData(-1) ; + + // Create a shaper pulse object + AliPHOSPulseGenerator *pulse = new AliPHOSPulseGenerator(); + + //Set Time step of sample + pulse->SetTimeStep(TMath::Abs(fgCalibData->GetSampleTimeStep())) ; + + Int_t *adcValuesLow = new Int_t[pulse->GetRawFormatTimeBins()]; + Int_t *adcValuesHigh= new Int_t[pulse->GetRawFormatTimeBins()]; + + const Int_t maxDDL = 20; + AliAltroBuffer *buffer[maxDDL]; + AliAltroMapping *mapping[maxDDL]; + + for(Int_t jDDL=0; jDDLGetEntries(); iDigit++) { - AliPHOSDigit* digit = dynamic_cast(digits->At(iDigit)) ; - if (digit->GetAmp() < kThreshold) + AliPHOSDigit* digit = static_cast(digits->At(iDigit)) ; + + // Skip small energy below treshold + if (digit->GetEnergy() < kThreshold) continue; + Int_t relId[4]; geom->AbsToRelNumbering(digit->GetId(), relId); - Int_t module = relId[0]; + Int_t module = 5-relId[0]; - // Begin FIXME + // Begin FIXME if (relId[1] != 0) continue; // ignore digits from CPV - // End FIXME - - // PHOS EMCA has 4 DDL per module. Splitting is done based on the row number - Int_t iDDL = 4 * (module - 1) + (4 * (relId[2] - 1)) / geom->GetNPhi(); + // End FIXME + + Int_t row = relId[2]-1; + Int_t col = relId[3]-1; + + Int_t iRCU = -111; + + //RCU0 + if(0<=row&&row<16 && 0<=col&&col<56) iRCU=0; + + //RCU1 + if(16<=row&&row<32 && 0<=col&&col<56) iRCU=1; + + //RCU2 + if(32<=row&&row<48 && 0<=col&&col<56) iRCU=2; + + //RCU3 + if(48<=row&&row<64 && 0<=col&&col<56) iRCU=3; + + + // PHOS EMCA has 4 DDL per module. Splitting is based on the (row,column) numbers. + // here module already in PHOS online convention: 0Flush(); - buffer->WriteDataHeader(kFALSE, kFALSE); - delete buffer; - } - - // open new file and write dummy header - TString fileName("PHOS_") ; - fileName += (iDDL + kDDLOffset) ; - fileName += ".ddl" ; - buffer = new AliAltroBuffer(fileName.Data(), 1); - buffer->WriteDataHeader(kTRUE, kFALSE); //Dummy; + if (buffer[iDDL] == 0) { + // open new file and write dummy header + TString fileName = AliDAQ::DdlFileName("PHOS",iDDL); + mapping[iDDL] = (AliAltroMapping*)maps->At(iDDL); + buffer[iDDL] = new AliAltroBuffer(fileName.Data(),mapping[iDDL]); + buffer[iDDL]->WriteDataHeader(kTRUE, kFALSE); //Dummy; + } prevDDL = iDDL; } - // 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, relId[3], relId[2], module); // trailer + AliDebug(2,Form("digit E=%.4f GeV, t=%g s, (mod,col,row)=(%d,%d,%d)\n", + digit->GetEnergy(),digit->GetTimeR(), + relId[0]-1,relId[3]-1,relId[2]-1)); + // if a signal is out of time range, write only trailer + if (digit->GetTimeR() > pulse->GetRawFormatTimeMax()*0.5 ) { + AliDebug(2,"Signal is out of time range.\n"); + buffer[iDDL]->FillBuffer(0); + buffer[iDDL]->FillBuffer(pulse->GetRawFormatTimeBins() ); // time bin + buffer[iDDL]->FillBuffer(3); // bunch length + buffer[iDDL]->WriteTrailer(3, relId[3]-1, relId[2]-1, 0); // trailer // calculate the time response function } else { - Double_t energy = 0 ; - if ( digit->GetId() <= geom->GetNModules() * geom->GetNCristalsInModule()) - energy = digit->GetAmp() * digitizer->GetEMCchannel() + digitizer->GetEMCpedestal() ; - else - energy = digit->GetAmp() * digitizer->GetCPVchannel() + digitizer->GetCPVpedestal() ; - - Bool_t lowgain = RawSampledResponse(digit->GetTimeR(), energy, adcValuesHigh, adcValuesLow) ; - - if (lowgain) - buffer->WriteChannel(relId[3], relId[2], module + fLowGainOffset, - GetRawFormatTimeBins(), adcValuesLow, kThreshold); - else - buffer->WriteChannel(relId[3], relId[2], module, - GetRawFormatTimeBins(), adcValuesHigh, kThreshold); - + Double_t energy = 0 ; + if (digit->GetId() <= geom->GetNModules() * geom->GetNCristalsInModule()) { + energy=digit->GetEnergy(); + if(energy>eMax) {eMax=energy; modMax=relId[0]; colMax=col; rowMax=row;} + } + else { + energy = 0; // CPV raw data format is now know yet + } + pulse->SetAmplitude(energy); + pulse->SetTZero(digit->GetTimeR()); + Double_t r =fgCalibData->GetHighLowRatioEmc(relId[0],relId[3],relId[2]) ; + pulse->SetHG2LGRatio(r) ; + pulse->MakeSamples(); + pulse->GetSamples(adcValuesHigh, adcValuesLow) ; + + buffer[iDDL]->WriteChannel(relId[3]-1, relId[2]-1, 0, + pulse->GetRawFormatTimeBins(), adcValuesLow , kAdcThreshold); + buffer[iDDL]->WriteChannel(relId[3]-1, relId[2]-1, 1, + pulse->GetRawFormatTimeBins(), adcValuesHigh, kAdcThreshold); } } - + delete [] adcValuesLow; + delete [] adcValuesHigh; + // write real header and close last file - if (buffer) { - buffer->Flush(); - buffer->WriteDataHeader(kFALSE, kFALSE); - delete buffer; + for (Int_t iDDL=0; iDDLFlush(); + buffer[iDDL]->WriteDataHeader(kFALSE, kFALSE); + delete buffer[iDDL]; + //if (mapping[iDDL]) delete mapping[iDDL]; + } } + AliDebug(1,Form("Digit with max. energy: modMax %d colMax %d rowMax %d eMax %f\n", + modMax,colMax,rowMax,eMax)); + + delete pulse; loader->UnloadDigits(); } @@ -508,9 +560,11 @@ void AliPHOS::Hits2SDigits() AliPHOSSDigitizer phosDigitizer(fLoader->GetRunLoader()->GetFileName().Data()) ; phosDigitizer.SetEventRange(0, -1) ; // do all the events + phosDigitizer.ExecuteTask("all") ; } + //____________________________________________________________________________ AliLoader* AliPHOS::MakeLoader(const char* topfoldername) { @@ -520,82 +574,15 @@ AliLoader* AliPHOS::MakeLoader(const char* topfoldername) return fLoader; } -//__________________________________________________________________ -Double_t AliPHOS::RawResponseFunction(Double_t *x, Double_t *par) -{ - // 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 - - 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)) ; - } - return signal ; -} - -//__________________________________________________________________ -Double_t AliPHOS::RawResponseFunctionMax(Double_t charge, Double_t gain) -{ - return ( charge * TMath::Power(fgOrder, fgOrder) * TMath::Power(gain, fgOrder) - / ( fgCapa * TMath::Exp(fgOrder) ) ); - -} - -//__________________________________________________________________ -Bool_t AliPHOS::RawSampledResponse(Double_t dtime, 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 AliPHOS::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 ; - } - 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 ; -} - //____________________________________________________________________________ void AliPHOS::SetTreeAddress() { // Links Hits in the Tree to Hits array TBranch *branch; char branchname[20]; - sprintf(branchname,"%s",GetName()); + snprintf(branchname,20,"%s",GetName()); // Branch address for hit tree - TTree *treeH = TreeH(); + TTree *treeH = fLoader->TreeH(); if (treeH) { branch = treeH->GetBranch(branchname); if (branch) @@ -607,3 +594,52 @@ void AliPHOS::SetTreeAddress() } } +//____________________________________________________________________________ +Bool_t AliPHOS::Raw2SDigits(AliRawReader* rawReader) +{ + + AliPHOSLoader * loader = static_cast(fLoader) ; + + TTree * tree = 0 ; + tree = loader->TreeS() ; + if ( !tree ) { + loader->MakeTree("S"); + tree = loader->TreeS() ; + } + + TClonesArray * sdigits = loader->SDigits() ; + if(!sdigits) { + loader->MakeSDigitsArray(); + sdigits = loader->SDigits(); + } + sdigits->Clear(); + + rawReader->Reset() ; + + const TObjArray* maps = AliPHOSRecoParam::GetMappings(); + if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!"); + + AliAltroMapping *mapping[20]; + for(Int_t i = 0; i < 20; i++) { + mapping[i] = (AliAltroMapping*)maps->At(i); + } + + AliPHOSRawFitterv0 fitter; + + fitter.SubtractPedestals(AliPHOSSimParam::GetInstance()->EMCSubtractPedestals()); + fitter.SetAmpOffset(AliPHOSSimParam::GetInstance()->GetGlobalAltroOffset()); + fitter.SetAmpThreshold(AliPHOSSimParam::GetInstance()->GetGlobalAltroThreshold()); + + AliPHOSRawDigiProducer pr(rawReader,mapping); + pr.SetSampleQualityCut(AliPHOSSimParam::GetInstance()->GetEMCSampleQualityCut()); + pr.MakeDigits(sdigits,&fitter); + + Int_t bufferSize = 32000 ; + // TBranch * sdigitsBranch = tree->Branch("PHOS",&sdigits,bufferSize); + tree->Branch("PHOS",&sdigits,bufferSize); + tree->Fill(); + + fLoader->WriteSDigits("OVERWRITE"); + return kTRUE; + +}