#include "AliTRDSimParam.h"
#include "AliTRDCommonParam.h"
+#include "Cal/AliTRDCalROC.h"
+#include "Cal/AliTRDCalDet.h"
+
ClassImp(AliTRDdigitizer)
//_____________________________________________________________________________
if (fMasks) {
delete [] fMasks;
- fMasks = 0;
+ fMasks = 0;
}
if (fTimeStruct1) {
fTimeStruct2 = 0;
}
+ if (fGeo) {
+ delete fGeo;
+ fGeo = 0;
+ }
+
}
//_____________________________________________________________________________
// Do not copy timestructs, just invalidate lastvdrift.
// Next time they are requested, they get recalculated
if (target.fTimeStruct1) {
- delete[] target.fTimeStruct1;
+ delete [] target.fTimeStruct1;
target.fTimeStruct1 = 0;
}
if (target.fTimeStruct2) {
- delete[] target.fTimeStruct2;
+ delete [] target.fTimeStruct2;
target.fTimeStruct2 = 0;
}
target.fTimeLastVdrift = -1;
TString optionString = option;
if (optionString.Contains("deb")) {
AliLog::SetClassDebugLevel("AliTRDdigitizer",1);
- AliInfo("Called with debug option\n");
+ AliInfo("Called with debug option");
}
// The AliRoot file is already connected by the manager
AliRunLoader *inrl;
if (gAlice) {
- AliDebug(1,"AliRun object found on file.\n");
+ AliDebug(1,"AliRun object found on file.");
}
else {
inrl = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(0));
inrl->LoadgAlice();
gAlice = inrl->GetAliRun();
if (!gAlice) {
- AliError("Could not find AliRun object.\n")
+ AliError("Could not find AliRun object.")
return;
}
}
// Initialization
//
- AliRunLoader* orl = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
+ AliRunLoader *orl = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
+
if (InitDetector()) {
- AliLoader* ogime = orl->GetLoader("TRDLoader");
+ AliLoader *ogime = orl->GetLoader("TRDLoader");
- TTree* tree = 0;
+ TTree *tree = 0;
if (fSDigits) {
// If we produce SDigits
tree = ogime->TreeS();
for (iInput = 0; iInput < nInput; iInput++) {
- AliDebug(1,Form("Add input stream %d\n",iInput));
+ AliDebug(1,Form("Add input stream %d",iInput));
// Check if the input tree exists
inrl = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(iInput));
}
if (treees == 0x0) {
- AliError(Form("Input stream %d does not exist\n",iInput));
+ AliError(Form("Input stream %d does not exist",iInput));
return;
}
}
// Convert the s-digits to normal digits
- AliDebug(1,"Do the conversion\n");
+ AliDebug(1,"Do the conversion");
SDigits2Digits();
// Store the digits
- AliDebug(1,"Write the digits\n");
+ AliDebug(1,"Write the digits");
fDigitsManager->WriteDigits();
// Write parameters
orl->CdGAFile();
- AliDebug(1,"Done\n");
+ AliDebug(1,"Done");
DeleteSDigitsManager();
fRunLoader = AliRunLoader::Open(file,evfoldname,"UPDATE");
}
if (!fRunLoader) {
- AliError(Form("Can not open session for file %s.",file));
- return kFALSE;
+ AliError(Form("Can not open session for file %s.",file));
+ return kFALSE;
}
if (!fRunLoader->GetAliRun()) {
gAlice = fRunLoader->GetAliRun();
if (gAlice) {
- AliDebug(1,"AliRun object found on file.\n");
+ AliDebug(1,"AliRun object found on file.");
}
else {
- AliError("Could not find AliRun object.\n");
+ AliError("Could not find AliRun object.");
return kFALSE;
}
}
else {
// If we produce Digits
+ tree = loader->TreeD();
+ if (!tree) {
+ loader->MakeTree("D");
+ tree = loader->TreeD();
+ }
+ }
+ return MakeBranch(tree);
+ }
+ else {
+ return kFALSE;
+ }
+
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDdigitizer::Open(AliRunLoader *runLoader, Int_t nEvent)
+{
+ //
+ // Opens a ROOT-file with TRD-hits and reads in the hit-tree
+ //
+ // Connect the AliRoot file containing Geometry, Kine, and Hits
+ //
+
+ fRunLoader = runLoader;
+ if (!fRunLoader) {
+ AliError("RunLoader does not exist");
+ return kFALSE;
+ }
+
+ if (!fRunLoader->GetAliRun()) {
+ fRunLoader->LoadgAlice();
+ }
+ gAlice = fRunLoader->GetAliRun();
+
+ if (gAlice) {
+ AliDebug(1,"AliRun object found on file.");
+ }
+ else {
+ AliError("Could not find AliRun object.");
+ return kFALSE;
+ }
+
+ fEvent = nEvent;
+
+ AliLoader *loader = fRunLoader->GetLoader("TRDLoader");
+ if (!loader) {
+ AliError("Can not get TRD loader from Run Loader");
+ return kFALSE;
+ }
+
+ if (InitDetector()) {
+ TTree *tree = 0;
+ if (fSDigits) {
+ // If we produce SDigits
+ tree = loader->TreeS();
+ if (!tree) {
+ loader->MakeTree("S");
+ tree = loader->TreeS();
+ }
+ }
+ else {
+ // If we produce Digits
+ tree = loader->TreeD();
if (!tree) {
loader->MakeTree("D");
tree = loader->TreeD();
// Get the pointer to the detector class and check for version 1
fTRD = (AliTRD *) gAlice->GetDetector("TRD");
if (!fTRD) {
- AliFatal("No TRD module found\n");
+ AliFatal("No TRD module found");
exit(1);
}
if (fTRD->IsVersion() != 1) {
- AliFatal("TRD must be version 1 (slow simulator)\n");
+ AliFatal("TRD must be version 1 (slow simulator)");
exit(1);
}
// Get the geometry
- fGeo = fTRD->GetGeometry();
- AliDebug(1,Form("Geometry version %d\n",fGeo->IsVersion()));
+ fGeo = new AliTRDgeometry();
// Create a digits manager
delete fDigitsManager;
}
//_____________________________________________________________________________
-Bool_t AliTRDdigitizer::MakeBranch(TTree* tree) const
+Bool_t AliTRDdigitizer::MakeBranch(TTree *tree) const
{
//
// Create the branches for the digits array
// Number of track dictionary arrays
const Int_t kNDict = AliTRDdigitsManager::kNDict;
- // Half the width of the amplification region
- const Float_t kAmWidth = AliTRDgeometry::AmThick() / 2.0;
+ // Width of the amplification region
+ const Float_t kAmWidth = AliTRDgeometry::AmThick();
// Width of the drift region
const Float_t kDrWidth = AliTRDgeometry::DrThick();
+ // Drift + amplification region
+ const Float_t kDrMin = - 0.5 * kAmWidth;
+ const Float_t kDrMax = kDrWidth + 0.5 * kAmWidth;
Int_t iRow;
Int_t iCol;
Int_t timeBinTRFend = 1;
Double_t pos[3];
- Double_t rot[3];
- Double_t xyz[3];
+ Double_t loc[3];
Double_t padSignal[kNpad];
Double_t signalOld[kNpad];
AliTRDpadPlane *padPlane = 0;
+ AliTRDCalROC *calVdriftROC = 0;
+ Float_t calVdriftDetValue = 0.0;
+ AliTRDCalROC *calT0ROC = 0;
+ Float_t calT0DetValue = 0.0;
+ AliTRDCalROC *calGainFactorROC = 0;
+ Float_t calGainFactorDetValue = 0.0;
+
if (!gGeoManager) {
AliFatal("No geometry manager!");
}
+ if (!fGeo) {
+ AliError("No geometry defined");
+ return kFALSE;
+ }
+ fGeo->ReadGeoMatrices();
+
AliTRDSimParam *simParam = AliTRDSimParam::Instance();
if (!simParam) {
- AliError("Could not get simulation parameters\n");
+ AliFatal("Could not get simulation parameters");
return kFALSE;
}
AliTRDCommonParam *commonParam = AliTRDCommonParam::Instance();
if (!commonParam) {
- AliError("Could not get common parameterss\n");
+ AliFatal("Could not get common parameterss");
return kFALSE;
}
-
- // Create a container for the amplitudes
- AliTRDsegmentArray *signalsArray = new AliTRDsegmentArray("AliTRDdataArrayF"
- ,AliTRDgeometry::Ndet());
- AliTRDcalibDB* calibration = AliTRDcalibDB::Instance();
+ AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
if (!calibration) {
- AliError("Could not get calibration object\n");
- return kFALSE;
- }
-
- if (simParam->TRFOn()) {
- timeBinTRFend = ((Int_t) (simParam->GetTRFhi()
- * calibration->GetSamplingFrequency())) - 1;
- AliDebug(1,Form("Sample the TRF up to bin %d\n",timeBinTRFend));
- }
-
- Float_t elAttachProp = simParam->GetElAttachProp() / 100.0;
-
- if (!fGeo) {
- AliError("No geometry defined\n");
+ AliFatal("Could not get calibration object");
return kFALSE;
}
- AliDebug(1,"Start creating digits.\n");
+ AliDebug(1,"Start creating digits");
+
+ // Create a container for the amplitudes
+ AliTRDsegmentArray *signalsArray = new AliTRDsegmentArray("AliTRDdataArrayF"
+ ,AliTRDgeometry::Ndet());
AliLoader *gimme = fRunLoader->GetLoader("TRDLoader");
if (!gimme->TreeH()) {
return kFALSE;
}
fTRD->SetTreeAddress();
-
+
// Get the number of entries in the hit tree
// (Number of primary particles creating a hit somewhere)
Int_t nTrack = (Int_t) hitTree->GetEntries();
- AliDebug(1,Form("Found %d primary particles\n",nTrack));
- AliDebug(1,Form("Sampling = %.0fMHz\n" ,calibration->GetSamplingFrequency()));
- AliDebug(1,Form("Gain = %d\n" ,((Int_t) simParam->GetGasGain())));
- AliDebug(1,Form("Noise = %d\n" ,((Int_t) simParam->GetNoise())));
+ AliDebug(1,Form("Found %d primary particles",nTrack));
+ AliDebug(1,Form("Sampling = %.0fMHz" ,commonParam->GetSamplingFrequency()));
+ AliDebug(1,Form("Gain = %d" ,((Int_t) simParam->GetGasGain())));
+ AliDebug(1,Form("Noise = %d" ,((Int_t) simParam->GetNoise())));
if (simParam->TimeStructOn()) {
- AliDebug(1,"Time Structure of drift cells implemented.\n");
+ AliDebug(1,"Time Structure of drift cells implemented.");
}
else {
- AliDebug(1,"Constant drift velocity in drift cells.\n");
+ AliDebug(1,"Constant drift velocity in drift cells.");
}
-
+ if (simParam->TRFOn()) {
+ timeBinTRFend = ((Int_t) (simParam->GetTRFhi()
+ * commonParam->GetSamplingFrequency())) - 1;
+ AliDebug(1,Form("Sample the TRF up to bin %d",timeBinTRFend));
+ }
+
+ // Get the detector wise calibration objects
+ const AliTRDCalDet *calVdriftDet = calibration->GetVdriftDet();
+ const AliTRDCalDet *calT0Det = calibration->GetT0Det();
+ const AliTRDCalDet *calGainFactorDet = calibration->GetGainFactorDet();
+
Int_t detectorOld = -1;
Int_t countHits = 0;
Int_t nTimeTotal = calibration->GetNumberOfTimeBins();
- Float_t samplingRate = calibration->GetSamplingFrequency();
+ Float_t samplingRate = commonParam->GetSamplingFrequency();
+ Float_t elAttachProp = simParam->GetElAttachProp() / 100.0;
// Loop through all entries in the tree
for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
pos[2] = hit->Z();
Int_t track = hit->Track();
Int_t detector = hit->GetDetector();
+ Float_t hittime = hit->GetTime();
Int_t plane = fGeo->GetPlane(detector);
Int_t chamber = fGeo->GetChamber(detector);
- Float_t time0 = AliTRDgeometry::GetTime0(plane);
- padPlane = commonParam->GetPadPlane(plane,chamber);
- Float_t row0 = padPlane->GetRow0();
+ padPlane = fGeo->GetPadPlane(plane,chamber);
+ Float_t row0 = padPlane->GetRow0ROC();
Int_t nRowMax = padPlane->GetNrows();
Int_t nColMax = padPlane->GetNcols();
Int_t inDrift = 1;
// Find the current volume with the geo manager
- gGeoManager->SetCurrentPoint(pos);
- gGeoManager->FindNode();
- if (strstr(gGeoManager->GetPath(),"/UK")) {
- inDrift = 0;
- }
+ gGeoManager->SetCurrentPoint(pos);
+ gGeoManager->FindNode();
+ if (strstr(gGeoManager->GetPath(),"/UK")) {
+ inDrift = 0;
+ }
if (detector != detectorOld) {
if (fCompress) dictionary[iDict]->Expand();
}
}
+
+ // Get the calibration objects
+ calVdriftROC = calibration->GetVdriftROC(detector);
+ calVdriftDetValue = calVdriftDet->GetValue(detector);
+ calT0ROC = calibration->GetT0ROC(detector);
+ calT0DetValue = calT0Det->GetValue(detector);
+
detectorOld = detector;
+
}
- // Rotate the sectors on top of each other
- // by using the geoManager
- Double_t aaa[3];
- gGeoManager->MasterToLocal(pos,aaa);
+ // Go to the local coordinate system:
+ // loc[0] - col direction in amplification or driftvolume
+ // loc[1] - row direction in amplification or driftvolume
+ // loc[2] - time direction in amplification or driftvolume
+ gGeoManager->MasterToLocal(pos,loc);
if (inDrift) {
- aaa[2] = time0 - (kDrWidth / 2.0 + kAmWidth) + aaa[2];
+ // Relative to middle of amplification region
+ loc[2] = loc[2] - kDrWidth/2.0 - kAmWidth/2.0;
}
- else {
- aaa[2] = time0 + aaa[2];
+
+ // The driftlength [cm] (w/o diffusion yet !).
+ // It is negative if the hit is between pad plane and anode wires.
+ Double_t driftlength = -1.0 * loc[2];
+
+ // Stupid patch to take care of TR photons that are absorbed
+ // outside the chamber volume. A real fix would actually need
+ // a more clever implementation of the TR hit generation
+ if (q < 0.0) {
+ if ((loc[1] < padPlane->GetRowEndROC()) ||
+ (loc[1] > padPlane->GetRow0ROC())) {
+ hit = (AliTRDhit *) fTRD->NextHit();
+ continue;
+ }
+ if ((driftlength < kDrMin) ||
+ (driftlength > kDrMax)) {
+ hit = (AliTRDhit *) fTRD->NextHit();
+ continue;
+ }
+ }
+
+ // Get row and col of unsmeared electron to retrieve drift velocity
+ // The pad row (z-direction)
+ Int_t rowE = padPlane->GetPadRowNumberROC(loc[1]);
+ if (rowE < 0) {
+ hit = (AliTRDhit *) fTRD->NextHit();
+ continue;
+ }
+ Double_t rowOffset = padPlane->GetPadRowOffsetROC(rowE,loc[1]);
+
+ // The pad column (rphi-direction)
+ Double_t offsetTilt = padPlane->GetTiltOffset(rowOffset);
+ Int_t colE = padPlane->GetPadColNumber(loc[0]+offsetTilt);
+ if (colE < 0) {
+ hit = (AliTRDhit *) fTRD->NextHit();
+ continue;
}
- aaa[1] = row0 + padPlane->GetLengthRim() + fGeo->RpadW()
- - 0.5 * fGeo->GetChamberLength(plane,chamber)
- + aaa[1];
- rot[0] = aaa[2];
- rot[1] = aaa[0];
- rot[2] = aaa[1];
+ Double_t colOffset = padPlane->GetPadColOffset(colE,loc[0]+offsetTilt);
- // The driftlength. It is negative if the hit is between pad plane and anode wires.
- Double_t driftlength = time0 - rot[0];
+ Float_t driftvelocity = calVdriftDetValue * calVdriftROC->GetValue(colE,rowE);
+
+ // Normalized drift length
+ Double_t absdriftlength = TMath::Abs(driftlength);
+ if (commonParam->ExBOn()) {
+ absdriftlength /= TMath::Sqrt(GetLorentzFactor(driftvelocity));
+ }
// Loop over all electrons of this hit
// TR photons produce hits with negative charge
Int_t nEl = ((Int_t) TMath::Abs(q));
for (Int_t iEl = 0; iEl < nEl; iEl++) {
- xyz[0] = rot[0];
- xyz[1] = rot[1];
- xyz[2] = rot[2];
-
- // Stupid patch to take care of TR photons that are absorbed
- // outside the chamber volume. A real fix would actually need
- // a more clever implementation of the TR hit generation
- if (q < 0.0) {
- if ((xyz[2] < padPlane->GetRowEnd()) ||
- (xyz[2] > padPlane->GetRow0())) {
- if (iEl == 0) {
- AliDebug(2,Form("Hit outside of sensitive volume, row (z=%f, row0=%f, rowE=%f)\n"
- ,xyz[2],padPlane->GetRow0(),padPlane->GetRowEnd()));
- }
- continue;
- }
- Float_t tt = driftlength + kAmWidth;
- if ((tt < 0.0) ||
- (tt > kDrWidth + 2.0*kAmWidth)) {
- if (iEl == 0) {
- AliDebug(2,Form("Hit outside of sensitive volume, time (Q = %d)\n"
- ,((Int_t) q)));
- }
- continue;
- }
- }
-
- // Get row and col of unsmeared electron to retrieve drift velocity
- // The pad row (z-direction)
- Int_t rowE = padPlane->GetPadRowNumber(xyz[2]);
- if (rowE < 0) continue;
- Double_t rowOffset = padPlane->GetPadRowOffset(rowE,xyz[2]);
-
- // The pad column (rphi-direction)
- Double_t offsetTilt = padPlane->GetTiltOffset(rowOffset);
- Int_t colE = padPlane->GetPadColNumber(xyz[1]+offsetTilt,rowOffset);
- if (colE < 0) continue;
- Double_t colOffset = padPlane->GetPadColOffset(colE,xyz[1]+offsetTilt);
-
- Float_t driftvelocity = calibration->GetVdrift(detector,colE,rowE);
-
- // Normalised drift length
- Double_t absdriftlength = TMath::Abs(driftlength);
- if (commonParam->ExBOn()) {
- absdriftlength /= TMath::Sqrt(GetLorentzFactor(driftvelocity));
- }
+ // Now the real local coordinate system of the ROC
+ // column direction: locC
+ // row direction: locR
+ // time direction: locT
+ // locR and locC are identical to the coordinates of the corresponding
+ // volumina of the drift or amplification region.
+ // locT is defined relative to the wire plane (i.e. middle of amplification
+ // region), meaming locT = 0, and is negative for hits coming from the
+ // drift region.
+ Double_t locC = loc[0];
+ Double_t locR = loc[1];
+ Double_t locT = loc[2];
// Electron attachment
if (simParam->ElAttachOn()) {
- if (gRandom->Rndm() < (absdriftlength * elAttachProp)) continue;
+ if (gRandom->Rndm() < (absdriftlength * elAttachProp)) {
+ continue;
+ }
}
// Apply the diffusion smearing
if (simParam->DiffusionOn()) {
- if (!(Diffusion(driftvelocity,absdriftlength,xyz))) continue;
+ if (!(Diffusion(driftvelocity,absdriftlength,locR,locC,locT))) {
+ continue;
+ }
}
// Apply E x B effects (depends on drift direction)
if (commonParam->ExBOn()) {
- if (!(ExB(driftvelocity,driftlength,xyz))) continue;
+ if (!(ExB(driftvelocity,driftlength,locC))) {
+ continue;
+ }
}
// The electron position after diffusion and ExB in pad coordinates.
// The pad row (z-direction)
- rowE = padPlane->GetPadRowNumber(xyz[2]);
+ rowE = padPlane->GetPadRowNumberROC(locR);
if (rowE < 0) continue;
- rowOffset = padPlane->GetPadRowOffset(rowE,xyz[2]);
+ rowOffset = padPlane->GetPadRowOffsetROC(rowE,locR);
// The pad column (rphi-direction)
offsetTilt = padPlane->GetTiltOffset(rowOffset);
- colE = padPlane->GetPadColNumber(xyz[1]+offsetTilt,rowOffset);
+ colE = padPlane->GetPadColNumber(locC+offsetTilt);
if (colE < 0) continue;
- colOffset = padPlane->GetPadColOffset(colE,xyz[1]+offsetTilt);
+ colOffset = padPlane->GetPadColOffset(colE,locC+offsetTilt);
// Also re-retrieve drift velocity because col and row may have changed
- driftvelocity = calibration->GetVdrift(detector,colE,rowE);
- Float_t t0 = calibration->GetT0(detector,colE,rowE);
-
- // Convert the position to drift time, using either constant drift velocity or
+ driftvelocity = calVdriftDetValue * calVdriftROC->GetValue(colE,rowE);
+ Float_t t0 = calT0DetValue + calT0ROC->GetValue(colE,rowE);
+
+ // Convert the position to drift time [mus], using either constant drift velocity or
// time structure of drift cells (non-isochronity, GARFIELD calculation).
+ // Also add absolute time of hits to take pile-up events into account properly
Double_t drifttime;
if (simParam->TimeStructOn()) {
- // Get z-position with respect to anode wire:
- Double_t Z = row0 - xyz[2] + simParam->GetAnodeWireOffset();
- Z -= ((Int_t)(2 * Z)) / 2.0;
- if (Z > 0.25) {
- Z = 0.5 - Z;
+ // Get z-position with respect to anode wire
+ Double_t zz = row0 - locR + simParam->GetAnodeWireOffset();
+ zz -= ((Int_t)(2 * zz)) / 2.0;
+ if (zz > 0.25) {
+ zz = 0.5 - zz;
}
// Use drift time map (GARFIELD)
- drifttime = TimeStruct(driftvelocity,time0-xyz[0]+kAmWidth,Z);
+ drifttime = TimeStruct(driftvelocity,0.5*kAmWidth-1.0*locT,zz)
+ + hittime;
}
else {
// Use constant drift velocity
- drifttime = TMath::Abs(time0 - xyz[0]) / driftvelocity;
+ drifttime = -1.0 * locT / driftvelocity
+ + hittime;
}
// Apply the gas gain including fluctuations
// in units of pad width
Double_t dist = (colOffset - 0.5*padPlane->GetColSize(colE))
/ padPlane->GetColSize(colE);
+ // This is still the fixed parametrization, i.e. not dependent on
+ // calibration values !!!!
if (!(calibration->PadResponse(signal,dist,plane,padSignal))) continue;
}
else {
padSignal[2] = 0.0;
}
- // The time bin (always positive), with t0 correction
+ // The time bin (always positive), with t0 distortion
Double_t timeBinIdeal = drifttime * samplingRate + t0;
// Protection
if (TMath::Abs(timeBinIdeal) > 2*nTimeTotal) {
// Sample the time response inside the drift region
// + additional time bins before and after.
// The sampling is done always in the middle of the time bin
- for (Int_t iTimeBin = TMath::Max(timeBinTruncated, 0);
- iTimeBin < TMath::Min(timeBinTruncated+timeBinTRFend,nTimeTotal);
- iTimeBin++) {
+ for (Int_t iTimeBin = TMath::Max(timeBinTruncated,0)
+ ;iTimeBin < TMath::Min(timeBinTruncated+timeBinTRFend,nTimeTotal)
+ ;iTimeBin++) {
// Apply the time response
Double_t timeResponse = 1.0;
} // Loop: primary tracks
- AliDebug(1,Form("Finished analyzing %d hits\n",countHits));
+ AliDebug(1,Form("Finished analyzing %d hits",countHits));
+
+ //____________________________________________________________________
+ //
+ // Create (s)digits
+ //____________________________________________________________________
+ //
// The coupling factor
Double_t coupling = simParam->GetPadCoupling()
Int_t plane = fGeo->GetPlane(iDet);
Int_t sector = fGeo->GetSector(iDet);
Int_t chamber = fGeo->GetChamber(iDet);
- Int_t nRowMax = commonParam->GetRowMax(plane,chamber,sector);
- Int_t nColMax = commonParam->GetColMax(plane);
+ Int_t nRowMax = fGeo->GetRowMax(plane,chamber,sector);
+ Int_t nColMax = fGeo->GetColMax(plane);
Double_t *inADC = new Double_t[nTimeTotal];
Double_t *outADC = new Double_t[nTimeTotal];
Int_t nDigits = 0;
// Don't create noise in detectors that are switched off / not installed, etc.
- if (calibration->GetChamberStatus(iDet)) {
+ if (( calibration->IsChamberInstalled(iDet)) &&
+ (!calibration->IsChamberMasked(iDet)) &&
+ ( fGeo->GetSMstatus(sector))) {
+
+ // Get the calibration objects
+ calGainFactorROC = calibration->GetGainFactorROC(iDet);
+ calGainFactorDetValue = calGainFactorDet->GetValue(iDet);
// Create the digits for this chamber
for (iRow = 0; iRow < nRowMax; iRow++ ) {
for (iCol = 0; iCol < nColMax; iCol++ ) {
- // Check whether pad is active / installed / whatever ...
- if (calibration->GetPadStatus(iDet,iCol,iRow)) continue;
- // Check whether MCM is active / installed / whatever ...
- if (calibration->GetMCMStatus(iDet,iCol,iRow)) continue;
+ // Check whether pad is masked
+ // Bridged pads are not considered yet!!!
+ if (calibration->IsPadMasked(iDet,iCol,iRow)) {
+ continue;
+ }
// Create summable digits
if (fSDigits) {
// Create normal digits
else {
+ Float_t padgain = calGainFactorDetValue
+ * calGainFactorROC->GetValue(iCol,iRow);
+ if (padgain <= 0) {
+ AliError(Form("Not a valid gain %f, %d %d %d",padgain,iDet,iCol,iRow));
+ }
+
for (iTime = 0; iTime < nTimeTotal; iTime++) {
Float_t signalAmp = signals->GetDataUnchecked(iRow,iCol,iTime);
// Pad and time coupling
signalAmp *= coupling;
-
- Float_t padgain = calibration->GetGainFactor(iDet,iCol,iRow);
- if (padgain <= 0) {
- AliError(Form("Not a valid gain %f, %d %d %d\n",padgain,iDet,iCol,iRow));
- }
+ // Gain factors
signalAmp *= padgain;
// Add the noise, starting from minus ADC baseline in electrons
- Double_t baselineEl = simParam->GetADCbaseline() * (simParam->GetADCinRange()
- / simParam->GetADCoutRange())
- / convert;
+ Double_t baselineEl = simParam->GetADCbaseline()
+ * (simParam->GetADCinRange() / simParam->GetADCoutRange())
+ / convert;
signalAmp = TMath::Max((Double_t) gRandom->Gaus(signalAmp,simParam->GetNoise())
,-baselineEl);
// Convert to mV
signalAmp *= convert;
// Add ADC baseline in mV
- signalAmp += simParam->GetADCbaseline() * (simParam->GetADCinRange()
- / simParam->GetADCoutRange());
+ signalAmp += simParam->GetADCbaseline()
+ * (simParam->GetADCinRange() / simParam->GetADCoutRange());
// Convert to ADC counts. Set the overflow-bit fADCoutRange if the
// signal is larger than fADCinRange
Int_t adc = 0;
for (iTime = 0; iTime < nTimeTotal; iTime++) {
// Store the amplitude of the digit if above threshold
- if (outADC[iTime] > simParam->GetADCthreshold()) {
+ if (outADC[iTime] > (simParam->GetADCbaseline() + simParam->GetADCthreshold())) {
nDigits++;
digits->SetDataUnchecked(iRow,iCol,iTime,((Int_t) outADC[iTime]));
}
if (nDigits > 0) {
Float_t nPixel = nRowMax * nColMax * nTimeTotal;
- AliDebug(1,Form("Found %d digits in detector %d (%3.0f).\n"
+ AliDebug(1,Form("Found %d digits in detector %d (%3.0f)."
,nDigits,iDet
,100.0 * ((Float_t) nDigits) / nPixel));
}
signalsArray = 0;
}
- AliDebug(1,Form("Total number of analyzed hits = %d\n",countHits));
- AliDebug(1,Form("Total digits data size = %d, %d, %d, %d\n",totalSizeDigits
- ,totalSizeDict0
- ,totalSizeDict1
- ,totalSizeDict2));
+ AliDebug(1,Form("Total number of analyzed hits = %d",countHits));
+ AliDebug(1,Form("Total digits data size = %d, %d, %d, %d",totalSizeDigits
+ ,totalSizeDict0
+ ,totalSizeDict1
+ ,totalSizeDict2));
return kTRUE;
Int_t iCol;
Int_t iTime;
- AliTRDSimParam *simParam = AliTRDSimParam::Instance();
+ AliTRDCalROC *calGainFactorROC = 0;
+ Float_t calGainFactorDetValue = 0.0;
+
+ AliTRDSimParam *simParam = AliTRDSimParam::Instance();
if (!simParam) {
- AliError("Could not get simulation parameters\n");
+ AliFatal("Could not get simulation parameters");
return kFALSE;
}
- AliTRDCommonParam *commonParam = AliTRDCommonParam::Instance();
- if (!commonParam) {
- AliError("Could not get common parameters\n");
- return kFALSE;
- }
-
- AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
+ AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
if (!calibration) {
- AliError("Could not get calibration object\n");
+ AliFatal("Could not get calibration object");
return kFALSE;
}
AliTRDdataArrayI *digitsOut;
AliTRDdataArrayI *dictionaryIn[kNDict];
AliTRDdataArrayI *dictionaryOut[kNDict];
+
+ // Get the detector wise calibration objects
+ const AliTRDCalDet *calGainFactorDet = calibration->GetGainFactorDet();
// Loop through the detectors
for (Int_t iDet = 0; iDet < AliTRDgeometry::Ndet(); iDet++) {
Int_t plane = fGeo->GetPlane(iDet);
Int_t sector = fGeo->GetSector(iDet);
Int_t chamber = fGeo->GetChamber(iDet);
- Int_t nRowMax = commonParam->GetRowMax(plane,chamber,sector);
- Int_t nColMax = commonParam->GetColMax(plane);
+ Int_t nRowMax = fGeo->GetRowMax(plane,chamber,sector);
+ Int_t nColMax = fGeo->GetColMax(plane);
Double_t *inADC = new Double_t[nTimeTotal];
Double_t *outADC = new Double_t[nTimeTotal];
dictionaryOut[iDict]->Allocate(nRowMax,nColMax,nTimeTotal);
}
- for (iRow = 0; iRow < nRowMax; iRow++ ) {
- for (iCol = 0; iCol < nColMax; iCol++ ) {
+ // Don't create noise in detectors that are switched off / not installed, etc.
+ if (( calibration->IsChamberInstalled(iDet)) &&
+ (!calibration->IsChamberMasked(iDet)) &&
+ ( fGeo->GetSMstatus(sector))) {
- for (iTime = 0; iTime < nTimeTotal; iTime++) {
+ // Get the calibration objects
+ calGainFactorROC = calibration->GetGainFactorROC(iDet);
+ calGainFactorDetValue = calGainFactorDet->GetValue(iDet);
- // Scale s-digits to normal digits
- Double_t signal = (Double_t) digitsIn->GetDataUnchecked(iRow,iCol,iTime);
- signal *= sDigitsScale;
- // Apply the pad-by-pad gain factors
- Float_t padgain = calibration->GetGainFactor(iDet,iCol,iRow);
- if (padgain <= 0.0) {
- AliError(Form("Not a valid gain %f, %d %d %d\n",padgain,iDet,iCol,iRow));
- }
- signal *= padgain;
- // Pad and time coupling
- signal *= coupling;
- // Add the noise, starting from minus ADC baseline in electrons
- Double_t baselineEl = adcBaseline * (adcInRange / adcOutRange) / convert;
- signal = TMath::Max((Double_t) gRandom->Gaus(signal,noise),-baselineEl);
- // Convert to mV
- signal *= convert;
- // add ADC baseline in mV
- signal += adcBaseline * (adcInRange / adcOutRange);
- // Convert to ADC counts. Set the overflow-bit adcOutRange if the
- // signal is larger than adcInRange
- Int_t adc = 0;
- if (signal >= adcInRange) {
- adc = ((Int_t) adcOutRange);
- }
- else {
- adc = TMath::Nint(signal * (adcOutRange / adcInRange));
+ for (iRow = 0; iRow < nRowMax; iRow++ ) {
+ for (iCol = 0; iCol < nColMax; iCol++ ) {
+
+ // Check whether pad is masked
+ // Bridged pads are not considered yet!!!
+ if (calibration->IsPadMasked(iDet,iCol,iRow)) {
+ continue;
}
- inADC[iTime] = adc;
- outADC[iTime] = adc;
- }
+ Float_t padgain = calGainFactorDetValue
+ * calGainFactorROC->GetValue(iCol,iRow);
+ if (padgain <= 0) {
+ AliError(Form("Not a valid gain %f, %d %d %d",padgain,iDet,iCol,iRow));
+ }
- for (iTime = 0; iTime < nTimeTotal; iTime++) {
- // Store the amplitude of the digit if above threshold
- if (outADC[iTime] > adcThreshold) {
- digitsOut->SetDataUnchecked(iRow,iCol,iTime,((Int_t) outADC[iTime]));
- // Copy the dictionary
- for (iDict = 0; iDict < kNDict; iDict++) {
- Int_t track = dictionaryIn[iDict]->GetDataUnchecked(iRow,iCol,iTime);
- dictionaryOut[iDict]->SetDataUnchecked(iRow,iCol,iTime,track);
+ for (iTime = 0; iTime < nTimeTotal; iTime++) {
+
+ // Scale s-digits to normal digits
+ Double_t signal = (Double_t) digitsIn->GetDataUnchecked(iRow,iCol,iTime);
+ signal *= sDigitsScale;
+ // Pad and time coupling
+ signal *= coupling;
+ // Gain factors
+ signal *= padgain;
+ // Add the noise, starting from minus ADC baseline in electrons
+ Double_t baselineEl = adcBaseline * (adcInRange / adcOutRange) / convert;
+ signal = TMath::Max((Double_t) gRandom->Gaus(signal,noise),-baselineEl);
+ // Convert to mV
+ signal *= convert;
+ // add ADC baseline in mV
+ signal += adcBaseline * (adcInRange / adcOutRange);
+ // Convert to ADC counts. Set the overflow-bit adcOutRange if the
+ // signal is larger than adcInRange
+ Int_t adc = 0;
+ if (signal >= adcInRange) {
+ adc = ((Int_t) adcOutRange);
+ }
+ else {
+ adc = TMath::Nint(signal * (adcOutRange / adcInRange));
+ }
+ inADC[iTime] = adc;
+ outADC[iTime] = adc;
+
+ }
+
+ for (iTime = 0; iTime < nTimeTotal; iTime++) {
+ // Store the amplitude of the digit if above threshold
+ if (outADC[iTime] > (adcBaseline + adcThreshold)) {
+ digitsOut->SetDataUnchecked(iRow,iCol,iTime,((Int_t) outADC[iTime]));
+ // Copy the dictionary
+ for (iDict = 0; iDict < kNDict; iDict++) {
+ Int_t track = dictionaryIn[iDict]->GetDataUnchecked(iRow,iCol,iTime);
+ dictionaryOut[iDict]->SetDataUnchecked(iRow,iCol,iTime,track);
+ }
}
}
- }
+ }
}
+
}
if (fCompress) {
AliTRDSimParam *simParam = AliTRDSimParam::Instance();
if (!simParam) {
- AliError("Could not get simulation parameters\n");
- return kFALSE;
- }
-
- AliTRDCommonParam *commonParam = AliTRDCommonParam::Instance();
- if (!commonParam) {
- AliError("Could not get common parameters\n");
+ AliFatal("Could not get simulation parameters");
return kFALSE;
}
AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
if (!calibration) {
- AliError("Could not get calibration object\n");
+ AliFatal("Could not get calibration object");
return kFALSE;
}
// Get the first s-digits
fSDigitsManager = (AliTRDdigitsManager *) fSDigitsManagerList->First();
if (!fSDigitsManager) {
- AliError("No SDigits manager\n");
+ AliError("No SDigits manager");
return kFALSE;
}
fSDigitsManagerList->After(fSDigitsManager);
if (mergeSDigitsManager) {
- AliDebug(1,Form("Merge %d input files.\n",fSDigitsManagerList->GetSize()));
+ AliDebug(1,Form("Merge %d input files.",fSDigitsManagerList->GetSize()));
}
else {
- AliDebug(1,"Only one input file.\n");
+ AliDebug(1,"Only one input file.");
}
Int_t nTimeTotal = calibration->GetNumberOfTimeBins();
Int_t plane = fGeo->GetPlane(iDet);
Int_t sector = fGeo->GetSector(iDet);
Int_t chamber = fGeo->GetChamber(iDet);
- Int_t nRowMax = commonParam->GetRowMax(plane,chamber,sector);
- Int_t nColMax = commonParam->GetColMax(plane);
+ Int_t nRowMax = fGeo->GetRowMax(plane,chamber,sector);
+ Int_t nColMax = fGeo->GetColMax(plane);
// Loop through the pixels of one detector and add the signals
digitsA = fSDigitsManager->GetDigits(iDet);
if (doMerge) {
- AliDebug(1,Form("Merge detector %d of input no.%d\n",iDet,iMerge+1));
+ AliDebug(1,Form("Merge detector %d of input no.%d",iDet,iMerge+1));
for (Int_t iRow = 0; iRow < nRowMax; iRow++ ) {
for (Int_t iCol = 0; iCol < nColMax; iCol++ ) {
return;
}
- AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
+ AliLoader *loader = fRunLoader->GetLoader("TRDLoader");
if (!loader) {
AliError("Can not get TRD loader from Run Loader");
return;
(r1 > 37) ||
(kz1 < 0) ||
(kz1 > 10)) {
- AliWarning(Form("Indices out of range: dist=%.2f, z=%.2f, r1=%d, kz1=%d\n"
+ AliWarning(Form("Indices out of range: dist=%.2f, z=%.2f, r1=%d, kz1=%d"
,dist,z,r1,kz1));
}
// Drift Time data calculated with Garfield (by C.Lippmann)
//
- // TODO make caching proper, if same timing structure is selected: do not update timestructs!
-
- // Noting to do
+ // Nothing to do
if (vdrift == fTimeLastVdrift) {
return;
}
-
fTimeLastVdrift = vdrift;
// Drift time maps are saved for some drift velocity values (in drift region):
fVDsmp[7] = 2.134;
if (vdrift < fVDsmp[0]) {
- AliWarning(Form("Drift Velocity too small (%.3f<%.3f)\n",vdrift,fVDsmp[0]));
+ AliWarning(Form("Drift Velocity too small (%.3f<%.3f)",vdrift,fVDsmp[0]));
vdrift = fVDsmp[0];
}
else if (vdrift > fVDsmp[7]) {
- AliWarning(Form("Drift Velocity too large (%.3f>%.3f)\n",vdrift,fVDsmp[6]));
+ AliWarning(Form("Drift Velocity too large (%.3f>%.3f)",vdrift,fVDsmp[6]));
vdrift = fVDsmp[7];
}
//
// Recalculates the diffusion parameters
//
+ // The B=0 case is not really included here.
+ // Should be revisited!
+ //
if (vdrift == fDiffLastVdrift) {
return;
AliTRDSimParam *simParam = AliTRDSimParam::Instance();
if (!simParam) {
- AliError("Could not get simulation parameters\n");
+ AliFatal("Could not get simulation parameters");
return;
}
AliTRDCommonParam *commonParam = AliTRDCommonParam::Instance();
if (!commonParam) {
- AliError("Could not get common parameters\n");
+ AliFatal("Could not get common parameters");
return;
}
AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
if (!calibration) {
- AliError("Could not get calibration object\n");
+ AliFatal("Could not get calibration object");
return;
}
- Float_t field = commonParam->GetField();
+ // The magnetic field strength
+ Double_t x[3] = { 0.0, 0.0, 0.0 };
+ Double_t b[3];
+ gAlice->Field(x,b); // b[] is in kilo Gauss
+ Float_t field = b[2] * 0.1; // Tesla
+
fDiffLastVdrift = vdrift;
// DiffusionL
+ p3T[ibT] * vdrift*vdrift*vdrift;
// OmegaTau
- fOmegaTau = calibration->GetOmegaTau(vdrift);
-
+ fOmegaTau = calibration->GetOmegaTau(vdrift,field);
+
// Lorentzfactor
if (commonParam->ExBOn()) {
fLorentzFactor = 1.0 / (1.0 + fOmegaTau*fOmegaTau);
}
//_____________________________________________________________________________
-Int_t AliTRDdigitizer::Diffusion(Float_t vdrift, Double_t driftlength, Double_t *xyz)
+Int_t AliTRDdigitizer::Diffusion(Float_t vdrift, Double_t absdriftlength
+ , Double_t &lRow, Double_t &lCol, Double_t &lTime)
{
//
- // Applies the diffusion smearing to the position of a single electron
+ // Applies the diffusion smearing to the position of a single electron.
+ // Depends on absolute drift length.
//
RecalcDiffusion(vdrift);
- Float_t driftSqrt = TMath::Sqrt(driftlength);
+ Float_t driftSqrt = TMath::Sqrt(absdriftlength);
Float_t sigmaT = driftSqrt * fDiffusionT;
Float_t sigmaL = driftSqrt * fDiffusionL;
- xyz[0] = gRandom->Gaus(xyz[0],sigmaL * GetLorentzFactor(vdrift));
- xyz[1] = gRandom->Gaus(xyz[1],sigmaT * GetLorentzFactor(vdrift));
- xyz[2] = gRandom->Gaus(xyz[2],sigmaT);
+ lRow = gRandom->Gaus(lRow ,sigmaT);
+ lCol = gRandom->Gaus(lCol ,sigmaT * GetLorentzFactor(vdrift));
+ lTime = gRandom->Gaus(lTime,sigmaL * GetLorentzFactor(vdrift));
return 1;
}
//_____________________________________________________________________________
-Int_t AliTRDdigitizer::ExB(Float_t vdrift, Double_t driftlength, Double_t *xyz)
+Int_t AliTRDdigitizer::ExB(Float_t vdrift, Double_t driftlength, Double_t &lCol)
{
//
- // Applies E x B effects to the position of a single electron
+ // Applies E x B effects to the position of a single electron.
+ // Depends on signed drift length.
//
RecalcDiffusion(vdrift);
-
- xyz[0] = xyz[0];
- xyz[1] = xyz[1] + fOmegaTau * driftlength;
- xyz[2] = xyz[2];
+
+ lCol = lCol + fOmegaTau * driftlength;
return 1;
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff