From 868783814e46938da03e3c64581b108e77506651 Mon Sep 17 00:00:00 2001 From: hdalsgaa Date: Wed, 21 Oct 2009 14:14:45 +0000 Subject: [PATCH] Fixing bugs in FMD reconstruction. Everything should work now --- FMD/AliFMDDigit.h | 30 ++++++++ FMD/AliFMDESDRevertexer.cxx | 26 ++++--- FMD/AliFMDESDRevertexer.h | 136 ++---------------------------------- FMD/AliFMDRawReader.cxx | 63 ++++++++++++----- FMD/AliFMDReconstructor.cxx | 25 ++++--- FMD/FMDrecLinkDef.h | 1 + FMD/libFMDrec.pkg | 3 +- 7 files changed, 115 insertions(+), 169 deletions(-) diff --git a/FMD/AliFMDDigit.h b/FMD/AliFMDDigit.h index 53cec000589..466f2e99249 100644 --- a/FMD/AliFMDDigit.h +++ b/FMD/AliFMDDigit.h @@ -105,6 +105,22 @@ public: * @param c Counts */ void SetCount(UShort_t s, Short_t c); + /** + * Initialize all counts to appropriate values for this oversampling + * rate. That is + * + * @verbatim + * Rate | Sample 1 | Sample 2 | Sample 3 | Sample 4 + * -----+----------+----------+----------+---------- + * 1 | 0 | -1 | -1 | -1 + * 2 | 0 | 0 | -1 | -1 + * 3 | 0 | 0 | 0 | -1 + * 4 | 0 | 0 | 0 | 0 + * @endverbatim + * + * @param rate Oversampling rate + */ + void SetDefaultCounts(UShort_t rate); protected: UShort_t fCount1; // Digital signal Short_t fCount2; // Digital signal (-1 if not used) @@ -113,6 +129,20 @@ protected: ClassDef(AliFMDDigit,2) // Normal FMD digit }; +inline void +AliFMDDigit::SetDefaultCounts(UShort_t rate) +{ + switch (rate) { + case 4: fCount4 = 0; // Fall through + case 3: fCount3 = 0; // Fall through + case 2: fCount2 = 0; // Fall through + case 1: fCount1 = 0; + break; + default: + fCount4 = fCount3 = fCount2 = fCount1 = 0; + break; + } +} inline UShort_t AliFMDDigit::Counts() const { diff --git a/FMD/AliFMDESDRevertexer.cxx b/FMD/AliFMDESDRevertexer.cxx index a0624b67f8e..e993cbddd6a 100644 --- a/FMD/AliFMDESDRevertexer.cxx +++ b/FMD/AliFMDESDRevertexer.cxx @@ -1,18 +1,25 @@ #include #include #include +#include +#include + +ClassImp(AliFMDESDRevertexer) +#if 0 // for emacs +; +#endif //____________________________________________________________________ AliFMDESDRevertexer::AliFMDESDRevertexer() { AliFMDGeometry* geom = AliFMDGeometry::Instance(); geom->Init(); - geom->InitTransforms(); + geom->InitTransformations(); } //____________________________________________________________________ Bool_t -AliFMDESDRevertexer::Revertex(AliESDFMD* fmdEsd, Double_t vz) +AliFMDESDRevertexer::Revertex(AliESDFMD* fmdEsd, Double_t vz) const { if (!fmdEsd) return kFALSE; @@ -24,30 +31,33 @@ AliFMDESDRevertexer::Revertex(AliESDFMD* fmdEsd, Double_t vz) UShort_t nsec = (ir == 0 ? 20 : 40); UShort_t nstr = (ir == 0 ? 512 : 256); for (UShort_t str = 0; str < nstr; str++) { - Double_t phi, r, theta, oldTheta; + Double_t phi, r, theta; Double_t eta = AliESDFMD::kInvalidEta; - Double_t oldEta = fmdEsd->Eta(det, rng, sec, str); + Double_t oldEta = fmdEsd->Eta(det, rng, 0, str); + if (oldEta == AliESDFMD::kInvalidEta) continue; + Double_t oldTheta = Eta2Theta(oldEta); Bool_t ret1 = PhysicalCoordinates(det, rng, 0, str, vz, - eta, phi, r, theta)); - fmdEsd->SetEta(det, rng, sec, str, eta); + eta, phi, r, theta); + fmdEsd->SetEta(det, rng, 0, str, eta); if (!ret1) { // If the was an error, then there's no reason to go on with // this strip-ring. Note, that the eta is correctly set to // AliESDFMD::kInvalidMult. AliWarning(Form("Failed to calculate eta, phi for " - "FMD%d%c[%02d,%03d] with v_z=%9.4f" + "FMD%d%c[%02d,%03d] with v_z=%9.4f", det, rng, 0, str, vz)); ret = kFALSE; continue; } - + Double_t corr = TMath::Abs(TMath::Cos(theta)); if (fmdEsd->IsAngleCorrected()) corr /= TMath::Abs(TMath::Cos(oldTheta)); for (UShort_t sec = 0; sec < nsec; sec++) { Double_t mult = fmdEsd->Multiplicity(det, rng, sec, str); + if (mult == AliESDFMD::kInvalidMult) continue; fmdEsd->SetMultiplicity(det, rng, sec, str, corr * mult); } } diff --git a/FMD/AliFMDESDRevertexer.h b/FMD/AliFMDESDRevertexer.h index 67e8a1fc8f9..5d048e7df8f 100644 --- a/FMD/AliFMDESDRevertexer.h +++ b/FMD/AliFMDESDRevertexer.h @@ -1,7 +1,7 @@ #ifndef ALIFMDESDREVERTEXER_H #define ALIFMDESDREVERTEXER_H -# include -class ::AliESDFMD; +#include +class AliESDFMD; // // Class to recaculate quantities in an AliESDFMD object based on new @@ -11,7 +11,7 @@ class ::AliESDFMD; // z-coordinate of the primary interaction, to recalibrate the signals // in the FMD ESD object, without having to redo the reconstruction. // -class AliFMDESDRevertexer +class AliFMDESDRevertexer : public TObject { public: /** @@ -69,137 +69,9 @@ public: */ Double_t Eta2Theta(Double_t eta) const; protected: + ClassDef(AliFMDESDRevertexer,0) // Revertex and FMD ESD Object. }; -# ifndef __CINT__ -# ifndef ALIFMDGEOMETRY_H -# include -# endif -# ifndef ALIESDFMD_H -# include -# endif - -//____________________________________________________________________ -inline -AliFMDESDRevertexer::AliFMDESDRevertexer() -{ - // Constructor - AliFMDGeometry* geom = AliFMDGeometry::Instance(); - geom->Init(); - geom->InitTransformations(); -} - -//____________________________________________________________________ -inline Bool_t -AliFMDESDRevertexer::Revertex(AliESDFMD* fmdEsd, Double_t vz) const -{ - // Revertex the passed ESD. The passed ESD object will be modified - // directly. - // - // Parameters: - // fmdEsd ESD object to revertex. - // vz New Z coordinate of primary vertex. - // - // Return: - // kTRUE on success, kFALSE failure. - if (!fmdEsd) return kFALSE; - - Bool_t ret = kTRUE; - const UShort_t sec0 = 0; - for (UShort_t det = 1; det <= 3; det++) { - UShort_t nrng = (det == 1 ? 1 : 2); - for (UShort_t ir = 0; ir < nrng; ir++) { - Char_t rng = (ir == 0 ? 'I' : 'O'); - UShort_t nsec = (ir == 0 ? 20 : 40); - UShort_t nstr = (ir == 0 ? 512 : 256); - for (UShort_t str = 0; str < nstr; str++) { - Double_t phi, r, theta; - Double_t eta = AliESDFMD::kInvalidEta; - Double_t oldEta = fmdEsd->Eta(det, rng, 0u, str); - Double_t oldTheta = Eta2Theta(oldEta); - Bool_t ret1 = PhysicalCoordinates(det, rng, sec0, str, vz, - eta, phi, r, theta); - fmdEsd->SetEta(det, rng, sec0, str, eta); - - if (!ret1) { - // If the was an error, then there's no reason to go on with - // this strip-ring. Note, that the eta is correctly set to - // AliESDFMD::kInvalidMult. - ret = kFALSE; - continue; - } - - Double_t corr = 1; - if (fmdEsd->IsAngleCorrected()) - corr = TMath::Abs(TMath::Cos(theta) / TMath::Cos(oldTheta)); - for (UShort_t sec = 0; sec < nsec; sec++) { - Double_t mult = fmdEsd->Multiplicity(det, rng, sec, str); - fmdEsd->SetMultiplicity(det, rng, sec, str, corr * mult); - } - } - } - } - - return ret; -} - -//____________________________________________________________________ -inline Double_t -AliFMDESDRevertexer::Eta2Theta(Double_t eta) const -{ - // Calculate the polar angle @f$ \theta@f$ corresponding to the - // psuedo-rapidity @f$ \eta@f$ - // - // Parameters: - // eta Psuedo rapidity @f$ \eta=-\log[\tan(\theta/2)]@f$ - // Return: - // Polar angle @f$ \theta=2\atan[\exp(-\eta)]@f$ - return 2 * TMath::ATan(TMath::Exp(-eta)); -} - - -//____________________________________________________________________ -inline Bool_t -AliFMDESDRevertexer::PhysicalCoordinates(UShort_t det, - Char_t rng, - UShort_t sec, - UShort_t str, - Double_t vz, - Double_t& eta, - Double_t& phi, - Double_t& r, - Double_t& theta) const -{ - // Calculate the physical coordinates (@a eta, @a phi) corresponding - // to the detector coordinates (@a det, @a rng, @a sec, @a str). - // - // Parameters: - // det The detector identifier - // rng The ring identifier - // sec The sector identifier - // str The strip identifier - // vz The z coordinate of the current primary interation vertex - // eta On return, the psuedo-rapidity - // phi On return, the azimuthal angle - // r On return, the radius - // phi On return, the polar angle - AliFMDGeometry* geom = AliFMDGeometry::Instance(); - Double_t x=0, y=0, z=0; - geom->Detector2XYZ(det, rng, sec, str, x, y, z); - - // Check that the conversion succeeded - if (x == 0 && y == 0 && z == 0) return kFALSE; - - // Correct for vertex offset. - z -= vz; - phi = TMath::ATan2(y, x); - r = TMath::Sqrt(y * y + x * x); - theta = TMath::ATan2(r, z); - eta = -TMath::Log(TMath::Tan(theta / 2)); - - return kTRUE; -} -# endif // __CINT__ #endif // // Local Variables: diff --git a/FMD/AliFMDRawReader.cxx b/FMD/AliFMDRawReader.cxx index 0206986b411..08260c04c34 100644 --- a/FMD/AliFMDRawReader.cxx +++ b/FMD/AliFMDRawReader.cxx @@ -251,9 +251,9 @@ AliFMDRawReader::NewSample(AliAltroRawStreamV3& input, map->Timebin2Strip(sec, t, fPreSamp, fSampleRate[ddl], stroff, samp); str = strbase + stroff; - AliFMDDebug(20, ("0x%04x/0x%03x/%04d maps to strip %3d sample %d " + AliFMDDebug(20, ("0x%04x/0x%03x/%04d=%4d maps to strip %3d sample %d " "(pre: %d, min: %d, max: %d, rate: %d)", - ddl, hwa, t, str, samp, fPreSamp, + ddl, hwa, t, adc, str, samp, fPreSamp, fMinStrip, fMaxStrip, fSampleRate[ddl])); if (str < 0) { AliFMDDebug(10, ("Got presamples at timebin %d", i)); @@ -266,9 +266,9 @@ AliFMDRawReader::NewSample(AliAltroRawStreamV3& input, AliFMDDebug(15, ("Checking if strip %d (%d) in range [%d,%d]", lstrip, str, fMinStrip, fMaxStrip)); if (lstrip < fMinStrip || lstrip > fMaxStrip) { - AliFMDDebug(20, ("Strip %03d-%d (%d,%d) from t=%d out of range (%3d->%3d)", + AliFMDDebug(10, ("Strip %03d-%d (%d,%d) from t=%d out of range (%3d->%3d)", str, samp, lstrip, stroff, t, fMinStrip, fMaxStrip)); - return -1; + adc = -1; } // Possibly do pedestal subtraction of signal if (adc > 1023) @@ -450,16 +450,18 @@ AliFMDRawReader::ReadAdcs(TClonesArray* array) det, ring, sec, str, samp, counts)); // Check the cache of indicies Int_t idx = fSeen(det, ring, sec, str); + AliFMDDigit* digit = 0; if (idx == kUShortMax) { // We haven't seen this strip yet. fSeen(det, ring, sec, str) = idx = array->GetEntriesFast(); AliFMDDebug(7,("making digit for FMD%d%c[%2d,%3d]-%d " "from timebin %4d", det, ring, sec, str, samp, t)); - new ((*array)[idx]) AliFMDDigit(det, ring, sec, str); + digit = new ((*array)[idx]) AliFMDDigit(det, ring, sec, str); + digit->SetDefaultCounts(fSampleRate[ddl]); } - AliFMDBaseDigit* digit = - static_cast(array->At(idx)); + else + digit = static_cast(array->At(idx)); AliFMDDebug(10, ("Setting FMD%d%c[%2d,%3d]-%d " "from timebin %4d=%4d (%4d)", det, ring, sec, str, samp, t, counts, data[i])); @@ -490,8 +492,14 @@ Bool_t AliFMDRawReader::ReadSODevent(AliFMDCalibSampleRate* sampleRate, AliFMDParameters* param = AliFMDParameters::Instance(); AliFMDAltroMapping* map = param->GetAltroMap(); - while(fReader->ReadNextData(fData)) { - + AliAltroRawStreamV3 streamer(fReader); + streamer.Reset(); + streamer.SelectRawData("FMD"); + //fReader->GetDDLID(); + //Int_t detID = fReader->GetDetectorID(); + + // while(fReader->ReadNextData(fData)) { + /* Int_t ddl = fReader->GetDDLID(); Int_t detID = fReader->GetDetectorID(); if (detectors) detectors[map->DDL2Detector(ddl)-1] = kTRUE; @@ -516,10 +524,28 @@ Bool_t AliFMDRawReader::ReadSODevent(AliFMDCalibSampleRate* sampleRate, AliFMDDebug(20, (" # trailer words: %d, # payload words: %d", nTrailerWords, nPayloadWords)); + ULong_t nPayloadWords = streamer.GetSOD...(); + UInt_t payloadWords* = streamer.GetSOD...(); + */ + + + + while (streamer.NextDDL()) { + Int_t ddl = streamer.GetDDLNumber(); + Int_t detID = fReader->GetDetectorID(); + if (detectors) detectors[map->DDL2Detector(ddl)-1] = kTRUE; + AliFMDDebug(0, (" From reader: DDL number is %d , det ID is %d",ddl,detID)); + + ULong_t nPayloadWords = streamer.GetRCUPayloadSizeInSOD(); + UChar_t* payloadData = streamer.GetRCUPayloadInSOD(); + UInt_t* payloadWords = reinterpret_cast(payloadData); + //UInt_t* payloadWords = streamer.GetRCUPayloadInSOD(); + + std::cout<> 21) & 0xf); @@ -544,14 +570,15 @@ Bool_t AliFMDRawReader::ReadSODevent(AliFMDCalibSampleRate* sampleRate, readDataWord = kTRUE; case 0x1: // Fec cmd case 0x2: // Fec write - i++; + i++; + payloadWords++; break; case 0x4: // Loop case 0x5: // Wait break; case 0x6: // End sequence case 0x7: // End Mem - i = nPayloadWords + 1; + i = nPayloadWords + 1; break; default: break; @@ -560,7 +587,7 @@ Bool_t AliFMDRawReader::ReadSODevent(AliFMDCalibSampleRate* sampleRate, //Don't read unless we have a FEC_RD if(!readDataWord) continue; - UInt_t dataWord = Get32bitWord(i); + UInt_t dataWord = *payloadWords;//Get32bitWord(i); UInt_t data = (0xFFFFF & dataWord) ; //UInt_t data = (0xFFFF & dataWord) ; @@ -674,7 +701,7 @@ Bool_t AliFMDRawReader::ReadSODevent(AliFMDCalibSampleRate* sampleRate, } AliFMDDebug(50, ("instruction 0x%x, dataword 0x%x", instruction,dataWord)); - } + } // End of loop over Result memory event UShort_t det,sector; Short_t strip; @@ -747,9 +774,9 @@ Bool_t AliFMDRawReader::ReadSODevent(AliFMDCalibSampleRate* sampleRate, strip_low[boards[i]], strip_high[boards[i]], shift_clk[boards[i]], sample_clk[boards[i]], pulse_size[boards[i]],pulse_length[boards[i]])); - } + } - } + } AliFMDParameters::Instance()->SetSampleRate(sampleRate); AliFMDParameters::Instance()->SetStripRange(stripRange); diff --git a/FMD/AliFMDReconstructor.cxx b/FMD/AliFMDReconstructor.cxx index f6493fdab33..ec353a0f10c 100644 --- a/FMD/AliFMDReconstructor.cxx +++ b/FMD/AliFMDReconstructor.cxx @@ -53,10 +53,7 @@ #include #include #include -// Import revertexer into a private namespace (to prevent conflicts) -namespace { -# include "AliFMDESDRevertexer.h" -} +#include "AliFMDESDRevertexer.h" class AliRawReader; @@ -321,13 +318,14 @@ AliFMDReconstructor::Reconstruct(TTree* digitsTree, // AliFMDDebug(1, ("Reconstructing from digits in a tree")); GetVertex(fESD); - + + static TClonesArray* digits = new TClonesArray("AliFMDDigit"); TBranch *digitBranch = digitsTree->GetBranch("FMD"); if (!digitBranch) { Error("Exec", "No digit branch for the FMD found"); return; } - TClonesArray* digits = new TClonesArray("AliFMDDigit"); + // TClonesArray* digits = new TClonesArray("AliFMDDigit"); digitBranch->SetAddress(&digits); if (fMult) fMult->Clear(); @@ -348,7 +346,7 @@ AliFMDReconstructor::Reconstruct(TTree* digitsTree, Int_t written = clusterTree->Fill(); AliFMDDebug(10, ("Filled %d bytes into cluster tree", written)); digits->Delete(); - delete digits; + // delete digits; } @@ -383,7 +381,8 @@ AliFMDReconstructor::ProcessDigit(AliFMDDigit* digit) const UShort_t sec = digit->Sector(); UShort_t str = digit->Strip(); Short_t adc = digit->Counts(); - + if (adc < 0) + digit->Print(); ProcessSignal(det, rng, sec, str, adc); } @@ -407,7 +406,7 @@ AliFMDReconstructor::ProcessSignal(UShort_t det, AliFMDParameters* param = AliFMDParameters::Instance(); // Check that the strip is not marked as dead if (param->IsDead(det, rng, sec, str)) { - AliFMDDebug(3, ("FMD%d%c[%2d,%3d] is dead", det, rng, sec, str)); + AliFMDDebug(1, ("FMD%d%c[%2d,%3d] is dead", det, rng, sec, str)); return; } @@ -428,6 +427,11 @@ AliFMDReconstructor::ProcessSignal(UShort_t det, // Make rough multiplicity Double_t mult = Energy2Multiplicity(det, rng, sec, str, edep); // Get rid of nonsense mult + if (mult > 20) { + AliWarning(Form("The mutliplicity in FMD%d%c[%2d,%3d]=%f > 20 " + "(ADC: %d, Energy: %f)", det, rng, sec, str, mult, + counts, edep)); + } if (mult < 0) return; AliFMDDebug(10, ("FMD%d%c[%2d,%3d]: " "ADC: %d, Counts: %d, Energy: %f, Mult: %f", @@ -540,7 +544,8 @@ AliFMDReconstructor::SubtractPedestal(UShort_t det, if (counts < noise * nf) counts = 0; if (counts > 0) AliDebugClass(15, "Got a hit strip"); - return counts; + UShort_t ret = counts < 0 ? 0 : counts; + return ret; } diff --git a/FMD/FMDrecLinkDef.h b/FMD/FMDrecLinkDef.h index 00f56d0e823..bf6ee9dd138 100644 --- a/FMD/FMDrecLinkDef.h +++ b/FMD/FMDrecLinkDef.h @@ -26,6 +26,7 @@ #pragma link C++ class AliFMDRawReader+; #pragma link C++ class AliFMDQADataMakerRec+; #pragma link C++ class AliFMDOfflineTrigger+; +#pragma link C++ class AliFMDESDRevertexer+; #else # error Not for compilation diff --git a/FMD/libFMDrec.pkg b/FMD/libFMDrec.pkg index bab58420471..0505030b7b8 100644 --- a/FMD/libFMDrec.pkg +++ b/FMD/libFMDrec.pkg @@ -8,7 +8,8 @@ SRCS = AliFMDReconstructor.cxx \ AliFMDRawReader.cxx \ AliFMDRecPoint.cxx \ AliFMDQADataMakerRec.cxx \ - AliFMDOfflineTrigger.cxx + AliFMDOfflineTrigger.cxx \ + AliFMDESDRevertexer.cxx HDRS = $(SRCS:.cxx=.h) DHDR := FMDrecLinkDef.h -- 2.43.5