From 1298ee54fea92fafb5365f7bd5417f19327016b0 Mon Sep 17 00:00:00 2001 From: hdalsgaa Date: Wed, 19 Jan 2011 09:51:28 +0000 Subject: [PATCH] Added dead-channel (acceptance) stuff --- .../analysis2/AliFMDCorrDeadChannels.cxx | 321 ++++++++++++++++++ .../analysis2/AliFMDCorrDeadChannels.h | 210 ++++++++++++ .../scripts/MakeDeadChannelCorrection.C | 231 +++++++++++++ 3 files changed, 762 insertions(+) create mode 100644 PWG2/FORWARD/analysis2/AliFMDCorrDeadChannels.cxx create mode 100644 PWG2/FORWARD/analysis2/AliFMDCorrDeadChannels.h create mode 100644 PWG2/FORWARD/analysis2/scripts/MakeDeadChannelCorrection.C diff --git a/PWG2/FORWARD/analysis2/AliFMDCorrDeadChannels.cxx b/PWG2/FORWARD/analysis2/AliFMDCorrDeadChannels.cxx new file mode 100644 index 00000000000..0547e58e695 --- /dev/null +++ b/PWG2/FORWARD/analysis2/AliFMDCorrDeadChannels.cxx @@ -0,0 +1,321 @@ +// +// This class contains the dead channels correction +// +// +#include "AliFMDCorrDeadChannels.h" +#include +#include +#include +#include + +//____________________________________________________________________ +AliFMDCorrDeadChannels::AliFMDCorrDeadChannels() + : fRingArray(), + fVertexAxis(0,0,0) +{ + // + // Default constructor + // + fRingArray.SetOwner(kTRUE); + fRingArray.SetName("rings"); + fVertexAxis.SetName("vtxAxis"); + fVertexAxis.SetTitle("v_{z} [cm]"); + +} +//____________________________________________________________________ +AliFMDCorrDeadChannels::AliFMDCorrDeadChannels(const + AliFMDCorrDeadChannels& o) + : TObject(o), + fRingArray(o.fRingArray), + fVertexAxis(o.fVertexAxis.GetNbins(), o.fVertexAxis.GetXmin(), + o.fVertexAxis.GetXmax()) +{ + // + // Copy constructor + // + // Parameters: + // o Object to copy from + // + fVertexAxis.SetName("vtxAxis"); + fVertexAxis.SetTitle("v_{z} [cm]"); +} +//____________________________________________________________________ +AliFMDCorrDeadChannels::~AliFMDCorrDeadChannels() +{ + // + // Destructor + // + // + fRingArray.Clear(); +} +//____________________________________________________________________ +AliFMDCorrDeadChannels& +AliFMDCorrDeadChannels::operator=(const AliFMDCorrDeadChannels& o) +{ + // + // Assignment operator + // + // Parameters: + // o Object to assign from + // + // Return: + // Reference to this object + // + fRingArray = o.fRingArray; + SetVertexAxis(o.fVertexAxis); + + return *this; +} +//____________________________________________________________________ +TH2D* +AliFMDCorrDeadChannels::GetCorrection(UShort_t d, Char_t r, Double_t v) const +{ + // + // Get the dead channels correction @f$ c_{r,v}@f$ + // + // Parameters: + // d Detector number (1-3) + // r Ring identifier (I or O) + // v Primary interaction point @f$z@f$ coordinate + // + // Return: + // The correction @f$ c_{r,v}@f$ + // + Int_t b = FindVertexBin(v); + if (b <= 0) return 0; + return GetCorrection(d, r, UShort_t(b)); +} +//____________________________________________________________________ +TH2D* +AliFMDCorrDeadChannels::GetCorrection(UShort_t d, Char_t r, UShort_t b) const +{ + // + // Get the dead channels correction @f$ c_{r,v}@f$ + // + // Parameters: + // d Detector number (1-3) + // r Ring identifier (I or O) + // b Bin corresponding to the primary interaction point + // @f$z@f$ coordinate (1 based) + // + // Return: + // The correction @f$ c_{r,v}@f$ + // + TObjArray* ringArray = GetRingArray(d, r); + if (!ringArray) return 0; + + if (b <= 0 || b > ringArray->GetEntriesFast()) { + AliWarning(Form("vertex bin %d out of range [1,%d]", + b, ringArray->GetEntriesFast())); + return 0; + } + + TObject* o = ringArray->At(b-1); + if (!o) { + AliWarning(Form("No dead channels map found for FMD%d%c in vertex bin %d", + d,r,b)); + return 0; + } + return static_cast(o); +} + +//____________________________________________________________________ +Int_t +AliFMDCorrDeadChannels::FindVertexBin(Double_t v) const +{ + // + // Find the vertex bin that corresponds to the passed vertex + // + // Parameters: + // vertex The interaction points @f$z@f$-coordinate + // + // Return: + // Vertex bin in @f$[1,N_{\mbox{vertex}}]@f$ or negative if + // out of range + // + if (fVertexAxis.GetNbins() <= 0) { + AliWarning("No vertex array defined"); + return 0; + } + Int_t bin = const_cast(fVertexAxis).FindBin(v); + if (bin <= 0 || bin > fVertexAxis.GetNbins()) { + AliWarning(Form("vertex %+8.4f out of range [%+8.4f,%+8.4f]", + v, fVertexAxis.GetXmin(), fVertexAxis.GetXmax())); + return 0; + } + return bin; +} +//____________________________________________________________________ +Int_t +AliFMDCorrDeadChannels::GetRingIndex(UShort_t d, Char_t r) const +{ + // + // Get the index corresponding to the given ring + // + // Parameters: + // d Detector + // r Ring + // + // Return: + // Index (0 based) or negative in case of errors + // + switch (d) { + case 1: return 0; + case 2: return (r == 'I' || r == 'i' ? 1 : 2); break; + case 3: return (r == 'I' || r == 'i' ? 3 : 4); break; + } + AliWarning(Form("Index for FMD%d%c not found", d, r)); + return -1; +} +//____________________________________________________________________ +TObjArray* +AliFMDCorrDeadChannels::GetRingArray(UShort_t d, Char_t r) const +{ + // + // Get the ring array corresponding to the specified ring + // + // Parameters: + // d Detector + // r Ring + // + // Return: + // Pointer to ring array, or null in case of problems + // + Int_t idx = GetRingIndex(d,r); + if (idx < 0) return 0; + + TObject* o = fRingArray.At(idx); + if (!o) { + AliWarning(Form("No array found for FMD%d%c", d, r)); + return 0; + } + + return static_cast(o); +} +//____________________________________________________________________ +TObjArray* +AliFMDCorrDeadChannels::GetOrMakeRingArray(UShort_t d, Char_t r) +{ + // + // Get the ring array corresponding to the specified ring + // + // Parameters: + // d Detector + // r Ring + // + // Return: + // Pointer to ring array, or newly created container + // + Int_t idx = GetRingIndex(d,r); + if (idx < 0) return 0; + + TObject* o = fRingArray.At(idx); + if (!o) { + TObjArray* a = new TObjArray(fVertexAxis.GetNbins()); + a->SetName(Form("FMD%d%c", d, r)); + a->SetOwner(kTRUE); + fRingArray.AddAtAndExpand(a, idx); + return a; + } + + return static_cast(fRingArray.At(idx)); +} + +//____________________________________________________________________ +Bool_t +AliFMDCorrDeadChannels::SetCorrection(UShort_t d, Char_t r, + UShort_t b, TH2D* h) +{ + // + // Set the dead channels map correction @f$ m_{r,v}(\eta)@f$ + // Note, that the object takes ownership of the passed pointer. + // + // Parameters: + // d Detector number (1-3) + // r Ring identifier (I or O) + // b Bin corresponding to the primary interaction point + // @f$z@f$ coordinate (1 based) + // h @f$ m_{r,v}(\eta)@f$ + // + // Return: + // true if operation succeeded + // + TObjArray* ringArray = GetOrMakeRingArray(d, r); + if (!ringArray) return false; + + if (b <= 0 || b > fVertexAxis.GetNbins()) { + AliWarning(Form("Vertex bin %3d out of range [1,%3d]", + b, fVertexAxis.GetNbins())); + return false; + } + h->SetName(Form("FMD%d%c_vtxbin%03d", d, r, b)); + h->SetTitle(Form("Dead Channels correction for FMD%d%c " + "in vertex bin %d [%+8.4f,%+8.4f]", + d, r, b, fVertexAxis.GetBinLowEdge(b), + fVertexAxis.GetBinUpEdge(b))); + h->SetXTitle("#eta"); + h->SetYTitle("dN_{ch}/d#eta / sum_i N_{ch,i}"); + h->SetFillStyle(3001); + h->SetDirectory(0); + h->SetStats(0); + ringArray->AddAtAndExpand(h, b-1); + return kTRUE; +} +//____________________________________________________________________ +Bool_t +AliFMDCorrDeadChannels::SetCorrection(UShort_t d, Char_t r, + Double_t v, TH2D* h) +{ + // + // Set the dead channels map correction @f$ m_{r,v}(\eta)@f$. + // Note, that the object takes ownership of the passed pointer. + // + // Parameters: + // d Detector number (1-3) + // r Ring identifier (I or O) + // v Primary interaction point @f$z@f$ coordinate + // h @f$ m_{r,v}(\eta)@f$ + // + // Return: + // true if operation succeeded + // + Int_t b = FindVertexBin(v); + if (b <= 0 || b > fVertexAxis.GetNbins()) { + AliWarning(Form("Vertex %+8.4f out of range [%+8.4f,%+8.4f]", + v, fVertexAxis.GetXmin(), fVertexAxis.GetXmax())); + return false; + } + return SetCorrection(d, r, UShort_t(b), h); +} +//____________________________________________________________________ +void +AliFMDCorrDeadChannels::Browse(TBrowser* b) +{ + // + // Browse this object in the browser + // + // Parameters: + // b + // + b->Add(&fRingArray); + b->Add(&fVertexAxis); +} +//____________________________________________________________________ +void +AliFMDCorrDeadChannels::Print(Option_t* option) const +{ + // + // Print this object + // + // Parameters: + // option + // + std::cout << "Merging efficiency correction" << std::endl; + fRingArray.Print(option); + fVertexAxis.Print(option); +} + +//____________________________________________________________________ +// +// EOF +// diff --git a/PWG2/FORWARD/analysis2/AliFMDCorrDeadChannels.h b/PWG2/FORWARD/analysis2/AliFMDCorrDeadChannels.h new file mode 100644 index 00000000000..c97eb70dd77 --- /dev/null +++ b/PWG2/FORWARD/analysis2/AliFMDCorrDeadChannels.h @@ -0,0 +1,210 @@ +// +// This class contains the dead channel correction +// +// +#ifndef ALIFMDCORRDEADCHANNELS_H +#define ALIFMDCORRDEADCHANNELS_H +#include +#include +#include +class TH2D; + +/** + * This class contains the merging efficiency correction. + * + * The secondary correction is given by + * @f[ + * m_{r,v}(\eta) = + * @f] + * + * These are generated from Monte-Carlo truth and ESD information. + * + * @ingroup pwg2_forward_corr + */ +class AliFMDCorrDeadChannels : public TObject +{ +public: + /** + * Default constructor + */ + AliFMDCorrDeadChannels(); + /** + * Copy constructor + * + * @param o Object to copy from + */ + AliFMDCorrDeadChannels(const AliFMDCorrDeadChannels& o); + /** + * Destructor + * + */ + virtual ~AliFMDCorrDeadChannels(); + /** + * @{ + * @name Get corrections and parameters + */ + /** + * Assignment operator + * + * @param o Object to assign from + * + * @return Reference to this object + */ + AliFMDCorrDeadChannels& operator=(const AliFMDCorrDeadChannels& o); + /** + * Get the secondary correction @f$ c_{r,v}@f$ + * + * @param d Detector number (1-3) + * @param r Ring identifier (I or O) + * @param v Primary interaction point @f$z@f$ coordinate + * + * @return The correction @f$ c_{r,v}@f$ + */ + TH2D* GetCorrection(UShort_t d, Char_t r, Double_t v) const; + /** + * Get the secondary correction @f$ c_{r,v}@f$ + * + * @param d Detector number (1-3) + * @param r Ring identifier (I or O) + * @param b Bin corresponding to the primary interaction point + * @f$z@f$ coordinate (1 based) + * + * @return The correction @f$ c_{r,v}@f$ + */ + TH2D* GetCorrection(UShort_t d, Char_t r, UShort_t b) const; + /** + * Get the vertex axis used + * + * @return vertex axis + */ + const TAxis& GetVertexAxis() const { return fVertexAxis; } + /* @} */ + + /** + * @{ + * @name Set corrections and parameters + */ + /** + * Set the secondary map correction @f$ m_{r,v}(\eta)@f$. + * Note, that the object takes ownership of the passed pointer. + * + * @param d Detector number (1-3) + * @param r Ring identifier (I or O) + * @param v Primary interaction point @f$z@f$ coordinate + * @param h @f$ m_{r,v}(\eta)@f$ + * + * @return true if operation succeeded + */ + Bool_t SetCorrection(UShort_t d, Char_t r, Double_t v, TH2D* h); + /** + * Set the secondary map correction @f$ m_{r,v}(\eta)@f$ + * Note, that the object takes ownership of the passed pointer. + * + * @param d Detector number (1-3) + * @param r Ring identifier (I or O) + * @param b Bin corresponding to the primary interaction point + * @f$z@f$ coordinate (1 based) + * @param h @f$ m_{r,v}(\eta)@f$ + * + * @return true if operation succeeded + */ + Bool_t SetCorrection(UShort_t d, Char_t r, UShort_t b, TH2D* h); + /** + * Set the vertex axis to use + * + * @param axis Vertex axis + */ + void SetVertexAxis(const TAxis& axis); + /** + * Set the vertex axis to use + * + * @param nBins Number of bins + * @param min Minimum + * @param max Maximum + */ + void SetVertexAxis(Int_t nBins, Double_t min, Double_t max); + /* @} */ + + /** + * @{ + * @name Auxiliary member functions + */ + /** + * Declare this as a folder + * + * @return Always true + */ + Bool_t IsFolder() const { return true; } + /** + * Browse this object in the browser + * + * @param b + */ + void Browse(TBrowser* b); + /** + * Print this object + * + * @param option + */ + void Print(Option_t* option="R") const; //*MENU* + /* @} */ +protected: + /** + * Find the vertex bin that corresponds to the passed vertex + * + * @param vertex The interaction points @f$z@f$-coordinate + * + * @return Vertex bin in @f$[1,N_{\mbox{vertex}}]@f$ or negative if + * out of range + */ + Int_t FindVertexBin(Double_t vertex) const; + /** + * Get the index corresponding to the given ring + * + * @param d Detector + * @param r Ring + * + * @return Index (0 based) or negative in case of errors + */ + Int_t GetRingIndex(UShort_t d, Char_t r) const; + /** + * Get the ring array corresponding to the specified ring + * + * @param d Detector + * @param r Ring + * + * @return Pointer to ring array, or null in case of problems + */ + TObjArray* GetRingArray(UShort_t d, Char_t r) const; + /** + * Get the ring array corresponding to the specified ring + * + * @param d Detector + * @param r Ring + * + * @return Pointer to ring array, or newly created container + */ + TObjArray* GetOrMakeRingArray(UShort_t d, Char_t r); + + TObjArray fRingArray; // Array of per-ring, per-vertex 2nd map + TAxis fVertexAxis; // The vertex axis + ClassDef(AliFMDCorrDeadChannels,1); // +}; + +//____________________________________________________________________ +inline void +AliFMDCorrDeadChannels::SetVertexAxis(Int_t nBins, Double_t min, + Double_t max) +{ + fVertexAxis.Set(nBins, min, max); +} +//____________________________________________________________________ +inline void +AliFMDCorrDeadChannels::SetVertexAxis(const TAxis& e) +{ + fVertexAxis.Set(e.GetNbins(), e.GetXmin(), e.GetXmax()); +} +#endif +// Local Variables: +// mode: C++ +// End: diff --git a/PWG2/FORWARD/analysis2/scripts/MakeDeadChannelCorrection.C b/PWG2/FORWARD/analysis2/scripts/MakeDeadChannelCorrection.C new file mode 100644 index 00000000000..d89faafeec7 --- /dev/null +++ b/PWG2/FORWARD/analysis2/scripts/MakeDeadChannelCorrection.C @@ -0,0 +1,231 @@ +//_____________________________________________________________________ +Float_t GetMaxR(Char_t ring) const { + //Get max R of ring + Float_t radius = 0; + if(ring == 'I') + radius = 17.2; + else if(ring == 'O') + radius = 28.0; + else + AliWarning("Unknown ring - must be I or O!"); + + return radius; +} +//_____________________________________________________________________ +Float_t GetMinR(Char_t ring) const{ + //Get min R of ring + Float_t radius = 0; + if(ring == 'I') + radius = 4.5213; + else if(ring == 'O') + radius = 15.4; + else + AliWarning("Unknown ring - must be I or O!"); + + return radius; + +} +//_____________________________________________________________________ +Float_t GetEtaFromStrip(UShort_t det, Char_t ring, UShort_t sec, UShort_t strip, Float_t zvtx) +{ + //Calculate eta from strip with vertex (redundant with AliESDFMD::Eta) + Float_t rad = GetMaxR(ring)-GetMinR(ring); + Float_t nStrips = (ring == 'I' ? 512 : 256); + Float_t segment = rad / nStrips; + Float_t r = GetMinR(ring) + segment*strip; + Float_t z = 0; + Int_t hybrid = sec / 2; + + if(det == 1) { + if(!(hybrid%2)) z = 320.266; else z = 319.766; + } + if(det == 2 && ring == 'I' ) { + if(!(hybrid%2)) z = 83.666; else z = 83.166; + } + if(det == 2 && ring == 'O' ) { + if(!(hybrid%2)) z = 74.966; else z = 75.466; + } + if(det == 3 && ring == 'I' ) { + if(!(hybrid%2)) z = -63.066; else z = -62.566; + } + if(det == 3 && ring == 'O' ) { + if(!(hybrid%2)) z = -74.966; else z = -75.466; + } + + Float_t theta = TMath::ATan2(r,z-zvtx); + Float_t eta = -1*TMath::Log(TMath::Tan(0.5*theta)); + + return eta; +} +//_____________________________________________________________________ +Float_t GetPhiFromSector(UShort_t det, Char_t ring, UShort_t sec) +{ + //Get phi from sector + Int_t nsec = (ring == 'I' ? 20 : 40); + Float_t basephi = 0; + if(det == 1) + basephi = 1.72787594; + if(det == 2 && ring == 'I') + basephi = 0.15707963; + if(det == 2 && ring == 'O') + basephi = 0.078539818; + if(det == 3 && ring == 'I') + basephi = 2.984513044; + if(det == 3 && ring == 'O') + basephi = 3.06305289; + + Float_t step = 2*TMath::Pi() / nsec; + Float_t phi = 0; + if(det == 3) + phi = basephi - sec*step; + else + phi = basephi + sec*step; + + if(phi < 0) + phi = phi +2*TMath::Pi(); + if(phi > 2*TMath::Pi() ) + phi = phi - 2*TMath::Pi(); + + return phi; +} +//_____________________________________________________________________ +void MakeDeadChannelCorrection(Int_t runnumber, Int_t nVtxBins=10, Float_t vtxLow=-10, Float_t vtxHigh=10){ + + gSystem->Load("libANALYSIS"); + gSystem->Load("libANALYSISalice"); + gSystem->Load("libPWG2forward2"); + + Float_t delta = (vtxHigh - vtxLow) / (Float_t)nVtxBins; + + //TGrid::Connect("alien://",0,0,"t"); + AliCDBManager* cdb = AliCDBManager::Instance(); + //cdb->SetDefaultStorage("alien://Folder=/alice/data/2010/OCDB"); + cdb->SetDefaultStorage("local:///home/canute/ALICE/AliRoot/OCDB"); + cdb->SetRun(runnumber); + + TObjArray* fReadArray = new TObjArray(); + fReadArray->SetName("DeadChannels"); + TObjArray* fAllArray = new TObjArray(); + fAllArray->SetName("AllChannels"); + + TH2D* hRead = 0; + TH2D* hAll = 0; + for(Int_t det =1; det<=3;det++) + { + TObjArray* detReadArray = new TObjArray(); + detReadArray->SetName(Form("FMD%d_Read",det)); + fReadArray->AddAtAndExpand(detReadArray,det); + + TObjArray* detArray = new TObjArray(); + detArray->SetName(Form("FMD%d_All",det)); + fAllArray->AddAtAndExpand(detArray,det); + + UShort_t nRings = (det==1 ? 1 : 2); + for(Int_t ring = 0;ringAddAtAndExpand(vtxArray,ring); + + TObjArray* vtxReadArray = new TObjArray(); + detReadArray->AddAtAndExpand(vtxReadArray,ring); + + + Char_t ringChar = (ring == 0 ? 'I' : 'O'); + Int_t nSec = (ringChar == 'I' ? 20 : 40); + + for(Int_t v=0; vSumw2(); + hAll->Sumw2(); + vtxArray->AddAtAndExpand(hAll,v); + vtxReadArray->AddAtAndExpand(hRead,v); + + } + + + } + } + + + + + AliFMDParameters* pars = AliFMDParameters::Instance(); + pars->Init(); + Int_t nDead = 0; + + for(UShort_t d=1;d<=1;d++) { + UShort_t nRings = (d==1 ? 1 : 2); + for (UShort_t ir = 0; ir < nRings; ir++) { + + Char_t ringChar = (ir == 0 ? 'I' : 'O'); + UShort_t nsec = (ir == 0 ? 20 : 40); + UShort_t nstr = (ir == 0 ? 512 : 256); + + std::cout<At(d); + TObjArray* vtxReadArray = (TObjArray*)detReadArray->At(ir); + TH2D* hRead = (TH2D*)vtxReadArray->At(v); + + TObjArray* detArray = (TObjArray*)fAllArray->At(d); + TObjArray* vtxArray = (TObjArray*)detArray->At(ir); + TH2D* hAll = (TH2D*)vtxArray->At(v); + + Float_t vtxMean = (v+0.5)*delta-vtxHigh; + for(UShort_t sec =0; sec < nsec; sec++) { + + for(UShort_t strip = 0; strip < nstr; strip++) { + + + if(pars->IsDead(d,ringChar,sec,strip)) { + nDead++; + } + else hRead->Fill(GetEtaFromStrip(d,ringChar,sec,strip,vtxMean),GetPhiFromSector(d,ringChar,sec)); + + hAll->Fill(GetEtaFromStrip(d,ringChar,sec,strip,vtxMean),GetPhiFromSector(d,ringChar,sec)); + + + } + } + } + } + } + + Float_t reldead = (Float_t)nDead / 51200.; + + std::cout<SetVertexAxis(nVtxBins, vtxLow, vtxHigh); + + for(UShort_t d=1;d<=3;d++) { + UShort_t nRings = (d==1 ? 1 : 2); + for (UShort_t ir = 0; ir < nRings; ir++) { + + Char_t ringChar = (ir == 0 ? 'I' : 'O'); + UShort_t nsec = (ir == 0 ? 20 : 40); + UShort_t nstr = (ir == 0 ? 512 : 256); + + std::cout<At(d); + TObjArray* vtxReadArray = (TObjArray*)detReadArray->At(ir); + TH2D* hRead = (TH2D*)vtxReadArray->At(vv); + + TObjArray* detArray = (TObjArray*)fAllArray->At(d); + TObjArray* vtxArray = (TObjArray*)detArray->At(ir); + TH2D* hAll = (TH2D*)vtxArray->At(vv); + + hRead->Divide(hRead,hAll,1,1,"B"); + UShort_t vtxbin = vv+1; + deadObject->SetCorrection(d,ringChar,vtxbin,hRead); + + } + } + } + +} -- 2.43.0