From: cholm Date: Wed, 31 Oct 2012 15:55:09 +0000 (+0000) Subject: Many changes in one. X-Git-Url: http://git.uio.no/git/?a=commitdiff_plain;h=5ca83feea0af7a210c2541f22521be25d1328abf;p=u%2Fmrichter%2FAliRoot.git Many changes in one. AOD generation: - Fixed how invalid signals are treated by the poisson calculator. Previously strips marked with kInvalidMult where considered empty. That is wrong because we don't know at all whether it was hit or not. Now they are filtered out of the calculation. - On-the-fly calculation of the phi acceptance per eta bin from the information stored in the ESDs. The phi acceptance is put in the overflow bin of the AOD object. - Possibly get the phi acceptance from the acceptance maps if they are used, Again, the phi acceptance is stored in the overflow bins. - Correct the phi angle of each strip from the (x,y) off-set of the IP. This means we have to propagate that information from the event inspector and on. - Possibility to flag extra strips as 'dead' in the sharing filter. This is to avoid using the acceptance maps for flagging dead-channels off-line * dN/deta generation - Possibly use the phi acceptance stored in the overflow bins - Fixed a few calculations of the normalisation. - Better diagnostics histograms --- diff --git a/PWGLF/FORWARD/analysis2/AliBaseMCTrackDensity.cxx b/PWGLF/FORWARD/analysis2/AliBaseMCTrackDensity.cxx index f693a802b25..81d6447d3dd 100644 --- a/PWGLF/FORWARD/analysis2/AliBaseMCTrackDensity.cxx +++ b/PWGLF/FORWARD/analysis2/AliBaseMCTrackDensity.cxx @@ -31,7 +31,7 @@ AliBaseMCTrackDensity::AliBaseMCTrackDensity() fDebug(false) { // Default constructor - DGUARD(0,0,"MC track density default construction"); + DGUARD(fDebug, 3,"Default CTOR of AliBasMCTrackDensity"); } //____________________________________________________________________ @@ -50,7 +50,7 @@ AliBaseMCTrackDensity::AliBaseMCTrackDensity(const char* name) fDebug(false) { // Normal constructor constructor - DGUARD(0,0,"MC track density named construction: %s", name); + DGUARD(fDebug, 3,"Named CTOR of AliBasMCTrackDensity: %s", name); } //____________________________________________________________________ @@ -69,7 +69,7 @@ AliBaseMCTrackDensity::AliBaseMCTrackDensity(const AliBaseMCTrackDensity& o) fDebug(o.fDebug) { // Normal constructor constructor - DGUARD(0,0,"MC track density copy construction"); + DGUARD(fDebug, 3,"Copy CTOR of AliBasMCTrackDensity"); } //____________________________________________________________________ diff --git a/PWGLF/FORWARD/analysis2/AliBasedNdetaTask.cxx b/PWGLF/FORWARD/analysis2/AliBasedNdetaTask.cxx index e98cf990bce..f511e7b7c92 100644 --- a/PWGLF/FORWARD/analysis2/AliBasedNdetaTask.cxx +++ b/PWGLF/FORWARD/analysis2/AliBasedNdetaTask.cxx @@ -43,7 +43,7 @@ AliBasedNdetaTask::AliBasedNdetaTask() // // Constructor // - DGUARD(0,0,"Default construction of AliBasedNdetaTask"); + DGUARD(fDebug,3,"Default CTOR of AliBasedNdetaTask"); } //____________________________________________________________________ @@ -75,7 +75,7 @@ AliBasedNdetaTask::AliBasedNdetaTask(const char* name) // // Constructor // - DGUARD(0,0,"Named construction of AliBasedNdetaTask: %s", name); + DGUARD(fDebug, 3,"Named CTOR of AliBasedNdetaTask: %s", name); fListOfCentralities = new TObjArray(1); // Set the normalisation scheme @@ -115,7 +115,7 @@ AliBasedNdetaTask::AliBasedNdetaTask(const AliBasedNdetaTask& o) fTriggerString(o.fTriggerString), fFinalMCCorrFile(o.fFinalMCCorrFile) { - DGUARD(0,0,"Copy construction of AliBasedNdetaTask"); + DGUARD(fDebug, 3,"Copy CTOR of AliBasedNdetaTask"); } //____________________________________________________________________ @@ -496,6 +496,25 @@ AliBasedNdetaTask::ScaleToCoverage(TH2D* copy, const TH1D* norm) } } } +//________________________________________________________________________ +void +AliBasedNdetaTask::ScaleToCoverage(TH1D* copy, const TH1D* norm) +{ + // Normalize to the acceptance - + // dndeta->Divide(accNorm); + for (Int_t i = 1; i <= copy->GetNbinsX(); i++) { + Double_t a = norm->GetBinContent(i); + if (a <= 0) { + copy->SetBinContent(i,0); + copy->SetBinError(i,0); + continue; + } + Double_t c = copy->GetBinContent(i); + Double_t e = copy->GetBinError(i); + copy->SetBinContent(i, c / a); + copy->SetBinError(i, e / a); + } +} //________________________________________________________________________ TH1D* @@ -535,10 +554,10 @@ AliBasedNdetaTask::ProjectX(const TH2D* h, Int_t first = firstbin; Int_t last = lastbin; - if (first < 0) first = 0; - else if (first >= yaxis->GetNbins()+1) first = yaxis->GetNbins(); + if (first < 0) first = 1; + else if (first >= yaxis->GetNbins()+2) first = yaxis->GetNbins()+1; if (last < 0) last = yaxis->GetNbins(); - else if (last > yaxis->GetNbins()+1) last = yaxis->GetNbins(); + else if (last >= yaxis->GetNbins()+2) last = yaxis->GetNbins()+1; if (last-first < 0) { AliWarningGeneral("AliBasedNdetaTask", Form("Nothing to project [%d,%d]", first, last)); @@ -1182,9 +1201,11 @@ AliBasedNdetaTask::Sum::CalcSum(TList* output, retCopy->SetMarkerStyle(marker); retCopy->SetDirectory(0); - TH1D* norm = ProjectX(fSum, "norm", 0, 0, rootProj, corrEmpty, false); - TH1D* norm0 = ProjectX(fSum0, "norm0", 0, 0, rootProj, corrEmpty, false); - TH1D* normAll = ProjectX(ret, "normAll", 0, 0, rootProj, corrEmpty, false); + Int_t nY = fSum->GetNbinsY(); + Int_t o = 0; // nY+1; + TH1D* norm = ProjectX(fSum, "norm", o, o, rootProj, corrEmpty, false); + TH1D* norm0 = ProjectX(fSum0, "norm0", o, o, rootProj, corrEmpty, false); + TH1D* normAll = ProjectX(ret, "normAll", o, o, rootProj, corrEmpty, false); norm->SetDirectory(0); norm0->SetDirectory(0); normAll->SetDirectory(0); @@ -1193,7 +1214,6 @@ AliBasedNdetaTask::Sum::CalcSum(TList* output, ScaleToCoverage(sum0Copy, norm0); ScaleToCoverage(retCopy, normAll); - Int_t nY = fSum->GetNbinsY(); TH1D* sumCopyPx = ProjectX(sumCopy, "average", 1, nY,rootProj,corrEmpty); TH1D* sum0CopyPx = ProjectX(sum0Copy, "average0", 1, nY,rootProj,corrEmpty); TH1D* retCopyPx = ProjectX(retCopy, "averageAll", 1, nY,rootProj,corrEmpty); @@ -1201,6 +1221,13 @@ AliBasedNdetaTask::Sum::CalcSum(TList* output, sum0CopyPx->SetDirectory(0); retCopyPx->SetDirectory(0); + TH1D* phi = ProjectX(fSum, "phi", nY+1, nY+1,rootProj,corrEmpty); + TH1D* phi0 = ProjectX(fSum0, "phi0", nY+1, nY+1,rootProj,corrEmpty); + TH1D* phiAll = ProjectX(ret, "phiAll", nY+1, nY+1,rootProj,corrEmpty); + phi->SetDirectory(0); + phi0->SetDirectory(0); + phiAll->SetDirectory(0); + // Scale our 1D histograms sumCopyPx->Scale(1., "width"); sum0CopyPx->Scale(1., "width"); @@ -1223,6 +1250,9 @@ AliBasedNdetaTask::Sum::CalcSum(TList* output, out->Add(norm); out->Add(norm0); out->Add(normAll); + out->Add(phi); + out->Add(phi0); + out->Add(phiAll); AliInfoF("Returning (1/%f * %s + 1/%f * %s), " "1/%f * %d + 1/%f * %d = %d", @@ -1255,7 +1285,7 @@ AliBasedNdetaTask::CentralityBin::CentralityBin() // // Constructor // - DGUARD(0,0,"Centrality bin default construction"); + DGUARD(fDebug,3,"Default CTOR of AliBasedNdeta::CentralityBin"); } //____________________________________________________________________ AliBasedNdetaTask::CentralityBin::CentralityBin(const char* name, @@ -1279,7 +1309,8 @@ AliBasedNdetaTask::CentralityBin::CentralityBin(const char* name, // low Lower centrality cut in percent // high Upper centrality cut in percent // - DGUARD(0,0,"Named Centrality bin construction: %s [%3d,%3d]",name,low,high); + DGUARD(fDebug,3,"Named CTOR of AliBasedNdeta::CentralityBin: %s [%3d,%3d]", + name,low,high); if (low <= 0 && high <= 0) { fLow = 0; fHigh = 0; @@ -1310,7 +1341,7 @@ AliBasedNdetaTask::CentralityBin::CentralityBin(const CentralityBin& o) // Parameters: // other Object to copy from // - DGUARD(0,0,"Copy Centrality bin construction"); + DGUARD(fDebug,3,"Copy CTOR of AliBasedNdeta::CentralityBin"); } //____________________________________________________________________ AliBasedNdetaTask::CentralityBin::~CentralityBin() @@ -1318,9 +1349,10 @@ AliBasedNdetaTask::CentralityBin::~CentralityBin() // // Destructor // - DGUARD(fDebug,3,"Centrality bin desctruction"); - if (fSums) fSums->Delete(); - if (fOutput) fOutput->Delete(); + DGUARD(fDebug,3,"DTOR of AliBasedNdeta::CentralityBin"); + + // if (fSums) fSums->Delete(); + // if (fOutput) fOutput->Delete(); } //____________________________________________________________________ @@ -1731,7 +1763,9 @@ AliBasedNdetaTask::CentralityBin::MakeResult(const TH2D* sum, DGUARD(fDebug,1,"Make centrality bin result from %s", sum->GetName()); TH2D* copy = static_cast(sum->Clone(Form("d2Ndetadphi%s%s", GetName(), postfix))); - TH1D* accNorm = ProjectX(sum, Form("norm%s%s",GetName(), postfix), 0, 0, + Int_t nY = sum->GetNbinsY(); + Int_t o = (corrEmpty ? 0 : nY+1); + TH1D* accNorm = ProjectX(sum, Form("norm%s%s",GetName(), postfix), o, o, rootProj, corrEmpty, false); accNorm->SetDirectory(0); @@ -1740,16 +1774,21 @@ AliBasedNdetaTask::CentralityBin::MakeResult(const TH2D* sum, else AliInfo("No shape correction specified, or disabled"); // --- Normalize to the coverage ----------------------------------- - ScaleToCoverage(copy, accNorm); - - // --- Event-level normalization ----------------------------------- - copy->Scale(scaler); + if (corrEmpty) { + ScaleToCoverage(copy, accNorm); + // --- Event-level normalization --------------------------------- + copy->Scale(scaler); + } // --- Project on X axis ------------------------------------------- TH1D* dndeta = ProjectX(copy, Form("dndeta%s%s",GetName(), postfix), - 1, sum->GetNbinsY(), rootProj, corrEmpty); + 1, nY, rootProj, corrEmpty); dndeta->SetDirectory(0); // Event-level normalization + if (!corrEmpty) { + ScaleToCoverage(dndeta, accNorm); + dndeta->Scale(scaler); + } dndeta->Scale(1., "width"); copy->Scale(1., "width"); diff --git a/PWGLF/FORWARD/analysis2/AliBasedNdetaTask.h b/PWGLF/FORWARD/analysis2/AliBasedNdetaTask.h index 28bfcf6ab8d..8ca3801f94c 100644 --- a/PWGLF/FORWARD/analysis2/AliBasedNdetaTask.h +++ b/PWGLF/FORWARD/analysis2/AliBasedNdetaTask.h @@ -314,6 +314,13 @@ public: * @param norm Coverage histogram */ static void ScaleToCoverage(TH2D* copy, const TH1D* norm); + /** + * Scale the copy of the 1D histogram by coverage in supplied 1D histogram + * + * @param copy Data to scale + * @param norm Coverage histogram + */ + static void ScaleToCoverage(TH1D* copy, const TH1D* norm); /** * Set histogram graphical options, etc. * diff --git a/PWGLF/FORWARD/analysis2/AliCentralMCCorrectionsTask.cxx b/PWGLF/FORWARD/analysis2/AliCentralMCCorrectionsTask.cxx index c1457ca189c..5cdb512c3d5 100644 --- a/PWGLF/FORWARD/analysis2/AliCentralMCCorrectionsTask.cxx +++ b/PWGLF/FORWARD/analysis2/AliCentralMCCorrectionsTask.cxx @@ -65,7 +65,7 @@ AliCentralMCCorrectionsTask::AliCentralMCCorrectionsTask() // Parameters: // name Name of task // - DGUARD(fDebug,0,"Default construction of AliCentralMCCorrectionsTask"); + DGUARD(fDebug, 3,"Default CTOR of AliCentralMCCorrectionsTask"); } //____________________________________________________________________ @@ -90,7 +90,7 @@ AliCentralMCCorrectionsTask::AliCentralMCCorrectionsTask(const char* name) // Parameters: // name Name of task // - DGUARD(fDebug,0,"Named construction of AliCentralMCCorrectionsTask: %s",name); + DGUARD(fDebug, 3,"Named CTOR of AliCentralMCCorrectionsTask: %s",name); DefineOutput(1, TList::Class()); DefineOutput(2, TList::Class()); fBranchNames = @@ -120,7 +120,7 @@ AliCentralMCCorrectionsTask::AliCentralMCCorrectionsTask(const AliCentralMCCorre // Parameters: // o Object to copy from // - DGUARD(fDebug,0,"Copy construction of AliCentralMCCorrectionsTask"); + DGUARD(fDebug, 3,"Copy CTOR of AliCentralMCCorrectionsTask"); } //____________________________________________________________________ @@ -381,7 +381,7 @@ AliCentralMCCorrectionsTask::UserExec(Option_t*) UInt_t triggers = 0; // Trigger bits Bool_t lowFlux = true; // Low flux flag UShort_t iVz = 0; // Vertex bin from ESD - Double_t vZ = 1000; // Z coordinate from ESD + TVector3 ip; // Z coordinate from ESD Double_t cent = -1; // Centrality UShort_t iVzMc = 0; // Vertex bin from MC Double_t vZMc = 1000; // Z coordinate of IP vertex from MC @@ -392,7 +392,7 @@ AliCentralMCCorrectionsTask::UserExec(Option_t*) Double_t phiR = 1000; // Reaction plane from MC UShort_t nClusters = 0; // Number of SPD clusters // Process the data - UInt_t retESD = fInspector.Process(esd, triggers, lowFlux, iVz, vZ, + UInt_t retESD = fInspector.Process(esd, triggers, lowFlux, iVz, ip, cent, nClusters); fInspector.ProcessMC(mcEvent, triggers, iVzMc, vZMc, b, cMC, nPart, nBin, phiR); diff --git a/PWGLF/FORWARD/analysis2/AliCentralMultiplicityTask.cxx b/PWGLF/FORWARD/analysis2/AliCentralMultiplicityTask.cxx index f04085bb843..831c9f7cb2e 100644 --- a/PWGLF/FORWARD/analysis2/AliCentralMultiplicityTask.cxx +++ b/PWGLF/FORWARD/analysis2/AliCentralMultiplicityTask.cxx @@ -51,7 +51,7 @@ AliCentralMultiplicityTask::AliCentralMultiplicityTask(const char* name) // // Constructor // - DGUARD(fDebug,0,"Named CTOR of AliCentralMultiplicityTask: %s", name); + DGUARD(fDebug, 3,"Named CTOR of AliCentralMultiplicityTask: %s", name); DefineOutput(1, TList::Class()); fBranchNames = "ESD:AliESDRun.,AliESDHeader.,AliMultiplicity.," @@ -80,7 +80,7 @@ AliCentralMultiplicityTask::AliCentralMultiplicityTask() // // Constructor // - DGUARD(fDebug,0,"Default CTOR of AliCentralMultiplicityTask"); + DGUARD(fDebug, 3,"Default CTOR of AliCentralMultiplicityTask"); } //____________________________________________________________________ AliCentralMultiplicityTask::AliCentralMultiplicityTask(const AliCentralMultiplicityTask& o) @@ -105,7 +105,7 @@ AliCentralMultiplicityTask::AliCentralMultiplicityTask(const AliCentralMultiplic // // Copy constructor // - DGUARD(fDebug,0,"COPY CTOR of AliCentralMultiplicityTask"); + DGUARD(fDebug, 3,"COPY CTOR of AliCentralMultiplicityTask"); } //____________________________________________________________________ @@ -419,11 +419,11 @@ void AliCentralMultiplicityTask::UserExec(Option_t* /*option*/) Bool_t lowFlux = kFALSE; UInt_t triggers = 0; UShort_t ivz = 0; - Double_t vz = 0; + TVector3 ip; Double_t cent = -1; UShort_t nClusters = 0; UInt_t found = fInspector.Process(esd, triggers, lowFlux, - ivz, vz, cent, nClusters); + ivz, ip, cent, nClusters); // No event or no trigger if (found & AliFMDEventInspector::kNoEvent) return; @@ -497,36 +497,37 @@ AliCentralMultiplicityTask::CorrectData(TH2D& aodHist, UShort_t vtxbin) const if (!hAcceptance) AliFatal("No acceptance!"); if (fUseSecondary && hSecMap) aodHist.Divide(hSecMap); - - for(Int_t nx = 1; nx <= aodHist.GetNbinsX(); nx++) { - Float_t accCor = hAcceptance->GetBinContent(nx); - Float_t accErr = hAcceptance->GetBinError(nx); + + Int_t nY = aodHist.GetNbinsY(); + for(Int_t ix = 1; ix <= aodHist.GetNbinsX(); ix++) { + Float_t accCor = hAcceptance->GetBinContent(ix); + Float_t accErr = hAcceptance->GetBinError(ix); Bool_t fiducial = true; - if (nx < fEtaMin[vtxbin-1] || nx > fEtaMax[vtxbin-1]) + if (ix < fEtaMin[vtxbin-1] || ix > fEtaMax[vtxbin-1]) fiducial = false; // Bool_t etabinSeen = kFALSE; - for(Int_t ny = 1; ny <= aodHist.GetNbinsY(); ny++) { + for(Int_t iy = 1; iy <= nY; iy++) { #if 1 if (!fiducial) { - aodHist.SetBinContent(nx, ny, 0); - aodHist.SetBinError(nx, ny, 0); + aodHist.SetBinContent(ix, iy, 0); + aodHist.SetBinError(ix, iy, 0); continue; } #endif // Get currrent value - Float_t aodValue = aodHist.GetBinContent(nx,ny); - Float_t aodErr = aodHist.GetBinError(nx,ny); + Float_t aodValue = aodHist.GetBinContent(ix,iy); + Float_t aodErr = aodHist.GetBinError(ix,iy); #if 0 // This is done once in the FindEtaBins function // Set underflow bin Float_t secCor = 0; - if(hSecMap) secCor = hSecMap->GetBinContent(nx,ny); + if(hSecMap) secCor = hSecMap->GetBinContent(ix,iy); if (secCor > 0.5) etabinSeen = kTRUE; #endif if (aodValue < 0.000001) { - aodHist.SetBinContent(nx,ny, 0); - aodHist.SetBinError(nx,ny, 0); + aodHist.SetBinContent(ix,iy, 0); + aodHist.SetBinError(ix,iy, 0); continue; } if (!fUseAcceptance) continue; @@ -536,14 +537,17 @@ AliCentralMultiplicityTask::CorrectData(TH2D& aodHist, UShort_t vtxbin) const Float_t aodNew = aodValue / accCor ; Float_t error = aodNew*TMath::Sqrt(TMath::Power(aodErr/aodValue,2) + TMath::Power(accErr/accCor,2) ); - aodHist.SetBinContent(nx,ny, aodNew); + aodHist.SetBinContent(ix,iy, aodNew); //test - aodHist.SetBinError(nx,ny,error); - aodHist.SetBinError(nx,ny,aodErr); + aodHist.SetBinError(ix,iy,error); + aodHist.SetBinError(ix,iy,aodErr); } //Filling underflow bin if we eta bin is in range - if (fiducial) aodHist.SetBinContent(nx,0, 1.); - // if (etabinSeen) aodHist.SetBinContent(nx,0, 1.); + if (fiducial) { + aodHist.SetBinContent(ix,0, 1.); + aodHist.SetBinContent(ix,nY+1, 1.); + } + // if (etabinSeen) aodHist.SetBinContent(ix,0, 1.); } } @@ -587,9 +591,13 @@ AliCentralMultiplicityTask::Print(Option_t* option) const for (Int_t v = 1; v <= vaxis.GetNbins(); v++) { std::cout << " " << std::setw(2) << v << " " << std::setw(5) << vaxis.GetBinLowEdge(v) << "-" - << std::setw(5) << vaxis.GetBinUpEdge(v) << " | " - << std::setw(3) << fEtaMin[v-1] << "-" - << std::setw(3) << fEtaMax[v-1] << std::endl; + << std::setw(5) << vaxis.GetBinUpEdge(v) << " | "; + if (fEtaMin.GetSize() <= 0) + std::cout << " ? - ?"; + else + std::cout << std::setw(3) << fEtaMin[v-1] << "-" + << std::setw(3) << fEtaMax[v-1]; + std::cout << std::endl; } } @@ -836,9 +844,10 @@ AliCentralMultiplicityTask::Manager::WriteFile(UShort_t what, Form("Failed to write %s to disk (%d)", ofName.Data(),ret)); return false; } - output->ls(); + // output->ls(); output->Close(); +#if 0 TString cName(obj->IsA()->GetName()); AliInfoGeneral("Manager", Form("Wrote %s object %s to %s\n", @@ -846,7 +855,7 @@ AliCentralMultiplicityTask::Manager::WriteFile(UShort_t what, if (!full) { TString dName(GetFileDir(what)); AliInfoGeneral("Manager", - Form("%s should be copied to %s\n" + Form("\n %s should be copied to %s\n" "Do for example\n\t" "aliroot $ALICE_ROOT/PWGLF/FORWARD/analysis2/scripts/" "MoveCorrections.C\\(%d\\)\nor\n\t" @@ -857,6 +866,7 @@ AliCentralMultiplicityTask::Manager::WriteFile(UShort_t what, } +#endif return true; } diff --git a/PWGLF/FORWARD/analysis2/AliCentraldNdetaTask.cxx b/PWGLF/FORWARD/analysis2/AliCentraldNdetaTask.cxx index 0c2f980097c..27df30e42ea 100644 --- a/PWGLF/FORWARD/analysis2/AliCentraldNdetaTask.cxx +++ b/PWGLF/FORWARD/analysis2/AliCentraldNdetaTask.cxx @@ -22,14 +22,14 @@ ClassImp(AliCentraldNdetaTask) AliCentraldNdetaTask::AliCentraldNdetaTask() : AliBasedNdetaTask() { - DGUARD(fDebug,0,"Default construction of AliCentraldNdetaTask"); + DGUARD(fDebug,3,"Default CTOR of AliCentraldNdetaTask"); } //____________________________________________________________________ AliCentraldNdetaTask::AliCentraldNdetaTask(const char*) : AliBasedNdetaTask("Central") { - DGUARD(fDebug,0,"Named construction of AliCentraldNdetaTask"); + DGUARD(fDebug,3,"Named CTOR of AliCentraldNdetaTask"); fSymmetrice = false; fCorrEmpty = false; SetTitle("Central"); diff --git a/PWGLF/FORWARD/analysis2/AliFMDCorrAcceptance.h b/PWGLF/FORWARD/analysis2/AliFMDCorrAcceptance.h index 1f3a0e6a1ff..411ed635c9f 100644 --- a/PWGLF/FORWARD/analysis2/AliFMDCorrAcceptance.h +++ b/PWGLF/FORWARD/analysis2/AliFMDCorrAcceptance.h @@ -18,6 +18,7 @@ #include #include class TH2D; +class TH1D; /** * This class contains the acceptance correction due to dead channels @@ -77,12 +78,37 @@ public: * @return The correction @f$ a_{r,v}@f$ */ TH2D* GetCorrection(UShort_t d, Char_t r, UShort_t b) const; + /** + * Get the phi acceptance + * + * @param d + * @param r + * @param v + * + * @return + */ + TH1D* GetPhiAcceptance(UShort_t d, Char_t r, Double_t v) const; + /** + * Get the phi acceptance + * + * @param d + * @param r + * @param v + * + * @return + */ + TH1D* GetPhiAcceptance(UShort_t d, Char_t r, UShort_t b) const; /** * Get the vertex axis used * * @return vertex axis */ const TAxis& GetVertexAxis() const { return fVertexAxis; } + /** + * @return true if the overflow bins along eta are set to the ratio + * of OK strips to all strips for a given eta bin. + */ + Bool_t HasOverflow() const { return fHasOverflow; } /* @} */ /** @@ -128,6 +154,12 @@ public: * @param max Maximum */ void SetVertexAxis(Int_t nBins, Double_t min, Double_t max); + /** + * Set that we have (or not) the overflow bin present + * + * @param present If true, assume the overflow bins are filled + */ + void SetHasOverflow(Bool_t present=true) { fHasOverflow = present; } /* @} */ /** @@ -180,7 +212,7 @@ protected: * * @return Pointer to ring array, or null in case of problems */ - TObjArray* GetRingArray(UShort_t d, Char_t r) const; + TObjArray* GetRingArray(const TObjArray& m, UShort_t d, Char_t r) const; /** * Get the ring array corresponding to the specified ring * @@ -189,11 +221,17 @@ protected: * * @return Pointer to ring array, or newly created container */ - TObjArray* GetOrMakeRingArray(UShort_t d, Char_t r); + TObjArray* GetOrMakeRingArray(TObjArray& m, UShort_t d, Char_t r) const; + TObject* GetObject(const TObjArray& m, UShort_t d, + Char_t r, UShort_t b) const; + + void FillCache() const; - TObjArray fRingArray; // Array of per-ring, per-vertex 2nd map + TObjArray fRingArray; // Array of per-ring, per-vertex 2d map + mutable TObjArray* fCache; //! Array of per-ring, per-vertex 1d factors TAxis fVertexAxis; // The vertex axis - ClassDef(AliFMDCorrAcceptance,1); // Acceptance correction due to dead areas + Bool_t fHasOverflow; // Whether we have the overflow bin set + ClassDef(AliFMDCorrAcceptance,2); // Acceptance correction due to dead areas }; //____________________________________________________________________ diff --git a/PWGLF/FORWARD/analysis2/AliFMDCorrector.cxx b/PWGLF/FORWARD/analysis2/AliFMDCorrector.cxx index 1a3dbc6e889..28e335c3582 100644 --- a/PWGLF/FORWARD/analysis2/AliFMDCorrector.cxx +++ b/PWGLF/FORWARD/analysis2/AliFMDCorrector.cxx @@ -35,7 +35,7 @@ AliFMDCorrector::AliFMDCorrector() fDebug(0) { // Constructor - DGUARD(fDebug, 0, "Default CTOR of AliFMDCorrector"); + DGUARD(fDebug, 3, "Default CTOR of AliFMDCorrector"); } //____________________________________________________________________ @@ -52,7 +52,7 @@ AliFMDCorrector::AliFMDCorrector(const char* title) // // Parameters: // title Title - DGUARD(fDebug, 0, "Named CTOR of AliFMDCorrector: %s", title); + DGUARD(fDebug, 3, "Named CTOR of AliFMDCorrector: %s", title); fRingHistos.SetName(GetName()); fRingHistos.Add(new RingHistos(1, 'I')); fRingHistos.Add(new RingHistos(2, 'I')); @@ -75,7 +75,7 @@ AliFMDCorrector::AliFMDCorrector(const AliFMDCorrector& o) // // Parameters: // o Object to copy from - DGUARD(fDebug, 0, "Copy CTOR of AliFMDCorrector"); + DGUARD(fDebug, 3, "Copy CTOR of AliFMDCorrector"); TIter next(&o.fRingHistos); TObject* obj = 0; while ((obj = next())) fRingHistos.Add(obj); @@ -159,6 +159,43 @@ AliFMDCorrector::GetRingHistos(UShort_t d, Char_t r) const return static_cast(fRingHistos.At(idx)); } +//____________________________________________________________________ +void +AliFMDCorrector::DivideMap(TH2* num, const TH2* denom, + Bool_t alsoUnderOver) const +{ + // + // Implement TH1::Divide but + // - Assume compatible histograms + // - Unless third argument is true, do not divide over/under flow bins + // + if (!num || !denom) return; + + Int_t first = (alsoUnderOver ? 0 : 1); + Int_t lastX = num->GetNbinsX() + (alsoUnderOver ? 1 : 0); + Int_t lastY = num->GetNbinsY() + (alsoUnderOver ? 1 : 0); + + for (Int_t ix = first; ix <= lastX; ix++) { + for (Int_t iy = first; iy <= lastY; iy++) { + Int_t bin = num->GetBin(ix,iy); + Double_t c0 = num->GetBinContent(bin); + Double_t c1 = denom->GetBinContent(bin); + if (!c1) { + num->SetBinContent(bin,0); + num->SetBinError(bin, 0); + continue; + } + Double_t w = c0 / c1; + Double_t e0 = num->GetBinError(bin); + Double_t e1 = denom->GetBinError(bin); + Double_t c12 = c1*c1; + Double_t e2 = (e0*e0*c1*c1 + e1*e1*c0*c0)/(c12*c12); + + num->SetBinContent(bin, w); + num->SetBinError(bin, TMath::Sqrt(e2)); + } + } +} //____________________________________________________________________ Bool_t AliFMDCorrector::Correct(AliForwardUtil::Histos& hists, @@ -192,7 +229,7 @@ AliFMDCorrector::Correct(AliForwardUtil::Histos& hists, continue; } // Divide by primary/total ratio - h->Divide(bg); + DivideMap(h, bg, false); } if (fUseVertexBias) { TH2D* ef = fcm.GetVertexBias()->GetCorrection(r, uvb); @@ -202,7 +239,7 @@ AliFMDCorrector::Correct(AliForwardUtil::Histos& hists, continue; } // Divide by the event selection efficiency - h->Divide(ef); + DivideMap(h, ef, false); } if (fUseAcceptance) { TH2D* ac = fcm.GetAcceptance()->GetCorrection(d, r, uvb); @@ -211,8 +248,12 @@ AliFMDCorrector::Correct(AliForwardUtil::Histos& hists, "vertex bin %d", d, r, uvb)); continue; } - // Divide by the acceptance correction - h->Divide(ac); + // Fill overflow bin with ones + for (Int_t i = 1; i <= h->GetNbinsX(); i++) + h->SetBinContent(i, h->GetNbinsY()+1, 1); + + // Divide by the acceptance correction - + DivideMap(h, ac, fcm.GetAcceptance()->HasOverflow()); } if (fUseMergingEfficiency) { @@ -231,10 +272,14 @@ AliFMDCorrector::Correct(AliForwardUtil::Histos& hists, for (Int_t ieta = 1; ieta <= h->GetNbinsX(); ieta++) { Float_t c = sf->GetBinContent(ieta); Float_t ec = sf->GetBinError(ieta); - - if (c == 0) continue; - + for (Int_t iphi = 1; iphi <= h->GetNbinsY(); iphi++) { + if (c == 0) { + h->SetBinContent(ieta,iphi,0); + h->SetBinError(ieta,iphi,0); + continue; + } + Double_t m = h->GetBinContent(ieta, iphi) / c; Double_t em = h->GetBinError(ieta, iphi); @@ -245,49 +290,6 @@ AliFMDCorrector::Correct(AliForwardUtil::Histos& hists, } } } - //HHD - /* - TH2D* bg = fcm.GetSecondaryMap()->GetCorrection(d, r, uvb); - TH2D hRing("hring","hring",bg->GetNbinsX(), - bg->GetXaxis()->GetXmin(), - bg->GetXaxis()->GetXmax(), - bg->GetNbinsY(), - bg->GetYaxis()->GetXmin(), - bg->GetYaxis()->GetXmax()); - - Int_t edgebin[4] = {0,0,0,0}; - for(Int_t ii = 1; ii <=bg->GetNbinsX(); ii++) { - for(Int_t jj = 1; jj <=bg->GetNbinsY(); jj++) { - Float_t bgcor = bg->GetBinContent(ii,jj); - if(bgcor<0.1) continue; - if(edgebin[0] == 0) edgebin[0] = ii; - if(edgebin[0] == ii) continue; - if(edgebin[0] > 0 && edgebin[1] == 0) edgebin[1] = ii; - if(edgebin[0]>0 && edgebin[1]>0) break; - } - } - for(Int_t ii = bg->GetNbinsX(); ii >= 1; ii--) { - for(Int_t jj = 1; jj <=bg->GetNbinsY(); jj++) { - Float_t bgcor = bg->GetBinContent(ii,jj); - if(bgcor<0.1) continue; - if(edgebin[2] == 0) edgebin[2] = ii; - if(edgebin[2] == ii) continue; - if(edgebin[2] > 0 && edgebin[3] == 0) edgebin[3] = ii; - if(edgebin[2]>0 && edgebin[3]>0) break; - } - } - for(Int_t ii = 1; ii <=bg->GetNbinsX(); ii++) { - for(Int_t jj = 1; jj <=bg->GetNbinsY(); jj++) { - Float_t data = h->GetBinContent(ii,jj); - if(data <0.000001) continue; - if(edgebin[0] == ii || edgebin[1] == ii || edgebin[2] == ii || edgebin[3] == ii) continue; - hRing.SetBinContent(ii,jj,data); - hRing.SetBinError(ii,jj,h->GetBinError(ii,jj)); - } - } - - //std::cout<fDensity->Add(h); } diff --git a/PWGLF/FORWARD/analysis2/AliFMDCorrector.h b/PWGLF/FORWARD/analysis2/AliFMDCorrector.h index 25f06477f6b..adedf391cac 100644 --- a/PWGLF/FORWARD/analysis2/AliFMDCorrector.h +++ b/PWGLF/FORWARD/analysis2/AliFMDCorrector.h @@ -17,7 +17,7 @@ #include #include "AliForwardUtil.h" class TH2D; - +class TH2; /** * @defgroup pwglf_forward_algo Algorithms * @@ -230,7 +230,16 @@ protected: * @return Ring histogram container */ RingHistos* GetRingHistos(UShort_t d, Char_t r) const; - + /** + * Divide a map with another map. This is a reimplementation of + * TH1::Divide, but we assume compatible histograms, and the under- + * and overflow bins are only divided if the third argument is true. + * + * @param num Numerator. On return contains the result + * @param denom Denominator + * @param alsoUnderOver If true, also divide under/overflow bins + */ + void DivideMap(TH2* num, const TH2* denom, Bool_t alsoUnderOver=false) const; TList fRingHistos; // List of histogram containers Bool_t fUseSecondaryMap; // Whether to do correction for secondaries Bool_t fUseVertexBias; // Whether to do correction for vertex bias @@ -238,7 +247,7 @@ protected: Bool_t fUseMergingEfficiency; // Whether to use the merging efficiency Int_t fDebug; // Debug level - ClassDef(AliFMDCorrector,2); // Correct the inclusive d2N/detadphi + ClassDef(AliFMDCorrector,3); // Correct the inclusive d2N/detadphi }; #endif diff --git a/PWGLF/FORWARD/analysis2/AliFMDDensityCalculator.cxx b/PWGLF/FORWARD/analysis2/AliFMDDensityCalculator.cxx index 919c36c5748..f1438ca6a7b 100644 --- a/PWGLF/FORWARD/analysis2/AliFMDDensityCalculator.cxx +++ b/PWGLF/FORWARD/analysis2/AliFMDDensityCalculator.cxx @@ -14,6 +14,7 @@ #include #include #include +#include #include #include @@ -22,6 +23,9 @@ ClassImp(AliFMDDensityCalculator) ; // For Emacs #endif +//____________________________________________________________________ +const char* AliFMDDensityCalculator::fgkFolderName = "fmdDensityCalculator"; + //____________________________________________________________________ AliFMDDensityCalculator::AliFMDDensityCalculator() : TNamed(), @@ -45,17 +49,18 @@ AliFMDDensityCalculator::AliFMDDensityCalculator() fPhiLumping(4), fDebug(0), fCuts(), - fRecalculateEta(false) + fRecalculateEta(false), + fRecalculatePhi(false) { // // Constructor // - DGUARD(fDebug, 0, "Default CTOR of FMD density calculator"); + DGUARD(fDebug, 3, "Default CTOR of FMD density calculator"); } //____________________________________________________________________ AliFMDDensityCalculator::AliFMDDensityCalculator(const char* title) - : TNamed("fmdDensityCalculator", title), + : TNamed(fgkFolderName, title), fRingHistos(), fSumOfWeights(0), fWeightedSum(0), @@ -76,7 +81,8 @@ AliFMDDensityCalculator::AliFMDDensityCalculator(const char* title) fPhiLumping(4), fDebug(0), fCuts(), - fRecalculateEta(false) + fRecalculateEta(false), + fRecalculatePhi(false) { // // Constructor @@ -84,7 +90,7 @@ AliFMDDensityCalculator::AliFMDDensityCalculator(const char* title) // Parameters: // name Name of object // - DGUARD(fDebug, 0, "Named CTOR of FMD density calculator: %s", title); + DGUARD(fDebug, 3, "Named CTOR of FMD density calculator: %s", title); fRingHistos.SetName(GetName()); fRingHistos.SetOwner(); fRingHistos.Add(new RingHistos(1, 'I')); @@ -143,7 +149,8 @@ AliFMDDensityCalculator::AliFMDDensityCalculator(const fPhiLumping(o.fPhiLumping), fDebug(o.fDebug), fCuts(o.fCuts), - fRecalculateEta(o.fRecalculateEta) + fRecalculateEta(o.fRecalculateEta), + fRecalculatePhi(o.fRecalculatePhi) { // // Copy constructor @@ -151,7 +158,7 @@ AliFMDDensityCalculator::AliFMDDensityCalculator(const // Parameters: // o Object to copy from // - DGUARD(fDebug, 0, "Copy CTOR of FMD density calculator"); + DGUARD(fDebug, 3, "Copy CTOR of FMD density calculator"); TIter next(&o.fRingHistos); TObject* obj = 0; while ((obj = next())) fRingHistos.Add(obj); @@ -201,6 +208,7 @@ AliFMDDensityCalculator::operator=(const AliFMDDensityCalculator& o) fPhiLumping = o.fPhiLumping; fCuts = o.fCuts; fRecalculateEta = o.fRecalculateEta; + fRecalculatePhi = o.fRecalculatePhi; fRingHistos.Delete(); TIter next(&o.fRingHistos); @@ -217,7 +225,7 @@ AliFMDDensityCalculator::Init(const TAxis& axis) // Intialize this sub-algorithm // // Parameters: - // etaAxis Not used + // etaAxis Eta axis DGUARD(fDebug, 1, "Initialize FMD density calculator"); CacheMaxWeights(axis); @@ -294,10 +302,9 @@ AliFMDDensityCalculator::GetMultCut(UShort_t d, Char_t r, Double_t eta, Bool_t AliFMDDensityCalculator::Calculate(const AliESDFMD& fmd, AliForwardUtil::Histos& hists, - UShort_t vtxbin, Bool_t lowFlux, Double_t /*cent*/, - Double_t zvtx) + const TVector3& ip) { // // Do the calculations @@ -311,7 +318,8 @@ AliFMDDensityCalculator::Calculate(const AliESDFMD& fmd, // Return: // true on successs DGUARD(fDebug, 1, "Calculate density in FMD density calculator"); - + + // --- Loop over detectors ----------------------------------------- for (UShort_t d=1; d<=3; d++) { UShort_t nr = (d == 1 ? 1 : 2); for (UShort_t q=0; qfPoisson.SetObject(d,r,vtxbin,cent); rh->fPoisson.Reset(0); + rh->fTotal->Reset(); + rh->fGood->Reset(); // rh->ResetPoissonHistos(h, fEtaLumping, fPhiLumping); + // --- Loop over sectors and strips ---------------------------- for (UShort_t s=0; sfTotal->Fill(eta); if (mult == AliESDFMD::kInvalidMult || mult > 20) { - rh->fPoisson.Fill(t , s, false); - rh->fEvsM->Fill(mult,0); + // Do not count invalid stuff + // rh->fPoisson.Fill(t , s, false); + rh->fEvsN->Fill(mult,-1); + rh->fEvsM->Fill(mult,-1); continue; } + // --- Automatic calculation of acceptance ----------------- + rh->fGood->Fill(eta); + + // --- If we asked to re-calculate phi for (x,y) IP -------- + if (fRecalculatePhi) { + oldPhi = phi; + phi = AliForwardUtil::GetPhiFromStrip(r, t, phi, ip.X(), ip.Y()); + } + + // --- Apply phi corner correction to eloss ---------------- if (fUsePhiAcceptance == kPhiCorrectELoss) mult *= AcceptanceCorrection(r,t); + // --- Get the low multiplicity cut ------------------------ Double_t cut = 1024; if (eta != AliESDFMD::kInvalidEta) cut = GetMultCut(d, r, eta,false); + // --- Now caluculate Nch for this strip using fits -------- Double_t n = 0; - if (cut > 0 && mult > cut) - n = NParticles(mult,d,r,s,t,vtxbin,eta,lowFlux); - + if (cut > 0 && mult > cut) n = NParticles(mult,d,r,eta,lowFlux); rh->fELoss->Fill(mult); rh->fEvsN->Fill(mult,n); rh->fEtaVsN->Fill(eta, n); - Double_t c = Correction(d,r,s,t,vtxbin,eta,lowFlux); + // --- Calculate correction if needed ---------------------- + Double_t c = Correction(d,r,t,eta,lowFlux); fCorrections->Fill(c); if (c > 0) n /= c; rh->fEvsM->Fill(mult,n); rh->fEtaVsM->Fill(eta, n); rh->fCorr->Fill(eta, c); + // --- Accumulate Poisson statistics ----------------------- Bool_t hit = (n > 0.9 && c > 0); - if (hit) rh->fELossUsed->Fill(mult); + if (hit) { + rh->fELossUsed->Fill(mult); + if (fRecalculatePhi) { + rh->fPhiBefore->Fill(oldPhi); + rh->fPhiAfter->Fill(phi); + } + } rh->fPoisson.Fill(t,s,hit,1./c); h->Fill(eta,phi,n); + + // --- If we use ELoss fits, apply now --------------------- if (!fUsePoisson) rh->fDensity->Fill(eta,phi,n); } // for t } // for s - + + // --- Automatic acceptance - Calculate as an efficiency ------- + rh->fGood->Divide(rh->fGood, rh->fTotal, 1, 1, "B"); + + // --- Make a copy and reset as needed ------------------------- TH2D* hclone = static_cast(h->Clone("hclone")); if (!fUsePoisson) hclone->Reset(); if ( fUsePoisson) h->Reset(); + // --- Store Poisson result ------------------------------------ TH2D* poisson = rh->fPoisson.Result(); for (Int_t t=0; t < poisson->GetNbinsX(); t++) { for (Int_t s=0; s < poisson->GetNbinsY(); s++) { @@ -386,7 +427,7 @@ AliFMDDensityCalculator::Calculate(const AliESDFMD& fmd, Double_t phi = fmd.Phi(d,r,s,t) / 180 * TMath::Pi(); Double_t eta = fmd.Eta(d,r,s,t); if(fRecalculateEta) - eta = AliForwardUtil::GetEtaFromStrip(d,r,s,t,zvtx); + eta = AliForwardUtil::GetEtaFromStrip(d,r,s,t,ip.Z()); if (fUsePoisson) h->Fill(eta,phi,poissonV); else @@ -395,8 +436,23 @@ AliFMDDensityCalculator::Calculate(const AliESDFMD& fmd, } } + // --- Make diagnostics - eloss vs poisson --------------------- + Int_t nY = h->GetNbinsY(); for (Int_t ieta=1; ieta <= h->GetNbinsX(); ieta++) { - for (Int_t iphi=1; iphi<= h->GetNbinsY(); iphi++) { + // Set the overflow bin to contain the phi acceptance + Double_t phiAcc = rh->fGood->GetBinContent(ieta); + Double_t phiAccE = rh->fGood->GetBinError(ieta); + h->SetBinContent(ieta, nY+1, phiAcc); + h->SetBinError(ieta, nY+1, phiAccE); + Double_t eta = h->GetXaxis()->GetBinCenter(ieta); + rh->fPhiAcc->Fill(eta, ip.Z(), phiAcc); +#if 0 + if (phiAcc > 1e-12) { + Info("", "FMD%d%c, eta=%3d phi acceptance: %f (%f)", + d, r, ieta, phiAcc, h->GetBinContent(ieta, nY+1)); + } +#endif + for (Int_t iphi=1; iphi<= nY; iphi++) { Double_t poissonV = 0; //h->GetBinContent(,s+1); Double_t eLossV = 0; @@ -410,6 +466,9 @@ AliFMDDensityCalculator::Calculate(const AliESDFMD& fmd, } rh->fELossVsPoisson->Fill(eLossV, poissonV); + rh->fDiffELossPoisson->Fill(poissonV < 1e-12 ? 0 : + (eLossV - poissonV) / poissonV); + } } delete hclone; @@ -575,9 +634,6 @@ Float_t AliFMDDensityCalculator::NParticles(Float_t mult, UShort_t d, Char_t r, - UShort_t /*s*/, - UShort_t /*t*/, - UShort_t /*v*/, Float_t eta, Bool_t lowFlux) const { @@ -631,9 +687,7 @@ AliFMDDensityCalculator::NParticles(Float_t mult, Float_t AliFMDDensityCalculator::Correction(UShort_t d, Char_t r, - UShort_t /*s*/, UShort_t t, - UShort_t /*v*/, Float_t eta, Bool_t lowFlux) const { @@ -832,33 +886,25 @@ AliFMDDensityCalculator::DefineOutput(TList* dir) // TNamed* sigma = new TNamed("sigma", // (fIncludeSigma ? "included" : "excluded")); -#if 0 - TNamed* maxP = new TNamed("maxParticle", Form("%d", fMaxParticles)); - TNamed* method = new TNamed("method", - (fUsePoisson ? "Poisson" : "Energy loss")); - TNamed* phiA = new TNamed("phiAcceptance", - (fUsePhiAcceptance == 0 ? "disabled" : - fUsePhiAcceptance == 1 ? "particles" : - "energy loss")); - TNamed* etaL = new TNamed("etaLumping", Form("%d", fEtaLumping)); - TNamed* phiL = new TNamed("phiLumping", Form("%d", fPhiLumping)); - // TParameter* nxi = new TParameter("nXi", fNXi); -#else TObject* maxP = AliForwardUtil::MakeParameter("maxParticle", fMaxParticles); TObject* method = AliForwardUtil::MakeParameter("method", fUsePoisson); TObject* phiA = AliForwardUtil::MakeParameter("phiAcceptance", fUsePhiAcceptance); TObject* etaL = AliForwardUtil::MakeParameter("etaLumping", fEtaLumping); TObject* phiL = AliForwardUtil::MakeParameter("phiLumping", fPhiLumping); -#endif + TObject* reEt = AliForwardUtil::MakeParameter("recalcEta", fRecalculateEta); + TObject* rePh = AliForwardUtil::MakeParameter("recalcPhi", fRecalculatePhi); + // d->Add(sigma); d->Add(maxP); d->Add(method); d->Add(phiA); d->Add(etaL); d->Add(phiL); + d->Add(reEt); + d->Add(rePh); // d->Add(nxi); - fCuts.Output(d,0); + fCuts.Output(d,"lCuts"); TIter next(&fRingHistos); RingHistos* o = 0; @@ -884,11 +930,18 @@ AliFMDDensityCalculator::Print(Option_t* option) const << std::boolalpha << ind << " Max(particles): " << fMaxParticles << '\n' << ind << " Poisson method: " << fUsePoisson << '\n' - << ind << " Use phi acceptance: " << fUsePhiAcceptance << '\n' << ind << " Eta lumping: " << fEtaLumping << '\n' << ind << " Phi lumping: " << fPhiLumping << '\n' + << ind << " Recalculate eta: " << fRecalculateEta << '\n' + << ind << " Recalculate phi: " << fRecalculatePhi << '\n' + << ind << " Use phi acceptance: " << std::noboolalpha << std::flush; + switch (fUsePhiAcceptance) { + case kPhiNoCorrect: std::cout << "none\n"; break; + case kPhiCorrectNch: std::cout << "correct Nch\n"; break; + case kPhiCorrectELoss: std::cout << "correct energy loss\n"; break; + } std::cout << ind << " Lower cut:" << std::endl; fCuts.Print(); TString opt(option); @@ -925,6 +978,7 @@ AliFMDDensityCalculator::Print(Option_t* option) const //==================================================================== AliFMDDensityCalculator::RingHistos::RingHistos() : AliForwardUtil::RingHistos(), + fList(0), fEvsN(0), fEvsM(0), fEtaVsN(0), @@ -932,10 +986,16 @@ AliFMDDensityCalculator::RingHistos::RingHistos() fCorr(0), fDensity(0), fELossVsPoisson(0), + fDiffELossPoisson(0), fPoisson(), fELoss(0), fELossUsed(0), - fMultCut(0) + fMultCut(0), + fTotal(0), + fGood(0), + fPhiAcc(0), + fPhiBefore(0), + fPhiAfter(0) { // // Default CTOR @@ -945,6 +1005,7 @@ AliFMDDensityCalculator::RingHistos::RingHistos() //____________________________________________________________________ AliFMDDensityCalculator::RingHistos::RingHistos(UShort_t d, Char_t r) : AliForwardUtil::RingHistos(d,r), + fList(0), fEvsN(0), fEvsM(0), fEtaVsN(0), @@ -952,10 +1013,16 @@ AliFMDDensityCalculator::RingHistos::RingHistos(UShort_t d, Char_t r) fCorr(0), fDensity(0), fELossVsPoisson(0), + fDiffELossPoisson(0), fPoisson("ignored"), fELoss(0), fELossUsed(0), - fMultCut(0) + fMultCut(0), + fTotal(0), + fGood(0), + fPhiAcc(0), + fPhiBefore(0), + fPhiAfter(0) { // // Constructor @@ -966,7 +1033,7 @@ AliFMDDensityCalculator::RingHistos::RingHistos(UShort_t d, Char_t r) // fEvsN = new TH2D("elossVsNnocorr", "#Delta E/#Delta E_{mip} vs uncorrected inclusive N_{ch}", - 250, -.5, 24.5, 250, -.5, 24.5); + 250, -.5, 24.5, 251, -1.5, 24.5); fEvsN->SetXTitle("#Delta E/#Delta E_{mip}"); fEvsN->SetYTitle("Inclusive N_{ch} (uncorrected)"); fEvsN->Sumw2(); @@ -1002,8 +1069,7 @@ AliFMDDensityCalculator::RingHistos::RingHistos(UShort_t d, Char_t r) 200, -4, 6, (r == 'I' || r == 'i' ? 20 : 40), 0, 2*TMath::Pi()); fDensity->SetDirectory(0); - fDensity->Sumw2(); - fDensity->SetMarkerColor(Color()); + fDensity->Sumw2(); fDensity->SetMarkerColor(Color()); fDensity->SetXTitle("#eta"); fDensity->SetYTitle("#phi [radians]"); fDensity->SetZTitle("Inclusive N_{ch} density"); @@ -1016,6 +1082,16 @@ AliFMDDensityCalculator::RingHistos::RingHistos(UShort_t d, Char_t r) fELossVsPoisson->SetYTitle("N_{ch} from Poisson"); fELossVsPoisson->SetZTitle("Correlation"); + fDiffELossPoisson = new TH1D("diffElossPoisson", + "(N_{ch,#DeltaE}-N_{ch,Poisson})/N_{ch,Poisson}", + 100, -1, 1); + fDiffELossPoisson->SetDirectory(0); + fDiffELossPoisson->SetXTitle(fDiffELossPoisson->GetTitle()); + fDiffELossPoisson->SetYTitle("Frequency"); + fDiffELossPoisson->SetMarkerColor(Color()); + fDiffELossPoisson->SetFillColor(Color()); + fDiffELossPoisson->Sumw2(); + fELoss = new TH1D("eloss", "#Delta/#Delta_{mip} in all strips", 600, 0, 15); fELoss->SetXTitle("#Delta/#Delta_{mip} (selected)"); @@ -1032,11 +1108,26 @@ AliFMDDensityCalculator::RingHistos::RingHistos(UShort_t d, Char_t r) fELossUsed->SetFillStyle(3002); fELossUsed->SetLineStyle(1); fELossUsed->SetDirectory(0); + + fPhiBefore = new TH1D("phiBefore", "#phi distribution (before recalc)", + (r == 'I' || r == 'i' ? 20 : 40), 0, 2*TMath::Pi()); + fPhiBefore->SetDirectory(0); + fPhiBefore->SetXTitle("#phi"); + fPhiBefore->SetYTitle("Events"); + fPhiBefore->SetMarkerColor(Color()); + fPhiBefore->SetLineColor(Color()); + fPhiBefore->SetFillColor(Color()); + fPhiBefore->SetMarkerStyle(20); + + fPhiAfter = static_cast(fPhiBefore->Clone("phiAfter")); + fPhiAfter->SetTitle("#phi distribution (after re-calc)"); + fPhiAfter->SetDirectory(0); } //____________________________________________________________________ AliFMDDensityCalculator::RingHistos::RingHistos(const RingHistos& o) - : AliForwardUtil::RingHistos(o), + : AliForwardUtil::RingHistos(o), + fList(o.fList), fEvsN(o.fEvsN), fEvsM(o.fEvsM), fEtaVsN(o.fEtaVsN), @@ -1044,10 +1135,16 @@ AliFMDDensityCalculator::RingHistos::RingHistos(const RingHistos& o) fCorr(o.fCorr), fDensity(o.fDensity), fELossVsPoisson(o.fELossVsPoisson), + fDiffELossPoisson(o.fDiffELossPoisson), fPoisson(o.fPoisson), fELoss(o.fELoss), fELossUsed(o.fELossUsed), - fMultCut(o.fMultCut) + fMultCut(o.fMultCut), + fTotal(o.fTotal), + fGood(o.fGood), + fPhiAcc(o.fPhiAcc), + fPhiBefore(o.fPhiBefore), + fPhiAfter(o.fPhiAfter) { // // Copy constructor @@ -1073,25 +1170,36 @@ AliFMDDensityCalculator::RingHistos::operator=(const RingHistos& o) if (&o == this) return *this; AliForwardUtil::RingHistos::operator=(o); - if (fEvsN) delete fEvsN; - if (fEvsM) delete fEvsM; - if (fEtaVsN) delete fEtaVsN; - if (fEtaVsM) delete fEtaVsM; - if (fCorr) delete fCorr; - if (fDensity) delete fDensity; - if (fELossVsPoisson) delete fELossVsPoisson; - - fEvsN = static_cast(o.fEvsN->Clone()); - fEvsM = static_cast(o.fEvsM->Clone()); - fEtaVsN = static_cast(o.fEtaVsN->Clone()); - fEtaVsM = static_cast(o.fEtaVsM->Clone()); - fCorr = static_cast(o.fCorr->Clone()); - fDensity = static_cast(o.fDensity->Clone()); - fELossVsPoisson = static_cast(o.fELossVsPoisson->Clone()); - fPoisson = o.fPoisson; - fELoss = static_cast(o.fELoss->Clone()); - fELossUsed = static_cast(o.fELossUsed->Clone()); - + if (fEvsN) delete fEvsN; + if (fEvsM) delete fEvsM; + if (fEtaVsN) delete fEtaVsN; + if (fEtaVsM) delete fEtaVsM; + if (fCorr) delete fCorr; + if (fDensity) delete fDensity; + if (fELossVsPoisson) delete fELossVsPoisson; + if (fDiffELossPoisson) delete fDiffELossPoisson; + if (fTotal) delete fTotal; + if (fGood) delete fGood; + if (fPhiAcc) delete fPhiAcc; + if (fPhiBefore) delete fPhiBefore; + if (fPhiAfter) delete fPhiAfter; + + fEvsN = static_cast(o.fEvsN->Clone()); + fEvsM = static_cast(o.fEvsM->Clone()); + fEtaVsN = static_cast(o.fEtaVsN->Clone()); + fEtaVsM = static_cast(o.fEtaVsM->Clone()); + fCorr = static_cast(o.fCorr->Clone()); + fDensity = static_cast(o.fDensity->Clone()); + fELossVsPoisson = static_cast(o.fELossVsPoisson->Clone()); + fDiffELossPoisson = static_cast(o.fDiffELossPoisson->Clone()); + fPoisson = o.fPoisson; + fELoss = static_cast(o.fELoss->Clone()); + fELossUsed = static_cast(o.fELossUsed->Clone()); + fTotal = static_cast(o.fTotal->Clone()); + fGood = static_cast(o.fGood->Clone()); + fPhiAcc = static_cast(o.fPhiAcc->Clone()); + fPhiBefore = static_cast(o.fPhiBefore->Clone()); + fPhiAfter = static_cast(o.fPhiAfter->Clone()); return *this; } //____________________________________________________________________ @@ -1105,9 +1213,29 @@ AliFMDDensityCalculator::RingHistos::~RingHistos() //____________________________________________________________________ void -AliFMDDensityCalculator::RingHistos::Init(const TAxis& /*eAxis*/) +AliFMDDensityCalculator::RingHistos::Init(const TAxis& eAxis) { + // Initialize + // This is called on first event fPoisson.Init(-1,-1); + fTotal = new TH1D("total", "Total # of strips per #eta", + eAxis.GetNbins(), eAxis.GetXmin(), eAxis.GetXmax()); + fTotal->SetDirectory(0); + fTotal->SetXTitle("#eta"); + fTotal->SetYTitle("# of strips"); + fGood = static_cast(fTotal->Clone("good")); + fGood->SetTitle("# of good strips per #eta"); + fGood->SetDirectory(0); + + fPhiAcc = new TH2D("phiAcc", "#phi acceptance vs Ip_{z}", + eAxis.GetNbins(), eAxis.GetXmin(), eAxis.GetXmax(), + 10, -10, 10); + fPhiAcc->SetXTitle("#eta"); + fPhiAcc->SetYTitle("v_{z} [cm]"); + fPhiAcc->SetZTitle("#phi acceptance"); + fPhiAcc->SetDirectory(0); + + if (fList) fList->Add(fPhiAcc); } //____________________________________________________________________ @@ -1115,7 +1243,7 @@ void AliFMDDensityCalculator::RingHistos::Output(TList* dir) { // - // Make output + // Make output. This is called as part of SlaveBegin // // Parameters: // dir Where to put it @@ -1128,12 +1256,15 @@ AliFMDDensityCalculator::RingHistos::Output(TList* dir) d->Add(fCorr); d->Add(fDensity); d->Add(fELossVsPoisson); + d->Add(fDiffELossPoisson); fPoisson.Output(d); fPoisson.GetOccupancy()->SetFillColor(Color()); fPoisson.GetMean()->SetFillColor(Color()); fPoisson.GetOccupancy()->SetFillColor(Color()); d->Add(fELoss); d->Add(fELossUsed); + d->Add(fPhiBefore); + d->Add(fPhiAfter); Bool_t inner = (fRing == 'I' || fRing == 'i'); Int_t nStr = inner ? 512 : 256; @@ -1145,6 +1276,7 @@ AliFMDDensityCalculator::RingHistos::Output(TList* dir) fPoisson.Define(x, y); d->Add(AliForwardUtil::MakeParameter("cut", fMultCut)); + fList = d; } //____________________________________________________________________ diff --git a/PWGLF/FORWARD/analysis2/AliFMDDensityCalculator.h b/PWGLF/FORWARD/analysis2/AliFMDDensityCalculator.h index 574eba11ad0..4168ee33ac9 100644 --- a/PWGLF/FORWARD/analysis2/AliFMDDensityCalculator.h +++ b/PWGLF/FORWARD/analysis2/AliFMDDensityCalculator.h @@ -16,6 +16,7 @@ #include #include #include +#include #include "AliForwardUtil.h" #include "AliFMDMultCuts.h" #include "AliPoissonCalculator.h" @@ -60,6 +61,10 @@ public: /** Correct the energy loss */ kPhiCorrectELoss }; + /** + * Folder name + */ + static const char* fgkFolderName; /** * Constructor */ @@ -99,7 +104,6 @@ public: * * @param fmd AliESDFMD object (possibly) corrected for sharing * @param hists Histogram cache - * @param vtxBin Vertex bin * @param lowFlux Low flux flag. * @param cent Centrality * @param vz Vertex Z position @@ -108,10 +112,9 @@ public: */ virtual Bool_t Calculate(const AliESDFMD& fmd, AliForwardUtil::Histos& hists, - UShort_t vtxBin, Bool_t lowFlux, Double_t cent=-1, - Double_t vz=0); + const TVector3& ip=TVector3(1024,1024,0)); /** * Scale the histograms to the total number of events * @@ -159,6 +162,13 @@ public: * */ void SetRecalculateEta(Bool_t use) { fRecalculateEta = use; } + /** + * In case of a displaced vertices recalculate eta and angle correction + * + * @param use recalculate or not + * + */ + void SetRecalculatePhi(Bool_t use) { fRecalculatePhi = use; } /** * Set whether to use the phi acceptance correction. * @@ -277,17 +287,16 @@ protected: * @param mult Signal * @param d Detector * @param r Ring - * @param s Sector - * @param t Strip (not used) - * @param v Vertex bin * @param eta Pseudo-rapidity * @param lowFlux Low-flux flag * * @return The number of particles */ - virtual Float_t NParticles(Float_t mult, - UShort_t d, Char_t r, UShort_t s, UShort_t t, - UShort_t v, Float_t eta, Bool_t lowFlux) const; + virtual Float_t NParticles(Float_t mult, + UShort_t d, + Char_t r, + Float_t eta, + Bool_t lowFlux) const; /** * Get the inverse correction factor. This consist of * @@ -305,8 +314,8 @@ protected: * * @return */ - virtual Float_t Correction(UShort_t d, Char_t r, UShort_t s, UShort_t t, - UShort_t v, Float_t eta, Bool_t lowFlux) const; + virtual Float_t Correction(UShort_t d, Char_t r, UShort_t t, + Float_t eta, Bool_t lowFlux) const; /** * Get the acceptance correction for strip @a t in an ring of type @a r * @@ -377,6 +386,7 @@ protected: * @param nEvents Number of events */ void ScaleHistograms(TList* dir, Int_t nEvents); + TList* fList; TH2D* fEvsN; // Correlation of Eloss vs uncorrected Nch TH2D* fEvsM; // Correlation of Eloss vs corrected Nch TProfile* fEtaVsN; // Average uncorrected Nch vs eta @@ -384,12 +394,17 @@ protected: TProfile* fCorr; // Average correction vs eta TH2D* fDensity; // Distribution inclusive Nch TH2D* fELossVsPoisson; // Correlation of energy loss vs Poisson N_ch + TH1D* fDiffELossPoisson;// Relative difference to Poisson AliPoissonCalculator fPoisson; // Calculate density using Poisson method TH1D* fELoss; // Energy loss as seen by this TH1D* fELossUsed; // Energy loss in strips with signal Double_t fMultCut; // If set, use this - - ClassDef(RingHistos,6); + TH1D* fTotal; // Total number of strips per eta + TH1D* fGood; // Number of good strips per eta + TH2D* fPhiAcc; // Phi acceptance vs IpZ + TH1D* fPhiBefore; // Phi before re-calce + TH1D* fPhiAfter; // Phi after re-calc + ClassDef(RingHistos,8); }; /** * Get the ring histogram container @@ -420,7 +435,8 @@ protected: Int_t fPhiLumping; // How to lump phi bins for Poisson Int_t fDebug; // Debug level AliFMDMultCuts fCuts; // Cuts - Bool_t fRecalculateEta; // //Whether to recalculate eta and angle correction (disp vtx) + Bool_t fRecalculateEta; // Whether to recalc eta and angle correction (disp vtx) + Bool_t fRecalculatePhi; // Whether to correct for (X,Y) offset ClassDef(AliFMDDensityCalculator,7); // Calculate Nch density }; diff --git a/PWGLF/FORWARD/analysis2/AliFMDEnergyFitterTask.cxx b/PWGLF/FORWARD/analysis2/AliFMDEnergyFitterTask.cxx index 223fb753fee..91fcc8d24a6 100644 --- a/PWGLF/FORWARD/analysis2/AliFMDEnergyFitterTask.cxx +++ b/PWGLF/FORWARD/analysis2/AliFMDEnergyFitterTask.cxx @@ -44,7 +44,7 @@ AliFMDEnergyFitterTask::AliFMDEnergyFitterTask() // // Constructor // - DGUARD(fDebug,0,"Default CTOR of AliFMDEnergyFitterTask"); + DGUARD(fDebug, 3,"Default CTOR of AliFMDEnergyFitterTask"); } //____________________________________________________________________ @@ -63,7 +63,7 @@ AliFMDEnergyFitterTask::AliFMDEnergyFitterTask(const char* name) // Parameters: // name Name of task // - DGUARD(fDebug,0,"Named CTOR of AliFMDEnergyFitterTask: %s", name); + DGUARD(fDebug, 3,"Named CTOR of AliFMDEnergyFitterTask: %s", name); DefineOutput(1, TList::Class()); DefineOutput(2, TList::Class()); fBranchNames = @@ -87,7 +87,7 @@ AliFMDEnergyFitterTask::AliFMDEnergyFitterTask(const AliFMDEnergyFitterTask& o) // Parameters: // o Object to copy from // - DGUARD(fDebug,0,"COPY CTOR of AliFMDEnergyFitterTask"); + DGUARD(fDebug, 3,"COPY CTOR of AliFMDEnergyFitterTask"); DefineOutput(1, TList::Class()); DefineOutput(2, TList::Class()); } @@ -250,11 +250,11 @@ AliFMDEnergyFitterTask::UserExec(Option_t*) Bool_t lowFlux = kFALSE; UInt_t triggers = 0; UShort_t ivz = 0; - Double_t vz = 0; + TVector3 ip; Double_t cent = 0; UShort_t nClusters = 0; UInt_t found = fEventInspector.Process(esd, triggers, lowFlux, - ivz, vz, cent, nClusters); + ivz, ip, cent, nClusters); if (found & AliFMDEventInspector::kNoEvent) return; if (found & AliFMDEventInspector::kNoTriggers) return; if (found & AliFMDEventInspector::kNoSPD) return; diff --git a/PWGLF/FORWARD/analysis2/AliFMDEventInspector.cxx b/PWGLF/FORWARD/analysis2/AliFMDEventInspector.cxx index bc199bc9c23..28cf603b12a 100644 --- a/PWGLF/FORWARD/analysis2/AliFMDEventInspector.cxx +++ b/PWGLF/FORWARD/analysis2/AliFMDEventInspector.cxx @@ -41,6 +41,9 @@ #include "AliVVZERO.h" //==================================================================== +const char* AliFMDEventInspector::fgkFolderName = "fmdEventInspector"; + +//____________________________________________________________________ AliFMDEventInspector::AliFMDEventInspector() : TNamed(), fHEventsTr(0), @@ -50,23 +53,23 @@ AliFMDEventInspector::AliFMDEventInspector() fHTriggers(0), fHTriggerCorr(0), fHType(0), - fHWords(0), - fHCent(0), - fHCentVsQual(0), - fHStatus(0), - fLowFluxCut(1000), - fMaxVzErr(0.2), - fList(0), - fEnergy(0), + fHWords(0), + fHCent(0), + fHCentVsQual(0), + fHStatus(0), + fLowFluxCut(1000), + fMaxVzErr(0.2), + fList(0), + fEnergy(0), fField(999), - fCollisionSystem(kUnknown), + fCollisionSystem(kUnknown), fDebug(0), fCentAxis(0), fVtxAxis(10,-10,10), - fUseFirstPhysicsVertex(true), - fUseV0AND(false), - fMinPileupContrib(3), - fMinPileupDistance(0.8), + fUseFirstPhysicsVertex(false), + fUseV0AND(false), + fMinPileupContrib(3), + fMinPileupDistance(0.8), fUseDisplacedVertices(false), fDisplacedVertex(), fCollWords(), @@ -80,7 +83,7 @@ AliFMDEventInspector::AliFMDEventInspector() //____________________________________________________________________ AliFMDEventInspector::AliFMDEventInspector(const char* name) - : TNamed("fmdEventInspector", name), + : TNamed(fgkFolderName, name), fHEventsTr(0), fHEventsTrVtx(0), fHEventsAccepted(0), @@ -101,7 +104,7 @@ AliFMDEventInspector::AliFMDEventInspector(const char* name) fDebug(0), fCentAxis(0), fVtxAxis(10,-10,10), - fUseFirstPhysicsVertex(true), + fUseFirstPhysicsVertex(false), fUseV0AND(false), fMinPileupContrib(3), fMinPileupDistance(0.8), @@ -232,6 +235,7 @@ Bool_t AliFMDEventInspector::FetchHistograms(const TList* d, TH1I*& hEventsTr, TH1I*& hEventsTrVtx, + TH1I*& hEventsAcc, TH1I*& hTriggers) const { // @@ -249,15 +253,20 @@ AliFMDEventInspector::FetchHistograms(const TList* d, DGUARD(fDebug,3,"Fetch histograms in AliFMDEventInspector"); hEventsTr = 0; hEventsTrVtx = 0; - hTriggers = 0; + hEventsAcc = 0; + hTriggers = 0; TList* dd = dynamic_cast(d->FindObject(GetName())); if (!dd) return kFALSE; hEventsTr = dynamic_cast(dd->FindObject("nEventsTr")); hEventsTrVtx = dynamic_cast(dd->FindObject("nEventsTrVtx")); + hEventsAcc = dynamic_cast(dd->FindObject("nEventsAccepted")); hTriggers = dynamic_cast(dd->FindObject("triggers")); - if (!hEventsTr || !hEventsTrVtx || !hTriggers) return kFALSE; + if (!hEventsTr || + !hEventsTrVtx || + !hEventsAcc || + !hTriggers) return kFALSE; return kTRUE; } //____________________________________________________________________ @@ -559,7 +568,7 @@ AliFMDEventInspector::Process(const AliESDEvent* event, UInt_t& triggers, Bool_t& lowFlux, UShort_t& ivz, - Double_t& vz, + TVector3& ip, Double_t& cent, UShort_t& nClusters) { @@ -627,35 +636,30 @@ AliFMDEventInspector::Process(const AliESDEvent* event, if (qual & (1 << i)) fHCentVsQual->Fill(Double_t(i+1), cent); } - // --- Get the vertex information ---------------------------------- - Double_t vx = 0; - Double_t vy = 0; - vz = 0; - - Bool_t vzOk = ReadVertex(*event, vz,vx,vy); - - fHEventsTr->Fill(vz); + // --- Get the interaction point ----------------------------------- + Bool_t vzOk = ReadVertex(*event, ip); + fHEventsTr->Fill(ip.Z()); if (!vzOk) { if (fDebug > 3) { AliWarning("Failed to read vertex from ESD"); } fHStatus->Fill(6); return kNoVertex; } - fHEventsTrVtx->Fill(vz); + fHEventsTrVtx->Fill(ip.Z()); // --- Get the vertex bin ------------------------------------------ - ivz = fVtxAxis.FindBin(vz); + ivz = fVtxAxis.FindBin(ip.Z()); if (ivz <= 0 || ivz > fVtxAxis.GetNbins()) { if (fDebug > 3) { AliWarning(Form("Vertex @ %f outside of range [%f,%f]", - vz, fVtxAxis.GetXmin(), fVtxAxis.GetXmax())); + ip.Z(), fVtxAxis.GetXmin(), fVtxAxis.GetXmax())); } ivz = 0; fHStatus->Fill(7); return kBadVertex; } - fHEventsAccepted->Fill(vz); - fHEventsAcceptedXY->Fill(vx,vy); + fHEventsAccepted->Fill(ip.Z()); + fHEventsAcceptedXY->Fill(ip.X(),ip.Y()); // --- Check the FMD ESD data -------------------------------------- if (!event->GetFMDData()) { @@ -999,10 +1003,7 @@ AliFMDEventInspector::CheckWords(const AliESDEvent& esd, UInt_t& triggers) const //____________________________________________________________________ Bool_t -AliFMDEventInspector::ReadVertex(const AliESDEvent& esd, - Double_t& vz, - Double_t& vx, - Double_t& vy) +AliFMDEventInspector::ReadVertex(const AliESDEvent& esd, TVector3& ip) { // // Read the vertex information from the ESD event @@ -1015,32 +1016,28 @@ AliFMDEventInspector::ReadVertex(const AliESDEvent& esd, // @c true on success, @c false otherwise // DGUARD(fDebug,2,"Read the vertex in AliFMDEventInspector"); - vz = 0; - vx = 1024; - vy = 1024; + ip.SetXYZ(1024, 1024, 0); if(fUseDisplacedVertices) { Double_t zvtx = fDisplacedVertex.GetVertexZ(); if(TMath::Abs(zvtx) < 999) { - vz = zvtx; + ip.SetZ(zvtx); return true; } return false; } - if(fUseFirstPhysicsVertex) return CheckPWGUDVertex(esd, vz, vx, vy); + if(fUseFirstPhysicsVertex) return CheckPWGUDVertex(esd, ip); - return CheckVertex(esd, vz, vx, vy); + return CheckVertex(esd, ip); } //____________________________________________________________________ Bool_t AliFMDEventInspector::CheckPWGUDVertex(const AliESDEvent& esd, - Double_t& vz, - Double_t& vx, - Double_t& vy) const + TVector3& ip) const { // This is the code used by the 1st physics people const AliESDVertex* vertex = esd.GetPrimaryVertex(); @@ -1066,20 +1063,18 @@ AliFMDEventInspector::CheckPWGUDVertex(const AliESDEvent& esd, return false; } } - vz = vertex->GetZ(); + ip.SetZ(vertex->GetZ()); if(!vertex->IsFromVertexerZ()) { - vx = vertex->GetX(); - vy = vertex->GetY(); + ip.SetX(vertex->GetX()); + ip.SetY(vertex->GetY()); } return true; } //____________________________________________________________________ Bool_t AliFMDEventInspector::CheckVertex(const AliESDEvent& esd, - Double_t& vz, - Double_t& vx, - Double_t& vy) const + TVector3& ip) const { // Use standard SPD vertex (perhaps preferable for Pb+Pb) // Get the vertex @@ -1094,7 +1089,7 @@ AliFMDEventInspector::CheckVertex(const AliESDEvent& esd, if(vertex->GetNContributors() <= 0) { DMSG(fDebug,2,"Number of contributors to vertex is %d<=0", vertex->GetNContributors()); - vz = 0; + ip.SetZ(0); return false; } // Check that the uncertainty isn't too large @@ -1105,12 +1100,13 @@ AliFMDEventInspector::CheckVertex(const AliESDEvent& esd, } // Get the z coordiante - vz = vertex->GetZ(); + ip.SetZ(vertex->GetZ()); const AliESDVertex* vertexXY = esd.GetPrimaryVertex(); + if(!vertexXY->IsFromVertexerZ()) { - vx = vertexXY->GetX(); - vy = vertexXY->GetY(); + ip.SetX(vertexXY->GetX()); + ip.SetY(vertexXY->GetY()); } return true; } diff --git a/PWGLF/FORWARD/analysis2/AliFMDEventInspector.h b/PWGLF/FORWARD/analysis2/AliFMDEventInspector.h index f26f3a6accf..fc035d6f906 100644 --- a/PWGLF/FORWARD/analysis2/AliFMDEventInspector.h +++ b/PWGLF/FORWARD/analysis2/AliFMDEventInspector.h @@ -26,6 +26,7 @@ class TH1F; class TH2F; class TH2I; class TAxis; +class TVector3; // class TList; /** @@ -94,6 +95,10 @@ public: kPP, kPbPb }; + /** + * Folder name + */ + static const char* fgkFolderName; /** * Constructor */ @@ -138,7 +143,9 @@ public: * event (according to the setting of fLowFluxCut) * @param ivz On return, the found vertex bin (1-based). A zero * means outside of the defined vertex range - * @param vz On return, the z position of the interaction + * @param vx On return, the x position of the interaction point + * @param vy On return, the y position of the interaction point + * @param vz On return, the z position of the interaction point * @param cent On return, the centrality (in percent) or < 0 * if not found * @param nClusters On return, number of SPD clusters in @f$ |\eta|<1@f$ @@ -149,7 +156,7 @@ public: UInt_t& triggers, Bool_t& lowFlux, UShort_t& ivz, - Double_t& vz, + TVector3& ip, Double_t& cent, UShort_t& nClusters); /** @@ -225,6 +232,7 @@ public: * @param d Input * @param hEventsTr On return, pointer to histogram, or null * @param hEventsTrVtx On return, pointer to histogram, or null + * @param hEventsAcc On return, pointer to histogram, or null * @param hTriggers On return, pointer to histogram, or null * * @return true on success, false otherwise @@ -232,6 +240,7 @@ public: Bool_t FetchHistograms(const TList* d, TH1I*& hEventsTr, TH1I*& hEventsTrVtx, + TH1I*& hEventsAcc, TH1I*& hTriggers) const; /** * Read the collision system, collision energy, and L3 field setting @@ -381,16 +390,11 @@ protected: * Read the vertex information from the ESD event * * @param esd ESD event - * @param vz On return, the vertex Z position - * @param vx On return, the vertex X position - * @param vy On return, the vertex Y position + * @param ip On return, the coordinates of the IP * * @return @c true on success, @c false otherwise */ - Bool_t ReadVertex(const AliESDEvent& esd, - Double_t& vz, - Double_t& vx, - Double_t& vy); + Bool_t ReadVertex(const AliESDEvent& esd, TVector3& ip); /** * Check the vertex using the method used in PWG-UD. That is * @@ -400,16 +404,12 @@ protected: * - Check that the dispersion and resolution is OK * * @param esd Data - * @param vz On return, the Z coordinate of the IP - * @param vx On return, possibly the X coordinate of the IP - * @param vy On return, possibly the Y coordinate of the IP + * @param ip On return, the coordinates of the IP * * @return true if the vertex was found and met the requirements */ virtual Bool_t CheckPWGUDVertex(const AliESDEvent& esd, - Double_t& vz, - Double_t& vx, - Double_t& vy) const; + TVector3& ip) const; /** * Check the vertex. That is * @@ -418,16 +418,11 @@ protected: * - Check that the reslution is OK * * @param esd Data - * @param vz On return, the Z coordinate of the IP - * @param vx On return, possibly the X coordinate of the IP - * @param vy On return, possibly the Y coordinate of the IP + * @param ip On return, the coordinates of the IP * * @return true if the vertex was found and met the requirements */ - virtual Bool_t CheckVertex(const AliESDEvent& esd, - Double_t& vz, - Double_t& vx, - Double_t& vy) const; + virtual Bool_t CheckVertex(const AliESDEvent& esd, TVector3& ip) const; /** * Read centrality from event * diff --git a/PWGLF/FORWARD/analysis2/AliFMDHistCollector.cxx b/PWGLF/FORWARD/analysis2/AliFMDHistCollector.cxx index 1b72027a8c3..cc5172766f0 100644 --- a/PWGLF/FORWARD/analysis2/AliFMDHistCollector.cxx +++ b/PWGLF/FORWARD/analysis2/AliFMDHistCollector.cxx @@ -48,7 +48,7 @@ AliFMDHistCollector::AliFMDHistCollector() fSkipFMDRings(0), fBgAndHitMaps(false) { - DGUARD(fDebug, 0, "Default CTOR of AliFMDHistCollector"); + DGUARD(fDebug, 3, "Default CTOR of AliFMDHistCollector"); } //____________________________________________________________________ @@ -67,7 +67,7 @@ AliFMDHistCollector::AliFMDHistCollector(const char* title) fSkipFMDRings(0), fBgAndHitMaps(false) { - DGUARD(fDebug, 0, "Named CTOR of AliFMDHistCollector: %s", title); + DGUARD(fDebug, 3, "Named CTOR of AliFMDHistCollector: %s", title); } //____________________________________________________________________ AliFMDHistCollector::AliFMDHistCollector(const AliFMDHistCollector& o) @@ -85,7 +85,7 @@ AliFMDHistCollector::AliFMDHistCollector(const AliFMDHistCollector& o) fSkipFMDRings(o.fSkipFMDRings), fBgAndHitMaps(o.fBgAndHitMaps) { - DGUARD(fDebug, 0, "Copy CTOR of AliFMDHistCollector"); + DGUARD(fDebug, 3, "Copy CTOR of AliFMDHistCollector"); } //____________________________________________________________________ @@ -555,6 +555,10 @@ AliFMDHistCollector::MergeBins(Double_t c, Double_t e, re = oe; } break; + case kSum: + rc = c + oc; + re = TMath::Sqrt(oe * oe + e * e);//Add in quadarature + break; default: AliError("No method for defining content of overlapping bins defined"); return; @@ -613,7 +617,8 @@ AliFMDHistCollector::Collect(const AliForwardUtil::Histos& hists, GetFirstAndLast(d, r, vtxbin, first, last); // Zero outside valid range - for (Int_t iPhi = 0; iPhi <= t->GetNbinsY()+1; iPhi++) { + Int_t nY = t->GetNbinsY(); + for (Int_t iPhi = 0; iPhi <= nY+1; iPhi++) { // Lower range for (Int_t iEta = 1; iEta < first; iEta++) { t->SetBinContent(iEta,iPhi,0); @@ -624,27 +629,43 @@ AliFMDHistCollector::Collect(const AliForwardUtil::Histos& hists, t->SetBinError(iEta,iPhi,0); } } - for (Int_t iEta = first; iEta <= last; iEta++) + for (Int_t iEta = first; iEta <= last; iEta++) { + // Double_t phiAcc = t->GetBinContent(iEta, nY+1); + // if (phiAcc > 1e-12 && phiAcc < 1) + // Info("", "FMD%d%c %3d phi acceptance: %f",d,r,iEta,phiAcc); t->SetBinContent(iEta,0,1); + } // Add to our per-ring sum o->Add(t); // Outer rings have better phi segmentation - rebin to same as inner. if (q == 1) t->RebinY(2); + nY = t->GetNbinsY(); // Now update profile output for (Int_t iEta = first; iEta <= last; iEta++) { // Get the possibly overlapping histogram Int_t overlap = GetOverlap(d,r,iEta,vtxbin); - // Fill eta acceptance for this event into the phi underlow bin + // Get factor + Float_t fac = (fMergeMethod != kSum && overlap >= 0 ? .5 : 1); + + // Fill eta acceptance for this event into the phi underflow bin Float_t ooc = out.GetBinContent(iEta,0); - Float_t noc = overlap >= 0? 0.5 : 1; - out.SetBinContent(iEta, 0, ooc + noc); + out.SetBinContent(iEta, 0, ooc + fac); + + // Fill phi acceptance for this event into the phi overflow bin + Float_t oop = out.GetBinContent(iEta,nY+1); + Float_t nop = t->GetBinContent(iEta,nY+1); +#if 0 + Info("", "etaBin=%3d Setting phi acceptance to %f(%f+%f)=%f", + iEta, fac, oop, nop, fac*(oop+nop)); +#endif + out.SetBinContent(iEta, nY+1, fac * nop + oop); // Should we loop over h or t Y bins - I think it's t - for (Int_t iPhi = 1; iPhi <= t->GetNbinsY(); iPhi++) { + for (Int_t iPhi = 1; iPhi <= nY; iPhi++) { Double_t c = t->GetBinContent(iEta,iPhi); Double_t e = t->GetBinError(iEta,iPhi); Double_t ee = t->GetXaxis()->GetBinCenter(iEta); @@ -700,13 +721,16 @@ AliFMDHistCollector::Print(Option_t* /* option */) const ind[gROOT->GetDirLevel()] = '\0'; std::cout << ind << ClassName() << ": " << GetName() << '\n' << ind << " # of cut bins: " << fNCutBins << '\n' - << ind << " Correction cut: " << fCorrectionCut << '\n' + << ind << " Fiducal method: " + << (fFiducialMethod == kByCut ? "cut" : "distance") << "\n" + << ind << " Fiducial cut: " << fCorrectionCut << "\n" << ind << " Merge method: "; switch (fMergeMethod) { case kStraightMean: std::cout << "straight mean\n"; break; case kStraightMeanNoZero: std::cout << "straight mean (no zeros)\n"; break; case kWeightedMean: std::cout << "weighted mean\n"; break; case kLeastError: std::cout << "least error\n"; break; + case kSum: std::cout << "straight sum\n"; break; } std::cout << ind << " Bin ranges:\n" << ind << " rings "; diff --git a/PWGLF/FORWARD/analysis2/AliFMDHistCollector.h b/PWGLF/FORWARD/analysis2/AliFMDHistCollector.h index c91abb2b448..c212a9f99a2 100644 --- a/PWGLF/FORWARD/analysis2/AliFMDHistCollector.h +++ b/PWGLF/FORWARD/analysis2/AliFMDHistCollector.h @@ -77,7 +77,11 @@ public: * c_2 & \mbox{otherwise}\end{array}\right. * @f] */ - kLeastError + kLeastError, + /** + * Just sum the signals + */ + kSum }; /** * How to obtain the fiducial cuts diff --git a/PWGLF/FORWARD/analysis2/AliFMDMCDensityCalculator.cxx b/PWGLF/FORWARD/analysis2/AliFMDMCDensityCalculator.cxx index f670404b46d..348ff61c2f0 100644 --- a/PWGLF/FORWARD/analysis2/AliFMDMCDensityCalculator.cxx +++ b/PWGLF/FORWARD/analysis2/AliFMDMCDensityCalculator.cxx @@ -133,7 +133,7 @@ AliFMDMCDensityCalculator::CalculateMC(const AliESDFMD& fmd, Float_t phi = fmd.Phi(d,r,s,t) / 180 * TMath::Pi(); Float_t eta = fmd.Eta(d,r,s,t); - Float_t corr = Correction(d, r, s, t, 0, eta, false); + Float_t corr = Correction(d, r, t, eta, false); Float_t sig = (corr <= 0 ? 0 : mult / corr); h->Fill(eta,phi,sig); } diff --git a/PWGLF/FORWARD/analysis2/AliFMDMCTrackDensity.h b/PWGLF/FORWARD/analysis2/AliFMDMCTrackDensity.h index e6695257c03..d59493f6e54 100644 --- a/PWGLF/FORWARD/analysis2/AliFMDMCTrackDensity.h +++ b/PWGLF/FORWARD/analysis2/AliFMDMCTrackDensity.h @@ -173,7 +173,7 @@ protected: * @return Reference to this object */ State& operator=(const State& o); - } fState; // State + } fState; //! State UShort_t fMaxConsequtiveStrips; // Max 'cluster' size @@ -183,7 +183,7 @@ protected: TH2D* fNcr; // Number of clusters per track AliESDFMD* fOutput; //! Output ESD object - ClassDef(AliFMDMCTrackDensity,3); // Calculate track-ref density + ClassDef(AliFMDMCTrackDensity,4); // Calculate track-ref density }; #endif diff --git a/PWGLF/FORWARD/analysis2/AliFMDSharingFilter.cxx b/PWGLF/FORWARD/analysis2/AliFMDSharingFilter.cxx index 825b5611ff9..de6d67ff92a 100644 --- a/PWGLF/FORWARD/analysis2/AliFMDSharingFilter.cxx +++ b/PWGLF/FORWARD/analysis2/AliFMDSharingFilter.cxx @@ -22,6 +22,7 @@ // // #include "AliFMDSharingFilter.h" +#include "AliFMDStripIndex.h" #include #include #include @@ -63,7 +64,8 @@ AliFMDSharingFilter::AliFMDSharingFilter() fHCuts(), fUseSimpleMerging(false), fThreeStripSharing(true), - fRecalculateEta(false) + fRecalculateEta(false), + fExtraDead(0) { // // Default Constructor - do not use @@ -86,7 +88,8 @@ AliFMDSharingFilter::AliFMDSharingFilter(const char* title) fHCuts(), fUseSimpleMerging(false), fThreeStripSharing(true), - fRecalculateEta(false) + fRecalculateEta(false), + fExtraDead(51200) { // // Constructor @@ -107,6 +110,8 @@ AliFMDSharingFilter::AliFMDSharingFilter(const char* title) fHCuts.SetNXi(1); fHCuts.SetIncludeSigma(1); fLCuts.SetMultCuts(.15); + + fExtraDead.Reset(-1); } //____________________________________________________________________ @@ -124,7 +129,8 @@ AliFMDSharingFilter::AliFMDSharingFilter(const AliFMDSharingFilter& o) fHCuts(o.fHCuts), fUseSimpleMerging(o.fUseSimpleMerging), fThreeStripSharing(o.fThreeStripSharing), - fRecalculateEta(o.fRecalculateEta) + fRecalculateEta(o.fRecalculateEta), + fExtraDead(o.fExtraDead) { // // Copy constructor @@ -211,6 +217,69 @@ AliFMDSharingFilter::GetRingHistos(UShort_t d, Char_t r) const return static_cast(fRingHistos.At(idx)); } +//____________________________________________________________________ +void +AliFMDSharingFilter::AddDead(UShort_t d, Char_t r, UShort_t s, UShort_t t) +{ + if (d < 1 || d > 3) { + Warning("AddDead", "Invalid detector FMD%d", d); + return; + } + Bool_t inner = (r == 'I' || r == 'i'); + if (d == 1 && !inner) { + Warning("AddDead", "Invalid ring FMD%d%c", d, r); + return; + } + if ((inner && s >= 20) || (!inner && s >= 40)) { + Warning("AddDead", "Invalid sector FMD%d%c[%02d]", d, r, s); + return; + } + if ((inner && t >= 512) || (!inner && t >= 256)) { + Warning("AddDead", "Invalid strip FMD%d%c[%02d,%03d]", d, r, s, t); + return; + } + + Int_t id = AliFMDStripIndex::Pack(d, r, s, t); + Int_t i = 0; + for (i = 0; i < fExtraDead.GetSize(); i++) { + Int_t j = fExtraDead.At(i); + if (j == id) return; // Already there + if (j < 0) break; // Free slot + } + if (i >= fExtraDead.GetSize()) { + Warning("AddDead", "No free slot to add FMD%d%c[%02d,%03d] at", + d, r, s, t); + return; + } + fExtraDead[i] = id; +} +//____________________________________________________________________ +void +AliFMDSharingFilter::AddDeadRegion(UShort_t d, Char_t r, + UShort_t s1, UShort_t s2, + UShort_t t1, UShort_t t2) +{ + // Add a dead region spanning from FMD[,] to + // FMD[,] (both inclusive) + for (Int_t s = s1; s <= s2; s++) + for (Int_t t = t1; t <= t2; t++) + AddDead(d, r, s, t); +} +//____________________________________________________________________ +Bool_t +AliFMDSharingFilter::IsDead(UShort_t d, Char_t r, UShort_t s, UShort_t t) const +{ + Int_t id = AliFMDStripIndex::Pack(d, r, s, t); + for (Int_t i = 0; i < fExtraDead.GetSize(); i++) { + Int_t j = fExtraDead.At(i); + if (j == id) { + //Info("IsDead", "FMD%d%c[%02d,%03d] marked as dead here", d, r, s, t); + return true; + } + if (j < 0) break; // High water mark + } + return false; +} //____________________________________________________________________ void AliFMDSharingFilter::Init(const TAxis& axis) @@ -320,10 +389,6 @@ AliFMDSharingFilter::Filter(const AliESDFMD& input, for (UShort_t t = 0; t < nstr; t++) status[t] = kCandidate; -#ifdef USE_OLDER_MERGING - Bool_t usedThis = kFALSE; - Bool_t usedPrev = kFALSE; -#endif //For simple merging Bool_t used = kFALSE; Double_t eTotal = -1; @@ -363,7 +428,7 @@ AliFMDSharingFilter::Filter(const AliESDFMD& input, } // Keep dead-channel information. - if(mult == AliESDFMD::kInvalidMult) + if(mult == AliESDFMD::kInvalidMult || IsDead(d,r,s,t)) output.SetMultiplicity(d,r,s,t,AliESDFMD::kInvalidMult); // If no signal or dead strip, go on. @@ -462,24 +527,15 @@ AliFMDSharingFilter::Filter(const AliESDFMD& input, } if (t != 0) histos->fNeighborsBefore->Fill(prevE, mult); -#ifdef USE_OLDER_MERGING - /*Double_t*/ mergedEnergy = MultiplicityOfStrip(mult,eta,prevE,nextE, - lowFlux,d,r,s,t, - usedPrev,usedThis); - status[t] = (usedPrev ? kMergedWithOther : kNone); - if (t != nstr - 1) status[t] = (usedThis ? kMergedWithOther : kNone); -#else - /*Double_t*/ mergedEnergy = MultiplicityOfStrip(mult, prevE, nextE, - eta, lowFlux, - d, r, s, t, - prevStatus, - thisStatus, - nextStatus); + mergedEnergy = MultiplicityOfStrip(mult, prevE, nextE, + eta, lowFlux, + d, r, s, t, + prevStatus, + thisStatus, + nextStatus); if (t != 0) status[t-1] = prevStatus; if (t != nstr-1) status[t+1] = nextStatus; - status[t] = thisStatus; - -#endif + status[t] = thisStatus; // If we're processing on non-angle corrected data, we // should do the angle correction here } // End of non-simple @@ -956,6 +1012,20 @@ AliFMDSharingFilter::DefineOutput(TList* dir) d->Add(AliForwardUtil::MakeParameter("lowSignal", fZeroSharedHitsBelowThreshold)); d->Add(AliForwardUtil::MakeParameter("simple", fUseSimpleMerging)); + + TObjArray* extraDead = new TObjArray; + extraDead->SetOwner(); + extraDead->SetName("extraDead"); + for (Int_t i = 0; i < fExtraDead.GetSize(); i++) { + if (fExtraDead.At(i) < 0) break; + UShort_t dd, s, t; + Char_t r; + Int_t id = fExtraDead.At(i); + AliFMDStripIndex::Unpack(id, dd, r, s, t); + extraDead->Add(AliForwardUtil::MakeParameter(Form("FMD%d%c[%02d,%03d]", + dd, r, s, t), id)); + } + d->Add(extraDead); fLCuts.Output(d,"lCuts"); fHCuts.Output(d,"hCuts"); diff --git a/PWGLF/FORWARD/analysis2/AliFMDSharingFilter.h b/PWGLF/FORWARD/analysis2/AliFMDSharingFilter.h index 85cc3bbdbf4..080d465ea57 100644 --- a/PWGLF/FORWARD/analysis2/AliFMDSharingFilter.h +++ b/PWGLF/FORWARD/analysis2/AliFMDSharingFilter.h @@ -220,6 +220,10 @@ public: * @param c Cuts object */ void SetHCuts(const AliFMDMultCuts& c) { fHCuts = c; } + + void AddDead(UShort_t d, Char_t r, UShort_t s, UShort_t t); + void AddDeadRegion(UShort_t d, Char_t r, UShort_t s1, UShort_t s2, + UShort_t t1, UShort_t t2); protected: /** * Internal data structure to keep track of the histograms @@ -368,7 +372,8 @@ protected: * @param nextStatus Next status * * @return The filtered signal in the strip - */ Double_t MultiplicityOfStrip(Double_t thisE, + */ + Double_t MultiplicityOfStrip(Double_t thisE, Double_t prevE, Double_t nextE, Double_t eta, @@ -426,6 +431,7 @@ protected: */ virtual Double_t GetLowCut(UShort_t d, Char_t r, Double_t eta) const; + virtual Bool_t IsDead(UShort_t d, Char_t r, UShort_t s, UShort_t t) const; TList fRingHistos; // List of histogram containers // Double_t fLowCut; // Low cut on sharing Bool_t fCorrectAngles; // Whether to work on angle corrected signals @@ -440,7 +446,8 @@ protected: Bool_t fUseSimpleMerging; //enable simple sharing by HHD Bool_t fThreeStripSharing; //In case of simple sharing allow 3 strips Bool_t fRecalculateEta; //Whether to recalculate eta and angle correction (disp vtx) - ClassDef(AliFMDSharingFilter,4); // + TArrayI fExtraDead; // List of extra dead channels + ClassDef(AliFMDSharingFilter,5); // }; #endif diff --git a/PWGLF/FORWARD/analysis2/AliForwardCorrectionManager.cxx b/PWGLF/FORWARD/analysis2/AliForwardCorrectionManager.cxx index c63e314dab9..f399e127fba 100644 --- a/PWGLF/FORWARD/analysis2/AliForwardCorrectionManager.cxx +++ b/PWGLF/FORWARD/analysis2/AliForwardCorrectionManager.cxx @@ -961,23 +961,25 @@ AliForwardCorrectionManager::WriteFile(ECorrection what, AliError(Form("Failed to write %s to disk (%d)", ofName.Data(), ret)); return false; } - output->ls(); + // output->ls(); output->Close(); +#if 0 TString cName(obj->IsA()->GetName()); AliInfo(Form("Wrote %s object %s to %s", cName.Data(), oName.Data(), ofName.Data())); if (!full) { TString dName(GetFileDir(what)); - AliInfo(Form("%s should be copied to %s" - "Do for example\n\n\t" - "aliroot $ALICE_ROOT/PWGLF/FORWARD/analysis2/scripts/" - "MoveCorrections.C\\(%d\\)\nor\n\t" - "cp %s %s/\n", - ofName.Data(),dName.Data(), - what, ofName.Data(), - gSystem->ExpandPathName(dName.Data()))); + AliInfoF("\n %s should be copied to %s" + "Do for example\n\n\t" + "aliroot $ALICE_ROOT/PWGLF/FORWARD/analysis2/scripts/" + "MoveCorrections.C\\(%d\\)\nor\n\t" + "cp %s %s/\n", + ofName.Data(),dName.Data(), + what, ofName.Data(), + gSystem->ExpandPathName(dName.Data())); } +#endif return true; } diff --git a/PWGLF/FORWARD/analysis2/AliForwardMCCorrectionsTask.cxx b/PWGLF/FORWARD/analysis2/AliForwardMCCorrectionsTask.cxx index 53da417636e..f40451f6da7 100644 --- a/PWGLF/FORWARD/analysis2/AliForwardMCCorrectionsTask.cxx +++ b/PWGLF/FORWARD/analysis2/AliForwardMCCorrectionsTask.cxx @@ -382,7 +382,7 @@ AliForwardMCCorrectionsTask::UserExec(Option_t*) UInt_t triggers = 0; // Trigger bits Bool_t lowFlux = true; // Low flux flag UShort_t iVz = 0; // Vertex bin from ESD - Double_t vZ = 1000; // Z coordinate from ESD + TVector3 ip(1024,1024,1000); Double_t cent = -1; // Centrality UShort_t iVzMc = 0; // Vertex bin from MC Double_t vZMc = 1000; // Z coordinate of IP vertex from MC @@ -393,7 +393,7 @@ AliForwardMCCorrectionsTask::UserExec(Option_t*) Double_t phiR = 100; // Reaction plane from MC UShort_t nClusters = 0; // Number of SPD clusters // Process the data - UInt_t retESD = fInspector.Process(esd, triggers, lowFlux, iVz, vZ, + UInt_t retESD = fInspector.Process(esd, triggers, lowFlux, iVz, ip, cent, nClusters); fInspector.ProcessMC(mcEvent, triggers, iVzMc, vZMc, b, cMC, nPart, nBin, phiR); diff --git a/PWGLF/FORWARD/analysis2/AliForwardMCMultiplicityTask.cxx b/PWGLF/FORWARD/analysis2/AliForwardMCMultiplicityTask.cxx index 3825b2a5cc7..a7a6756c87b 100644 --- a/PWGLF/FORWARD/analysis2/AliForwardMCMultiplicityTask.cxx +++ b/PWGLF/FORWARD/analysis2/AliForwardMCMultiplicityTask.cxx @@ -327,11 +327,11 @@ AliForwardMCMultiplicityTask::UserExec(Option_t*) Bool_t lowFlux = kFALSE; UInt_t triggers = 0; UShort_t ivz = 0; - Double_t vz = 0; + TVector3 ip(1024, 1024, 0); Double_t cent = -1; UShort_t nClusters = 0; UInt_t found = fEventInspector.Process(esd, triggers, lowFlux, - ivz, vz, cent, nClusters); + ivz, ip, cent, nClusters); UShort_t ivzMC = 0; Double_t vzMC = 0; Double_t phiR = 0; @@ -342,7 +342,7 @@ AliForwardMCMultiplicityTask::UserExec(Option_t*) // UInt_t foundMC = fEventInspector.ProcessMC(mcEvent, triggers, ivzMC, vzMC, b, cMC, npart, nbin, phiR); - fEventInspector.CompareResults(vz, vzMC, cent, b, npart, nbin); + fEventInspector.CompareResults(ip.Z(), vzMC, cent, b, npart, nbin); //Store all events MarkEventForStore(); @@ -375,8 +375,8 @@ AliForwardMCMultiplicityTask::UserExec(Option_t*) if (found & AliFMDEventInspector::kNoVertex) isAccepted = false; // return; if (isAccepted) { - fAODFMD.SetIpZ(vz); - fMCAODFMD.SetIpZ(vz); + fAODFMD.SetIpZ(ip.Z()); + fMCAODFMD.SetIpZ(ip.Z()); } if (found & AliFMDEventInspector::kBadVertex) isAccepted = false; // return; @@ -389,11 +389,11 @@ AliForwardMCMultiplicityTask::UserExec(Option_t*) AliESDFMD* esdFMD = esd->GetFMDData(); // Apply the sharing filter (or hit merging or clustering if you like) - if (isAccepted && !fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD, vz)) { + if (isAccepted && !fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD,ip.Z())) { AliWarning("Sharing filter failed!"); return; } - if (!fSharingFilter.FilterMC(*esdFMD, *mcEvent, vz, fMCESDFMD, fPrimary)) { + if (!fSharingFilter.FilterMC(*esdFMD, *mcEvent, ip.Z(),fMCESDFMD,fPrimary)) { AliWarning("MC Sharing filter failed!"); return; } @@ -404,7 +404,7 @@ AliForwardMCMultiplicityTask::UserExec(Option_t*) fSharingFilter.CompareResults(fESDFMD, fMCESDFMD); // Calculate the inclusive charged particle density - if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux, cent, vz)) { + if (!fDensityCalculator.Calculate(fESDFMD, fHistos, lowFlux, cent, ip.Z())) { AliWarning("Density calculator failed!"); return; } @@ -463,46 +463,15 @@ AliForwardMCMultiplicityTask::Terminate(Option_t*) if (GetOutputData(1)) GetOutputData(1)->Print(); return; } - - // Get our histograms from the container - TH1I* hEventsTr = 0; - TH1I* hEventsTrVtx = 0; - TH1I* hTriggers = 0; - if (!fEventInspector.FetchHistograms(list, hEventsTr, - hEventsTrVtx, hTriggers)) { - AliError(Form("Didn't get histograms from event selector " - "(hEventsTr=%p,hEventsTrVtx=%p)", - hEventsTr, hEventsTrVtx)); - list->ls(); - return; - } - - TH2D* hData = static_cast(list->FindObject("d2Ndetadphi")); - if (!hData) { - AliError(Form("Couldn't get our summed histogram from output " - "list %s (d2Ndetadphi=%p)", list->GetName(), hData)); - list->ls(); - return; - } - - // TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e"); - TH1D* dNdeta = hData->ProjectionX("dNdeta", 1, -1, "e"); - TH1D* norm = hData->ProjectionX("norm", 0, 1, ""); - dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions"); - dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}"); - dNdeta->Divide(norm); - dNdeta->SetStats(0); - dNdeta->Scale(Double_t(hEventsTrVtx->GetEntries())/hEventsTr->GetEntries(), - "width"); - list->Add(dNdeta); - list->Add(norm); + Double_t nTr = 0, nTrVtx = 0, nAcc = 0; + MakeSimpledNdeta(list, list, nTr, nTrVtx, nAcc); MakeRingdNdeta(list, "ringSums", list, "ringResults"); MakeRingdNdeta(list, "mcRingSums", list, "mcRingResults", 24); - fSharingFilter.ScaleHistograms(list,Int_t(hEventsTr->Integral())); - fDensityCalculator.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral())); - fCorrections.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral())); + fSharingFilter.ScaleHistograms(list,Int_t(nTr)); + fDensityCalculator.ScaleHistograms(list,Int_t(nTrVtx)); + fCorrections.ScaleHistograms(list,Int_t(nTrVtx)); } diff --git a/PWGLF/FORWARD/analysis2/AliForwardMultiplicityBase.cxx b/PWGLF/FORWARD/analysis2/AliForwardMultiplicityBase.cxx index 137300169f5..71658caed23 100644 --- a/PWGLF/FORWARD/analysis2/AliForwardMultiplicityBase.cxx +++ b/PWGLF/FORWARD/analysis2/AliForwardMultiplicityBase.cxx @@ -40,7 +40,7 @@ AliForwardMultiplicityBase::AliForwardMultiplicityBase(const char* name) fFirstEvent(true), fCorrManager(0) { - DGUARD(0,0,"Named Construction of AliForwardMultiplicityBase %s",name); + DGUARD(fDebug, 3,"Named CTOR of AliForwardMultiplicityBase %s",name); // Set our persistent pointer fCorrManager = &AliForwardCorrectionManager::Instance(); fBranchNames = @@ -251,6 +251,94 @@ AliForwardMultiplicityBase::MarkEventForStore() const ah->SetFillAOD(kTRUE); } +//____________________________________________________________________ +Bool_t +AliForwardMultiplicityBase::MakeSimpledNdeta(const TList* input, + TList* output, + Double_t& nTr, + Double_t& nTrVtx, + Double_t& nAcc) +{ + // Get our histograms from the container + TH1I* hEventsTr = 0; + TH1I* hEventsTrVtx = 0; + TH1I* hEventsAcc = 0; + TH1I* hTriggers = 0; + if (!GetEventInspector().FetchHistograms(input, + hEventsTr, + hEventsTrVtx, + hEventsAcc, + hTriggers)) { + AliError(Form("Didn't get histograms from event selector " + "(hEventsTr=%p,hEventsTrVtx=%p,hEventsAcc=%p,hTriggers=%p)", + hEventsTr, hEventsTrVtx, hEventsAcc, hTriggers)); + input->ls(); + return false; + } + nTr = hEventsTr->Integral(); + nTrVtx = hEventsTrVtx->Integral(); + nAcc = hEventsAcc->Integral(); + Double_t vtxEff = nTrVtx / nTr; + TH2D* hData = static_cast(input->FindObject("d2Ndetadphi")); + if (!hData) { + AliError(Form("Couldn't get our summed histogram from output " + "list %s (d2Ndetadphi=%p)", input->GetName(), hData)); + input->ls(); + return false; + } + + Int_t nY = hData->GetNbinsY(); + TH1D* dNdeta = hData->ProjectionX("dNdeta", 1, nY, "e"); + TH1D* dNdeta_ = hData->ProjectionX("dNdeta_", 1, nY, "e"); + TH1D* norm = hData->ProjectionX("norm", 0, 0, ""); + TH1D* phi = hData->ProjectionX("phi", nY+1, nY+1, ""); + dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions"); + dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}"); + dNdeta->SetMarkerColor(kRed+1); + dNdeta->SetMarkerStyle(20); + dNdeta->SetDirectory(0); + + dNdeta_->SetTitle("dN_{ch}/d#eta in the forward regions"); + dNdeta_->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}"); + dNdeta_->SetMarkerColor(kMagenta+1); + dNdeta_->SetMarkerStyle(21); + dNdeta_->SetDirectory(0); + + norm->SetTitle("Normalization to #eta coverage"); + norm->SetYTitle("#eta coverage"); + norm->SetMarkerColor(kBlue+1); + norm->SetMarkerStyle(21); + norm->SetFillColor(kBlue+1); + norm->SetFillStyle(3005); + norm->SetDirectory(0); + + phi->SetTitle("Normalization to #phi acceptance"); + phi->SetYTitle("#phi acceptance"); + phi->SetMarkerColor(kGreen+1); + phi->SetMarkerStyle(20); + phi->SetFillColor(kGreen+1); + phi->SetFillStyle(3004); + // phi->Scale(1. / nAcc); + phi->SetDirectory(0); + + // dNdeta->Divide(norm); + dNdeta->Divide(phi); + dNdeta->SetStats(0); + dNdeta->Scale(vtxEff, "width"); + + dNdeta_->Divide(norm); + dNdeta_->SetStats(0); + dNdeta_->Scale(vtxEff, "width"); + + output->Add(dNdeta); + output->Add(dNdeta_); + output->Add(norm); + output->Add(phi); + + return true; +} + + //____________________________________________________________________ void AliForwardMultiplicityBase::MakeRingdNdeta(const TList* input, @@ -295,46 +383,71 @@ AliForwardMultiplicityBase::MakeRingdNdeta(const TList* input, ptr++; continue; } - TH2D* copy = static_cast(h->Clone("sum")); - copy->SetDirectory(0); - thisList->Add(copy); + TH2D* sumPhi = static_cast(h->Clone("sum_phi")); + sumPhi->SetDirectory(0); + thisList->Add(sumPhi); + + TH2D* sumEta = static_cast(h->Clone("sum_eta")); + sumEta->SetDirectory(0); + thisList->Add(sumEta); - TH1D* norm =static_cast(h->ProjectionX("norm", 0, 0, "")); - for (Int_t i = 1; i <= copy->GetNbinsX(); i++) { - for (Int_t j = 1; j <= copy->GetNbinsY(); j++) { - Double_t c = copy->GetBinContent(i, j); - Double_t e = copy->GetBinError(i, j); - Double_t a = norm->GetBinContent(i); - copy->SetBinContent(i, j, a <= 0 ? 0 : c / a); - copy->SetBinError(i, j, a <= 0 ? 0 : e / a); - } - } + Int_t nY = sumEta->GetNbinsY(); + TH1D* etaCov =static_cast(h->ProjectionX("etaCov", 0, 0, "")); + TH1D* phiAcc =static_cast(h->ProjectionX("phiAcc", nY+1, nY+1, "")); - TH1D* res =static_cast(copy->ProjectionX("dndeta",1, - h->GetNbinsY(),"e")); - TH1D* proj =static_cast(h->ProjectionX("proj",1,h->GetNbinsY(),"e")); - res->SetTitle(*ptr); - res->Scale(1., "width"); - copy->Scale(1., "width"); + etaCov->SetTitle("Normalization to #eta coverage"); + etaCov->SetYTitle("#eta coverage"); + etaCov->SetMarkerColor(kBlue+1); + etaCov->SetFillColor(kBlue+1); + etaCov->SetFillStyle(3005); + etaCov->SetDirectory(0); - if(norm->GetMaximum() > 0) { - proj->Scale(1. / norm->GetMaximum(), "width"); - norm->Scale(1. / norm->GetMaximum()); + phiAcc->SetTitle("Normalization to #phi acceptance"); + phiAcc->SetYTitle("#phi acceptance"); + phiAcc->SetMarkerColor(kGreen+1); + phiAcc->SetFillColor(kGreen+1); + phiAcc->SetFillStyle(3004); + // phiAcc->Scale(1. / nAcc); + phiAcc->SetDirectory(0); + + // Double_t s = (etaCov->GetMaximum() > 0 ? 1. / etaCov->GetMaximum() : 1); + for (Int_t i = 1; i <= sumEta->GetNbinsX(); i++) { + for (Int_t j = 1; j <= nY; j++) { + Double_t c = sumEta->GetBinContent(i, j); + Double_t e = sumEta->GetBinError(i, j); + Double_t a = etaCov->GetBinContent(i); + Double_t p = phiAcc->GetBinContent(i); + // Double_t t = p; // * a + sumEta->SetBinContent(i, j, a <= 0 ? 0 : c / a); + sumEta->SetBinError( i, j, a <= 0 ? 0 : e / a); + sumPhi->SetBinContent(i, j, p <= 0 ? 0 : c / p); + sumPhi->SetBinError( i, j, p <= 0 ? 0 : e / p); + } } - - res->SetMarkerStyle(style); - norm->SetDirectory(0); - res->SetDirectory(0); - proj->SetDirectory(0); - thisList->Add(norm); - thisList->Add(res); - thisList->Add(proj); - dndetaRings->Add(res); + // etaCov->Scale(s); + // phiAcc->Scale(s); + + TH1D* resPhi =static_cast(sumPhi->ProjectionX("dndeta_phi", + 1,nY,"e")); + resPhi->SetMarkerStyle(style); + resPhi->SetDirectory(0); + resPhi->Scale(1, "width"); + + TH1D* resEta =static_cast(sumEta->ProjectionX("dndeta_eta", + 1,nY,"e")); + resEta->SetMarkerStyle(style); + resEta->SetDirectory(0); + resEta->Scale(1, "width"); + + thisList->Add(resEta); + thisList->Add(etaCov); + thisList->Add(resPhi); + thisList->Add(phiAcc); + dndetaRings->Add(resPhi); ptr++; } out->Add(dndetaRings); } - //____________________________________________________________________ void AliForwardMultiplicityBase::Print(Option_t* option) const diff --git a/PWGLF/FORWARD/analysis2/AliForwardMultiplicityBase.h b/PWGLF/FORWARD/analysis2/AliForwardMultiplicityBase.h index 59b4a3c3bdb..0dbba2850c1 100644 --- a/PWGLF/FORWARD/analysis2/AliForwardMultiplicityBase.h +++ b/PWGLF/FORWARD/analysis2/AliForwardMultiplicityBase.h @@ -295,6 +295,22 @@ protected: * */ virtual void MarkEventForStore() const; + /** + * Calculate a simple dN/deta from all accepted events + * + * @param input Input list + * @param output Output list + * @param nTr On return, number of triggers + * @param nTrVtx On return, number of trigger+vertex events + * @param nAcc On return, number of accepted events + * + * @return true on success + */ + virtual Bool_t MakeSimpledNdeta(const TList* input, + TList* output, + Double_t& nTr, + Double_t& nTrVtx, + Double_t& nAcc); /** * Make Ring @f$ dN/d\eta @f$ histogram and a stack * diff --git a/PWGLF/FORWARD/analysis2/AliForwardMultiplicityTask.cxx b/PWGLF/FORWARD/analysis2/AliForwardMultiplicityTask.cxx index 57241ed21be..4ab265271af 100644 --- a/PWGLF/FORWARD/analysis2/AliForwardMultiplicityTask.cxx +++ b/PWGLF/FORWARD/analysis2/AliForwardMultiplicityTask.cxx @@ -47,7 +47,7 @@ AliForwardMultiplicityTask::AliForwardMultiplicityTask() // // Constructor // - DGUARD(0,0,"Default construction of AliForwardMultiplicityTask"); + DGUARD(fDebug, 3,"Default CTOR of AliForwardMultiplicityTask"); } //____________________________________________________________________ @@ -73,7 +73,7 @@ AliForwardMultiplicityTask::AliForwardMultiplicityTask(const char* name) // Parameters: // name Name of task // - DGUARD(0,0,"named construction of AliForwardMultiplicityTask: %s", name); + DGUARD(fDebug, 3,"named CTOR of AliForwardMultiplicityTask: %s", name); DefineOutput(1, TList::Class()); } @@ -100,7 +100,7 @@ AliForwardMultiplicityTask::AliForwardMultiplicityTask(const AliForwardMultiplic // Parameters: // o Object to copy from // - DGUARD(0,0,"Copy construction of AliForwardMultiplicityTask"); + DGUARD(fDebug, 3,"Copy CTOR of AliForwardMultiplicityTask"); DefineOutput(1, TList::Class()); } @@ -265,11 +265,11 @@ AliForwardMultiplicityTask::UserExec(Option_t*) Bool_t lowFlux = kFALSE; UInt_t triggers = 0; UShort_t ivz = 0; - Double_t vz = 0; + TVector3 ip; Double_t cent = -1; UShort_t nClusters = 0; UInt_t found = fEventInspector.Process(esd, triggers, lowFlux, - ivz, vz, cent, nClusters); + ivz, ip, cent, nClusters); if (found & AliFMDEventInspector::kNoEvent) return; if (found & AliFMDEventInspector::kNoTriggers) return; @@ -288,7 +288,7 @@ AliForwardMultiplicityTask::UserExec(Option_t*) if (triggers & AliAODForwardMult::kPileUp) return; - fAODFMD.SetIpZ(vz); + fAODFMD.SetIpZ(ip.Z()); if (found & AliFMDEventInspector::kBadVertex) return; @@ -298,20 +298,21 @@ AliForwardMultiplicityTask::UserExec(Option_t*) // Get FMD data AliESDFMD* esdFMD = esd->GetFMDData(); // // Apply the sharing filter (or hit merging or clustering if you like) - if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD, vz)) { + if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD, ip.Z())) { AliWarning("Sharing filter failed!"); return; } // Calculate the inclusive charged particle density - if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux, cent,vz)) { + if (!fDensityCalculator.Calculate(fESDFMD, fHistos, lowFlux, cent, ip)) { // if (!fDensityCalculator.Calculate(*esdFMD, fHistos, ivz, lowFlux)) { AliWarning("Density calculator failed!"); return; } if (fEventInspector.GetCollisionSystem() == AliFMDEventInspector::kPbPb) { - if (!fEventPlaneFinder.FindEventplane(esd, fAODEP, &(fAODFMD.GetHistogram()), &fHistos)) + if (!fEventPlaneFinder.FindEventplane(esd, fAODEP, + &(fAODFMD.GetHistogram()), &fHistos)) AliWarning("Eventplane finder failed!"); } @@ -364,45 +365,14 @@ AliForwardMultiplicityTask::Terminate(Option_t*) if (GetOutputData(1)) GetOutputData(1)->Print(); return; } - - // Get our histograms from the container - TH1I* hEventsTr = 0;//static_cast(list->FindObject("nEventsTr")); - TH1I* hEventsTrVtx = 0;//static_cast(list->FindObject("nEventsTrVtx")); - TH1I* hTriggers = 0; - if (!fEventInspector.FetchHistograms(list, hEventsTr, - hEventsTrVtx, hTriggers)) { - AliError(Form("Didn't get histograms from event selector " - "(hEventsTr=%p,hEventsTrVtx=%p)", - hEventsTr, hEventsTrVtx)); - list->ls(); - return; - } - - TH2D* hData = static_cast(list->FindObject("d2Ndetadphi")); - if (!hData) { - AliError(Form("Couldn't get our summed histogram from output " - "list %s (d2Ndetadphi=%p)", list->GetName(), hData)); - list->ls(); - return; - } - - // TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e"); - TH1D* dNdeta = hData->ProjectionX("dNdeta", 1, -1, "e"); - TH1D* norm = hData->ProjectionX("norm", 0, 0, ""); - dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions"); - dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}"); - dNdeta->Divide(norm); - dNdeta->SetStats(0); - dNdeta->Scale(Double_t(hEventsTrVtx->GetEntries())/hEventsTr->GetEntries(), - "width"); - list->Add(dNdeta); - list->Add(norm); + Double_t nTr = 0, nTrVtx = 0, nAcc = 0; + MakeSimpledNdeta(list, list, nTr, nTrVtx, nAcc); MakeRingdNdeta(list, "ringSums", list, "ringResults"); - fSharingFilter.ScaleHistograms(list,Int_t(hEventsTr->Integral())); - fDensityCalculator.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral())); - fCorrections.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral())); + fSharingFilter.ScaleHistograms(list,Int_t(nTr)); + fDensityCalculator.ScaleHistograms(list,Int_t(nTrVtx)); + fCorrections.ScaleHistograms(list,Int_t(nTrVtx)); } // diff --git a/PWGLF/FORWARD/analysis2/AliForwardQATask.cxx b/PWGLF/FORWARD/analysis2/AliForwardQATask.cxx index 8112079c7ed..394e81a830d 100644 --- a/PWGLF/FORWARD/analysis2/AliForwardQATask.cxx +++ b/PWGLF/FORWARD/analysis2/AliForwardQATask.cxx @@ -352,11 +352,11 @@ AliForwardQATask::UserExec(Option_t*) Bool_t lowFlux = kFALSE; UInt_t triggers = 0; UShort_t ivz = 0; - Double_t vz = 0; + TVector3 ip; Double_t cent = -1; UShort_t nClusters = 0; UInt_t found = fEventInspector.Process(esd, triggers, lowFlux, - ivz, vz, cent, nClusters); + ivz, ip, cent, nClusters); Bool_t ok = true; if (found & AliFMDEventInspector::kNoEvent) ok = false; @@ -371,7 +371,8 @@ AliForwardQATask::UserExec(Option_t*) AliFMDEventInspector::CodeString(found)); return; } - DMSG(fDebug,2,"Event triggers: %s", AliAODForwardMult::GetTriggerString(triggers)); + DMSG(fDebug,2,"Event triggers: %s", + AliAODForwardMult::GetTriggerString(triggers)); // We we do not want to use low flux specific code, we disable it here. if (!fEnableLowFlux) lowFlux = false; @@ -387,7 +388,7 @@ AliForwardQATask::UserExec(Option_t*) } // // Apply the sharing filter (or hit merging or clustering if you like) - if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD, vz)) { + if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD, ip.Z())) { AliWarning("Sharing filter failed!"); return; } @@ -426,12 +427,16 @@ AliForwardQATask::Terminate(Option_t*) // Get our histograms from the container TH1I* hEventsTr = 0;//static_cast(list->FindObject("nEventsTr")); TH1I* hEventsTrVtx = 0;//static_cast(list->FindObject("nEventsTrVtx")); + TH1I* hEventsAcc = 0; TH1I* hTriggers = 0; - if (!fEventInspector.FetchHistograms(list, hEventsTr, - hEventsTrVtx, hTriggers)) { + if (!fEventInspector.FetchHistograms(list, + hEventsTr, + hEventsTrVtx, + hEventsAcc, + hTriggers)) { AliError(Form("Didn't get histograms from event selector " - "(hEventsTr=%p,hEventsTrVtx=%p)", - hEventsTr, hEventsTrVtx)); + "(hEventsTr=%p,hEventsTrVtx=%p,hEventsAcc=%p)", + hEventsTr, hEventsTrVtx,hEventsAcc)); return; } diff --git a/PWGLF/FORWARD/analysis2/AliForwardUtil.cxx b/PWGLF/FORWARD/analysis2/AliForwardUtil.cxx index f7c50a91726..85e229298fa 100644 --- a/PWGLF/FORWARD/analysis2/AliForwardUtil.cxx +++ b/PWGLF/FORWARD/analysis2/AliForwardUtil.cxx @@ -319,31 +319,40 @@ void AliForwardUtil::GetParameter(TObject* o, Bool_t& value) } //_____________________________________________________________________ -Double_t AliForwardUtil::GetEtaFromStrip(UShort_t det, Char_t ring, UShort_t sec, UShort_t strip, Double_t zvtx) +Double_t AliForwardUtil::GetStripR(Char_t ring, UShort_t strip) { - //Calculate eta from strip with vertex (redundant with AliESDFMD::Eta but support displaced vertices) - //Get max R of ring - Double_t maxR = 0; - Double_t minR = 0; - Bool_t inner = false; - switch (ring) { - case 'i': case 'I': maxR = 17.2; minR = 4.5213; inner = true; break; - case 'o': case 'O': maxR = 28.0; minR = 15.4; inner = false; break; - default: - return -99999; - } + const Double_t iMinR = 4.5213; + const Double_t iMaxR = 17.2; + const Double_t oMinR = 15.4; + const Double_t oMaxR = 28.0; + + Double_t minR = (ring == 'I' || ring == 'i') ? iMinR : oMinR; + Double_t maxR = (ring == 'I' || ring == 'i') ? iMaxR : oMaxR; + Double_t nStrips = (ring == 'I' || ring == 'i') ? 512 : 256; + Double_t rad = maxR - minR; + Double_t segment = rad / nStrips; + Double_t r = minR + segment*strip; - Double_t rad = maxR- minR; - Double_t nStrips = (ring == 'I' ? 512 : 256); - Double_t segment = rad / nStrips; - Double_t r = minR + segment*strip; - Int_t hybrid = sec / 2; + return r; +} + +//_____________________________________________________________________ +Double_t AliForwardUtil::GetEtaFromStrip(UShort_t det, Char_t ring, + UShort_t sec, UShort_t strip, + Double_t zvtx) +{ + // Calculate eta from strip with vertex (redundant with + // AliESDFMD::Eta but support displaced vertices) - Double_t z = 0; + //Get max R of ring + Double_t r = GetStripR(ring, strip); + Int_t hybrid = sec / 2; + Bool_t inner = (ring == 'I' || ring == 'i'); + Double_t z = 0; switch (det) { - case 1: z = 320.266; break; - case 2: z = (inner ? 83.666 : 74.966); break; + case 1: z = 320.266; break; + case 2: z = (inner ? 83.666 : 74.966); break; case 3: z = (inner ? -63.066 : -74.966); break; default: return -999999; } @@ -355,6 +364,27 @@ Double_t AliForwardUtil::GetEtaFromStrip(UShort_t det, Char_t ring, UShort_t sec return eta; } +//_____________________________________________________________________ +Double_t AliForwardUtil::GetPhiFromStrip(Char_t ring, UShort_t strip, + Double_t phi, + Double_t xvtx, Double_t yvtx) +{ + // Calculate eta from strip with vertex (redundant with + // AliESDFMD::Eta but support displaced vertices) + + // Unknown x,y -> no change + if (yvtx > 999 || xvtx > 999) return phi; + + //Get max R of ring + Double_t r = GetStripR(ring, strip); + Double_t amp = TMath::Sqrt(xvtx*xvtx+yvtx*yvtx) / r; + Double_t pha = (TMath::Abs(yvtx) < 1e-12 ? 0 : TMath::ATan2(xvtx, yvtx)); + Double_t cha = amp * TMath::Cos(phi+pha); + phi += cha; + if (phi < 0) phi += TMath::TwoPi(); + if (phi > TMath::TwoPi()) phi -= TMath::TwoPi(); + return phi; +} //==================================================================== Int_t AliForwardUtil::fgConvolutionSteps = 100; diff --git a/PWGLF/FORWARD/analysis2/AliForwardUtil.h b/PWGLF/FORWARD/analysis2/AliForwardUtil.h index 86056dc3efd..86bcbcb22cf 100644 --- a/PWGLF/FORWARD/analysis2/AliForwardUtil.h +++ b/PWGLF/FORWARD/analysis2/AliForwardUtil.h @@ -153,6 +153,15 @@ public: * @return Short integer value of magnetic field in kG */ static Short_t ParseMagneticField(Float_t field); + /** + * Get the radius of a strip. + * + * @param ring Ring identifier 'I' or 'O' + * @param strip Strip number + * + * @return Radial distance from beam of the strip + */ + static Double_t GetStripR(Char_t ring, UShort_t strip); /** * Get eta from strip * @@ -160,7 +169,21 @@ public: * * @return eta */ - static Double_t GetEtaFromStrip(UShort_t det, Char_t ring, UShort_t sec, UShort_t strip, Double_t zvtx) ; + static Double_t GetEtaFromStrip(UShort_t det, Char_t ring, + UShort_t sec, UShort_t strip, Double_t zvtx); + /** + * Get the azimuthal angle of a strip + * + * @param ring Ring identifier 'I' or 'O' + * @param strip Strip number + * @param phi Straight forward strip phi + * @param xvtx Ip X coordinate + * @param yvtx Ip Y coordinate + * + * @return The phi angle correctef for (X,Y) off set. + */ + static Double_t GetPhiFromStrip(Char_t ring, UShort_t strip, + Double_t phi, Double_t xvtx, Double_t yvtx); /** * Get a string representation of the magnetic field * diff --git a/PWGLF/FORWARD/analysis2/AliPoissonCalculator.cxx b/PWGLF/FORWARD/analysis2/AliPoissonCalculator.cxx index 9d74048b191..54e59211205 100644 --- a/PWGLF/FORWARD/analysis2/AliPoissonCalculator.cxx +++ b/PWGLF/FORWARD/analysis2/AliPoissonCalculator.cxx @@ -202,7 +202,7 @@ void AliPoissonCalculator::MakeOutput() { fMean->SetDirectory(0); fOcc = new TH1D("occupancy", "Occupancy = #int_{1}^{#infty}dN P(N)", - 1000, 0, 100); + 101, -.5, 100.5); fOcc->SetXTitle("#int_{1}^{#infty}dN P(N) [%]"); fOcc->SetYTitle("Events"); fOcc->SetFillColor(kBlue+1);