From faf80567e50c3b411106dcc26a9d3e95b6a1fa1a Mon Sep 17 00:00:00 2001 From: cholm Date: Thu, 22 Jan 2009 16:07:40 +0000 Subject: [PATCH] A lot of changes after detector review: * AliFMDBaseDigit (and AliFMDDigit and AliFMDSDigit by extension) derive from silly AliDigit. * Up to 3 track labels propegated from hits to digits and sdigits. * Implemented AliFMD::Raw2SDigits (see also scripts/TestRaw2SDigits.C). * Added survey data of FMD1 and 2. * Added FMD2 survey markers to geomtry. * Survey to alignment is almost done. The code is in AliFMDSurveyToAlignObjs. Note, FMD2 survey can only give FMD2 z-offset, while FMD1 survey can in priniple also give x,y offsets as well as rotations. FMD2 z-offset is determined as a best-fit of a plane to the available data. * FMD2 and 3 are still roughly 4mm off from design. I need to talk to ITS people about this. --- FMD/AliFMD.cxx | 67 ++++++++-- FMD/AliFMD.h | 45 ++++--- FMD/AliFMDBackgroundCorrection.cxx | 180 ++++++++++++------------- FMD/AliFMDBaseDigit.cxx | 57 +++++++- FMD/AliFMDBaseDigit.h | 127 ++++++++++++++---- FMD/AliFMDBaseDigitizer.cxx | 27 ++-- FMD/AliFMDBaseDigitizer.h | 3 +- FMD/AliFMDDigit.cxx | 36 +++-- FMD/AliFMDDigit.h | 108 ++++++++++----- FMD/AliFMDDigitizer.cxx | 16 ++- FMD/AliFMDGeometryBuilder.cxx | 49 +++++++ FMD/AliFMDHitDigitizer.cxx | 11 +- FMD/AliFMDInput.cxx | 111 +++++++++++----- FMD/AliFMDInput.h | 3 +- FMD/AliFMDRawReader.cxx | 130 ++++++++++++------ FMD/AliFMDRawReader.h | 77 +++++++++-- FMD/AliFMDReconstructor.cxx | 203 +++++++++++++++++++++++++---- FMD/AliFMDReconstructor.h | 76 +++++++++++ FMD/AliFMDSDigit.cxx | 28 +++- FMD/AliFMDSDigit.h | 111 ++++++++++++---- FMD/AliFMDStripIndex.h | 53 ++++++-- FMD/AliFMDSurveyToAlignObjs.cxx | 178 ++++++++++++++++++++----- FMD/AliFMDSurveyToAlignObjs.h | 13 +- FMD/Config.C | 89 +------------ FMD/Reconstruct.C | 14 +- FMD/Simulate.C | 4 +- FMD/Survey_943928_FMD.txt | 49 +++++++ FMD/Survey_976326_FMD.txt | 51 ++++++++ FMD/scripts/DrawBothDigits.C | 110 ++++++++++++++++ FMD/scripts/DrawDigits.C | 2 + FMD/scripts/DrawESD.C | 83 +++++++----- FMD/scripts/DrawHits.C | 43 +++--- FMD/scripts/DrawHitsDigits.C | 7 + FMD/scripts/DrawHitsSDigits.C | 121 +++++++++++++++++ FMD/scripts/DrawSDigits.C | 3 +- FMD/scripts/TestRaw2SDigits.C | 10 ++ 36 files changed, 1753 insertions(+), 542 deletions(-) create mode 100644 FMD/Survey_943928_FMD.txt create mode 100644 FMD/Survey_976326_FMD.txt create mode 100644 FMD/scripts/DrawBothDigits.C create mode 100644 FMD/scripts/DrawHitsSDigits.C create mode 100644 FMD/scripts/TestRaw2SDigits.C diff --git a/FMD/AliFMD.cxx b/FMD/AliFMD.cxx index 0f404539c22..8e778889eba 100644 --- a/FMD/AliFMD.cxx +++ b/FMD/AliFMD.cxx @@ -115,8 +115,11 @@ //#endif // #include "AliFMDGeometryBuilder.h" #include "AliFMDRawWriter.h" // ALIFMDRAWWRITER_H +#include "AliFMDRawReader.h" // ALIFMDRAWREADER_H #include "AliTrackReference.h" #include "AliFMDStripIndex.h" +#include "AliFMDParameters.h" +#include "AliFMDReconstructor.h" //____________________________________________________________________ ClassImp(AliFMD) @@ -696,14 +699,15 @@ AliFMD::AddDigit(Int_t* digits, Int_t*) //____________________________________________________________________ void -AliFMD::AddDigitByFields(UShort_t detector, - Char_t ring, - UShort_t sector, - UShort_t strip, - UShort_t count1, - Short_t count2, - Short_t count3, - Short_t count4) +AliFMD::AddDigitByFields(UShort_t detector, + Char_t ring, + UShort_t sector, + UShort_t strip, + UShort_t count1, + Short_t count2, + Short_t count3, + Short_t count4, + const TArrayI& refs) { // add a real digit - as coming from data // @@ -719,11 +723,13 @@ AliFMD::AddDigitByFields(UShort_t detector, TClonesArray& a = *(DigitsArray()); new (a[fNdigits++]) - AliFMDDigit(detector, ring, sector, strip, count1, count2, count3, count4); - AliFMDDebug(15, ("Adding digit # %5d/%5d for FMD%d%c[%2d,%3d]=(%d,%d,%d,%d)", + AliFMDDigit(detector, ring, sector, strip, + count1, count2, count3, count4, refs); + AliFMDDebug(15, ("Adding digit # %5d/%5d for FMD%d%c[%2d,%3d]" + "=(%d,%d,%d,%d) with %d tracks", fNdigits-1, a.GetEntriesFast(), detector, ring, sector, strip, - count1, count2, count3, count4)); + count1, count2, count3, count4, refs.fN)); } @@ -907,6 +913,45 @@ AliFMD::Digits2Raw() writer.Exec(); } +//==================================================================== +// +// Raw data reading +// +//__________________________________________________________________ +Bool_t +AliFMD::Raw2SDigits(AliRawReader* reader) +{ + // Turn digits into raw data. + // + // This uses the class AliFMDRawWriter to do the job. Please refer + // to that class for more information. + AliFMDParameters::Instance()->Init(); + MakeTree("S"); + MakeBranch("S"); + + TClonesArray* sdigits = SDigits(); + AliFMDReconstructor rec; + + // The two boolean arguments + // Make sdigits instead of digits + // Subtract the pedestal off the signal + rec.Digitize(reader, sdigits); + // + // Bool_t ret = fmdReader.ReadAdcs(sdigits, kTRUE, kTRUE); + // sdigits->ls(); + UShort_t ns = sdigits->GetEntriesFast(); + for (UShort_t i = 0; i < ns; i++) + sdigits->At(i)->Print("pl"); + + AliFMDDebug(1, ("Got a total of %d SDigits", ns)); + + fLoader->TreeS()->Fill(); + ResetSDigits(); + fLoader->WriteSDigits("OVERWRITE"); + + return kTRUE; +} + //==================================================================== // diff --git a/FMD/AliFMD.h b/FMD/AliFMD.h index b1c16407adb..0e3c0372819 100644 --- a/FMD/AliFMD.h +++ b/FMD/AliFMD.h @@ -433,14 +433,15 @@ public: @param count1 ADC count (a 10-bit word) @param count2 ADC count (a 10-bit word), or -1 if not used @param count3 ADC count (a 10-bit word), or -1 if not used */ - virtual void AddDigitByFields(UShort_t detector=0, - Char_t ring='\0', - UShort_t sector=0, - UShort_t strip=0, - UShort_t count1=0, - Short_t count2=-1, - Short_t count3=-1, - Short_t count4=-1); + virtual void AddDigitByFields(UShort_t detector=0, + Char_t ring='\0', + UShort_t sector=0, + UShort_t strip=0, + UShort_t count1=0, + Short_t count2=-1, + Short_t count3=-1, + Short_t count4=-1, + const TArrayI& refs=TArrayI()); /** Add a digit to the Digit tree @param digits - digits[0] [UShort_t] Detector # @@ -495,16 +496,30 @@ public: /** @{ */ /** @name Raw data */ - /** Turn digits into raw data. This uses the class AliFMDRawWriter - to do the job. Please refer to that class for more - information. */ - virtual void Digits2Raw(); + /** + * Turn digits into raw data. This uses the class AliFMDRawWriter to + * do the job. Please refer to that class for more information. + */ + virtual void Digits2Raw(); + /** + * Convert raw data to sdigits + * + * @param reader Raw reader + * + * @return @c true on success + */ + virtual Bool_t Raw2SDigits(AliRawReader* reader); /** @}*/ /** @{ */ - /** @name Utility */ - /** Browse this object - @param b Browser to show this object in */ + /** + * @name Utility + */ + /** + * Browse this object + * + * @param b Browser to show this object in + */ void Browse(TBrowser* b); /** @}*/ protected: diff --git a/FMD/AliFMDBackgroundCorrection.cxx b/FMD/AliFMDBackgroundCorrection.cxx index 94daa71a60d..17a5770cf60 100644 --- a/FMD/AliFMDBackgroundCorrection.cxx +++ b/FMD/AliFMDBackgroundCorrection.cxx @@ -303,37 +303,48 @@ void AliFMDBackgroundCorrection::Simulate(Int_t nEvents) { } //_____________________________________________________________________ -Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::ProcessHit( AliFMDHit* h, TParticle* p) { +Bool_t +AliFMDBackgroundCorrection::AliFMDInputBG::ProcessHit(AliFMDHit* h, + TParticle* /*p*/) +{ if(!h) return kTRUE; - Bool_t retval = ProcessEvent(h->Detector(),h->Ring(),h->Sector(),h->Strip(),h->Track(),h->Q()); + Bool_t retval = ProcessEvent(h->Detector(), + h->Ring(), + h->Sector() + ,h->Strip(), + h->Track(), + h->Q()); return retval; } //_____________________________________________________________________ -Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::ProcessTrackRef( AliTrackReference* tr, TParticle* p) { - +Bool_t +AliFMDBackgroundCorrection::AliFMDInputBG::ProcessTrackRef(AliTrackReference* tr, + TParticle* p) +{ if(!tr) return kTRUE; UShort_t det,sec,strip; Char_t ring; AliFMDStripIndex::Unpack(tr->UserId(),det,ring,sec,strip); - Int_t nTrack = tr->GetTrack(); - TDatabasePDG* pdgDB = TDatabasePDG::Instance(); + Int_t nTrack = tr->GetTrack(); + TDatabasePDG* pdgDB = TDatabasePDG::Instance(); TParticlePDG* pdgPart = pdgDB->GetParticle(p->GetPdgCode()); - Float_t charge = (pdgPart ? pdgPart->Charge() : 0); - Bool_t retval = ProcessEvent(det,ring,sec,strip,nTrack,charge); + Float_t charge = (pdgPart ? pdgPart->Charge() : 0); + Bool_t retval = ProcessEvent(det,ring,sec,strip,nTrack,charge); return retval; } //_____________________________________________________________________ -Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::ProcessEvent(UShort_t det, - Char_t ring, - UShort_t sec, - UShort_t strip, - Int_t nTrack, - Float_t charge) +Bool_t +AliFMDBackgroundCorrection::AliFMDInputBG::ProcessEvent(UShort_t det, + Char_t ring, + UShort_t sec, + UShort_t strip, + Int_t nTrack, + Float_t charge) { AliGenEventHeader* genHeader = fLoader->GetHeader()->GenEventHeader(); @@ -343,12 +354,14 @@ Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::ProcessEvent(UShort_t det, if(TMath::Abs(vertex.At(2)) > fZvtxCut) return kTRUE; - Double_t delta = 2*fZvtxCut/fNvtxBins; + Double_t delta = 2 * fZvtxCut / fNvtxBins; Double_t vertexBinDouble = (vertex.At(2) + fZvtxCut) / delta; - Int_t vertexBin = (Int_t)vertexBinDouble; + Int_t vertexBin = Int_t(vertexBinDouble); - if(charge != 0 && (nTrack != fPrevTrack || det != fPrevDetector || ring != fPrevRing)) { - + if(charge != 0 && + (nTrack != fPrevTrack || + det != fPrevDetector || + ring != fPrevRing)) { Double_t x,y,z; AliFMDGeometry* fmdgeo = AliFMDGeometry::Instance(); fmdgeo->Detector2XYZ(det,ring,sec,strip,x,y,z); @@ -359,11 +372,11 @@ Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::ProcessEvent(UShort_t det, TObjArray* vtxArray = (TObjArray*)detArray->At(iring); TH2F* hHits = (TH2F*)vtxArray->At(vertexBin); - Float_t phi = TMath::ATan2(y,x); - if(phi<0) - phi = phi+2*TMath::Pi(); - Float_t theta = TMath::ATan2(TMath::Sqrt(TMath::Power(x,2)+TMath::Power(y,2)),z+vertex.At(2)); - Float_t eta = -1*TMath::Log(TMath::Tan(0.5*theta)); + Float_t phi = TMath::ATan2(y,x); + if(phi<0) phi = phi+2*TMath::Pi(); + Float_t r = TMath::Sqrt(TMath::Power(x,2)+TMath::Power(y,2)); + Float_t theta = TMath::ATan2(r,z+vertex.At(2)); + Float_t eta = -1*TMath::Log(TMath::Tan(0.5*theta)); hHits->Fill(eta,phi); fHits++; fPrevDetector = det; @@ -376,8 +389,8 @@ Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::ProcessEvent(UShort_t det, } //_____________________________________________________________________ -Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::Init() { - +Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::Init() +{ fPrimaryArray.SetOwner(); fPrimaryArray.SetName("FMD_primary"); @@ -401,49 +414,40 @@ Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::Init() { fHitArray.SetOwner(); fHitArray.SetName("FMD_hits"); - for(Int_t det =1; det<=3;det++) - { - TObjArray* detArrayHits = new TObjArray(); - detArrayHits->SetName(Form("FMD%d",det)); - fHitArray.AddAtAndExpand(detArrayHits,det); - Int_t nRings = (det==1 ? 1 : 2); - for(Int_t ring = 0;ringSetName(Form("FMD%d%c",det,ringChar)); - detArrayHits->AddAtAndExpand(vtxArrayHits,ring); - for(Int_t v=0; vSumw2(); - vtxArrayHits->AddAtAndExpand(hHits,v); - - } - } - + for(Int_t det =1; det<=3;det++) { + TObjArray* detArrayHits = new TObjArray(); + detArrayHits->SetName(Form("FMD%d",det)); + fHitArray.AddAtAndExpand(detArrayHits,det); + Int_t nRings = (det==1 ? 1 : 2); + for(Int_t ring = 0;ringSetName(Form("FMD%d%c",det,ringChar)); + detArrayHits->AddAtAndExpand(vtxArrayHits,ring); + for(Int_t v=0; vSumw2(); + vtxArrayHits->AddAtAndExpand(hHits,v); + + } } + } - - AliFMDInput::Init(); - return kTRUE; } //_____________________________________________________________________ -Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::Begin(Int_t event ) { +Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::Begin(Int_t event ) +{ - Bool_t retVal = AliFMDInput::Begin(event); - - AliStack* partStack=fLoader->Stack(); - - Int_t nTracks=partStack->GetNtrack(); + Bool_t retVal = AliFMDInput::Begin(event); + AliStack* partStack = fLoader->Stack(); + Int_t nTracks = partStack->GetNtrack(); AliGenEventHeader* genHeader = fLoader->GetHeader()->GenEventHeader(); TArrayF vertex; genHeader->PrimaryVertex(vertex); @@ -451,37 +455,35 @@ Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::Begin(Int_t event ) { if(TMath::Abs(vertex.At(2)) > fZvtxCut) return kTRUE; - Double_t delta = 2*fZvtxCut/fNvtxBins; + Double_t delta = 2*fZvtxCut/fNvtxBins; Double_t vertexBinDouble = (vertex.At(2) + fZvtxCut) / delta; - Int_t vertexBin = (Int_t)vertexBinDouble; - - for(Int_t j=0;jParticle(j); - TDatabasePDG* pdgDB = TDatabasePDG::Instance(); - TParticlePDG* pdgPart = pdgDB->GetParticle(p->GetPdgCode()); - Float_t charge = (pdgPart ? pdgPart->Charge() : 0); - Float_t phi = TMath::ATan2(p->Py(),p->Px()); - - if(phi<0) - phi = phi+2*TMath::Pi(); - if(p->Theta() == 0) continue; - Float_t eta = -1*TMath::Log(TMath::Tan(0.5*p->Theta())); - - Bool_t primary = partStack->IsPhysicalPrimary(j); - //(charge!=0)&&(TMath::Abs(p->Vx() - vertex.At(0))<0.01)&&(TMath::Abs(p->Vy() - vertex.At(1))<0.01)&&(TMath::Abs(p->Vz() - vertex.At(2))<0.01); - if(primary && charge !=0) { - - fPrim++; - TObjArray* innerArray = (TObjArray*)fPrimaryArray.At(0); - TObjArray* outerArray = (TObjArray*)fPrimaryArray.At(1); - - TH2F* hPrimaryInner = (TH2F*)innerArray->At(vertexBin); - TH2F* hPrimaryOuter = (TH2F*)outerArray->At(vertexBin); - hPrimaryInner->Fill(eta,phi); - hPrimaryOuter->Fill(eta,phi); - } - } + Int_t vertexBin = (Int_t)vertexBinDouble; + + for(Int_t j=0;jParticle(j); + TDatabasePDG* pdgDB = TDatabasePDG::Instance(); + TParticlePDG* pdgPart = pdgDB->GetParticle(p->GetPdgCode()); + Float_t charge = (pdgPart ? pdgPart->Charge() : 0); + Float_t phi = TMath::ATan2(p->Py(),p->Px()); + + if(phi<0) phi = phi+2*TMath::Pi(); + if(p->Theta() == 0) continue; + Float_t eta = -1*TMath::Log(TMath::Tan(0.5*p->Theta())); + + Bool_t primary = partStack->IsPhysicalPrimary(j); + //(charge!=0)&&(TMath::Abs(p->Vx() - vertex.At(0))<0.01)&&(TMath::Abs(p->Vy() - vertex.At(1))<0.01)&&(TMath::Abs(p->Vz() - vertex.At(2))<0.01); + if(primary && charge !=0) { + + fPrim++; + TObjArray* innerArray = (TObjArray*)fPrimaryArray.At(0); + TObjArray* outerArray = (TObjArray*)fPrimaryArray.At(1); + + TH2F* hPrimaryInner = (TH2F*)innerArray->At(vertexBin); + TH2F* hPrimaryOuter = (TH2F*)outerArray->At(vertexBin); + hPrimaryInner->Fill(eta,phi); + hPrimaryOuter->Fill(eta,phi); + } + } return retVal; } diff --git a/FMD/AliFMDBaseDigit.cxx b/FMD/AliFMDBaseDigit.cxx index 5ac0b767e61..41d69200d48 100644 --- a/FMD/AliFMDBaseDigit.cxx +++ b/FMD/AliFMDBaseDigit.cxx @@ -64,6 +64,7 @@ ////////////////////////////////////////////////////////////////////// #include "AliFMDBaseDigit.h" // ALIFMDDIGIT_H +#include "AliFMDStripIndex.h" #include "Riostream.h" // ROOT_Riostream // #include // #include @@ -89,7 +90,8 @@ AliFMDBaseDigit::AliFMDBaseDigit(UShort_t detector, Char_t ring, UShort_t sector, UShort_t strip) - : fDetector(detector), + : AliDigit(), + fDetector(detector), fRing(ring), fSector(sector), fStrip(strip), @@ -106,12 +108,37 @@ AliFMDBaseDigit::AliFMDBaseDigit(UShort_t detector, // strip Strip # (For inner/outer rings: 0-511/0-255) } +//____________________________________________________________________ +AliFMDBaseDigit::AliFMDBaseDigit(Int_t* tracks, + UShort_t detector, + Char_t ring, + UShort_t sector, + UShort_t strip) + : AliDigit(tracks), + fDetector(detector), + fRing(ring), + fSector(sector), + fStrip(strip), + fName("") +{ + // + // Creates a base data digit object + // + // Parameters + // + // tracks Array of 3 track labels + // detector Detector # (1, 2, or 3) + // ring Ring ID ('I' or 'O') + // sector Sector # (For inner/outer rings: 0-19/0-39) + // strip Strip # (For inner/outer rings: 0-511/0-255) +} + //____________________________________________________________________ void AliFMDBaseDigit::Print(Option_t* /* option*/) const { // Print digit to standard out - cout << ClassName() << ": FMD" << GetName() << flush; + cout << ClassName() << ": " << GetName() << flush; } //____________________________________________________________________ @@ -132,9 +159,13 @@ ULong_t AliFMDBaseDigit::Hash() const { // Calculate a hash value based on the detector coordinates. +#if 1 + return AliFMDStripIndex::Pack(fDetector, fRing, fSector, fStrip); +#else size_t ringi = (fRing == 'I' || fRing == 'i' ? 0 : 1); return fStrip + fMaxStrips * (fSector + fMaxSectors * (ringi + fMaxRings * (fDetector - 1))); +#endif } @@ -164,6 +195,28 @@ AliFMDBaseDigit::Compare(const TObject* o) const return 1; } +//____________________________________________________________________ +void +AliFMDBaseDigit::AddTrack(Int_t track) +{ + if (fTracks[0] == -1) fTracks[0] = track; + else if (fTracks[1] == -1) fTracks[1] = track; + else if (fTracks[2] == -1) fTracks[2] = track; + else + AliWarning(Form("While adding track label to %s for %s: " + "All 3 track labels used, cannot add reference to track %d", + ClassName(), GetName(), track)); +} + +//____________________________________________________________________ +UShort_t +AliFMDBaseDigit::GetNTrack() const +{ + for (Int_t i = 3; i > 0; i--) + if (fTracks[i-1] != -1) return i; + return 0; +} + //____________________________________________________________________ // diff --git a/FMD/AliFMDBaseDigit.h b/FMD/AliFMDBaseDigit.h index fe9b7e70d2d..a6d88a7b440 100644 --- a/FMD/AliFMDBaseDigit.h +++ b/FMD/AliFMDBaseDigit.h @@ -13,61 +13,138 @@ // AliFMDDigit - Normal (smeared) digit // AliFMDSDigit - Summable (non-smeared) digit // -#ifndef ROOT_TObject -# include +#ifndef ALIDIGIT_H +# include #endif #ifndef ROOT_TString # include #endif //____________________________________________________________________ -/** @class AliFMDBaseDigit AliFMDDigit.h - @brief base class for digits - @ingroup FMD_base +/** + * @class AliFMDBaseDigit AliFMDDigit.h + * + * @brief base class for digits + * + * @ingroup FMD_base */ -class AliFMDBaseDigit : public TObject +class AliFMDBaseDigit : public AliDigit { public: - /** CTOR */ + /** + * CTOR + */ AliFMDBaseDigit(); - /** Constrctor - @param detector Detector - @param ring Ring - @param sector Sector - @param strip Strip */ + /** + * Constrctor + * + * @param detector Detector + * @param ring Ring + * @param sector Sector + * @param strip Strip + */ AliFMDBaseDigit(UShort_t detector, Char_t ring='\0', UShort_t sector=0, UShort_t strip=0); - /** DTOR */ + /** + * Constrctor + * + * @param tracks Array of 3 track indicies + * @param detector Detector + * @param ring Ring + * @param sector Sector + * @param strip Strip + */ + AliFMDBaseDigit(Int_t* tracks, + UShort_t detector, + Char_t ring='\0', + UShort_t sector=0, + UShort_t strip=0); + /** + * DTOR + */ virtual ~AliFMDBaseDigit() {} - /** @return Detector # */ + /** + * + * @return Detector # + */ UShort_t Detector() const { return fDetector; } - /** @return Ring ID */ + /** + * + * @return Ring ID + */ Char_t Ring() const { return fRing; } - /** @return sector # */ + /** + * + * @return sector # + */ UShort_t Sector() const { return fSector; } - /** @return strip # */ + /** + * + * @return strip # + */ UShort_t Strip() const { return fStrip; } - /** Print information - @param opt Not used */ + /** + * Print information + * + * @param opt Not used + */ virtual void Print(Option_t* opt="") const; - /** @return Name */ + /** + * + * @return Name + */ const char* GetName() const; - /** @param rhs Other digit to compare to - @return -1 if this is less than @a rhs, 0 if the refer to the - same, and 1 if @a rhs is larger than this */ + /** + * @param rhs Other digit to compare to + * + * @return -1 if this is less than @a rhs, 0 if the refer to the + * same, and 1 if @a rhs is larger than this + */ Int_t Compare(const TObject* o) const; - /** @return Always true */ + /** + * + * @return Always true + */ Bool_t IsSortable() const { return kTRUE; } + + /** + * Add a track referenc + * + * @param trackno The track number + */ + void AddTrack(Int_t trackno); + + /** + * Get the number of track references (max 3) + * + * + * @return Number of valid track references. + */ + UShort_t GetNTrack() const; + + /** + * Set the count value + * + * @param s Sample number + * @param c Counts + */ + virtual void SetCount(UShort_t s, Short_t c) = 0; protected: + /** + * Calculate the hash value + * + * + * @return Hash value + */ ULong_t Hash() const; UShort_t fDetector; // (Sub) Detector # (1,2, or 3) Char_t fRing; // Ring ID ('I' or 'O') UShort_t fSector; // Sector # (phi division) UShort_t fStrip; // Strip # (radial division) mutable TString fName; //! Name (cached, but not stored) - ClassDef(AliFMDBaseDigit, 2) // Base class for FMD digits + ClassDef(AliFMDBaseDigit, 3) // Base class for FMD digits }; #endif diff --git a/FMD/AliFMDBaseDigitizer.cxx b/FMD/AliFMDBaseDigitizer.cxx index 780b27c3f1d..36bf25549b7 100644 --- a/FMD/AliFMDBaseDigitizer.cxx +++ b/FMD/AliFMDBaseDigitizer.cxx @@ -231,7 +231,7 @@ AliFMDBaseDigitizer::AliFMDBaseDigitizer() AliFMDMap::kMaxSectors, AliFMDMap::kMaxStrips), fShapingTime(6), - fStoreTrackRefs(kFALSE) + fStoreTrackRefs(kTRUE) { AliFMDDebug(1, ("Constructed")); // Default ctor - don't use it @@ -247,7 +247,7 @@ AliFMDBaseDigitizer::AliFMDBaseDigitizer(AliRunDigitizer* manager) AliFMDMap::kMaxSectors, AliFMDMap::kMaxStrips), fShapingTime(6), - fStoreTrackRefs(kFALSE) + fStoreTrackRefs(kTRUE) { // Normal CTOR AliFMDDebug(1, ("Constructed")); @@ -265,7 +265,7 @@ AliFMDBaseDigitizer::AliFMDBaseDigitizer(const Char_t* name, AliFMDMap::kMaxSectors, AliFMDMap::kMaxStrips), fShapingTime(6), - fStoreTrackRefs(kFALSE) + fStoreTrackRefs(kTRUE) { // Normal CTOR AliFMDDebug(1, (" Constructed")); @@ -314,7 +314,8 @@ AliFMDBaseDigitizer::AddContribution(UShort_t detector, UShort_t strip, Float_t edep, Bool_t isPrimary, - Int_t trackno) + Int_t nTrack, + Int_t* tracknos) { // Add edep contribution from (detector,ring,sector,strip) to cache AliFMDParameters* param = AliFMDParameters::Instance(); @@ -341,15 +342,17 @@ AliFMDBaseDigitizer::AddContribution(UShort_t detector, // Sum energy deposition entry.fEdep += edep; - entry.fN += 1; - if (isPrimary) entry.fNPrim += 1; + entry.fN += nTrack; + if (isPrimary) entry.fNPrim += nTrack; if (fStoreTrackRefs) { + Int_t oldN = entry.fLabels.fN; entry.fLabels.Set(entry.fN); - entry.fLabels[entry.fN-1] = trackno; + for (Int_t i = 0; i < nTrack; i++) + entry.fLabels[oldN + i] = tracknos[i]; } - AliFMDDebug(15, ("Adding contribution %f to FMD%d%c[%2d,%3d] (%f)", + AliFMDDebug(15, ("Adding contribution %f to FMD%d%c[%2d,%3d] (%f) track %d", edep, detector, ring, sector, strip, - entry.fEdep)); + entry.fEdep, (nTrack > 0 ? tracknos[0] : -1))); } @@ -428,7 +431,7 @@ AliFMDBaseDigitizer::DigitizeHits() const 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])); + 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 @@ -550,12 +553,12 @@ AliFMDBaseDigitizer::AddDigit(UShort_t detector, Short_t count4, UShort_t /* ntot */, UShort_t /* nprim */, - const TArrayI& /* refs */) const + const TArrayI& refs) const { // Add a digit or summable digit fFMD->AddDigitByFields(detector, ring, sector, strip, - count1, count2, count3, count4); + count1, count2, count3, count4, refs); } //____________________________________________________________________ diff --git a/FMD/AliFMDBaseDigitizer.h b/FMD/AliFMDBaseDigitizer.h index b5d77f8d9ad..1508eec3646 100644 --- a/FMD/AliFMDBaseDigitizer.h +++ b/FMD/AliFMDBaseDigitizer.h @@ -230,7 +230,8 @@ protected: UShort_t strip, Float_t edep, Bool_t isPrimary, - Int_t trackno); + Int_t nTrackno, + Int_t* tracknos); /** Add a digit to output */ virtual void AddDigit(UShort_t detector, Char_t ring, diff --git a/FMD/AliFMDDigit.cxx b/FMD/AliFMDDigit.cxx index 2de6a9a039a..5941b9c0e51 100644 --- a/FMD/AliFMDDigit.cxx +++ b/FMD/AliFMDDigit.cxx @@ -83,14 +83,15 @@ AliFMDDigit::AliFMDDigit() } //____________________________________________________________________ -AliFMDDigit::AliFMDDigit(UShort_t detector, - Char_t ring, - UShort_t sector, - UShort_t strip, - UShort_t count1, - Short_t count2, - Short_t count3, - Short_t count4) +AliFMDDigit::AliFMDDigit(UShort_t detector, + Char_t ring, + UShort_t sector, + UShort_t strip, + UShort_t count1, + Short_t count2, + Short_t count3, + Short_t count4, + const TArrayI& refs) : AliFMDBaseDigit(detector, ring, sector, strip), fCount1(count1), fCount2(count2), @@ -109,6 +110,7 @@ AliFMDDigit::AliFMDDigit(UShort_t detector, // count1 ADC count (a 10-bit word) // count2 ADC count (a 10-bit word) -1 if not used // count3 ADC count (a 10-bit word) -1 if not used + for (Int_t i = 0; i < refs.fN; i++) AddTrack(refs.fArray[i]); } //____________________________________________________________________ @@ -123,13 +125,23 @@ AliFMDDigit::GetTitle() const //____________________________________________________________________ void -AliFMDDigit::Print(Option_t* /* option*/) const +AliFMDDigit::Print(Option_t* option) const { // Print digit to standard out AliFMDBaseDigit::Print(); - cout << "\t" - << fCount1 << " (" << fCount2 << "," << fCount3 << "," << fCount4 - << ") = " << Counts() << endl; + std::cout << "\t" + << std::setw(4) << fCount1 << " (" + << std::setw(4) << fCount2 << "," + << std::setw(4) << fCount3 << "," + << std::setw(4) << fCount4 << ") = " + << std::setw(4) << Counts() << std::flush; + TString opt(option); + if (opt.Contains("l", TString::kIgnoreCase)) { + std::cout << " "; + for (Int_t i = 0; i < 3; i++) + std::cout << (i == 0 ? "" : ",") << std::setw(5) << fTracks[i]; + } + std::cout << std::endl; } //____________________________________________________________________ diff --git a/FMD/AliFMDDigit.h b/FMD/AliFMDDigit.h index da905fbda17..815f90e5d91 100644 --- a/FMD/AliFMDDigit.h +++ b/FMD/AliFMDDigit.h @@ -15,6 +15,10 @@ #ifndef ALIFMDBASEDIGIT_H # include #endif +#ifndef ROOT_TArrayI +# include +#endif + //____________________________________________________________________ /** @class AliFMDDigit AliFMDDigit.h @@ -26,45 +30,79 @@ class AliFMDDigit : public AliFMDBaseDigit public: /** CTOR */ AliFMDDigit(); - /** Constrctor - @param detector Detector - @param ring Ring - @param sector Sector - @param strip Strip - @param count ADC (first sample) - @param count2 ADC (second sample, or -1 if not used) - @param count3 ADC (third sample, or -1 if not used) */ - AliFMDDigit(UShort_t detector, - Char_t ring='\0', - UShort_t sector=0, - UShort_t strip=0, - UShort_t count=0, - Short_t count2=-1, - Short_t count3=-1, - Short_t count4=-1); - /** DTOR */ + /** + * Constrctor + * + * @param detector Detector + * @param ring Ring + * @param sector Sector + * @param strip Strip + * @param count ADC (first sample) + * @param count2 ADC (second sample, or -1 if not used) + * @param count3 ADC (third sample, or -1 if not used) + * @param refs Track references + */ + AliFMDDigit(UShort_t detector, + Char_t ring='\0', + UShort_t sector=0, + UShort_t strip=0, + UShort_t count=0, + Short_t count2=-1, + Short_t count3=-1, + Short_t count4=-1, + const TArrayI& refs=TArrayI()); + /** + * DTOR + */ virtual ~AliFMDDigit() {} - /** @param i # of sample to get - @return sample # @a i */ + /** + * @param i # of sample to get + * + * @return sample # @a i + */ Int_t Count(UShort_t i=0) const; - /** @return ADC count (first sample) */ - UShort_t Count1() const { return fCount1; } - /** @return ADC count (second sample, or -1 if not used) */ - Short_t Count2() const { return fCount2; } - /** @return ADC count (third sample, or -1 if not used) */ - Short_t Count3() const { return fCount3; } - /** @return ADC count (third sample, or -1 if not used) */ - Short_t Count4() const { return fCount4; } - /** @return Canonical ADC counts */ - UShort_t Counts() const; - /** Print info - @param opt Not used */ + /** + * + * @return ADC count (first sample) + */ + UShort_t Count1() const { return fCount1; } + /** + * + * @return ADC count (second sample, or -1 if not used) + */ + Short_t Count2() const { return fCount2; } + /** + * + * @return ADC count (third sample, or -1 if not used) + */ + Short_t Count3() const { return fCount3; } + /** + * + * @return ADC count (third sample, or -1 if not used) + */ + Short_t Count4() const { return fCount4; } + /** + * + * @return Canonical ADC counts + */ + UShort_t Counts() const; + /** + * Print info + * + * @param opt Not used + */ void Print(Option_t* opt="") const; - /** @return Title */ + /** + * + * @return Title + */ const char* GetTitle() const; - /** Set the count value - @param s Sample number - @param c Counts */ + /** + * Set the count value + * + * @param s Sample number + * @param c Counts + */ void SetCount(UShort_t s, Short_t c); protected: UShort_t fCount1; // Digital signal diff --git a/FMD/AliFMDDigitizer.cxx b/FMD/AliFMDDigitizer.cxx index 91e4ff8270a..c17a2650732 100644 --- a/FMD/AliFMDDigitizer.cxx +++ b/FMD/AliFMDDigitizer.cxx @@ -270,7 +270,7 @@ AliFMDDigitizer::Exec(Option_t*) folderName.Data())); return; } - runLoader->LoadgAlice(); + if (!runLoader->GetAliRun()) runLoader->LoadgAlice(); runLoader->GetAliRun(); } if (!gAlice) { @@ -339,6 +339,7 @@ AliFMDDigitizer::Exec(Option_t*) fFMD->SetTreeAddress(); // Sum contributions from the sdigits + AliFMDDebug(3, ("Will now sum contributions from SDigits")); SumContributions(sdigitsBranch); // Unload the sdigits @@ -398,6 +399,16 @@ AliFMDDigitizer::SumContributions(TBranch* sdigitsBranch) // Get the sdigit number `sdigit' AliFMDSDigit* fmdSDigit = static_cast(fmdSDigits->UncheckedAt(sdigit)); + + AliFMDDebug(5, ("Adding contribution of %d tracks", + fmdSDigit->GetNTrack())); + AliFMDDebug(15, ("Contrib from FMD%d%c[%2d,%3d] (%s) from track %d", + fmdSDigit->Detector(), + fmdSDigit->Ring(), + fmdSDigit->Sector(), + fmdSDigit->Strip(), + fmdSDigit->GetName(), + fmdSDigit->GetTrack(0))); // Extract parameters AddContribution(fmdSDigit->Detector(), @@ -406,7 +417,8 @@ AliFMDDigitizer::SumContributions(TBranch* sdigitsBranch) fmdSDigit->Strip(), fmdSDigit->Edep(), kTRUE, - -1); + fmdSDigit->GetNTrack(), + fmdSDigit->GetTracks()); } // sdigit loop } // event loop diff --git a/FMD/AliFMDGeometryBuilder.cxx b/FMD/AliFMDGeometryBuilder.cxx index ecb994f64f5..85a39f187af 100644 --- a/FMD/AliFMDGeometryBuilder.cxx +++ b/FMD/AliFMDGeometryBuilder.cxx @@ -654,9 +654,12 @@ AliFMDGeometryBuilder::FMD1Geometry(AliFMD1* fmd1, AliFMDRing* r = fmd1->GetInner(); Double_t z = fmd1->GetInnerZ(); + // `Top' or `Outside' master volume TGeoVolume* fmd1TopVolume = new TGeoVolumeAssembly(Form(fgkFMDName, fmd1->GetId(), 'T')); fmd1TopVolume->SetTitle("FMD1 top half"); + + // `Bottom' or `Inside' master volume TGeoVolume* fmd1BotVolume = new TGeoVolumeAssembly(Form(fgkFMDName, fmd1->GetId(), 'B')); fmd1BotVolume->SetTitle("FMD1 bottom half"); @@ -728,11 +731,57 @@ AliFMDGeometryBuilder::FMD1Geometry(AliFMD1* fmd1, TGeoRotation* rot = new TGeoRotation("FMD1 rotatation"); rot->RotateZ(90); TGeoMatrix* matrix = new TGeoCombiTrans("FMD1 trans", 0, 0, z, rot); + AliFMDDebug(5, ("Placing volumes %s and %s in ALIC at z=%f", fmd1TopVolume->GetName(), fmd1BotVolume->GetName(), z)); top->AddNode(fmd1TopVolume, fmd1->GetId(), matrix); top->AddNode(fmd1BotVolume, fmd1->GetId(), matrix); + + + // Survey points on V0A (screw holes for the FMD) + const Double_t icb[] = { +12.700, -21.997, 324.670 }; + const Double_t ict[] = { +12.700, +21.997, 324.670 }; + const Double_t ocb[] = { -12.700, -21.997, 324.670 }; + const Double_t oct[] = { -12.700, +21.997, 324.670 }; + + TGeoTube* surveyShape = new TGeoTube("FMD1_survey_marker", + 0, .2, .001); + + TGeoMatrix* outMat = matrix; +#if 0 + if (gGeoManager->cd("/ALIC_1/F1MT_1")) + outMat = gGeoManager->GetCurrentMatrix(); + else + AliWarning("Couldn't cd to /ALIC_1/F1MT_1"); +#endif + + Double_t loct[3], locb[3]; + outMat->MasterToLocal(oct, loct); + outMat->MasterToLocal(ocb, locb); + TGeoVolume* vOct = new TGeoVolume("V0L_OCT", surveyShape, fPlastic); + TGeoVolume* vOcb = new TGeoVolume("V0L_OCB", surveyShape, fPlastic); + + fmd1TopVolume->AddNode(vOct, 1, new TGeoTranslation(loct[0],loct[1],loct[2])); + fmd1TopVolume->AddNode(vOcb, 1, new TGeoTranslation(locb[0],locb[1],locb[2])); + + + TGeoMatrix* inMat = matrix; +#if 0 + if (gGeoManager->cd("/ALIC_1/F1MT_1")) + inMat = gGeoManager->GetCurrentMatrix(); + else + AliWarning("Couldn't cd to /ALIC_1/F1MT_1"); +#endif + + Double_t lict[3], licb[3]; + inMat->MasterToLocal(ict, lict); + inMat->MasterToLocal(icb, licb); + TGeoVolume* vIct = new TGeoVolume("V0L_ICT", surveyShape, fPlastic); + TGeoVolume* vIcb = new TGeoVolume("V0L_ICB", surveyShape, fPlastic); + fmd1BotVolume->AddNode(vIct, 1, new TGeoTranslation(lict[0],lict[1],lict[2])); + fmd1BotVolume->AddNode(vIcb, 1, new TGeoTranslation(licb[0],licb[1],licb[2])); + return 0; } diff --git a/FMD/AliFMDHitDigitizer.cxx b/FMD/AliFMDHitDigitizer.cxx index 20714d60f4f..50132d81b4e 100644 --- a/FMD/AliFMDHitDigitizer.cxx +++ b/FMD/AliFMDHitDigitizer.cxx @@ -430,7 +430,8 @@ AliFMDHitDigitizer::SumContributions(TBranch* hitsBranch) fmdHit->Strip(), fmdHit->Edep(), isPrimary, - trackno); + 1, + &trackno); } // hit loop } // track loop AliFMDDebug(5, ("Size of cache: %d bytes, read %d bytes", @@ -469,6 +470,9 @@ AliFMDHitDigitizer::AddDigit(UShort_t detector, { // Add a digit or summable digit if (fOutput == kDigits) { + AliFMDDebug(15,("Adding digit for FMD%d%c[%2d,%3d] = (%x,%x,%x,%x)", + detector, ring, sector, strip, + count1, count2, count3, count4)); AliFMDBaseDigitizer::AddDigit(detector, ring, sector, strip, 0, count1, count2, count3, count4, 0, 0, refs); return; @@ -484,9 +488,10 @@ AliFMDHitDigitizer::AddDigit(UShort_t detector, detector, ring, sector, strip)); return; } - AliFMDDebug(15, ("Adding digit for FMD%d%c[%2d,%3d] = (%x,%x,%x,%x) [%d/%d]", + AliFMDDebug(15, ("Adding sdigit for FMD%d%c[%2d,%3d] = " + "(%x,%x,%x,%x) [%d/%d] %d", detector, ring, sector, strip, - count1, count2, count3, count4, nprim, ntotal)); + count1, count2, count3, count4, nprim, ntotal, refs.fN)); fFMD->AddSDigitByFields(detector, ring, sector, strip, edep, count1, count2, count3, count4, ntotal, nprim, refs); diff --git a/FMD/AliFMDInput.cxx b/FMD/AliFMDInput.cxx index a819d920afb..a8056824c60 100644 --- a/FMD/AliFMDInput.cxx +++ b/FMD/AliFMDInput.cxx @@ -93,6 +93,7 @@ AliFMDInput::AliFMDInput() fChainE(0), fArrayE(0), fArrayH(0), + fArrayTR(0), fArrayD(0), fArrayS(0), fArrayR(0), @@ -136,6 +137,7 @@ AliFMDInput::AliFMDInput(const char* gAliceFile) fChainE(0), fArrayE(0), fArrayH(0), + fArrayTR(0), fArrayD(0), fArrayS(0), fArrayR(0), @@ -160,6 +162,7 @@ AliFMDInput::NEvents() const { // Get number of events if (TESTBIT(fTreeMask, kRaw)) return fReader->GetNumberOfEvents(); + if (fChainE) return fChainE->GetEntriesFast(); if (fTreeE) return fTreeE->GetEntries(); return -1; } @@ -175,38 +178,71 @@ AliFMDInput::Init() AliWarning("Already initialized"); return fIsInit; } + Info("Init","Initialising w/mask 0x%04x\n" + "\tHits: %d\n" + "\tKinematics: %d\n" + "\tDigits: %d\n" + "\tSDigits: %d\n" + "\tHeader: %d\n" + "\tRecPoints: %d\n" + "\tESD: %d\n" + "\tRaw: %d\n" + "\tGeometry: %d\n" + "\tTracks: %d\n" + "\tTracksRefs: %d", + fTreeMask, + TESTBIT(fTreeMask, kHits), + TESTBIT(fTreeMask, kKinematics), + TESTBIT(fTreeMask, kDigits), + TESTBIT(fTreeMask, kSDigits), + TESTBIT(fTreeMask, kHeader), + TESTBIT(fTreeMask, kRecPoints), + TESTBIT(fTreeMask, kESD), + TESTBIT(fTreeMask, kRaw), + TESTBIT(fTreeMask, kGeometry), + TESTBIT(fTreeMask, kTracks), + TESTBIT(fTreeMask, kTrackRefs)); // Get the run - if (!TESTBIT(fTreeMask, kRaw) && - !TESTBIT(fTreeMask, kESD)) { - if (fGAliceFile.IsNull()) fGAliceFile = "galice.root"; - // Get the loader - fLoader = AliRunLoader::Open(fGAliceFile.Data(), "Alice", "read"); - if (!fLoader) { - AliError(Form("Coulnd't read the file %s", fGAliceFile.Data())); - return kFALSE; + if (TESTBIT(fTreeMask, kDigits) || + TESTBIT(fTreeMask, kSDigits) || + TESTBIT(fTreeMask, kKinematics) || + TESTBIT(fTreeMask, kTrackRefs) || + TESTBIT(fTreeMask, kTracks) || + TESTBIT(fTreeMask, kHeader)) { + if (!gSystem->FindFile(".:/", fGAliceFile)) { + AliWarning(Form("Cannot find file %s in .:/", fGAliceFile.Data())); } - - if (fLoader->LoadgAlice()) return kFALSE; - fRun = fLoader->GetAliRun(); - - // Get the FMD - fFMD = static_cast(fRun->GetDetector("FMD")); - if (!fFMD) { - AliError("Failed to get detector FMD from loader"); - return kFALSE; - } - - // Get the FMD loader - fFMDLoader = fLoader->GetLoader("FMDLoader"); - if (!fFMDLoader) { - AliError("Failed to get detector FMD loader from loader"); - return kFALSE; - } - if (fLoader->LoadHeader()) { - AliError("Failed to get event header information from loader"); - return kFALSE; + else { + fLoader = AliRunLoader::Open(fGAliceFile.Data(), "Alice", "read"); + if (!fLoader) { + AliError(Form("Coulnd't read the file %s", fGAliceFile.Data())); + return kFALSE; + } + AliInfo(Form("Opened GAlice file %s", fGAliceFile.Data())); + + if (fLoader->LoadgAlice()) return kFALSE; + + fRun = fLoader->GetAliRun(); + + // Get the FMD + fFMD = static_cast(fRun->GetDetector("FMD")); + if (!fFMD) { + AliError("Failed to get detector FMD from loader"); + return kFALSE; + } + + // Get the FMD loader + fFMDLoader = fLoader->GetLoader("FMDLoader"); + if (!fFMDLoader) { + AliError("Failed to get detector FMD loader from loader"); + return kFALSE; + } + if (fLoader->LoadHeader()) { + AliError("Failed to get event header information from loader"); + return kFALSE; + } + fTreeE = fLoader->TreeE(); } - fTreeE = fLoader->TreeE(); } // Optionally, get the ESD files @@ -361,7 +397,10 @@ AliFMDInput::Begin(Int_t event) // Possibly load FMD Sdigit information if (TESTBIT(fTreeMask, kSDigits)) { // AliInfo("Getting FMD summable digits"); - if (!fFMDLoader || fFMDLoader->LoadSDigits("READ")) return kFALSE; + if (!fFMDLoader || fFMDLoader->LoadSDigits("READ")) { + AliWarning("Failed to load SDigits!"); + return kFALSE; + } fTreeS = fFMDLoader->TreeS(); if (!fArrayS) fArrayS = fFMD->SDigits(); } @@ -428,10 +467,10 @@ AliFMDInput::Event() if (!ProcessTrackRefs()) return kFALSE; if (TESTBIT(fTreeMask, kTracks)) if (!ProcessTracks()) return kFALSE; - if (TESTBIT(fTreeMask, kDigits)) - if (!ProcessDigits()) return kFALSE; if (TESTBIT(fTreeMask, kSDigits)) if (!ProcessSDigits()) return kFALSE; + if (TESTBIT(fTreeMask, kDigits)) + if (!ProcessDigits()) return kFALSE; if (TESTBIT(fTreeMask, kRaw)) if (!ProcessRawDigits()) return kFALSE; if (TESTBIT(fTreeMask, kRecPoints)) @@ -606,9 +645,13 @@ AliFMDInput::ProcessSDigits() Int_t nEv = fTreeS->GetEntries(); for (Int_t i = 0; i < nEv; i++) { Int_t sdigitRead = fTreeS->GetEntry(i); - if (sdigitRead <= 0) continue; - Int_t nSdigit = fArrayS->GetEntries(); + if (sdigitRead <= 0) { + AliInfo(Form("Read nothing from tree")); + continue; + } + Int_t nSdigit = fArrayS->GetEntriesFast(); AliFMDDebug(0, ("Got %5d digits for this event", nSdigit)); + AliInfo(Form("Got %5d digits for this event", nSdigit)); if (nSdigit <= 0) continue; for (Int_t j = 0; j < nSdigit; j++) { AliFMDSDigit* sdigit = static_cast(fArrayS->At(j)); diff --git a/FMD/AliFMDInput.h b/FMD/AliFMDInput.h index 355c26dd177..9fbe3dc2096 100644 --- a/FMD/AliFMDInput.h +++ b/FMD/AliFMDInput.h @@ -327,7 +327,8 @@ protected: }; inline Bool_t AliFMDInput::ProcessHit(AliFMDHit*,TParticle*) { return kTRUE; } -inline Bool_t AliFMDInput::ProcessTrackRef(AliTrackReference* trackRef, TParticle* track) { return kTRUE; } +inline Bool_t AliFMDInput::ProcessTrackRef(AliTrackReference*, + TParticle*) { return kTRUE; } inline Bool_t AliFMDInput::ProcessTrack(Int_t,TParticle*, AliFMDHit*) { return kTRUE; } inline Bool_t AliFMDInput::ProcessDigit(AliFMDDigit*) { return kTRUE; } diff --git a/FMD/AliFMDRawReader.cxx b/FMD/AliFMDRawReader.cxx index 332c31991f1..74c44e133a6 100644 --- a/FMD/AliFMDRawReader.cxx +++ b/FMD/AliFMDRawReader.cxx @@ -50,6 +50,7 @@ #include "AliFMDDebug.h" // Better debug macros #include "AliFMDParameters.h" // ALIFMDPARAMETERS_H #include "AliFMDDigit.h" // ALIFMDDIGIT_H +#include "AliFMDSDigit.h" // ALIFMDSDIGIT_H #include "AliFMDRawStream.h" // ALIFMDRAWSTREAM_H #include "AliRawReader.h" // ALIRAWREADER_H #include "AliFMDRawReader.h" // ALIFMDRAWREADER_H @@ -77,12 +78,17 @@ AliFMDRawReader::AliFMDRawReader(AliRawReader* reader, TTree* tree) : TTask("FMDRawReader", "Reader of Raw ADC values from the FMD"), fTree(tree), fReader(reader), - fSampleRate(1), + // fSampleRate(1), fData(0), fNbytes(0), fSeen() { // Default CTOR + for (Int_t i = 0; i < 3; i++) { + fSampleRate[i] = 0; + fZeroSuppress[i] = kFALSE; + fNoiseFactor[i] = 1; + } } //____________________________________________________________________ @@ -104,13 +110,11 @@ AliFMDRawReader::Exec(Option_t*) array->GetEntriesFast(), nWrite)); } - //____________________________________________________________________ Bool_t -AliFMDRawReader::NextSignal(UShort_t& det, Char_t& rng, - UShort_t& sec, UShort_t& str, - Short_t& adc, Bool_t& zs, - UShort_t& fac) +AliFMDRawReader::NextSample(UShort_t& det, Char_t& rng, UShort_t& sec, + UShort_t& str, UShort_t& sam, UShort_t& rat, + Short_t& adc, Bool_t& zs, UShort_t& fac) { // Scan current event for next signal. It returns kFALSE when // there's no more data in the event. @@ -118,12 +122,14 @@ AliFMDRawReader::NextSignal(UShort_t& det, Char_t& rng, static AliFMDParameters* pars = 0; static AliFMDAltroMapping* map = 0; static Int_t ddl = -1; - static UInt_t rate = 0; + // static UInt_t rate = 0; static UShort_t tdet = 0; static Char_t trng = '\0'; static UShort_t tsec = 0; static Short_t tstr = 0; static Short_t bstr = -1; + static Short_t tsam = -1; + static UInt_t trate = 0; static Int_t hwaddr = -1; static UShort_t stripMin = 0; static UShort_t stripMax = 0; // 127; @@ -146,11 +152,12 @@ AliFMDRawReader::NextSignal(UShort_t& det, Char_t& rng, // Reset variables ddl = -1; - rate = 0; + trate = 0; tdet = 0; trng = '\0'; tsec = 0; tstr = 0; + tsam = -1; hwaddr = -1; } do { @@ -163,14 +170,14 @@ AliFMDRawReader::NextSignal(UShort_t& det, Char_t& rng, Int_t thisDDL = stream.GetDDLNumber(); AliFMDDebug(10, ("RCU @ DDL %d", thisDDL)); if (thisDDL != ddl) { - ddl = thisDDL; - zs = stream.GetZeroSupp(); - fac = stream.GetNPostsamples(); - tdet = map->DDL2Detector(ddl); - rate = 0; // stream.GetNPresamples(); + ddl = thisDDL; + fZeroSuppress[ddl] = zs = stream.GetZeroSupp(); + fNoiseFactor[ddl] = fac = stream.GetNPostsamples(); + fSampleRate[ddl] = trate = 0; // stream.GetNPresamples(); + tdet = map->DDL2Detector(ddl); AliFMDDebug(10, ("RCU @ DDL %d zero suppression: %s",ddl, zs?"yes":"no")); AliFMDDebug(10, ("RCU @ DDL %d noise factor: %d", ddl,fac)); - AliFMDDebug(10, ("RCU @ DDL %d sample rate: %d", ddl, rate)); + AliFMDDebug(10, ("RCU @ DDL %d sample rate: %d", ddl, trate)); } Int_t thisAddr = stream.GetHWAddress(); AliFMDDebug(10, ("RCU @ DDL %d, Address 0x%03x", ddl, thisAddr)); @@ -185,9 +192,10 @@ AliFMDRawReader::NextSignal(UShort_t& det, Char_t& rng, stripMin = pars->GetMinStrip(tdet, trng, tsec, bstr); stripMax = pars->GetMaxStrip(tdet, trng, tsec, bstr); preSamp = pars->GetPreSamples(tdet, trng, tsec, bstr); - if (rate == 0) rate = pars->GetSampleRate(tdet, trng, tsec, bstr); + if (trate == 0) + fSampleRate[ddl] = trate = pars->GetSampleRate(tdet, trng, tsec, bstr); AliFMDDebug(10, ("RCU @ DDL %d, Address 0x%03x sample rate: %d", - ddl, hwaddr, rate)); + ddl, hwaddr, trate)); Int_t nChAddrMismatch = stream.GetNChAddrMismatch(); Int_t nChLenMismatch = stream.GetNChLengthMismatch(); @@ -220,22 +228,12 @@ AliFMDRawReader::NextSignal(UShort_t& det, Char_t& rng, Short_t strOff = 0; UShort_t samp = 0; - map->Timebin2Strip(tsec, t, preSamp, rate, strOff, samp); + map->Timebin2Strip(tsec, t, preSamp, trate, strOff, samp); tstr = bstr + strOff; + tsam = samp; AliFMDDebug(20, ("0x%04x/0x%03x/%04d maps to FMD%d%c[%2d,%3d]-%d", ddl, hwaddr, t, tdet, trng, tsec, tstr, samp)); - Bool_t take = kFALSE; - switch (rate) { - case 1: take = kTRUE; break; - case 2: if (samp == 1) take = kTRUE; break; - case 3: if (samp == 1) take = kTRUE; break; - case 4: if (samp == 2) take = kTRUE; break; - default: if (samp == rate-2) take = kTRUE; break; - } - if (!take) continue; - - // Local strip number to channel Short_t l = (tstr > bstr ? tstr - bstr : bstr - tstr); AliFMDDebug(10, ("Checking if strip %d in range [%d,%d]", @@ -252,6 +250,8 @@ AliFMDRawReader::NextSignal(UShort_t& det, Char_t& rng, rng = trng; sec = tsec; str = tstr; + sam = tsam; + rat = trate; // adc = stream.GetSignal(); break; @@ -259,6 +259,49 @@ AliFMDRawReader::NextSignal(UShort_t& det, Char_t& rng, return kTRUE; } +//____________________________________________________________________ +Bool_t +AliFMDRawReader::NextSignal(UShort_t& det, Char_t& rng, + UShort_t& sec, UShort_t& str, + Short_t& adc, Bool_t& zs, + UShort_t& fac) +{ + + do { + UShort_t samp, rate; + if (!NextSample(det, rng, sec, str, samp, rate, adc, zs, fac)) + return kFALSE; + + Bool_t take = kFALSE; + switch (rate) { + case 1: take = kTRUE; break; + case 2: if (samp == 1) take = kTRUE; break; + case 3: if (samp == 1) take = kTRUE; break; + case 4: if (samp == 2) take = kTRUE; break; + default: if (samp == rate-2) take = kTRUE; break; + } + if (!take) continue; + break; + } while (true); + return kTRUE; +} + +//____________________________________________________________________ +Bool_t +AliFMDRawReader::SelectSample(UShort_t samp, UShort_t rate) +{ + Bool_t take = kFALSE; + switch (rate) { + case 1: take = kTRUE; break; + case 2: if (samp == 1) take = kTRUE; break; + case 3: if (samp == 1) take = kTRUE; break; + case 4: if (samp == 2) take = kTRUE; break; + default: if (samp == rate-2) take = kTRUE; break; + } + + return take; +} + #if 1 //____________________________________________________________________ Bool_t @@ -291,7 +334,7 @@ AliFMDRawReader::ReadAdcs(TClonesArray* array) Int_t oldddl = -1; UInt_t ddl = 0; - UInt_t rate = 0; + // UInt_t rate = 0; UInt_t last = 0; UInt_t hwaddr = 0; // Data array is approx twice the size needed. @@ -318,6 +361,12 @@ AliFMDRawReader::ReadAdcs(TClonesArray* array) fNoiseFactor[ddl] = input.GetNPostsamples(); AliFMDDebug(20, ("RCU @ DDL %d noise factor: %d", ddl,fNoiseFactor[ddl])); + // WARNING: We store the noise factor in the 2nd baseline + // filters excluded post samples, since we'll never use that + // mode. + fSampleRate[ddl] = input.GetNPretriggerSamples(); + AliFMDDebug(20, ("RCU @ DDL %d Sample rate: %d", ddl,fNoiseFactor[ddl])); + Int_t nChAddrMismatch = input.GetNChAddrMismatch(); Int_t nChLenMismatch = input.GetNChLengthMismatch(); if (nChAddrMismatch != 0) @@ -343,16 +392,14 @@ AliFMDRawReader::ReadAdcs(TClonesArray* array) "hardware address 0x%03x", ddl, hwaddr)); continue; } - AliFMDDebug(1, ("Board: 0x%02x, Altro: 0x%x, Channel: 0x%x, Length: %4d", + AliFMDDebug(5, ("Board: 0x%02x, Altro: 0x%x, Channel: 0x%x, Length: %4d", board, chip, channel, last)); stripMin = pars->GetMinStrip(det, ring, sec, strbase); stripMax = pars->GetMaxStrip(det, ring, sec, strbase); preSamp = pars->GetPreSamples(det, ring, sec, strbase); - // WARNING: We use the number of pre-samples to store the - // oversampling rate in. - rate = input.GetNPretriggerSamples(); - if (rate == 0) rate = pars->GetSampleRate(det, ring, sec, strbase); + if (fSampleRate[ddl] == 0) + fSampleRate[ddl] = pars->GetSampleRate(det, ring, sec, strbase); // Loop over the `timebins', and make the digits for (size_t i = 0; i < last; i++) { @@ -360,7 +407,7 @@ AliFMDRawReader::ReadAdcs(TClonesArray* array) AliFMDDebug(15, ("0x%04x/0x%03x/%04d %4d", ddl, hwaddr, i, data[i])); Short_t stroff = 0; - map->Timebin2Strip(sec, i, preSamp, rate, stroff, samp); + map->Timebin2Strip(sec, i, preSamp, fSampleRate[ddl], stroff, samp); Short_t str = strbase + stroff; AliFMDDebug(10, ("0x%04x/0x%03x/%04d maps to FMD%d%c[%2d,%3d]-%d", @@ -371,7 +418,7 @@ AliFMDRawReader::ReadAdcs(TClonesArray* array) continue; } - Short_t lstrip = (i - preSamp) / rate + stripMin; + Short_t lstrip = (i - preSamp) / fSampleRate[ddl] + stripMin; AliFMDDebug(15, ("Checking if strip %d (%d) in range [%d,%d]", lstrip, str, stripMin, stripMax)); @@ -381,6 +428,9 @@ AliFMDRawReader::ReadAdcs(TClonesArray* array) data[i] = 0; // Reset cache continue; } + // Possibly do pedestal subtraction of signal + Int_t counts = data[i]; + // Check the cache of indicies Int_t idx = fSeen(det, ring, sec, str); @@ -391,11 +441,11 @@ AliFMDRawReader::ReadAdcs(TClonesArray* array) det, ring, sec, str, samp, i)); new ((*array)[idx]) AliFMDDigit(det, ring, sec, str); } - AliFMDDigit* digit = static_cast(array->At(idx)); + AliFMDBaseDigit* digit = static_cast(array->At(idx)); AliFMDDebug(10, - ("Setting from FMD%d%c[%2d,%3d]-%d from timebin %4d = %4d", - det, ring, sec, str, samp, i, data[i])); - digit->SetCount(samp, data[i]); + ("Setting FMD%d%c[%2d,%3d]-%d from timebin %4d=%4d (%4d)", + det, ring, sec, str, samp, i, counts, data[i])); + digit->SetCount(samp, counts); data[i] = 0; // Reset cache } } diff --git a/FMD/AliFMDRawReader.h b/FMD/AliFMDRawReader.h index b1e2a537fdd..2499ff8326c 100644 --- a/FMD/AliFMDRawReader.h +++ b/FMD/AliFMDRawReader.h @@ -47,18 +47,34 @@ class AliFMDCalibStripRange; class AliFMDRawReader : public TTask { public: - /** CTOR - @param reader Raw reader - @param array Output tree */ + /** + * CTOR + * + * @param reader Raw reader + * @param array Output tree + */ AliFMDRawReader(AliRawReader* reader, TTree* array); - /** DTOR */ + /** + * DTOR + */ virtual ~AliFMDRawReader() {} - /** Read in, and store in output tree - @param option Not used */ + /** + * Read in, and store in output tree + * + * @param option Not used + */ virtual void Exec(Option_t* option=""); - /** Read ADC's into a TClonesArray of AliFMDDigit objects. - @param array Array to read into - @return @c true on success */ + /** + * Read ADC's into a TClonesArray of AliFMDDigit objects. + * + * @param array Array to read into + * @param summable Create SDigits rather than digits + * @param pedSub Whether to do pedestal subtraction. + * @param noiseFactor If we do pedestal subtraction, then this is + * the number we use to suppress remenants of the noise. + * + * @return @c true on success + */ virtual Bool_t ReadAdcs(TClonesArray* array); /** * Read SOD event into passed objects. @@ -92,16 +108,55 @@ public: */ UShort_t NoiseFactor(UShort_t ddl) const { return fNoiseFactor[ddl]; } + /** + * Get the next signal + * + * @param det On return, the detector + * @param rng On return, the ring + * @param sec On return, the sector + * @param str On return, the strip + * @param sam On return, the sample + * @param rat On return, the sample rate + * @param adc On return, the ADC value + * @param zs On return, whether zero-supp. is enabled + * @param fac On return, the usd noise factor + * + * @return true if valid data is returned + */ + Bool_t NextSample(UShort_t& det, Char_t& rng, UShort_t& sec, UShort_t& str, + UShort_t& sam, UShort_t& rat, Short_t& adc, + Bool_t& zs, UShort_t& fac); + /** + * Get the next signal + * + * @param det On return, the detector + * @param rng On return, the ring + * @param sec On return, the sector + * @param str On return, the strip + * @param adc On return, the ADC value + * @param zs On return, whether zero-supp. is enabled + * @param fac On return, the usd noise factor + * + * @return true if valid data is returned + */ Bool_t NextSignal(UShort_t& det, Char_t& rng, UShort_t& sec, UShort_t& str, Short_t& adc, Bool_t& zs, UShort_t& fac); + /** + * Whether to keep a sample based on the rate used. + * + * @param samp Sample number + * @param rate Over sampling rate + * + * @return Whether to keep the sample or not + */ + static Bool_t SelectSample(UShort_t samp, UShort_t rate); protected: AliFMDRawReader(const AliFMDRawReader& o) : TTask(o), fTree(0), fReader(0), - fSampleRate(0), fData(0), fNbytes(0), fSeen() @@ -112,7 +167,7 @@ protected: Int_t GetHalfringIndex(UShort_t det, Char_t ring, UShort_t board); TTree* fTree; //! Pointer to tree to read into AliRawReader* fReader; //! Pointer to raw reader - UShort_t fSampleRate; // The sample rate (if 0, inferred from data) + UShort_t fSampleRate[3]; // The sample rate (if 0, inferred from data) UChar_t* fData; ULong_t fNbytes; Bool_t fZeroSuppress[3]; diff --git a/FMD/AliFMDReconstructor.cxx b/FMD/AliFMDReconstructor.cxx index 0bd09864de5..5e9cab6d9ea 100644 --- a/FMD/AliFMDReconstructor.cxx +++ b/FMD/AliFMDReconstructor.cxx @@ -39,6 +39,7 @@ #include "AliFMDParameters.h" // ALIFMDPARAMETERS_H #include "AliFMDAltroMapping.h" // ALIFMDALTROMAPPING_H #include "AliFMDDigit.h" // ALIFMDDIGIT_H +#include "AliFMDSDigit.h" // ALIFMDDIGIT_H #include "AliFMDReconstructor.h" // ALIFMDRECONSTRUCTOR_H #include "AliFMDRawReader.h" // ALIFMDRAWREADER_H #include "AliFMDRecPoint.h" // ALIFMDMULTNAIIVE_H @@ -205,7 +206,7 @@ AliFMDReconstructor::ConvertDigits(AliRawReader* reader, TTree* digitsTree) const { // Convert Raw digits to AliFMDDigit's in a tree - AliFMDDebug(2, ("Reading raw data into digits tree")); + AliFMDDebug(1, ("Reading raw data into digits tree")); AliFMDRawReader rawRead(reader, digitsTree); // rawRead.SetSampleRate(fFMD->GetSampleRate()); rawRead.Exec(); @@ -276,6 +277,34 @@ AliFMDReconstructor::Reconstruct(AliRawReader* reader, TTree*) const } } +//____________________________________________________________________ +void +AliFMDReconstructor::Digitize(AliRawReader* reader, TClonesArray* sdigits) const +{ + // Reconstruct directly from raw data (no intermediate output on + // digit tree or rec point tree). + // + // Parameters: + // reader Raw event reader + // ctree Not used. + AliFMDRawReader rawReader(reader, 0); + + UShort_t det, sec, str, sam, rat, fac; + Short_t adc, oldDet = -1; + Bool_t zs; + Char_t rng; + + while (rawReader.NextSample(det, rng, sec, str, sam, rat, adc, zs, fac)) { + if (!rawReader.SelectSample(sam, rat)) continue; + if (det != oldDet) { + fZS[det-1] = zs; + fZSFactor[det-1] = fac; + oldDet = det; + } + DigitizeSignal(sdigits, det, rng, sec, str, sam, adc); + } +} + //____________________________________________________________________ void AliFMDReconstructor::Reconstruct(TTree* digitsTree, @@ -311,7 +340,7 @@ AliFMDReconstructor::Reconstruct(TTree* digitsTree, AliFMDDebug(5, ("Getting entry 0 from digit branch")); digitBranch->GetEntry(0); - AliFMDDebug(5, ("Processing digits")); + AliFMDDebug(1, ("Processing digits")); ProcessDigits(digits); Int_t written = clusterTree->Fill(); @@ -418,6 +447,99 @@ AliFMDReconstructor::ProcessSignal(UShort_t det, } +//____________________________________________________________________ +void +AliFMDReconstructor::DigitizeSignal(TClonesArray* sdigits, + UShort_t det, + Char_t rng, + UShort_t sec, + UShort_t str, + UShort_t /* sam */, + Short_t adc) const +{ + // Process the signal from a single strip + // + // Parameters: + // det Detector ID + // rng Ring ID + // sec Sector ID + // rng Strip ID + // adc ADC counts + // + AliFMDParameters* param = AliFMDParameters::Instance(); + // Check that the strip is not marked as dead + if (param->IsDead(det, rng, sec, str)) { + AliFMDDebug(10, ("FMD%d%c[%2d,%3d] is dead", det, rng, sec, str)); + return; + } + + // Substract pedestal. + UShort_t counts = SubtractPedestal(det, rng, sec, str, adc); + if(counts == USHRT_MAX || counts == 0) return; + + // Gain match digits. + Double_t edep = Adc2Energy(det, rng, sec, str, counts); + // Get rid of nonsense energy + if(edep < 0) return; + + Int_t n = sdigits->GetEntriesFast(); + // AliFMDSDigit* sdigit = + new ((*sdigits)[n]) + AliFMDSDigit(det, rng, sec, str, edep, counts, counts, counts, counts); + // sdigit->SetCount(sam, counts); +} + +//____________________________________________________________________ +UShort_t +AliFMDReconstructor::SubtractPedestal(UShort_t det, + Char_t rng, + UShort_t sec, + UShort_t str, + UShort_t adc, + Float_t noiseFactor, + Bool_t zsEnabled, + UShort_t zsNoiseFactor) const +{ + AliFMDParameters* param = AliFMDParameters::Instance(); + Float_t ped = (zsEnabled ? 0 : + param->GetPedestal(det, rng, sec, str)); + Float_t noise = param->GetPedestalWidth(det, rng, sec, str); + if(ped < 0 || noise < 0) { + AliWarningClass(Form("Invalid pedestal (%f) or noise (%f) " + "for FMD%d%c[%02d,%03d]", + ped, noise, det, rng, sec, str)); + return USHRT_MAX; + } + AliDebugClass(15, Form("Subtracting pedestal for FMD%d%c[%2d,%3d]=%4d " + "(%s w/factor %d, noise factor %f, " + "pedestal %8.2f+/-%8.2f)", + det, rng, sec, str, adc, + (zsEnabled ? "zs'ed" : "straight"), + zsNoiseFactor, noiseFactor, ped, noise)); + + Int_t counts = adc + Int_t(zsEnabled ? zsNoiseFactor * noise : - ped); + counts = TMath::Max(Int_t(counts), 0); + // Calculate the noise factor for suppressing remenants of the noise + // peak. If we have done on-line zero suppression, we only check + // for noise signals that are larger than the suppressed noise. If + // the noise factor used on line is larger than the factor used + // here, we do not do this check at all. + // + // For example: + // Online factor | Read factor | Result + // ---------------+--------------+------------------------------- + // 2 | 3 | Check if signal > 1 * noise + // 3 | 3 | Check if signal > 0 + // 3 | 2 | Check if signal > 0 + // + // In this way, we make sure that we do not suppress away too much + // data, and that the read-factor is the most stringent cut. + Float_t nf = TMath::Max(0.F, noiseFactor - (zsEnabled ? zsNoiseFactor : 0)); + if (counts < noise * nf) counts = 0; + if (counts > 0) AliDebugClass(15, "Got a hit strip"); + + return counts; +} //____________________________________________________________________ @@ -439,29 +561,12 @@ AliFMDReconstructor::SubtractPedestal(UShort_t det, // Return: // Pedestal subtracted signal or USHRT_MAX in case of problems // - AliFMDParameters* param = AliFMDParameters::Instance(); - Bool_t zs = fZS[det-1]; - UShort_t fac = fZSFactor[det-1]; - Float_t ped = (zs ? 0 : - param->GetPedestal(det, rng, sec, str)); - Float_t noise = param->GetPedestalWidth(det, rng, sec, str); - if(ped < 0 || noise < 0) { - AliWarning(Form("Invalid pedestal (%f) or noise (%f) " - "for FMD%d%c[%02d,%03d]", ped, noise, det, rng, sec, str)); - return USHRT_MAX; - } - - AliFMDDebug(5, ("Subtracting pedestal %f from signal %d", ped, adc)); - // if (digit->Count3() > 0) adc = digit->Count3(); - // else if (digit->Count2() > 0) adc = digit->Count2(); - // else adc = digit->Count1(); - Int_t counts = adc + Int_t(zs ? fac * noise : - ped); - counts = TMath::Max(Int_t(counts), 0); - if (counts < noise * fNoiseFactor) counts = 0; - if (counts > 0) AliFMDDebug(15, ("Got a hit strip")); + UShort_t counts = SubtractPedestal(det, rng, sec, str, adc, + fNoiseFactor, fZS[det-1], + fZSFactor[det-1]); if (fDiagStep1) fDiagStep1->Fill(adc, counts); - return UShort_t(counts); + return counts; } //____________________________________________________________________ @@ -470,7 +575,6 @@ AliFMDReconstructor::Adc2Energy(UShort_t det, Char_t rng, UShort_t sec, UShort_t str, - Float_t eta, UShort_t count) const { // Converts number of ADC counts to energy deposited. @@ -504,7 +608,6 @@ AliFMDReconstructor::Adc2Energy(UShort_t det, // rng Ring ID // sec Sector ID // rng Strip ID - // eta Psuedo-rapidity // counts Number of ADC counts over pedestal // Return // The energy deposited in a single strip, or -1 in case of problems @@ -523,6 +626,56 @@ AliFMDReconstructor::Adc2Energy(UShort_t det, Double_t edep = ((count * param->GetEdepMip()) / (gain * param->GetDACPerMIP())); + return edep; +} + +//____________________________________________________________________ +Float_t +AliFMDReconstructor::Adc2Energy(UShort_t det, + Char_t rng, + UShort_t sec, + UShort_t str, + Float_t eta, + UShort_t count) const +{ + // Converts number of ADC counts to energy deposited. + // Note, that this member function can be overloaded by derived + // classes to do strip-specific look-ups in databases or the like, + // to find the proper gain for a strip. + // + // In the first simple version, we calculate the energy deposited as + // + // EnergyDeposited = cos(theta) * gain * count + // + // where + // + // Pre_amp_MIP_Range + // gain = ----------------- * Energy_deposited_per_MIP + // ADC_channel_size + // + // is constant and the same for all strips. + // + // For the production we use the conversion measured in the NBI lab. + // The total conversion is then: + // + // gain = ADC / DAC + // + // EdepMip * count + // => energy = ---------------- + // gain * DACPerADC + // + // Parameters: + // det Detector ID + // rng Ring ID + // sec Sector ID + // rng Strip ID + // eta Psuedo-rapidity + // counts Number of ADC counts over pedestal + // Return + // The energy deposited in a single strip, or -1 in case of problems + // + Double_t edep = Adc2Energy(det, rng, sec, str, count); + if (fDiagStep2) fDiagStep2->Fill(count, edep); if (fAngleCorrect) { Double_t theta = 2 * TMath::ATan(TMath::Exp(-eta)); diff --git a/FMD/AliFMDReconstructor.h b/FMD/AliFMDReconstructor.h index 446b32e3c52..a04989d0c16 100644 --- a/FMD/AliFMDReconstructor.h +++ b/FMD/AliFMDReconstructor.h @@ -117,6 +117,16 @@ public: */ virtual void FillESD(AliRawReader*, TTree* clusterTree, AliESDEvent* esd) const; + + /** + * Create SDigits from raw data + * + * @param reader The raw reader + * @param sdigits Array to fill with AliFMDSDigit objects. + */ + virtual void Digitize(AliRawReader* reader, + TClonesArray* sdigits) const; + /** * Not used */ @@ -199,6 +209,49 @@ protected: UShort_t sec, UShort_t str, Short_t adc) const; + /** + * Process the signal from a single strip. + * + * @param sdigits Array to fill + * @param det Detector number + * @param rng Ring identifier + * @param sec Sector number + * @param str Strip number + * @param sam Sample number + * @param adc Number of ADC counts for this strip + */ + virtual void DigitizeSignal(TClonesArray* sdigits, + UShort_t det, + Char_t rng, + UShort_t sec, + UShort_t str, + UShort_t sam, + Short_t adc) const; + /** + * Subtract the pedestal off the ADC counts. + * + * @param det Detector number + * @param rng Ring identifier + * @param sec Sector number + * @param str Strip number + * @param adc ADC counts + * @param noiseFactor If pedestal substracted pedestal is less then + * this times the noise, then consider this to be 0. + * @param zsEnabled Whether zero-suppression is on. + * @param zsNoiseFactor Noise factor used in on-line pedestal + * subtraction. + * + * @return The pedestal subtracted ADC counts (possibly 0), or @c + * USHRT_MAX in case of problems. + */ + virtual UShort_t SubtractPedestal(UShort_t det, + Char_t rng, + UShort_t sec, + UShort_t str, + UShort_t adc, + Float_t noiseFactor, + Bool_t zsEnabled, + UShort_t zsNoiseFactor) const; /** * Substract pedestals from raw ADC in @a digit * @@ -233,6 +286,29 @@ protected: * * @return Energy deposited @f$ E_i@f$ */ + virtual Float_t Adc2Energy(UShort_t det, + Char_t rng, + UShort_t sec, + UShort_t str, + UShort_t count) const; + /** + * Converts number of ADC counts to energy deposited. This is + * done by + * @f[ + * E_i = A_i g_i + * @f] + * where @f$ A_i@f$ is the pedestal subtracted ADC counts, and @f$ + * g_i@f$ is the gain for the @f$ i^{\mbox{th}}@f$ strip. + * + * @param det Detector number + * @param rng Ring identifier + * @param sec Sector number + * @param str Strip number + * @param eta Psuedo-rapidity of digit. + * @param count Pedestal subtracted ADC counts + * + * @return Energy deposited @f$ E_i@f$ + */ virtual Float_t Adc2Energy(UShort_t det, Char_t rng, UShort_t sec, diff --git a/FMD/AliFMDSDigit.cxx b/FMD/AliFMDSDigit.cxx index 9077fc5cec7..e4aeabdf2e5 100644 --- a/FMD/AliFMDSDigit.cxx +++ b/FMD/AliFMDSDigit.cxx @@ -79,7 +79,9 @@ AliFMDSDigit::AliFMDSDigit() fCount2(-1), fCount3(-1), fCount4(-1), - fLabels(0) + fNParticles(0), + fNPrimaries(0) + // , fLabels(0) { // cTOR } @@ -104,8 +106,8 @@ AliFMDSDigit::AliFMDSDigit(UShort_t detector, fCount3(count3), fCount4(count4), fNParticles(npart), - fNPrimaries(nprim), - fLabels(refs) + fNPrimaries(nprim) + // , fLabels(refs) { // // Creates a real data digit object @@ -120,6 +122,7 @@ AliFMDSDigit::AliFMDSDigit(UShort_t detector, // count1 ADC count (a 10-bit word) // count2 ADC count (a 10-bit word) -1 if not used // count3 ADC count (a 10-bit word) -1 if not used + for (Int_t i = 0; i < refs.fN; i++) AddTrack(refs.fArray[i]); } //____________________________________________________________________ @@ -128,17 +131,28 @@ AliFMDSDigit::Print(Option_t* option) const { // Print digit to standard out AliFMDBaseDigit::Print(); - std::cout << "\t" << fEdep << " -> " - << fCount1 << " (" << fCount2 << "," << fCount3 << "," - << fCount4 << ") = " << Counts() << std::flush; + std::cout << "\t" + << std::setw(10) << fEdep << " -> " + << std::setw(4) << fCount1 << " (" + << std::setw(4) << fCount2 << "," + << std::setw(4) << fCount3 << "," + << std::setw(4) << fCount4 << ") = " + << std::setw(4) << Counts() << std::flush; + TString opt(option); if (opt.Contains("p", TString::kIgnoreCase)) - std::cout << " [" << fNPrimaries << "/" << fNParticles << "]" + std::cout << " [" + << std::setw(2) << fNPrimaries << "/" + << std::setw(2) << fNParticles << "]" << std::flush; if (opt.Contains("l", TString::kIgnoreCase)) { std::cout << " "; + for (Int_t i = 0; i < GetNTrack(); i++) + std::cout << (i == 0 ? "" : ",") << std::setw(5) << fTracks[i]; +#if 0 for (Int_t i = 0; i < fLabels.fN; i++) std::cout << (i == 0 ? "" : ",") << fLabels.fArray[i]; +#endif } std::cout << std::endl; } diff --git a/FMD/AliFMDSDigit.h b/FMD/AliFMDSDigit.h index 18598eeeddf..c786cc26401 100644 --- a/FMD/AliFMDSDigit.h +++ b/FMD/AliFMDSDigit.h @@ -27,17 +27,22 @@ class AliFMDSDigit : public AliFMDBaseDigit { public: - /** CTOR */ + /** + * CTOR + */ AliFMDSDigit(); - /** Constrctor - @param detector Detector - @param ring Ring - @param sector Sector - @param strip Strip - @param edep Energy deposited - @param count ADC (first sample) - @param count2 ADC (second sample, or -1 if not used) - @param count3 ADC (third sample, or -1 if not used) */ + /** + * Constrctor + * + * @param detector Detector + * @param ring Ring + * @param sector Sector + * @param strip Strip + * @param edep Energy deposited + * @param count ADC (first sample) + * @param count2 ADC (second sample, or -1 if not used) + * @param count3 ADC (third sample, or -1 if not used) + */ AliFMDSDigit(UShort_t detector, Char_t ring='\0', UShort_t sector=0, @@ -50,29 +55,67 @@ public: UShort_t npart=0, UShort_t nprim=0, const TArrayI& refs=TArrayI()); - /** DTOR */ + /** + * DTOR + */ virtual ~AliFMDSDigit() {} - /** @return ADC count (first sample) */ - UShort_t Count1() const { return fCount1; } - /** @return ADC count (second sample, or -1 if not used) */ - Short_t Count2() const { return fCount2; } - /** @return ADC count (third sample, or -1 if not used) */ - Short_t Count3() const { return fCount3; } - /** @return ADC count (third sample, or -1 if not used) */ - Short_t Count4() const { return fCount4; } - /** @return Canonical ADC counts */ - UShort_t Counts() const; - /** @return Energy deposited */ - Float_t Edep() const { return fEdep; } - /** @return Number of particles that hit this strip */ + /** + * + * @return ADC count (first sample) + */ + UShort_t Count1() const { return fCount1; } + /** + * + * @return ADC count (second sample, or -1 if not used) + */ + Short_t Count2() const { return fCount2; } + /** + * + * @return ADC count (third sample, or -1 if not used) + */ + Short_t Count3() const { return fCount3; } + /** + * + * @return ADC count (third sample, or -1 if not used) + */ + Short_t Count4() const { return fCount4; } + /** + * + * @return Canonical ADC counts + */ + UShort_t Counts() const; + /** + * + * @return Energy deposited + */ + Float_t Edep() const { return fEdep; } + /** + * + * @return Number of particles that hit this strip + */ UShort_t NParticles() const { return fNParticles; } - /** @return Number of primary particles that hit this strip */ + /** + * + * @return Number of primary particles that hit this strip + */ UShort_t NPrimaries() const { return fNPrimaries; } +#if 0 /** @return the track labels */ const TArrayI& TrackLabels() const { return fLabels; } - /** Print info - @param opt Not used */ +#endif + /** + * Print info + * + * @param opt Not used + */ void Print(Option_t* opt="") const; + /** + * Set the count value + * + * @param s Sample number + * @param c Counts + */ + void SetCount(UShort_t s, Short_t c); protected: Float_t fEdep; // Energy deposited UShort_t fCount1; // Digital signal @@ -81,8 +124,10 @@ protected: Short_t fCount4; // Digital signal (-1 if not used) UShort_t fNParticles; // Total number of particles that hit this strip UShort_t fNPrimaries; // Number of primary particles that his this strip +#if 0 TArrayI fLabels; // MC-truth track labels - ClassDef(AliFMDSDigit,4) // Summable FMD digit +#endif + ClassDef(AliFMDSDigit,5) // Summable FMD digit }; inline UShort_t @@ -94,6 +139,16 @@ AliFMDSDigit::Counts() const return fCount1; } +inline void +AliFMDSDigit::SetCount(UShort_t i, Short_t c) +{ + switch (i) { + case 0: fCount1 = c; break; + case 1: fCount2 = c; break; + case 2: fCount3 = c; break; + case 3: fCount4 = c; break; + } +} #endif //____________________________________________________________________ diff --git a/FMD/AliFMDStripIndex.h b/FMD/AliFMDStripIndex.h index 3ca4fbca759..b6d524cff99 100644 --- a/FMD/AliFMDStripIndex.h +++ b/FMD/AliFMDStripIndex.h @@ -15,16 +15,39 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -//Struct to encode a strip address into one integer -//developed by Christian Holm Christensen (cholm@nbi.dk). -// The functions are static -// to ensure applicability from anywhere. This is needed to smoothly store -//strip addresses in track references. +// Struct to encode a strip address into one integer +// developed by Christian Holm Christensen (cholm@nbi.dk). +// +// The functions are static to ensure applicability from +// anywhere. This is needed to smoothly store strip addresses in track +// references. +// // Added by Hans H. Dalsgaard (hans.dalsgaard@cern.ch) -struct AliFMDStripIndex +class AliFMDStripIndex { +public: + /** + * Constructor + * + */ + AliFMDStripIndex() {} + /** + * Destructor + * + */ + virtual ~AliFMDStripIndex() {} + /** + * Pack an identifier from detector coordinates + * + * @param det Detector + * @param rng Ring + * @param sec Sector + * @param str Strip + * + * @return Packed identifier + */ static UInt_t Pack(UShort_t det, Char_t rng, UShort_t sec, UShort_t str) { UInt_t irg = (rng == 'I' || rng == 'i' ? 0 : 1); @@ -34,8 +57,17 @@ struct AliFMDStripIndex ((det & 0x003) << 17)); return id; } + /** + * Unpack an identifier to detector coordinates + * + * @param id Identifier to unpack + * @param det On return, the detector + * @param rng On return, the ring + * @param sec On return, the sector + * @param str On return, the strip + */ static void Unpack(UInt_t id, - UShort_t& det, Char_t& rng, UShort_t& sec, UShort_t& str) + UShort_t& det, Char_t& rng, UShort_t& sec, UShort_t& str) { str = ((id >> 0) & 0x1FF); sec = ((id >> 9) & 0x03F); @@ -43,6 +75,11 @@ struct AliFMDStripIndex det = ((id >> 17) & 0x003); } - ClassDef(AliFMDStripIndex,0) + ClassDef(AliFMDStripIndex,1) }; #endif +// +// Local Variables: +// mode: C++ +// End: +// diff --git a/FMD/AliFMDSurveyToAlignObjs.cxx b/FMD/AliFMDSurveyToAlignObjs.cxx index bac9fd7dc7f..3c43d8a2653 100644 --- a/FMD/AliFMDSurveyToAlignObjs.cxx +++ b/FMD/AliFMDSurveyToAlignObjs.cxx @@ -12,43 +12,157 @@ #include #include "AliFMDGeometry.h" +//____________________________________________________________________ +Double_t +AliFMDSurveyToAlignObjs::GetUnitFactor() const +{ + // Returns the conversion factor from the measured values to + // centimeters. + TString units(fSurveyObj->GetUnits()); + if (units.CompareTo("mm", TString::kIgnoreCase) == 0) return .1; + else if (units.CompareTo("cm", TString::kIgnoreCase) == 0) return 1.; + else if (units.CompareTo("m", TString::kIgnoreCase) == 0) return 100.; + return 1; +} + //____________________________________________________________________ Bool_t -AliFMDSurveyToAlignObjs::GetFMD2Plane(Double_t* rot, - Double_t* trans) +AliFMDSurveyToAlignObjs::GetPoint(const char* name, + TVector3& point, + TVector3& error) const { + // Get named point. On return, point will contain the point + // coordinates in centimeters, and error will contain the + // meassurement errors in centimeters too. If no point is found, + // returns false, otherwise true. + Double_t unit = GetUnitFactor(); + TObject* obj = fSurveyPoints->FindObject(name); + if (!obj) return kFALSE; + + AliSurveyPoint* p = static_cast(obj); + point.SetXYZ(unit * p->GetX(), + unit * p->GetY(), + unit * p->GetZ()); + error.SetXYZ(unit * p->GetPrecisionX(), + unit * p->GetPrecisionY(), + unit * p->GetPrecisionZ()); + return kTRUE; +} - // The possile survey points - const char* names[] = { "FMD2_ITOP", "FMD2_OTOP", - "FMD2_IBOTM", "FMD2_OBOTM", - "FMD2_IBOT", "FMD2_OBOT", - 0 }; - const char** name = names; - Int_t i = 0; - TGraph2DErrors g; +//____________________________________________________________________ +Bool_t +AliFMDSurveyToAlignObjs::CalculatePlane(const TVector3& a, + const TVector3& b, + const TVector3& c, + Double_t* trans, + Double_t* rot) const +{ + // Vector a->b, b->c, and normal to plane defined by these two + // vectors. + TVector3 ab(b-a), bc(c-b); - Double_t unit = 1.; - TString units(fSurveyObj->GetUnits()); - if (units.CompareTo("mm", TString::kIgnoreCase) == 0) unit = .1; - else if (units.CompareTo("cm", TString::kIgnoreCase) == 0) unit = 1.; - else if (units.CompareTo("m", TString::kIgnoreCase) == 0) unit = 100.; + // Normal vector to the plane of the fiducial marks obtained + // as cross product of the two vectors on the plane d0^d1 + TVector3 nn(ab.Cross(bc)); + if (nn.Mag() < 1e-8) { + Info("DoIt", "Normal vector is null vector"); + return kFALSE; + } + // We express the plane in Hessian normal form. + // + // n x = -p, + // + // where n is the normalised normal vector given by + // + // n_x = a / l, n_y = b / l, n_z = c / l, p = d / l + // + // with l = sqrt(a^2+b^2+c^2) and a, b, c, and d are from the + // normal plane equation + // + // ax + by + cz + d = 0 + // + // Normalize + TVector3 n(nn.Unit()); + Double_t p = - (n * a); + + // The center of the square with the fiducial marks as the + // corners. The mid-point of one diagonal - md. Used to get the + // center of the surveyd box. + TVector3 md(a + c); + md *= 1/2.; + + // The center of the box. + TVector3 orig(md); + trans[0] = orig[0]; + trans[1] = orig[1]; + trans[2] = orig[2]; + + // Normalize the spanning vectors + TVector3 uab(ab.Unit()); + TVector3 ubc(bc.Unit()); + + for (size_t i = 0; i < 3; i++) { + rot[i * 3 + 0] = uab[i]; + rot[i * 3 + 1] = ubc[i]; + rot[i * 3 + 2] = n[i]; + } + return kTRUE; +} + + +//____________________________________________________________________ +Bool_t +AliFMDSurveyToAlignObjs::GetFMD1Plane(Double_t* rot, Double_t* trans) const +{ + + // The possile survey points + TVector3 icb, ict, ocb, oct, dummy; + Int_t missing = 0; + if (!GetPoint("V0L_ICB", icb, dummy)) missing++; + if (!GetPoint("V0L_ICT", icb, dummy)) missing++; + if (!GetPoint("V0L_OCB", icb, dummy)) missing++; + if (!GetPoint("V0L_OCT", icb, dummy)) missing++; + + // Check that we have enough points + if (missing > 1) { + Error("GetFMD1Plane", "Only got %d survey points - no good for FMD1 plane", + 4-missing); + return kFALSE; + } + + if (!CalculatePlane(icb, ict, oct, rot, trans)) return kFALSE; + + return kTRUE; +} + +//____________________________________________________________________ +Bool_t +AliFMDSurveyToAlignObjs::DoFMD1() +{ + return kTRUE; +} + +//____________________________________________________________________ +Bool_t +AliFMDSurveyToAlignObjs::GetFMD2Plane(Double_t* rot, Double_t* trans) const +{ + + // The possile survey points + const char* names[] = { "FMD2_ITOP", "FMD2_OTOP", + "FMD2_IBOTM", "FMD2_OBOTM", + "FMD2_IBOT", "FMD2_OBOT", + 0 }; + const char** name = names; + Int_t i = 0; + TGraph2DErrors g; // Loop and fill graph while (*name) { - TObject* obj = fSurveyPoints->FindObject(*name); - name++; - if (!obj) continue; - - AliSurveyPoint* p = static_cast(obj); - Double_t x = unit * p->GetX(); - Double_t y = unit * p->GetY(); - Double_t z = unit * p->GetZ(); - Double_t ex = unit * p->GetPrecisionX(); - Double_t ey = unit * p->GetPrecisionY(); - Double_t ez = unit * p->GetPrecisionZ(); - - g.SetPoint(i, x, y, z); - g.SetPointError(i, ex, ey, ez); + TVector3 p, e; + if (!GetPoint(*name++, p, e)) continue; + + g.SetPoint(i, p.X(), p.Y(), p.Z()); + g.SetPointError(i, e.X(), e.Y(), e.Z()); i++; } @@ -83,9 +197,9 @@ AliFMDSurveyToAlignObjs::GetFMD2Plane(Double_t* rot, TVector3 ua(a.Unit()); TVector3 ub(b.Unit()); Double_t angAb = ua.Angle(ub); - PrintVector("ua: ", ua); - PrintVector("ub: ", ub); - std::cout << "Angle: " << angAb * 180 / TMath::Pi() << std::endl; + // PrintVector("ua: ", ua); + // PrintVector("ub: ", ub); + // std::cout << "Angle: " << angAb * 180 / TMath::Pi() << std::endl; for (size_t i = 0; i < 3; i++) { rot[i * 3 + 0] = ua[i]; diff --git a/FMD/AliFMDSurveyToAlignObjs.h b/FMD/AliFMDSurveyToAlignObjs.h index e510dbdfb28..02b793cf47a 100644 --- a/FMD/AliFMDSurveyToAlignObjs.h +++ b/FMD/AliFMDSurveyToAlignObjs.h @@ -13,9 +13,18 @@ public: void Run(); Bool_t CreateAlignObjs() { return kTRUE; } protected: + Bool_t DoFMD1(); + Bool_t GetFMD1Plane(Double_t* rot, Double_t* trans) const; Bool_t DoFMD2(); - Bool_t GetFMD2Plane(Double_t* rot, Double_t* trans); - + Bool_t GetFMD2Plane(Double_t* rot, Double_t* trans) const; + + Double_t GetUnitFactor() const; + Bool_t GetPoint(const char* name, TVector3& p, TVector3& e) const; + Bool_t CalculatePlane(const TVector3& a, + const TVector3& b, + const TVector3& c, + Double_t* trans, + Double_t* rot) const; static void PrintVector(const char* text, const Double_t* v); static void PrintVector(const char* text, const TVector3& v); static void PrintRotation(const char* text, const Double_t* rot); diff --git a/FMD/Config.C b/FMD/Config.C index a5bc57bb915..8e4105e328e 100644 --- a/FMD/Config.C +++ b/FMD/Config.C @@ -170,7 +170,7 @@ enum MC_t { // Functions Float_t EtaToTheta(Float_t eta); EG_t LookupEG(const Char_t* name); -AliGenerator* GeneratorFactory(EG_t eg, Rad_t rad, TString& comment); +AliGenerator* GeneratorFactory(EG_t eg, Rad_t rad); AliGenHijing* HijingStandard(); void ProcessEnvironmentVars(EG_t& eg, Int_t& seed); @@ -192,10 +192,6 @@ Config() Int_t seed = 12345; //Set 0 to use the current time MC_t mc = kGEANT3TGEO; - //____________________________________________________________________ - // Comment line - static TString comment; - //____________________________________________________________________ // Get settings from environment variables ProcessEnvironmentVars(eg, seed); @@ -367,7 +363,7 @@ Config() //__________________________________________________________________ // Generator Configuration - AliGenerator* gener = GeneratorFactory(eg, rad, comment); + AliGenerator* gener = GeneratorFactory(eg, rad); gener->SetOrigin(0, 0, 0); // vertex position gener->SetSigma(0, 0, 5.3); // Sigma in (X,Y,Z) (cm) on IP position gener->SetCutVertexZ(1.); // Truncate at 1 sigma @@ -375,30 +371,6 @@ Config() gener->SetTrackingFlag(1); gener->Init(); - //__________________________________________________________________ - // - // Comments - // - switch (mag) { - case k2kG: comment = comment.Append(" | L3 field 0.2 T"); break; - case k4kG: comment = comment.Append(" | L3 field 0.4 T"); break; - case k5kG: comment = comment.Append(" | L3 field 0.5 T"); break; - } - - switch (rad) { - case kGluonRadiation: - comment = comment.Append(" | Gluon Radiation On"); break; - default: - comment = comment.Append(" | Gluon Radiation Off"); break; - } - - switch(geo) { - case kHoles: comment = comment.Append(" | Holes for PHOS/HMPID"); break; - default: comment = comment.Append(" | No holes for PHOS/HMPID"); break; - } - - std::cout << "\n\n Comment: " << comment << "\n" << std::endl; - //__________________________________________________________________ // Field (L3 0.4 T) AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 2, 1., 10., mag); @@ -488,7 +460,7 @@ LookupEG(const Char_t* name) //____________________________________________________________________ AliGenerator* -GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) +GeneratorFactory(EG_t eg, Rad_t rad) { Int_t isw = 3; if (rad == kNoGluonRadiation) isw = 0; @@ -520,7 +492,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) switch (eg) { case test50: { - comment = comment.Append(":HIJINGparam test 50 particles"); AliGenHIJINGpara *gener = new AliGenHIJINGpara(50); gener->SetMomentumRange(0, 999999.); gener->SetPhiRange(0., 360.); @@ -533,7 +504,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kParam_8000: { - comment = comment.Append(":HIJINGparam N=8000"); AliGenHIJINGpara *gener = new AliGenHIJINGpara(86030); gener->SetMomentumRange(0, 999999.); gener->SetPhiRange(0., 360.); @@ -546,7 +516,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kParam_4000: { - comment = comment.Append("HIJINGparam N=4000"); AliGenHIJINGpara *gener = new AliGenHIJINGpara(43015); gener->SetMomentumRange(0, 999999.); gener->SetPhiRange(0., 360.); @@ -559,7 +528,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kParam_2000: { - comment = comment.Append("HIJINGparam N=2000"); AliGenHIJINGpara *gener = new AliGenHIJINGpara(21507); gener->SetMomentumRange(0, 999999.); gener->SetPhiRange(0., 360.); @@ -572,7 +540,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kParam_fmd: { - comment = comment.Append("HIJINGparam N=100"); AliGenHIJINGpara *gener = new AliGenHIJINGpara(500); gener->SetMomentumRange(0, 999999.); gener->SetPhiRange(0., 360.); @@ -588,7 +555,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) // case kHijing_cent1: { - comment = comment.Append("HIJING cent1"); AliGenHijing *gener = HijingStandard(); // impact parameter range gener->SetImpactParameterRange(0., 5.); @@ -597,7 +563,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kHijing_cent2: { - comment = comment.Append("HIJING cent2"); AliGenHijing *gener = HijingStandard(); // impact parameter range gener->SetImpactParameterRange(0., 2.); @@ -609,7 +574,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) // case kHijing_per1: { - comment = comment.Append("HIJING per1"); AliGenHijing *gener = HijingStandard(); // impact parameter range gener->SetImpactParameterRange(5., 8.6); @@ -618,7 +582,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kHijing_per2: { - comment = comment.Append("HIJING per2"); AliGenHijing *gener = HijingStandard(); // impact parameter range gener->SetImpactParameterRange(8.6, 11.2); @@ -627,7 +590,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kHijing_per3: { - comment = comment.Append("HIJING per3"); AliGenHijing *gener = HijingStandard(); // impact parameter range gener->SetImpactParameterRange(11.2, 13.2); @@ -636,7 +598,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kHijing_per4: { - comment = comment.Append("HIJING per4"); AliGenHijing *gener = HijingStandard(); // impact parameter range gener->SetImpactParameterRange(13.2, 15.); @@ -645,7 +606,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kHijing_per5: { - comment = comment.Append("HIJING per5"); AliGenHijing *gener = HijingStandard(); // impact parameter range gener->SetImpactParameterRange(15., 100.); @@ -657,7 +617,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) // case kHijing_jj25: { - comment = comment.Append("HIJING Jet 25 GeV"); AliGenHijing *gener = HijingStandard(); // impact parameter range gener->SetImpactParameterRange(0., 5.); @@ -674,7 +633,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) case kHijing_jj50: { - comment = comment.Append("HIJING Jet 50 GeV"); AliGenHijing *gener = HijingStandard(); // impact parameter range gener->SetImpactParameterRange(0., 5.); @@ -691,7 +649,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) case kHijing_jj75: { - comment = comment.Append("HIJING Jet 75 GeV"); AliGenHijing *gener = HijingStandard(); // impact parameter range gener->SetImpactParameterRange(0., 5.); @@ -708,7 +665,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) case kHijing_jj100: { - comment = comment.Append("HIJING Jet 100 GeV"); AliGenHijing *gener = HijingStandard(); // impact parameter range gener->SetImpactParameterRange(0., 5.); @@ -725,7 +681,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) case kHijing_jj200: { - comment = comment.Append("HIJING Jet 200 GeV"); AliGenHijing *gener = HijingStandard(); // impact parameter range gener->SetImpactParameterRange(0., 5.); @@ -744,7 +699,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) // case kHijing_gj25: { - comment = comment.Append("HIJING Gamma 25 GeV"); AliGenHijing *gener = HijingStandard(); // impact parameter range gener->SetImpactParameterRange(0., 5.); @@ -761,7 +715,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) case kHijing_gj50: { - comment = comment.Append("HIJING Gamma 50 GeV"); AliGenHijing *gener = HijingStandard(); // impact parameter range gener->SetImpactParameterRange(0., 5.); @@ -778,7 +731,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) case kHijing_gj75: { - comment = comment.Append("HIJING Gamma 75 GeV"); AliGenHijing *gener = HijingStandard(); // impact parameter range gener->SetImpactParameterRange(0., 5.); @@ -795,7 +747,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) case kHijing_gj100: { - comment = comment.Append("HIJING Gamma 100 GeV"); AliGenHijing *gener = HijingStandard(); // impact parameter range gener->SetImpactParameterRange(0., 5.); @@ -812,7 +763,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) case kHijing_gj200: { - comment = comment.Append("HIJING Gamma 200 GeV"); AliGenHijing *gener = HijingStandard(); // impact parameter range gener->SetImpactParameterRange(0., 5.); @@ -828,7 +778,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kHijing_pA: { - comment = comment.Append("HIJING pA"); AliGenCocktail *gener = new AliGenCocktail(); @@ -865,7 +814,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kPythia6: { - comment = comment.Append(":Pythia p-p @ 14 TeV"); AliGenPythia *gener = new AliGenPythia(-1); gener->SetMomentumRange(0,999999); gener->SetThetaRange(0., 180.); @@ -878,7 +826,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kPythia6Jets20_24: { - comment = comment.Append(":Pythia jets 20-24 GeV @ 5.5 TeV"); AliGenPythia * gener = new AliGenPythia(-1); gener->SetEnergyCMS(5500.);// Centre of mass energy gener->SetProcess(kPyJets);// Process type @@ -897,7 +844,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kPythia6Jets24_29: { - comment = comment.Append(":Pythia jets 24-29 GeV @ 5.5 TeV"); AliGenPythia * gener = new AliGenPythia(-1); gener->SetEnergyCMS(5500.);// Centre of mass energy gener->SetProcess(kPyJets);// Process type @@ -916,7 +862,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kPythia6Jets29_35: { - comment = comment.Append(":Pythia jets 29-35 GeV @ 5.5 TeV"); AliGenPythia * gener = new AliGenPythia(-1); gener->SetEnergyCMS(5500.);// Centre of mass energy gener->SetProcess(kPyJets);// Process type @@ -935,7 +880,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kPythia6Jets35_42: { - comment = comment.Append(":Pythia jets 35-42 GeV @ 5.5 TeV"); AliGenPythia * gener = new AliGenPythia(-1); gener->SetEnergyCMS(5500.);// Centre of mass energy gener->SetProcess(kPyJets);// Process type @@ -954,7 +898,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kPythia6Jets42_50: { - comment = comment.Append(":Pythia jets 42-50 GeV @ 5.5 TeV"); AliGenPythia * gener = new AliGenPythia(-1); gener->SetEnergyCMS(5500.);// Centre of mass energy gener->SetProcess(kPyJets);// Process type @@ -973,7 +916,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kPythia6Jets50_60: { - comment = comment.Append(":Pythia jets 50-60 GeV @ 5.5 TeV"); AliGenPythia * gener = new AliGenPythia(-1); gener->SetEnergyCMS(5500.);// Centre of mass energy gener->SetProcess(kPyJets);// Process type @@ -992,7 +934,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kPythia6Jets60_72: { - comment = comment.Append(":Pythia jets 60-72 GeV @ 5.5 TeV"); AliGenPythia * gener = new AliGenPythia(-1); gener->SetEnergyCMS(5500.);// Centre of mass energy gener->SetProcess(kPyJets);// Process type @@ -1011,7 +952,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kPythia6Jets72_86: { - comment = comment.Append(":Pythia jets 72-86 GeV @ 5.5 TeV"); AliGenPythia * gener = new AliGenPythia(-1); gener->SetEnergyCMS(5500.);// Centre of mass energy gener->SetProcess(kPyJets);// Process type @@ -1030,7 +970,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kPythia6Jets86_104: { - comment = comment.Append(":Pythia jets 86-104 GeV @ 5.5 TeV"); AliGenPythia * gener = new AliGenPythia(-1); gener->SetEnergyCMS(5500.);// Centre of mass energy gener->SetProcess(kPyJets);// Process type @@ -1049,7 +988,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kPythia6Jets104_125: { - comment = comment.Append(":Pythia jets 105-125 GeV @ 5.5 TeV"); AliGenPythia * gener = new AliGenPythia(-1); gener->SetEnergyCMS(5500.);// Centre of mass energy gener->SetProcess(kPyJets);// Process type @@ -1068,7 +1006,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kPythia6Jets125_150: { - comment = comment.Append(":Pythia jets 125-150 GeV @ 5.5 TeV"); AliGenPythia * gener = new AliGenPythia(-1); gener->SetEnergyCMS(5500.);// Centre of mass energy gener->SetProcess(kPyJets);// Process type @@ -1087,7 +1024,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kPythia6Jets150_180: { - comment = comment.Append(":Pythia jets 150-180 GeV @ 5.5 TeV"); AliGenPythia * gener = new AliGenPythia(-1); gener->SetEnergyCMS(5500.);// Centre of mass energy gener->SetProcess(kPyJets);// Process type @@ -1106,7 +1042,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kD0PbPb5500: { - comment = comment.Append(" D0 in Pb-Pb at 5.5 TeV"); AliGenPythia * gener = new AliGenPythia(10); gener->SetProcess(kPyD0PbPbMNR); gener->SetStrucFunc(kCTEQ4L); @@ -1123,7 +1058,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kCharmSemiElPbPb5500: { - comment = comment.Append(" Charm in Pb-Pb at 5.5 TeV"); AliGenPythia * gener = new AliGenPythia(10); gener->SetProcess(kPyCharmPbPbMNR); gener->SetStrucFunc(kCTEQ4L); @@ -1139,7 +1073,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kBeautySemiElPbPb5500: { - comment = comment.Append(" Beauty in Pb-Pb at 5.5 TeV"); AliGenPythia *gener = new AliGenPythia(10); gener->SetProcess(kPyBeautyPbPbMNR); gener->SetStrucFunc(kCTEQ4L); @@ -1155,7 +1088,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kCocktailTRD: { - comment = comment.Append(" Cocktail for TRD at 5.5 TeV"); AliGenCocktail *gener = new AliGenCocktail(); AliGenParam *jpsi = new AliGenParam(10, @@ -1201,7 +1133,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kPyJJ: { - comment = comment.Append(" Jet-jet at 5.5 TeV"); AliGenPythia *gener = new AliGenPythia(-1); gener->SetEnergyCMS(5500.); gener->SetProcess(kPyJets); @@ -1215,7 +1146,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kPyGJ: { - comment = comment.Append(" Gamma-jet at 5.5 TeV"); AliGenPythia *gener = new AliGenPythia(-1); gener->SetEnergyCMS(5500.); gener->SetProcess(kPyDirectGamma); @@ -1230,7 +1160,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kMuonCocktailCent1: { - comment = comment.Append(" Muon Cocktail Cent1"); AliGenMUONCocktail * gener = new AliGenMUONCocktail(); gener->SetPtRange(1.0,100.); // Transverse momentum range gener->SetPhiRange(0.,360.); // Azimuthal angle range @@ -1245,7 +1174,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kMuonCocktailPer1: { - comment = comment.Append(" Muon Cocktail Per1"); AliGenMUONCocktail * gener = new AliGenMUONCocktail(); gener->SetPtRange(1.0,100.); // Transverse momentum range gener->SetPhiRange(0.,360.); // Azimuthal angle range @@ -1260,7 +1188,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kMuonCocktailPer4: { - comment = comment.Append(" Muon Cocktail Per4"); AliGenMUONCocktail * gener = new AliGenMUONCocktail(); gener->SetPtRange(1.0,100.); // Transverse momentum range gener->SetPhiRange(0.,360.); // Azimuthal angle range @@ -1275,7 +1202,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kMuonCocktailCent1HighPt: { - comment = comment.Append(" Muon Cocktail HighPt Cent1"); AliGenMUONCocktail * gener = new AliGenMUONCocktail(); gener->SetPtRange(1.0,100.); // Transverse momentum range gener->SetPhiRange(0.,360.); // Azimuthal angle range @@ -1290,7 +1216,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kMuonCocktailPer1HighPt : { - comment = comment.Append(" Muon Cocktail HighPt Per1"); AliGenMUONCocktail * gener = new AliGenMUONCocktail(); gener->SetPtRange(1.0,100.); // Transverse momentum range gener->SetPhiRange(0.,360.); // Azimuthal angle range @@ -1305,7 +1230,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kMuonCocktailPer4HighPt: { - comment = comment.Append(" Muon Cocktail HighPt Per4"); AliGenMUONCocktail * gener = new AliGenMUONCocktail(); gener->SetPtRange(1.0,100.); // Transverse momentum range gener->SetPhiRange(0.,360.); // Azimuthal angle range @@ -1320,7 +1244,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kMuonCocktailCent1Single: { - comment = comment.Append(" Muon Cocktail Single Cent1"); AliGenMUONCocktail * gener = new AliGenMUONCocktail(); gener->SetPtRange(1.0,100.); // Transverse momentum range gener->SetPhiRange(0.,360.); // Azimuthal angle range @@ -1335,7 +1258,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kMuonCocktailPer1Single : { - comment = comment.Append(" Muon Cocktail Single Per1"); AliGenMUONCocktail * gener = new AliGenMUONCocktail(); gener->SetPtRange(1.0,100.); // Transverse momentum range gener->SetPhiRange(0.,360.); // Azimuthal angle range @@ -1350,7 +1272,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kMuonCocktailPer4Single: { - comment = comment.Append(" Muon Cocktail Single Per4"); AliGenMUONCocktail * gener = new AliGenMUONCocktail(); gener->SetPtRange(1.0,100.); // Transverse momentum range gener->SetPhiRange(0.,360.); // Azimuthal angle range @@ -1365,7 +1286,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kFMD1Flat: { - comment = comment.Append(" Flat in FMD1 range"); AliGenBox* gener = new AliGenBox(2000); gener->SetPart(kPiPlus); gener->SetMomentumRange(3,4); @@ -1376,7 +1296,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kFMD2Flat: { - comment = comment.Append(" Flat in FMD2 range"); AliGenBox* gener = new AliGenBox(10); gener->SetPart(kPiPlus); gener->SetMomentumRange(3,4); @@ -1387,7 +1306,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kFMD3Flat: { - comment = comment.Append(" Flat in FMD3 range"); AliGenBox* gener = new AliGenBox(2000); gener->SetPart(kPiPlus); gener->SetMomentumRange(3,4); @@ -1398,7 +1316,6 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment) break; case kFMDFlat: { - comment = comment.Append(" Flat in FMD range"); AliGenCocktail* gener = new AliGenCocktail(); gener->SetMomentumRange(3,4); gener->SetPhiRange(0, 360); diff --git a/FMD/Reconstruct.C b/FMD/Reconstruct.C index a1a051bb4c2..5a73e1f4155 100644 --- a/FMD/Reconstruct.C +++ b/FMD/Reconstruct.C @@ -26,16 +26,16 @@ void Reconstruct() { // Debug the FMD. - // AliLog::SetModuleDebugLevel("FMD", 1); + AliLog::SetModuleDebugLevel("FMD", 1); // To reconstruct raw data from FDR-I, please enable below lines: // AliFMDParameters::Instance()->UseRcuTrailer(false); // AliFMDParameters::Instance()->UseCompleteHeader(false); - TFile* magF = TFile::Open("mag.root", "READ"); - AliMagF* mag = static_cast(magF->Get("mag")); - if (!mag) return; - AliTracker::SetFieldMap(mag, true); + // TFile* magF = TFile::Open("mag.root", "READ"); + // AliMagF* mag = static_cast(magF->Get("mag")); + // if (!mag) return; + // AliTracker::SetFieldMap(mag, true); AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT"); AliCDBManager::Instance()->SetRun(0); @@ -45,8 +45,8 @@ Reconstruct() rec.SetRunReconstruction("FMD"); rec.SetRunTracking(""); rec.SetFillESD("FMD"); - rec.SetRunQA(""); - rec.SetInput("."); + rec.SetRunQA(":"); + rec.SetInput("raw.root"); // rec.SetRecoParam("TOF", new AliTOFRecoParam()); rec.Run(); diff --git a/FMD/Simulate.C b/FMD/Simulate.C index da627ef8cda..7415b8b5986 100644 --- a/FMD/Simulate.C +++ b/FMD/Simulate.C @@ -28,9 +28,9 @@ Simulate(Int_t n=1) sim.SetConfigFile("$(ALICE_ROOT)/FMD/Config.C"); sim.SetMakeSDigits("FMD"); sim.SetMakeDigits("FMD"); - sim.SetWriteRawData("FMD"); + sim.SetWriteRawData("FMD", "raw.root"); // sim.SetMakeDigitsFromHits("FMD"); - sim.SetRunQA(":"); + sim.SetRunQA("FMD:ALL"); TStopwatch w; w.Start(); sim.Run(n); diff --git a/FMD/Survey_943928_FMD.txt b/FMD/Survey_943928_FMD.txt new file mode 100644 index 00000000000..4cb8b95d374 --- /dev/null +++ b/FMD/Survey_943928_FMD.txt @@ -0,0 +1,49 @@ +> Title: + +Measurement of targets on FMD2 + +> Date: + +12.08.2008 + +> Subdetector: + +FMD1 + +> Report URL: + +https://edms.cern.ch/document/943928 + +> Version: + +1 + +> General Observations: + +The V0A detector has been installed and first measured on 18.07.2008. +Afterwards it has been moved close to final position and measured on +25.07.2008. The final alignment of the V0A detector has been done in +three iterations (called sessions) on 12.08.2008. + +> Coordinate System: + +ALICEPH + +> Units: + +mm + +> Nr Columns: + +7 + +> Column Names: + +Point Name,XPH,YPH,ZPH,Point Type,Target Used,Precision(mm) + + +> Data: +V0L_ICB 127.73 -220.41 3246.58 M Y 1 +V0L_ICT 127.37 219.48 3247.24 M Y 1 +V0L_OCB -126.99 -220.40 3247.27 M Y 1 +V0L_OCT -126.99 219.51 3247.49 M Y 1 diff --git a/FMD/Survey_976326_FMD.txt b/FMD/Survey_976326_FMD.txt new file mode 100644 index 00000000000..618f2473e8f --- /dev/null +++ b/FMD/Survey_976326_FMD.txt @@ -0,0 +1,51 @@ +> Title: + +Measurement of targets on FMD2 + +> Date: + +05.06.2008 + +> Subdetector: + +FMD2 + +> Report URL: + +https://edms.cern.ch/document/976326 + +> Version: + +1 + +> General Observations: + +During the installation period of detectors in the Space Frame, +several targets on the FMD2 have been measured regularly in order to +detect relative movements of the TPC. + + +> Coordinate System: + +ALICEPH + +> Units: + +mm + +> Nr Columns: + +7 + +> Column Names: + +Point Name,XPH,YPH,ZPH,Point Type,Target Used,Precision(mm) + + +> Data: + +FMD2_IBOT 192.7 -257.4 892.4 M Y 1 +FMD2_OBOT -194.5 -257.2 894.3 M Y 1 +FMD2_OBOTM -271.7 -122.0 892.1 M Y 1 +FMD2_ITOP 246.9 100.4 892.5 M Y 1 +FMD2_OTOP -232.4 113.6 891.7 M Y 1 diff --git a/FMD/scripts/DrawBothDigits.C b/FMD/scripts/DrawBothDigits.C new file mode 100644 index 00000000000..3a018e56cff --- /dev/null +++ b/FMD/scripts/DrawBothDigits.C @@ -0,0 +1,110 @@ +//____________________________________________________________________ +// +// $Id: DrawBothDigits.C 20907 2007-09-25 08:44:03Z cholm $ +// +// Script that contains a class to draw eloss from hits, versus ADC +// counts from digits, using the AliFMDInputHits class in the util library. +// +// It draws the energy loss versus the p/(mq^2). It can be overlayed +// with the Bethe-Bloc curve to show how the simulation behaves +// relative to the expected. +// +// Use the script `Compile.C' to compile this class using ACLic. +// +#include +#include +#include +#include +#include +#include +#include +#include +#include + +/** @class DrawBothDigits + @brief Draw hit energy loss versus digit ADC + @code + Root> .L Compile.C + Root> Compile("DrawBothDigits.C") + Root> DrawBothDigits c + Root> c.Run(); + @endcode + @ingroup FMD_script + */ +class DrawBothDigits : public AliFMDInput +{ +private: + TH2F* fTrackNos; // Histogram + AliFMDEdepMap fCache; +public: + //__________________________________________________________________ + DrawBothDigits(Int_t max=300) + : AliFMDInput("galice.root") + { + AddLoad(kDigits); + AddLoad(kSDigits); + fTrackNos = new TH2F("trackNos", "Track numbers", + max+1, -1.5, max-.5, max+1, -1.5, max-.5); + fTrackNos->SetXTitle("Digit track"); + fTrackNos->SetYTitle("SDigit track"); + } + //__________________________________________________________________ + Bool_t Begin(Int_t evno) + { + fCache.Reset(); + return AliFMDInput::Begin(evno); + } + //__________________________________________________________________ + Bool_t ProcessSDigit(AliFMDSDigit* sdigit) + { + if (!sdigit) return kTRUE; + AliFMDEdepHitPair& entry = fCache(sdigit->Detector(), + sdigit->Ring(), + sdigit->Sector(), + sdigit->Strip()); + entry.fLabels.Set(sdigit->GetNTrack()); + Info("ProcessSDigit", "Got %d SDigit tracks", sdigit->GetNTrack()); + for (size_t i = 0; i < sdigit->GetNTrack(); i++) + entry.fLabels.fArray[i] = sdigit->GetTrack(i); + return kTRUE; + } + //__________________________________________________________________ + Bool_t ProcessDigit(AliFMDDigit* digit) + { + if (!digit) return kTRUE; + AliFMDEdepHitPair& entry = fCache(digit->Detector(), + digit->Ring(), + digit->Sector(), + digit->Strip()); + TArrayI& stracks = entry.fLabels; + Info("ProcessDigit", "Got %d SDigit tracks, and %d Digit tracks", + stracks.fN, digit->GetNTrack()); + for (Int_t i = 0; (i < stracks.fN || i < digit->GetNTrack()); i++) { + Int_t strack = (i < stracks.fN ? stracks.fArray[i] : -1); + Int_t dtrack = (i < digit->GetNTrack() ? digit->GetTrack(i) : -1); + fTrackNos->Fill(strack, dtrack); + } + return kTRUE; + } + //__________________________________________________________________ + Bool_t Finish() + { + gStyle->SetPalette(1); + gStyle->SetOptTitle(0); + gStyle->SetCanvasColor(0); + gStyle->SetCanvasBorderSize(0); + gStyle->SetPadColor(0); + gStyle->SetPadBorderSize(0); + fTrackNos->SetStats(kFALSE); + // fTrackNos->Scale(1. / fAdc->GetEntries()); + fTrackNos->Draw("colz"); + return kTRUE; + } + + ClassDef(DrawBothDigits,0); +}; + +//____________________________________________________________________ +// +// EOF +// diff --git a/FMD/scripts/DrawDigits.C b/FMD/scripts/DrawDigits.C index 107ba89cf10..62ddea5b4d5 100644 --- a/FMD/scripts/DrawDigits.C +++ b/FMD/scripts/DrawDigits.C @@ -38,6 +38,7 @@ private: public: //__________________________________________________________________ DrawDigits(Int_t m=1100, Double_t amin=-0.5, Double_t amax=1023.5) + : AliFMDInput("galice.root") { AddLoad(kDigits); fAdc = new TH1D("adc", "ADC", m, amin, amax); @@ -49,6 +50,7 @@ public: { if (!digit) return kTRUE; fAdc->Fill(digit->Counts()); + digit->Print("l"); return kTRUE; } //__________________________________________________________________ diff --git a/FMD/scripts/DrawESD.C b/FMD/scripts/DrawESD.C index e5e2b7cf6b8..94924fca07b 100644 --- a/FMD/scripts/DrawESD.C +++ b/FMD/scripts/DrawESD.C @@ -64,20 +64,24 @@ public: return AliFMDInput::Begin(ev); } //__________________________________________________________________ - Bool_t ProcessESD(UShort_t /* det */, - Char_t /* ring */, - UShort_t /* sec */, - UShort_t /* strip */, + Bool_t ProcessESD(UShort_t det, + Char_t rng, + UShort_t sec, + UShort_t str, Float_t /* eta */, Float_t mult) { // Cache the energy loss - if (mult/fCorr > 0.01) fMult->Fill(mult/fCorr); + if (mult > 0) + Info("ProcessESD", "FMD%d%c[%2d,%3d]=%f", det, rng, sec, str, mult); + if (mult/fCorr > 0.001) fMult->Fill(mult/fCorr); return kTRUE; } //__________________________________________________________________ TF1* FitPeak(Int_t n, TH1D* hist, Double_t min, Double_t& max) { + if (TMath::Abs(max-min) < .25) return 0; + std::cout << "Fit peack in range " << min << " to " << max << std::endl; TF1* l = new TF1(Form("l%02d", n), "landau", min, max); hist->GetXaxis()->SetRangeUser(0, 4); hist->Fit(l, "0Q", "", min, max); @@ -112,6 +116,7 @@ public: //__________________________________________________________________ Bool_t Finish() { + Info("Finish", "Will draw results"); gStyle->SetPalette(1); gStyle->SetOptTitle(0); gStyle->SetOptFit(1111111); @@ -120,16 +125,20 @@ public: gStyle->SetPadColor(0); gStyle->SetPadBorderSize(0); + if (fMult->GetEntries() <= 0) return kFALSE; + TCanvas* c = new TCanvas("c", "C"); c->cd(); - c->SetLogy(); + // c->SetLogy(); fMult->GetXaxis()->SetRangeUser(0,4); fMult->Scale(1. / fMult->GetEntries()); fMult->SetStats(kFALSE); fMult->SetFillColor(2); fMult->SetFillStyle(3001); fMult->Draw("hist e"); - + + return kTRUE; + Double_t mean, rms; MaxInRange(fMult, 0.2, mean, rms); Double_t x1 = mean-rms/2; // .75; // .8; // .65 / fCorr; @@ -139,49 +148,61 @@ public: MaxInRange(fMult, x2, mean, rms); Double_t x3 = mean + rms; TF1* l2 = FitPeak(2, fMult, x2, x3); - Double_t diff = l2->GetParameter(1)-l1->GetParameter(1); - TF1* f = new TF1("user", "landau(0)+landau(3)", x1, x3); - + TF1* f = 0; + Double_t diff = 0; + if (l2) { + diff = l2->GetParameter(1)-l1->GetParameter(1); + f = new TF1("user", "landau(0)+landau(3)", x1, x3); + f->SetParNames("A_{1}", "Mpv_{1}", "#sigma_{1}", + "A_{2}", "Mpv_{2}", "#sigma_{2}", + "A_{3}", "Mpv_{3}", "#sigma_{3}"); + f->SetParameters(l1->GetParameter(0), + l1->GetParameter(1), + l1->GetParameter(2), + l2 ? l2->GetParameter(0) : 0, + l2 ? l2->GetParameter(1) : 0, + l2 ? l2->GetParameter(2) : 0, + l2->GetParameter(0)/10, + l2->GetParameter(1) + diff, + l2->GetParameter(2)); + } + else { + x3 = x2; + f = new TF1("usr", "landau", x1, x3); + } + + std::cout << "Range: " << x1 << "-" << x3 << std::endl; + fMult->GetXaxis()->SetRangeUser(0, 4); - f->SetParNames("A_{1}", "Mpv_{1}", "#sigma_{1}", - "A_{2}", "Mpv_{2}", "#sigma_{2}", - "A_{3}", "Mpv_{3}", "#sigma_{3}"); - f->SetParameters(l1->GetParameter(0), - l1->GetParameter(1), - l1->GetParameter(2), - l2->GetParameter(0), - l2->GetParameter(1), - l2->GetParameter(2), - l2->GetParameter(0)/10, - l2->GetParameter(1) + diff, - l2->GetParameter(2)); - fMult->Fit(f, "0Q", "", x1, x3); - fMult->Fit(f, "ME0", "E1", x1, x3); + fMult->Fit(f, "0QR", "", x1, x3); + fMult->Fit(f, "ME0R", "E1", x1, x3); fMult->DrawClone("same hist"); l1->SetLineColor(3); l1->SetLineWidth(2); l1->SetRange(0, 4); l1->Draw("same"); - l2->SetLineColor(4); - l2->SetLineWidth(2); - l2->SetRange(0, 4); - l2->Draw("same"); + if (l2) { + l2->SetLineColor(4); + l2->SetLineWidth(2); + l2->SetRange(0, 4); + l2->Draw("same"); + } f->SetLineWidth(2); f->SetRange(0, 4); f->Draw("same"); TLegend* l = new TLegend(0.6, 0.6, .89, .89); l->AddEntry(l1, "1 particle Landau", "l"); - l->AddEntry(l2, "2 particle Landau", "l"); + if (l2) l->AddEntry(l2, "2 particle Landau", "l"); l->AddEntry(f, "1+2 particle Landau", "l"); l->SetFillColor(0); l->Draw("same"); -#if 1 +#if 0 c = new TCanvas("c2", "Landaus"); - c->SetLogy(); + // c->SetLogy(); fMult->DrawClone("axis"); f->Draw("same"); Double_t* p1 = f->GetParameters(); diff --git a/FMD/scripts/DrawHits.C b/FMD/scripts/DrawHits.C index 2d6fabb3c3e..c314038cf4e 100644 --- a/FMD/scripts/DrawHits.C +++ b/FMD/scripts/DrawHits.C @@ -90,7 +90,7 @@ public: std::cout << "No track" << std::endl; return kFALSE; } - if (!p->IsPrimary()) return kTRUE; + // if (!p->IsPrimary()) return kTRUE; if (hit->IsStop()) return kTRUE; Float_t x = hit->P(); @@ -103,7 +103,7 @@ public: y /= q * q; fElossVsPMQ->Fill(x, y); fEloss->Fill(y); - + fNHits++; return kTRUE; } //__________________________________________________________________ @@ -155,24 +155,27 @@ public: c = new TCanvas("eloss", "Energy loss per unit material"); // c->SetLogx(); - c->SetLogy(); - fEloss->Scale(1. / fEloss->GetEntries()); - fEloss->GetXaxis()->SetRangeUser(1, 10); - fEloss->Fit("landau", "", "", 1, 10); - TF1* land = fEloss->GetFunction("landau"); - land->SetLineWidth(2); - Double_t max = fEloss->GetMaximum(); - TGraph* resp = GetResp(); - Double_t* x = resp->GetX(); - Double_t* y = resp->GetY(); - TGraph* g = new TGraph(resp->GetN()); - TGraph* co = GetCorr(); - std::cout << "Correction factor: " << co->Eval(fBetaGammaMip) << std::endl; - Double_t xs = fRho; // * 1.19; // / - for (Int_t i = 0; i < g->GetN(); i++) - g->SetPoint(i, x[i] * xs, y[i] * max); - g->Draw("C same"); - + if (fEloss->GetEntries() != 0) { + c->SetLogy(); + fEloss->Scale(1. / fEloss->GetEntries()); + fEloss->GetXaxis()->SetRangeUser(1, 10); + fEloss->Fit("landau", "", "", 1, 10); + TF1* land = fEloss->GetFunction("landau"); + land->SetLineWidth(2); + Double_t max = fEloss->GetMaximum(); + TGraph* resp = GetResp(); + Double_t* x = resp->GetX(); + Double_t* y = resp->GetY(); + TGraph* g = new TGraph(resp->GetN()); + TGraph* co = GetCorr(); + std::cout << "Correction factor: " << co->Eval(fBetaGammaMip) + << std::endl; + Double_t xs = fRho; // * 1.19; // / + for (Int_t i = 0; i < g->GetN(); i++) + g->SetPoint(i, x[i] * xs, y[i] * max); + g->Draw("C same"); + } + l = new TLegend(.6, .6, .89, .89); l->AddEntry(fEloss, fEloss->GetTitle(), "lf"); l->AddEntry(land, "Landau fit", "l"); diff --git a/FMD/scripts/DrawHitsDigits.C b/FMD/scripts/DrawHitsDigits.C index 94c8df27685..e5bd31b8806 100644 --- a/FMD/scripts/DrawHitsDigits.C +++ b/FMD/scripts/DrawHitsDigits.C @@ -65,6 +65,7 @@ public: //__________________________________________________________________ Bool_t ProcessHit(AliFMDHit* hit, TParticle*) { + // Info("ProcessHit", "Processing hit %p", hit); // Cache the energy loss if (!hit) return kFALSE; UShort_t det = hit->Detector(); @@ -77,11 +78,13 @@ public: } fMap(det, rng, sec, str).fEdep += hit->Edep(); fMap(det, rng, sec, str).fN++; + // hit->Print(); return kTRUE; } //__________________________________________________________________ Bool_t ProcessDigit(AliFMDDigit* digit) { + // Info("ProcessDigit", "Processing digit %p", digit); if (!digit) return kFALSE; UShort_t det = digit->Detector(); Char_t rng = digit->Ring(); @@ -92,6 +95,10 @@ public: return kFALSE; } fElossVsAdc->Fill(fMap(det, rng, sec, str).fEdep, digit->Counts()); + if (fMap(det, rng, sec, str).fEdep > 0) + Info("ProcessDigit", "FMD%d%c[%2d,%3d] Edep=%8.5f -> ADC=0x%03x", + det, rng, sec, str, + fMap(det, rng, sec, str).fEdep, digit->Counts()); return kTRUE; } //__________________________________________________________________ diff --git a/FMD/scripts/DrawHitsSDigits.C b/FMD/scripts/DrawHitsSDigits.C new file mode 100644 index 00000000000..6a51ec098a7 --- /dev/null +++ b/FMD/scripts/DrawHitsSDigits.C @@ -0,0 +1,121 @@ +//____________________________________________________________________ +// +// $Id: DrawHitsSDigits.C 20907 2007-09-25 08:44:03Z cholm $ +// +// Script that contains a class to draw eloss from hits, versus ADC +// counts from sdigits, using the AliFMDInputHits class in the util library. +// +// It draws the energy loss versus the p/(mq^2). It can be overlayed +// with the Bethe-Bloc curve to show how the simulation behaves +// relative to the expected. +// +// Use the script `Compile.C' to compile this class using ACLic. +// +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +/** @class DrawHitsSDigits + @brief Draw hit energy loss versus sdigit ADC + @code + Root> .L Compile.C + Root> Compile("DrawHitsSDigits.C") + Root> DrawHitsSDigits c + Root> c.Run(); + @endcode + @ingroup FMD_script + */ +class DrawHitsSDigits : public AliFMDInput +{ +private: + TH2D* fElossVsAdc; // Histogram + AliFMDEdepMap fMap; +public: + //__________________________________________________________________ + DrawHitsSDigits(Int_t n=900, Double_t emin=1e-3, Double_t emax=10, + Int_t m=1100, Double_t amin=-0.5, Double_t amax=1099.5) + { + AddLoad(kSDigits); + AddLoad(kHits); + TArrayF eloss(MakeLogScale(n, emin, emax)); + TArrayF adcs(m+1); + adcs[0] = amin; + for (Int_t i = 1; i < m+1; i++) adcs[i] = adcs[i-1] + (amax-amin)/m; + fElossVsAdc = new TH2D("bad", "#Delta E vs. ADC", + eloss.fN-1, eloss.fArray, adcs.fN-1, adcs.fArray); + fElossVsAdc->SetXTitle("#Delta E/#Delta x [MeV/cm]"); + fElossVsAdc->SetYTitle("ADC value"); + } + //__________________________________________________________________ + /** Begining of event + @param ev Event number + @return @c false on error */ + Bool_t Begin(Int_t ev) + { + fMap.Reset(); + return AliFMDInput::Begin(ev); + } + //__________________________________________________________________ + Bool_t ProcessHit(AliFMDHit* hit, TParticle*) + { + Info("ProcessHit", "Processing hit %p", hit); + // Cache the energy loss + if (!hit) return kFALSE; + UShort_t det = hit->Detector(); + Char_t rng = hit->Ring(); + UShort_t sec = hit->Sector(); + UShort_t str = hit->Strip(); + if (str > 511) { + AliWarning(Form("Bad strip number %d in hit", str)); + return kTRUE; + } + fMap(det, rng, sec, str).fEdep += hit->Edep(); + fMap(det, rng, sec, str).fN++; + hit->Print(); + return kTRUE; + } + //__________________________________________________________________ + Bool_t ProcessSDigit(AliFMDSDigit* sdigit) + { + Info("ProcessSDigit", "Processing sdigit %p", sdigit); + if (!sdigit) return kFALSE; + UShort_t det = sdigit->Detector(); + Char_t rng = sdigit->Ring(); + UShort_t sec = sdigit->Sector(); + UShort_t str = sdigit->Strip(); + if (str > 511) { + AliWarning(Form("Bad strip number %d in sdigit", str)); + return kFALSE; + } + fElossVsAdc->Fill(fMap(det, rng, sec, str).fEdep, sdigit->Counts()); + sdigit->Print(); + return kTRUE; + } + //__________________________________________________________________ + Bool_t Finish() + { + gStyle->SetPalette(1); + gStyle->SetOptTitle(0); + gStyle->SetCanvasColor(0); + gStyle->SetCanvasBorderSize(0); + gStyle->SetPadColor(0); + gStyle->SetPadBorderSize(0); + fElossVsAdc->SetStats(kFALSE); + fElossVsAdc->Draw("COLZ"); + return kTRUE; + } + + ClassDef(DrawHitsSDigits,0); +}; + +//____________________________________________________________________ +// +// EOF +// diff --git a/FMD/scripts/DrawSDigits.C b/FMD/scripts/DrawSDigits.C index 53d404a1e76..9ba1a86a4f5 100644 --- a/FMD/scripts/DrawSDigits.C +++ b/FMD/scripts/DrawSDigits.C @@ -68,6 +68,7 @@ public: //__________________________________________________________________ DrawSDigits(Int_t m=1100, Double_t amin=-0.5, Double_t amax=1023.5) + : AliFMDInput("galice.root") { AddLoad(kSDigits); AddLoad(kGeometry); @@ -121,6 +122,7 @@ public: { if (!digit) return kTRUE; fAdc->Fill(digit->Counts()); + digit->Print("lp"); if (digit->NParticles() == 0) return kTRUE; @@ -143,7 +145,6 @@ public: Double_t ratio = digit->NPrimaries() / digit->NParticles(); if (phi < 0) phi += 360; fPrimRatio[primIdx]->Fill(eta, phi, ratio); - return kTRUE; } //__________________________________________________________________ diff --git a/FMD/scripts/TestRaw2SDigits.C b/FMD/scripts/TestRaw2SDigits.C new file mode 100644 index 00000000000..874c70c57bd --- /dev/null +++ b/FMD/scripts/TestRaw2SDigits.C @@ -0,0 +1,10 @@ +void +TestRaw2SDigits() +{ + AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT"); + AliCDBManager::Instance()->SetRun(0); + AliLog::SetModuleDebugLevel("FMD", 2); + AliSimulation sim; + sim.SetConfigFile("$ALICE_ROOT/FMD/Config.C"); + sim.ConvertRaw2SDigits("raw.root", 0); +} -- 2.39.3