X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=FMD%2FAliFMDBaseDigitizer.cxx;h=1970eb348a35e2758d13ff83f6ccc712fcfd4441;hb=2691ce9cc4770afad101bec5dcff0317b931c46d;hp=4bc8ecacf20b56cf423e3f420a66d406733766c6;hpb=8b26caab2fbd5941ca80598b84dc6133427bcb28;p=u%2Fmrichter%2FAliRoot.git diff --git a/FMD/AliFMDBaseDigitizer.cxx b/FMD/AliFMDBaseDigitizer.cxx index 4bc8ecacf20..1970eb348a3 100644 --- a/FMD/AliFMDBaseDigitizer.cxx +++ b/FMD/AliFMDBaseDigitizer.cxx @@ -212,7 +212,9 @@ // #include // ALIRUNDIGITIZER_H //#include // ALIRUN_H #include // ALILOADER_H +#include // ALILOADER_H #include // ALIRUNLOADER_H +#include //==================================================================== ClassImp(AliFMDBaseDigitizer) @@ -222,28 +224,30 @@ ClassImp(AliFMDBaseDigitizer) //____________________________________________________________________ AliFMDBaseDigitizer::AliFMDBaseDigitizer() - : fRunLoader(0), + : fFMD(0), + fRunLoader(0), fEdep(AliFMDMap::kMaxDetectors, AliFMDMap::kMaxRings, AliFMDMap::kMaxSectors, AliFMDMap::kMaxStrips), - fShapingTime(6) + fShapingTime(6), + fStoreTrackRefs(kTRUE) { + AliFMDDebug(1, ("Constructed")); // Default ctor - don't use it } //____________________________________________________________________ AliFMDBaseDigitizer::AliFMDBaseDigitizer(AliRunDigitizer* manager) : AliDigitizer(manager, "AliFMDBaseDigitizer", "FMD Digitizer base class"), + fFMD(0), fRunLoader(0), - fEdep(AliFMDMap::kMaxDetectors, - AliFMDMap::kMaxRings, - AliFMDMap::kMaxSectors, - AliFMDMap::kMaxStrips), - fShapingTime(6) + fEdep(0), // nDet==0 means 51200 slots + fShapingTime(6), + fStoreTrackRefs(kTRUE) { // Normal CTOR - AliFMDDebug(1, (" processed")); + AliFMDDebug(1, ("Constructed")); SetShapingTime(); } @@ -251,15 +255,14 @@ AliFMDBaseDigitizer::AliFMDBaseDigitizer(AliRunDigitizer* manager) AliFMDBaseDigitizer::AliFMDBaseDigitizer(const Char_t* name, const Char_t* title) : AliDigitizer(name, title), + fFMD(0), fRunLoader(0), - fEdep(AliFMDMap::kMaxDetectors, - AliFMDMap::kMaxRings, - AliFMDMap::kMaxSectors, - AliFMDMap::kMaxStrips), - fShapingTime(6) + fEdep(0), // nDet==0 means 51200 slots + fShapingTime(6), + fStoreTrackRefs(kTRUE) { // Normal CTOR - AliFMDDebug(1, (" processed")); + AliFMDDebug(1, (" Constructed")); SetShapingTime(); } @@ -273,132 +276,115 @@ AliFMDBaseDigitizer::~AliFMDBaseDigitizer() Bool_t AliFMDBaseDigitizer::Init() { - // Initialization + // Initialization. Get a pointer to the parameter manager, and + // initialize it. AliFMDParameters::Instance()->Init(); - if (AliLog::GetDebugLevel("FMD","") >= 10) - AliFMDParameters::Instance()->Print("ALL"); + if (AliLog::GetDebugLevel("FMD","") >= 15) + AliFMDParameters::Instance()->Print(""); return kTRUE; } - //____________________________________________________________________ UShort_t -AliFMDBaseDigitizer::MakePedestal(UShort_t, - Char_t, - UShort_t, - UShort_t) const +AliFMDBaseDigitizer::MakePedestal(UShort_t detector, + Char_t ring, + UShort_t sector, + UShort_t strip) const { - // Make a pedestal - return 0; + // Make a pedestal. The pedestal value is drawn from a Gaussian + // distribution. The mean of the distribution is the measured + // pedestal, and the width is the measured noise. + AliFMDParameters* param =AliFMDParameters::Instance(); + Float_t mean =param->GetPedestal(detector,ring,sector,strip); + Float_t width =param->GetPedestalWidth(detector,ring,sector,strip); + return UShort_t(TMath::Max(gRandom->Gaus(mean, width), 0.)); } //____________________________________________________________________ void -AliFMDBaseDigitizer::SumContributions(AliFMD* fmd) +AliFMDBaseDigitizer::AddContribution(UShort_t detector, + Char_t ring, + UShort_t sector, + UShort_t strip, + Float_t edep, + Bool_t isPrimary, + Int_t nTrack, + Int_t* tracknos) { - // Sum energy deposited contributions from each hit in a cache - // (fEdep). - if (!fRunLoader) - Fatal("SumContributions", "no run loader"); - - // Clear array of deposited energies - fEdep.Reset(); - - // Get the FMD loader - AliLoader* inFMD = fRunLoader->GetLoader("FMDLoader"); - // And load the hits - inFMD->LoadHits("READ"); - - // Get the tree of hits - TTree* hitsTree = inFMD->TreeH(); - if (!hitsTree) { - // Try again - inFMD->LoadHits("READ"); - hitsTree = inFMD->TreeH(); + // Add edep contribution from (detector,ring,sector,strip) to cache + AliFMDParameters* param = AliFMDParameters::Instance(); + AliFMDDebug(10, ("Adding contribution %7.5f for FMD%d%c[%2d,%3d] " + " from %d tracks (%s)", + edep, + detector, + ring, + sector, + strip, + nTrack, + (isPrimary ? "primary" : "secondary"))); + // Check if strip is `dead' + if (param->IsDead(detector, ring, sector, strip)) { + AliFMDDebug(5, ("FMD%d%c[%2d,%3d] is marked as dead", + detector, ring, sector, strip)); + return; } + // Check if strip is out-side read-out range + // if (strip < minstrip || strip > maxstrip) { + // AliFMDDebug(5, ("FMD%d%c[%2d,%3d] is outside range [%3d,%3d]", + // detector,ring,sector,strip,minstrip,maxstrip)); + // continue; + // } - // Get the FMD branch - TBranch* hitsBranch = hitsTree->GetBranch("FMD"); - if (hitsBranch) fmd->SetHitsAddressBranch(hitsBranch); - else AliFatal("Branch FMD hit not found"); - - // Get a list of hits from the FMD manager - TClonesArray *fmdHits = fmd->Hits(); - - // Get number of entries in the tree - Int_t ntracks = Int_t(hitsTree->GetEntries()); - - AliFMDParameters* param = AliFMDParameters::Instance(); - Int_t read = 0; - // Loop over the tracks in the - for (Int_t track = 0; track < ntracks; track++) { - // Read in entry number `track' - read += hitsBranch->GetEntry(track); - - // Get the number of hits - Int_t nhits = fmdHits->GetEntries (); - for (Int_t hit = 0; hit < nhits; hit++) { - // Get the hit number `hit' - AliFMDHit* fmdHit = - static_cast(fmdHits->UncheckedAt(hit)); - - // Extract parameters - UShort_t detector = fmdHit->Detector(); - Char_t ring = fmdHit->Ring(); - UShort_t sector = fmdHit->Sector(); - UShort_t strip = fmdHit->Strip(); - Float_t edep = fmdHit->Edep(); - // UShort_t minstrip = param->GetMinStrip(detector, ring, sector, strip); - // UShort_t maxstrip = param->GetMaxStrip(detector, ring, sector, strip); - // Check if strip is `dead' - AliFMDDebug(2, ("Hit in FMD%d%c[%2d,%3d]=%f", - detector, ring, sector, strip, edep)); - if (param->IsDead(detector, ring, sector, strip)) { - AliFMDDebug(1, ("FMD%d%c[%2d,%3d] is marked as dead", - detector, ring, sector, strip)); - continue; - } - // Check if strip is out-side read-out range - // if (strip < minstrip || strip > maxstrip) { - // AliFMDDebug(5, ("FMD%d%c[%2d,%3d] is outside range [%3d,%3d]", - // detector,ring,sector,strip,minstrip,maxstrip)); - // continue; - // } - - // Give warning in case of double hit - if (fEdep(detector, ring, sector, strip).fEdep != 0) - AliFMDDebug(5, ("Double hit in %d%c(%d,%d)", - detector, ring, sector, strip)); + AliFMDEdepHitPair& entry = fEdep(detector, ring, sector, strip); + + // Give warning in case of double sdigit + if (entry.fEdep != 0) + AliFMDDebug(5, ("Double digit in FMD%d%c[%2d,%3d]", + detector, ring, sector, strip)); - // Sum energy deposition - fEdep(detector, ring, sector, strip).fEdep += edep; - fEdep(detector, ring, sector, strip).fN += 1; - // Add this to the energy deposited for this strip - } // hit loop - } // track loop - AliFMDDebug(1, ("Size of cache: %d bytes, read %d bytes", - sizeof(fEdep), read)); + // Sum energy deposition + Int_t oldN = entry.fN; + entry.fEdep += edep; + entry.fN += nTrack; + if (isPrimary) entry.fNPrim += nTrack; + if (fStoreTrackRefs) { + if (entry.fLabels.fN < entry.fN) { + AliFMDDebug(15, ("== New label array size %d, was %d, added %d", + entry.fN, entry.fLabels.fN, nTrack)); + entry.fLabels.Set(entry.fN); + } + for (Int_t i = 0; i < nTrack; i++) { + AliFMDDebug(15, ("=> Setting track label # %d", oldN+i)); + entry.fLabels[oldN + i] = tracknos[i]; + AliFMDDebug(15, ("<= Setting track label # %d", oldN+i)); + } + } + AliFMDDebug(15,("Adding contribution %f to FMD%d%c[%2d,%3d] (%f) track %d", + edep, detector, ring, sector, strip, + entry.fEdep, (nTrack > 0 ? tracknos[0] : -1))); + } //____________________________________________________________________ void -AliFMDBaseDigitizer::DigitizeHits(AliFMD* fmd) const +AliFMDBaseDigitizer::DigitizeHits() const { // For the stored energy contributions in the cache (fEdep), convert // the energy signal to ADC counts, and store the created digit in // the digits array (AliFMD::fDigits) // + AliFMDDebug(5, ("Will now digitize all the summed signals")); AliFMDGeometry* geometry = AliFMDGeometry::Instance(); TArrayI counts(4); for (UShort_t detector=1; detector <= 3; detector++) { - AliFMDDebug(5, ("Processing hits in FMD%d", detector)); + AliFMDDebug(10, ("Processing hits in FMD%d", detector)); // Get pointer to subdetector AliFMDDetector* det = geometry->GetDetector(detector); if (!det) continue; for (UShort_t ringi = 0; ringi <= 1; ringi++) { Char_t ring = ringi == 0 ? 'I' : 'O'; - AliFMDDebug(5, (" Processing hits in FMD%d%c", detector,ring)); + AliFMDDebug(10, (" Processing hits in FMD%d%c", detector,ring)); // Get pointer to Ring AliFMDRing* r = det->GetRing(ring); if (!r) continue; @@ -407,7 +393,7 @@ AliFMDBaseDigitizer::DigitizeHits(AliFMD* fmd) const UShort_t nSectors = UShort_t(360. / r->GetTheta()); // Loop over the number of sectors for (UShort_t sector = 0; sector < nSectors; sector++) { - AliFMDDebug(5, (" Processing hits in FMD%d%c[%2d]", + AliFMDDebug(10, (" Processing hits in FMD%d%c[%2d]", detector,ring,sector)); // Get number of strips UShort_t nStrips = r->GetNStrips(); @@ -420,10 +406,18 @@ AliFMDBaseDigitizer::DigitizeHits(AliFMD* fmd) const // VA1_ALICE channel. if (strip % 128 == 0) last = 0; - Float_t edep = fEdep(detector, ring, sector, strip).fEdep; + const AliFMDEdepHitPair& entry = fEdep(detector,ring,sector,strip); + Float_t edep = entry.fEdep; + UShort_t ntot = entry.fN; + UShort_t nprim = entry.fNPrim; + const TArrayI& labels = entry.fLabels; + if (edep > 0) + AliFMDDebug(15, ("Edep = %f for FMD%d%c[%2d,%3d]", + edep, detector, ring, sector, strip)); ConvertToCount(edep, last, detector, ring, sector, strip, counts); last = edep; + // The following line was introduced - wrongly - by Peter // Hristov. It _will_ break the digitisation and the // following reconstruction. The behviour of the @@ -441,11 +435,12 @@ AliFMDBaseDigitizer::DigitizeHits(AliFMD* fmd) const // here. // // if (edep<=0) continue; - AddDigit(fmd, detector, ring, sector, strip, edep, + AddDigit(detector, ring, sector, strip, edep, UShort_t(counts[0]), Short_t(counts[1]), - Short_t(counts[2]), Short_t(counts[3])); - AliFMDDebug(10, (" Adding digit in FMD%d%c[%2d,%3d]=%d", - detector,ring,sector,strip,counts[0])); + Short_t(counts[2]), Short_t(counts[3]), + ntot, nprim, labels); + AliFMDDebug(15, (" Adding digit in FMD%d%c[%2d,%3d]=%d", + detector,ring,sector,strip,counts[0])); #if 0 // This checks if the digit created will give the `right' // number of particles when reconstructed, using a naiive @@ -508,7 +503,8 @@ AliFMDBaseDigitizer::ConvertToCount(Float_t edep, // = E + (l - E) * ext(-B * t) // AliFMDParameters* param = AliFMDParameters::Instance(); - Float_t convF = 1./param->GetPulseGain(detector,ring,sector,strip); + Float_t convF = (param->GetDACPerMIP() / param->GetEdepMip() * + param->GetPulseGain(detector,ring,sector,strip)); Int_t ped = MakePedestal(detector,ring,sector,strip); Int_t maxAdc = param->GetAltroChannelSize()-1; if (maxAdc < 0) { @@ -516,18 +512,25 @@ AliFMDBaseDigitizer::ConvertToCount(Float_t edep, maxAdc = 1023; } UShort_t rate = param->GetSampleRate(detector,ring,sector,strip); - if (rate < 1 || rate > 4) rate = 1; - + AliFMDDebug(15, ("Sample rate for FMD%d%c[%2d,%3d] = %d", + detector, ring, sector, strip, rate)); + if (rate < 1 || rate > 4) { + AliWarning(Form("Invalid sample rate for for FMD%d%c[%2d,%3d] = %d", + detector, ring, sector, strip, rate)); + rate = 1; + } + // In case we don't oversample, just return the end value. if (rate == 1) { Float_t a = edep * convF + ped; if (a < 0) a = 0; counts[0] = UShort_t(TMath::Min(a, Float_t(maxAdc))); - AliFMDDebug(2, ("FMD%d%c[%2d,%3d]: converting ELoss %f to " + AliFMDDebug(15, ("FMD%d%c[%2d,%3d]: converting ELoss %f to " "ADC %4d (%f,%d)", detector,ring,sector,strip,edep,counts[0],convF,ped)); return; } + // Create a pedestal Float_t b = fShapingTime; @@ -538,9 +541,68 @@ AliFMDBaseDigitizer::ConvertToCount(Float_t edep, if (a < 0) a = 0; counts[i] = UShort_t(TMath::Min(a, Float_t(maxAdc))); } + AliFMDDebug(15, ("Converted edep = %f to ADC (%x,%x,%x,%x) " + "[gain: %f=(%f/%f*%f), pedestal: %d, rate: %d]", + edep, counts[0], counts[1], counts[2], counts[3], + convF, param->GetDACPerMIP(),param->GetEdepMip(), + param->GetPulseGain(detector,ring,sector,strip), + ped, rate)); +} + +//____________________________________________________________________ +void +AliFMDBaseDigitizer::AddDigit(UShort_t detector, + Char_t ring, + UShort_t sector, + UShort_t strip, + Float_t /* edep */, + UShort_t count1, + Short_t count2, + Short_t count3, + Short_t count4, + UShort_t ntot, + UShort_t /* nprim */, + const TArrayI& refs) const +{ + // Add a digit or summable digit + fFMD->AddDigitByFields(detector, ring, sector, strip, + count1, count2, count3, count4, + ntot, fStoreTrackRefs ? refs.fArray : 0); } +//____________________________________________________________________ +TTree* +AliFMDBaseDigitizer::MakeOutputTree(AliLoader* loader) +{ + // Create output tree using loader. If the passed loader differs + // from the currently set loader in the FMD object, reset the FMD + // loader to be the passed loader. This is for the cases wher the + // output is different from the output. + AliFMDDebug(5, ("Making digits tree")); + loader->LoadDigits("UPDATE"); // "RECREATE"); + TTree* out = loader->TreeD(); + if (!out) loader->MakeTree("D"); + out = loader->TreeD(); + if (out) { + out->Reset(); + if (loader != fFMD->GetLoader()) + fFMD->SetLoader(loader); + fFMD->MakeBranch("D"); + } + return out; +} +//____________________________________________________________________ +void +AliFMDBaseDigitizer::StoreDigits(AliLoader* loader) +{ + // Write the digits to disk + AliFMDDebug(5, ("Storing %d digits", fFMD->Digits()->GetEntries())); + loader->WriteDigits("OVERWRITE"); + loader->UnloadDigits(); + // Reset the digits in the AliFMD object + fFMD->ResetDigits(); +} //____________________________________________________________________ //