//Float_t binLimitsVtx[] = {-20,-15,-10,-6,-3,0,3,6,10,15,20};
//Float_t binLimitsVtx[] = {-20,-19,-18,-17,-16,-15,-14,-13,-12,-11,-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20};
Float_t binLimitsVtx[] = {-20,-15,-10,-8,-6,-4,-2,0,2,4,6,8,10,15,20};
- Float_t binLimitsEta[] = {-2,-1.9,-1.8,-1.7,-1.6,-1.5,-1.4,-1.3,-1.2,-1.1,-1.0,-0.9,-0.8,-0.7,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,0.0,
- 0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0};
+ Float_t binLimitsEta[] = {-2,-1.9,-1.8,-1.7,-1.6,-1.5,-1.4,-1.3,-1.2,-1.1,-1.0,-0.9,-0.8,-0.7,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0};
+// Float_t binLimitsEta[] = { -7.0, -6.9, -6.8, -6.7, -6.6, -6.5, -6.4, -6.3, -6.2, -6.1, -6.0, -5.9, -5.8, -5.7, -5.6, -5.5, -5.4, -5.3, -5.2, -5.1, -5.0, -4.9, -4.8, -4.7, -4.6, -4.5, -4.4, -4.3, -4.2, -4.1, -4.0, -3.9, -3.8, -3.7, -3.6, -3.5, -3.4, -3.3, -3.2, -3.1, -3.0, -2.9, -2.8, -2.7, -2.6, -2.5, -2.4, -2.3, -2.2, -2.1, -2.0, -1.9, -1.8, -1.7, -1.6, -1.5, -1.4, -1.3, -1.2, -1.1, -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, -0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3.0, 3.1, 3.2, 3.3, 3.4, 3.5, 3.6, 3.7, 3.8, 3.9, 4.0, 4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9, 5.0, 5.1, 5.2, 5.3, 5.4, 5.5, 5.6, 5.7, 5.8, 5.9, 6.0, 6.1, 6.2, 6.3, 6.4, 6.5, 6.6, 6.7, 6.8, 6.9, 7.0 };
TH3F* dummyBinning = new TH3F("","",14, binLimitsVtx, 40, binLimitsEta , 28, binLimitsPt);
/* $Id$ */
#include <TList.h>
+#include <TParticle.h>
#include <AliPWG0depHelper.h>
#include <AliHeader.h>
+#include <AliStack.h>
+#include <AliLog.h>
#include <AliGenEventHeader.h>
#include <AliGenPythiaEventHeader.h>
ClassImp(AliPWG0depHelper)
//____________________________________________________________________
-const Int_t AliPWG0depHelper::GetPythiaEventProcessType(AliHeader* aHeader, Bool_t adebug) {
+Int_t AliPWG0depHelper::GetPythiaEventProcessType(AliHeader* aHeader, Bool_t adebug) {
//
// get the process type of the event.
//
return pythiaGenHeader->ProcessType();
}
+
+//____________________________________________________________________
+TParticle* AliPWG0depHelper::FindPrimaryMother(AliStack* stack, Int_t label)
+{
+ //
+ // Finds the first mother among the primary particles of the particle identified by <label>,
+ // i.e. the primary that "caused" this particle
+ //
+
+ Int_t nPrim = stack->GetNprimary();
+ TParticle* mother = stack->Particle(label);
+
+ if (!mother)
+ {
+ AliDebugGeneral("FindPrimaryMother", AliLog::kError, Form("UNEXPECTED: Could not retrieve particle with label %d.", label));
+ return 0;
+ }
+
+ while (label >= nPrim)
+ {
+ //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0));
+
+ if (mother->GetMother(0) == -1)
+ {
+ AliDebugGeneral("FindPrimaryMother", AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label));
+ mother = 0;
+ break;
+ }
+
+ label = mother->GetMother(0);
+
+ mother = stack->Particle(label);
+ if (!mother)
+ {
+ AliDebugGeneral("FindPrimaryMother", AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (find mother loop).", label));
+ break;
+ }
+ }
+
+ return mother;
+}
// static helper functions that depend on more than ESD
class AliHeader;
+class TParticle;
+class AliStack;
class AliPWG0depHelper : public TObject
{
public:
- static const Int_t GetPythiaEventProcessType(AliHeader* aHeader, Bool_t adebug = kFALSE);
+ static Int_t GetPythiaEventProcessType(AliHeader* aHeader, Bool_t adebug = kFALSE);
+ static TParticle* FindPrimaryMother(AliStack* stack, Int_t label);
protected:
ClassDef(AliPWG0depHelper, 0)
#pragma link C++ class AliROCESDAnalysisSelector+;
#pragma link C++ class AliROCRawAnalysisSelector+;
#pragma link C++ class AliROCClusterAnalysisSelector+;
+#pragma link C++ class AliHighMultiplicitySelector+;
#endif
#include <TF1.h>
#include <TMath.h>
#include <TCollection.h>
+#include <TLegend.h>
+#include <TLine.h>
ClassImp(AliMultiplicityCorrection)
-const Int_t AliMultiplicityCorrection::fgMaxParams = 200;
+const Int_t AliMultiplicityCorrection::fgMaxInput = 250; // bins in measured histogram
+const Int_t AliMultiplicityCorrection::fgMaxParams = 250; // bins in unfolded histogram = number of fit params
+
TH1* AliMultiplicityCorrection::fCurrentESD = 0;
TH1* AliMultiplicityCorrection::fCurrentCorrelation = 0;
TH1* AliMultiplicityCorrection::fCurrentEfficiency = 0;
-TMatrixF* AliMultiplicityCorrection::fCorrelationMatrix = 0;
-TMatrixF* AliMultiplicityCorrection::fCorrelationCovarianceMatrix = 0;
-TVectorF* AliMultiplicityCorrection::fCurrentESDVector = 0;
+TMatrixD* AliMultiplicityCorrection::fCorrelationMatrix = 0;
+TMatrixD* AliMultiplicityCorrection::fCorrelationCovarianceMatrix = 0;
+TVectorD* AliMultiplicityCorrection::fCurrentESDVector = 0;
+TVectorD* AliMultiplicityCorrection::fEntropyAPriori = 0;
AliMultiplicityCorrection::RegularizationType AliMultiplicityCorrection::fRegularizationType = AliMultiplicityCorrection::kPol1;
-Float_t AliMultiplicityCorrection::fRegularizationWeight = 1e4;
+Float_t AliMultiplicityCorrection::fRegularizationWeight = 5000;
TF1* AliMultiplicityCorrection::fNBD = 0;
//____________________________________________________________________
#define VTXBINNING 10, binLimitsVtx
#define NBINNING fgMaxParams, binLimitsN*/
- #define NBINNING 251, -0.5, 250.5
- #define VTXBINNING 10, -10, 10
+ #define NBINNING 501, -0.5, 500.5
+ #define VTXBINNING 1, -10, 10
for (Int_t i = 0; i < kESDHists; ++i)
fMultiplicityESD[i] = new TH2F(Form("fMultiplicityESD%d", i), "fMultiplicityESD;vtx-z;Ntracks;Count", VTXBINNING, NBINNING);
}
//____________________________________________________________________
-Double_t AliMultiplicityCorrection::RegularizationPol0(Double_t *params)
+Double_t AliMultiplicityCorrection::RegularizationPol0(TVectorD& params)
{
// homogenity term for minuit fitting
// pure function of the parameters
Double_t chi2 = 0;
- for (Int_t i=1; i<fgMaxParams; ++i)
+ // ignore the first bin here. on purpose...
+ for (Int_t i=2; i<fgMaxParams; ++i)
{
Double_t right = params[i];
Double_t left = params[i-1];
- Double_t diff = (right - left);
- chi2 += diff * diff;
+ if (right != 0)
+ {
+ Double_t diff = 1 - left / right;
+ chi2 += diff * diff;
+ }
}
- return chi2 * 100;
+ return chi2 / 100.0;
}
//____________________________________________________________________
-Double_t AliMultiplicityCorrection::RegularizationPol1(Double_t *params)
+Double_t AliMultiplicityCorrection::RegularizationPol1(TVectorD& params)
{
// homogenity term for minuit fitting
// pure function of the parameters
Double_t chi2 = 0;
- for (Int_t i=2; i<fgMaxParams; ++i)
+ // ignore the first bin here. on purpose...
+ for (Int_t i=2+1; i<fgMaxParams; ++i)
{
if (params[i-1] == 0)
continue;
}
//____________________________________________________________________
-Double_t AliMultiplicityCorrection::RegularizationTest(Double_t *params)
+Double_t AliMultiplicityCorrection::RegularizationTest(TVectorD& params)
{
// homogenity term for minuit fitting
// pure function of the parameters
}
//____________________________________________________________________
-Double_t AliMultiplicityCorrection::RegularizationTotalCurvature(Double_t *params)
+Double_t AliMultiplicityCorrection::RegularizationTotalCurvature(TVectorD& params)
{
// homogenity term for minuit fitting
// pure function of the parameters
Double_t chi2 = 0;
- for (Int_t i=2; i<fgMaxParams; ++i)
+ // ignore the first bin here. on purpose...
+ for (Int_t i=3; i<fgMaxParams; ++i)
{
Double_t right = params[i];
Double_t middle = params[i-1];
Double_t der1 = (right - middle);
Double_t der2 = (middle - left);
- Double_t secDer = (der1 - der2);
+ Double_t diff = (der1 - der2);
- // square and weight with the bin width
- chi2 += secDer * secDer;
+ chi2 += diff * diff;
}
- return chi2 * 100;
+ return chi2 * 1e4;
}
//____________________________________________________________________
-Double_t AliMultiplicityCorrection::RegularizationEntropy(Double_t *params)
+Double_t AliMultiplicityCorrection::RegularizationEntropy(TVectorD& params)
{
// homogenity term for minuit fitting
// pure function of the parameters
// The method of reduced cross-entropy (M. Schmelling 1993)
Double_t paramSum = 0;
- for (Int_t i=0; i<fgMaxParams; ++i)
+ // ignore the first bin here. on purpose...
+ for (Int_t i=1; i<fgMaxParams; ++i)
paramSum += params[i];
Double_t chi2 = 0;
- for (Int_t i=0; i<fgMaxParams; ++i)
+ for (Int_t i=1; i<fgMaxParams; ++i)
{
Double_t tmp = params[i] / paramSum;
- if (tmp > 0)
- chi2 += tmp * log(tmp);
+ if (tmp > 0 && (*fEntropyAPriori)[i] > 0)
+ chi2 += tmp * log(tmp / (*fEntropyAPriori)[i]);
}
- return chi2;
+ return 10.0 + chi2;
}
//____________________________________________________________________
// does: nbd
// func = new TF1("nbd", "[0] * TMath::Binomial([2]+TMath::Nint(x)-1, [2]-1) * pow([1] / ([1]+[2]), TMath::Nint(x)) * pow(1 + [1]/[2], -[2])", 0, 50);
// func->SetParNames("scaling", "averagen", "k");
+ // func->SetParLimits(0, 0.001, fCurrentESD->GetMaximum() * 1000);
+ // func->SetParLimits(1, 0.001, 1000);
+ // func->SetParLimits(2, 0.001, 1000);
//
fNBD->SetParameters(params[0], params[1], params[2]);
// does: (m - Ad)W(m - Ad) where m = measured, A correlation matrix, d = guess, W = covariance matrix
//
- static Int_t callCount = 0;
-
- Double_t chi2FromFit = 0;
-
// d
- TVectorF paramsVector(fgMaxParams);
+ TVectorD paramsVector(fgMaxParams);
for (Int_t i=0; i<fgMaxParams; ++i)
- paramsVector[i] = params[i];
+ paramsVector[i] = params[i] * params[i];
+
+ // calculate penalty factor
+ Double_t penaltyVal = 0;
+ switch (fRegularizationType)
+ {
+ case kNone: break;
+ case kPol0: penaltyVal = RegularizationPol0(paramsVector); break;
+ case kPol1: penaltyVal = RegularizationPol1(paramsVector); break;
+ case kCurvature: penaltyVal = RegularizationTotalCurvature(paramsVector); break;
+ case kEntropy: penaltyVal = RegularizationEntropy(paramsVector); break;
+ case kTest: penaltyVal = RegularizationTest(paramsVector); break;
+ }
+
+ //if (penaltyVal > 1e10)
+ // paramsVector2.Print();
+
+ penaltyVal *= fRegularizationWeight;
// Ad
- paramsVector = (*fCorrelationMatrix) * paramsVector;
+ TVectorD measGuessVector(fgMaxInput);
+ measGuessVector = (*fCorrelationMatrix) * paramsVector;
// Ad - m
- paramsVector -= (*fCurrentESDVector);
+ measGuessVector -= (*fCurrentESDVector);
- TVectorF copy(paramsVector);
+ TVectorD copy(measGuessVector);
// (Ad - m) W
- paramsVector *= (*fCorrelationCovarianceMatrix);
+ measGuessVector *= (*fCorrelationCovarianceMatrix);
- // (Ad - m) W (Ad - m)
- chi2FromFit = paramsVector * copy;
-
- Double_t penaltyVal = 0;
+ //measGuessVector.Print();
- switch (fRegularizationType)
- {
- case kNone: break;
- case kPol0: penaltyVal = RegularizationPol0(params); break;
- case kPol1: penaltyVal = RegularizationPol1(params); break;
- case kCurvature: penaltyVal = RegularizationTotalCurvature(params); break;
- case kEntropy: penaltyVal = RegularizationEntropy(params); break;
- case kTest: penaltyVal = RegularizationTest(params); break;
- }
-
- penaltyVal *= fRegularizationWeight;
+ // (Ad - m) W (Ad - m)
+ Double_t chi2FromFit = measGuessVector * copy * 1e6;
chi2 = chi2FromFit + penaltyVal;
- if ((callCount++ % 1000) == 0)
- printf("%f %f --> %f\n", chi2FromFit, penaltyVal, chi2);
+ static Int_t callCount = 0;
+ if ((callCount++ % 10000) == 0)
+ printf("%d) %f %f --> %f\n", callCount, chi2FromFit, penaltyVal, chi2);
}
//____________________________________________________________________
-void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType)
+void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t createBigBin)
{
//
// fills fCurrentESD, fCurrentCorrelation
}
fCurrentCorrelation = hist->Project3D("zy");
+ //((TH2*) fCurrentCorrelation)->Rebin2D(2, 1);
+ //fMultiplicityESDCorrected[correlationID]->Rebin(2);
fCurrentCorrelation->Sumw2();
+ if (createBigBin)
+ {
+ Int_t maxBin = 0;
+ for (Int_t i=1; i<=fCurrentESD->GetNbinsX(); ++i)
+ {
+ if (fCurrentESD->GetBinContent(i) <= 5)
+ {
+ maxBin = i;
+ break;
+ }
+ }
+
+ if (maxBin > 0)
+ {
+ TCanvas* canvas = new TCanvas("StatSolution", "StatSolution", 1000, 800);
+ canvas->Divide(2, 2);
+
+ canvas->cd(1);
+ fCurrentESD->DrawCopy();
+ gPad->SetLogy();
+
+ canvas->cd(2);
+ fCurrentCorrelation->DrawCopy("COLZ");
+
+ printf("Bin limit in measured spectrum is %d.\n", maxBin);
+ fCurrentESD->SetBinContent(maxBin, fCurrentESD->Integral(maxBin, fCurrentESD->GetNbinsX()));
+ for (Int_t i=maxBin+1; i<=fCurrentESD->GetNbinsX(); ++i)
+ {
+ fCurrentESD->SetBinContent(i, 0);
+ fCurrentESD->SetBinError(i, 0);
+ }
+ // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
+ fCurrentESD->SetBinError(maxBin, TMath::Sqrt(fCurrentESD->GetBinContent(maxBin)));
+
+ printf("This bin has now %f +- %f entries\n", fCurrentESD->GetBinContent(maxBin), fCurrentESD->GetBinError(maxBin));
+
+ for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
+ {
+ fCurrentCorrelation->SetBinContent(i, maxBin, fCurrentCorrelation->Integral(i, i, maxBin, fCurrentCorrelation->GetNbinsY()));
+ // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
+ fCurrentCorrelation->SetBinError(i, maxBin, TMath::Sqrt(fCurrentCorrelation->GetBinContent(i, maxBin)));
+
+ for (Int_t j=maxBin+1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
+ {
+ fCurrentCorrelation->SetBinContent(i, j, 0);
+ fCurrentCorrelation->SetBinError(i, j, 0);
+ }
+ }
+
+ printf("Adjusted correlation matrix!\n");
+
+ canvas->cd(3);
+ fCurrentESD->DrawCopy();
+ gPad->SetLogy();
+
+ canvas->cd(4);
+ fCurrentCorrelation->DrawCopy("COLZ");
+ }
+ }
+
+ //normalize ESD
+ fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
+
fCurrentEfficiency = fMultiplicityVtx[inputRange]->ProjectionY("CurrentEfficiency");
TH2* divisor = 0;
switch (eventType)
case kINEL: divisor = fMultiplicityINEL[inputRange]; break;
}
fCurrentEfficiency->Divide(divisor->ProjectionY());
+ //fCurrentEfficiency->Rebin(2);
+ //fCurrentEfficiency->Scale(0.5);
}
//____________________________________________________________________
-void AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check, TH1* inputDist)
+void AliMultiplicityCorrection::SetRegularizationParameters(RegularizationType type, Float_t weight)
+{
+ fRegularizationType = type;
+ fRegularizationWeight = weight;
+
+ printf("AliMultiplicityCorrection::SetRegularizationParameters --> Regularization set to %d with weight %f\n", (Int_t) type, weight);
+}
+
+//____________________________________________________________________
+Int_t AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check, TH1* inputDist)
{
//
// correct spectrum using minuit chi2 method
Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
- SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
+ SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE); // TODO FAKE kTRUE
- fCorrelationMatrix = new TMatrixF(fgMaxParams, fgMaxParams);
- fCorrelationCovarianceMatrix = new TMatrixF(fgMaxParams, fgMaxParams);
- fCurrentESDVector = new TVectorF(fgMaxParams);
+ fCorrelationMatrix = new TMatrixD(fgMaxInput, fgMaxParams);
+ fCorrelationCovarianceMatrix = new TMatrixD(fgMaxInput, fgMaxInput);
+ fCurrentESDVector = new TVectorD(fgMaxInput);
+ fEntropyAPriori = new TVectorD(fgMaxParams);
+
+ /*new TCanvas; fCurrentESD->DrawCopy();
+ fCurrentESD = ((TH2*) fCurrentCorrelation)->ProjectionY("check-proj2");
+ fCurrentESD->Sumw2();
+ fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
+ fCurrentESD->SetLineColor(2);
+ fCurrentESD->DrawCopy("SAME");*/
// normalize correction for given nPart
for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
- if (i <= fgMaxParams && j <= fgMaxParams)
+ if (i <= fgMaxParams && j <= fgMaxInput)
(*fCorrelationMatrix)(j-1, i-1) = fCurrentCorrelation->GetBinContent(i, j);
}
//printf("MPV for Ntrue = %f is %f\n", fCurrentCorrelation->GetXaxis()->GetBinCenter(i), fCurrentCorrelation->GetYaxis()->GetBinCenter(maxBin));
}
- // small hack to get around charge conservation for full phase space ;-)
- /*if (fullPhaseSpace)
- {
- for (Int_t i=2; i<=50; i+=2)
- for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
- fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i-1, j));
- }*/
-
// Initialize TMinuit via generic fitter interface
TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgMaxParams);
+ Double_t arglist[100];
+ arglist[0] = 0;
+ minuit->ExecuteCommand("SET PRINT", arglist, 1);
minuit->SetFCN(MinuitFitFunction);
- static Double_t results[fgMaxParams];
+ for (Int_t i=0; i<fgMaxParams; i++)
+ (*fEntropyAPriori)[i] = 1;
if (inputDist)
{
printf("Using different starting conditions...\n");
new TCanvas;
inputDist->DrawCopy();
+
+ inputDist->Scale(1.0 / inputDist->Integral());
+ for (Int_t i=0; i<fgMaxParams; i++)
+ if (inputDist->GetBinContent(i+1) > 0)
+ (*fEntropyAPriori)[i] = inputDist->GetBinContent(i+1);
}
else
inputDist = fCurrentESD;
- // normalize ESD
- fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
- inputDist->Scale(1.0 / inputDist->Integral());
- Float_t minStartValue = 1e-3;
+ //Float_t minStartValue = 0; //1e-3;
+ //new TCanvas; fMultiplicityVtx[mcTarget]->Draw("COLZ");
TH1* proj = fMultiplicityVtx[mcTarget]->ProjectionY("check-proj");
+ //proj->Rebin(2);
proj->Scale(1.0 / proj->Integral());
+
+ Double_t results[fgMaxParams];
for (Int_t i=0; i<fgMaxParams; ++i)
{
results[i] = inputDist->GetBinContent(i+1);
+
if (check)
results[i] = proj->GetBinContent(i+1);
- printf("%f %f %f\n", inputDist->GetBinContent(i+1), TMath::Sqrt(inputDist->GetBinContent(i+1)), inputDist->GetBinError(i+1));
+ // minimum value
+ results[i] = TMath::Max(results[i], 1e-3);
+
+ Float_t stepSize = 0.1;
+
+ // minuit sees squared values to prevent it from going negative...
+ results[i] = TMath::Sqrt(results[i]);
+ //stepSize /= results[i]; // keep relative step size
+
+ minuit->SetParameter(i, Form("param%d", i), results[i], stepSize, 0, 0);
+ }
+ // bin 0 is filled with value from bin 1 (otherwise it's 0)
+ //minuit->SetParameter(0, "param0", (results[1] > minStartValue) ? results[1] : minStartValue, 0.1, 0, 0);
+ //results[0] = 0;
+ //minuit->SetParameter(0, "param0", 0, 0, 0, 0);
+
+ for (Int_t i=0; i<fgMaxInput; ++i)
+ {
+ //printf("%f %f %f\n", inputDist->GetBinContent(i+1), TMath::Sqrt(inputDist->GetBinContent(i+1)), inputDist->GetBinError(i+1));
(*fCurrentESDVector)[i] = fCurrentESD->GetBinContent(i+1);
if (fCurrentESD->GetBinError(i+1) > 0)
- (*fCorrelationCovarianceMatrix)(i, i) = 1.0 / fCurrentESD->GetBinError(i+1);
+ (*fCorrelationCovarianceMatrix)(i, i) = (Double_t) 1e-6 / fCurrentESD->GetBinError(i+1) / fCurrentESD->GetBinError(i+1);
+
+ if ((*fCorrelationCovarianceMatrix)(i, i) > 1e7)
+ (*fCorrelationCovarianceMatrix)(i, i) = 0;
- minuit->SetParameter(i, Form("param%d", i), ((results[i] > minStartValue) ? results[i] : minStartValue), 0.1, 0, 1);
+ //printf("%d --> %e\n", i, (*fCorrelationCovarianceMatrix)(i, i));
}
- // bin 0 is filled with value from bin 1 (otherwise it's 0)
- minuit->SetParameter(0, "param0", (results[1] > minStartValue) ? results[1] : minStartValue, 0.1, 0, 1);
Int_t dummy;
Double_t chi2;
printf("Chi2 of initial parameters is = %f\n", chi2);
if (check)
- return;
-
- Double_t arglist[100];
- arglist[0] = 0;
- minuit->ExecuteCommand("SET PRINT", arglist, 1);
- minuit->ExecuteCommand("MIGRAD", arglist, 1);
+ return -1;
+
+ // first param is number of iterations, second is precision....
+ arglist[0] = 1e6;
+ //arglist[1] = 1e-5;
+ //minuit->ExecuteCommand("SCAN", arglist, 0);
+ Int_t status = minuit->ExecuteCommand("MIGRAD", arglist, 1);
+ printf("MINUIT status is %d\n", status);
+ //minuit->ExecuteCommand("MIGRAD", arglist, 1);
+ //minuit->ExecuteCommand("MIGRAD", arglist, 1);
//printf("!!!!!!!!!!!!!! MIGRAD finished: Starting MINOS !!!!!!!!!!!!!!");
//minuit->ExecuteCommand("MINOS", arglist, 0);
{
if (fCurrentEfficiency->GetBinContent(i+1) > 0)
{
- fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, minuit->GetParameter(i) / fCurrentEfficiency->GetBinContent(i+1));
- fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, minuit->GetParError(i) / fCurrentEfficiency->GetBinContent(i+1));
+ fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, minuit->GetParameter(i) * minuit->GetParameter(i) / fCurrentEfficiency->GetBinContent(i+1));
+ // error is : (relError) * (value) = (minuit->GetParError(i) / minuit->GetParameter(i)) * (minuit->GetParameter(i) * minuit->GetParameter(i))
+ fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, minuit->GetParError(i) * minuit->GetParameter(i) / fCurrentEfficiency->GetBinContent(i+1));
}
else
{
fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, 0);
}
}
+
+ return status;
}
//____________________________________________________________________
Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
- SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx);
+ SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx, kFALSE);
- fCorrelationMatrix = new TMatrixF(fgMaxParams, fgMaxParams);
- fCorrelationCovarianceMatrix = new TMatrixF(fgMaxParams, fgMaxParams);
- fCurrentESDVector = new TVectorF(fgMaxParams);
+ fCorrelationMatrix = new TMatrixD(fgMaxParams, fgMaxParams);
+ fCorrelationCovarianceMatrix = new TMatrixD(fgMaxParams, fgMaxParams);
+ fCurrentESDVector = new TVectorD(fgMaxParams);
// normalize correction for given nPart
for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
Double_t arglist[100];
arglist[0] = 0;
minuit->ExecuteCommand("SET PRINT", arglist, 1);
- minuit->ExecuteCommand("MIGRAD", arglist, 1);
+ minuit->ExecuteCommand("MIGRAD", arglist, 0);
//minuit->ExecuteCommand("MINOS", arglist, 0);
fNBD->SetParameters(minuit->GetParameter(0), minuit->GetParameter(1), minuit->GetParameter(2));
}
//____________________________________________________________________
-void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRange, Bool_t fullPhaseSpace, Bool_t normalizeESD, TH1* mcHist)
+void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRange, Bool_t fullPhaseSpace, Bool_t normalizeESD, TH1* mcHist, Bool_t simple)
{
+ //mcHist->Rebin(2);
+
Int_t esdCorrId = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
TString tmpStr;
tmpStr.Form("%s_DrawComparison_%d", name, esdCorrId);
+ if (fMultiplicityESDCorrected[esdCorrId]->Integral() == 0)
+ {
+ printf("ERROR. Unfolded histogram is empty\n");
+ return;
+ }
+
+ //regain measured distribution used for unfolding, because the bins at high mult. were modified in SetupCurrentHists
+ fCurrentESD = fMultiplicityESD[esdCorrId]->ProjectionY();
+ fCurrentESD->Sumw2();
+ fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
+
+ // normalize unfolded result to 1
+ fMultiplicityESDCorrected[esdCorrId]->Scale(1.0 / fMultiplicityESDCorrected[esdCorrId]->Integral());
+
+ //fCurrentESD->Scale(mcHist->Integral(2, 200));
+
+ //new TCanvas;
+ /*TH1* ratio = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("ratio");
+ ratio->Divide(mcHist);
+ ratio->Draw("HIST");
+ ratio->Fit("pol0", "W0", "", 20, 120);
+ Float_t scalingFactor = ratio->GetFunction("pol0")->GetParameter(0);
+ delete ratio;
+ mcHist->Scale(scalingFactor);*/
+
+ mcHist->Scale(1.0 / mcHist->Integral());
+
// calculate residual
// for that we convolute the response matrix with the gathered result
// if normalizeESD is not set, the histogram is already normalized, this needs to be passed to CalculateMultiplicityESD
TH1* tmpESDEfficiencyRecorrected = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("tmpESDEfficiencyRecorrected");
tmpESDEfficiencyRecorrected->Multiply(fCurrentEfficiency);
- TH2* convoluted = CalculateMultiplicityESD(tmpESDEfficiencyRecorrected, esdCorrId, !normalizeESD);
- TH1* proj2 = convoluted->ProjectionY("proj2", -1, -1, "e");
- NormalizeToBinWidth(proj2);
+ TH2* convoluted = CalculateMultiplicityESD(tmpESDEfficiencyRecorrected, esdCorrId);
+ TH1* convolutedProj = convoluted->ProjectionY("convolutedProj", -1, -1, "e");
+ if (convolutedProj->Integral() > 0)
+ {
+ convolutedProj->Scale(1.0 / convolutedProj->Integral());
+ }
+ else
+ printf("ERROR: convolutedProj is empty. Something went wrong calculating the convoluted histogram.\n");
+
+ //NormalizeToBinWidth(proj2);
- TH1* residual = (TH1*) proj2->Clone("residual");
- residual->SetTitle("Residuals;Ntracks;(folded unfolded measured - measured) / measured");
+ TH1* residual = (TH1*) convolutedProj->Clone("residual");
+ residual->SetTitle("Residuals;Ntracks;(folded unfolded measured - measured) / e");
residual->Add(fCurrentESD, -1);
- residual->Divide(residual, fCurrentESD, 1, 1, "B");
+ //residual->Divide(residual, fCurrentESD, 1, 1, "B");
+
+ TH1* residualHist = new TH1F("residualHist", "residualHist", 50, -5, 5);
// TODO fix errors
- for (Int_t i=1; i<residual->GetNbinsX(); ++i)
+ Float_t chi2 = 0;
+ for (Int_t i=1; i<=residual->GetNbinsX(); ++i)
{
- proj2->SetBinError(i, 0);
+ if (fCurrentESD->GetBinError(i) > 0)
+ {
+ Float_t value = residual->GetBinContent(i) / fCurrentESD->GetBinError(i);
+ if (i > 1)
+ chi2 += value * value;
+ residual->SetBinContent(i, value);
+ residualHist->Fill(value);
+ }
+ else
+ {
+ //printf("Residual bin %d set to 0\n", i);
+ residual->SetBinContent(i, 0);
+ }
+ convolutedProj->SetBinError(i, 0);
residual->SetBinError(i, 0);
- }
- // normalize mc to 1
- //mcHist->Scale(1.0 / mcHist->Integral(1, 200));
+ if (i == 200)
+ fLastChi2Residuals = chi2;
+ }
- // normalize unfolded esd to 1
- //fMultiplicityESDCorrected[esdCorrId]->Scale(1.0 / fMultiplicityESDCorrected[esdCorrId]->Integral(1, 200));
+ residualHist->Fit("gaus", "N");
+ delete residualHist;
- // normalize unfolded result to mc hist
- fMultiplicityESDCorrected[esdCorrId]->Scale(1.0 / fMultiplicityESDCorrected[esdCorrId]->Integral());
- fMultiplicityESDCorrected[esdCorrId]->Scale(mcHist->Integral(1, 200));
+ printf("Difference (Residuals) is %f for bin 2-200\n", chi2);
- TCanvas* canvas1 = new TCanvas(tmpStr, tmpStr, 1200, 800);
- canvas1->Divide(2, 2);
+ TCanvas* canvas1 = 0;
+ if (simple)
+ {
+ canvas1 = new TCanvas(tmpStr, tmpStr, 900, 400);
+ canvas1->Divide(2, 1);
+ }
+ else
+ {
+ canvas1 = new TCanvas(tmpStr, tmpStr, 1200, 1200);
+ canvas1->Divide(2, 3);
+ }
canvas1->cd(1);
TH1* proj = (TH1*) mcHist->Clone("proj");
NormalizeToBinWidth(fMultiplicityESDCorrected[esdCorrId]);
proj->GetXaxis()->SetRangeUser(0, 200);
+ proj->SetTitle(";true multiplicity;Entries");
+ proj->SetStats(kFALSE);
proj->DrawCopy("HIST");
gPad->SetLogy();
//fMultiplicityESDCorrected[esdCorrId]->SetMarkerStyle(3);
fMultiplicityESDCorrected[esdCorrId]->SetLineColor(2);
- fMultiplicityESDCorrected[esdCorrId]->DrawCopy("SAME HIST");
+ fMultiplicityESDCorrected[esdCorrId]->DrawCopy("SAME HIST E");
- Float_t chi2 = 0;
+ TLegend* legend = new TLegend(0.3, 0.8, 0.93, 0.93);
+ legend->AddEntry(proj, "true distribution");
+ legend->AddEntry(fMultiplicityESDCorrected[esdCorrId], "unfolded distribution");
+ legend->SetFillColor(0);
+ legend->Draw();
+ // unfortunately does not work. maybe a bug? --> legend->SetTextSizePixels(14);
+
+ canvas1->cd(2);
+
+ gPad->SetLogy();
+ fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
+ //fCurrentESD->SetLineColor(2);
+ fCurrentESD->SetTitle(";measured multiplicity;Entries");
+ fCurrentESD->SetStats(kFALSE);
+ fCurrentESD->DrawCopy("HIST E");
+
+ convolutedProj->SetLineColor(2);
+ //proj2->SetMarkerColor(2);
+ //proj2->SetMarkerStyle(5);
+ convolutedProj->DrawCopy("HIST SAME");
+
+ legend = new TLegend(0.3, 0.8, 0.93, 0.93);
+ legend->AddEntry(fCurrentESD, "measured distribution");
+ legend->AddEntry(convolutedProj, "R #otimes unfolded distribution");
+ legend->SetFillColor(0);
+ legend->Draw();
+
+ TH1* diffMCUnfolded = dynamic_cast<TH1*> (proj->Clone("diffMCUnfolded"));
+ diffMCUnfolded->Add(fMultiplicityESDCorrected[esdCorrId], -1);
+
+ /*Float_t chi2 = 0;
+ Float_t chi = 0;
fLastChi2MCLimit = 0;
+ Int_t limit = 0;
for (Int_t i=2; i<=200; ++i)
{
if (proj->GetBinContent(i) != 0)
{
Float_t value = (proj->GetBinContent(i) - fMultiplicityESDCorrected[esdCorrId]->GetBinContent(i)) / proj->GetBinContent(i);
chi2 += value * value;
+ chi += TMath::Abs(value);
- if (chi2 < 0.1)
- fLastChi2MCLimit = i;
-
- if (i == 100)
- fLastChi2MC = chi2;
- }
- }
+ //printf("%d %f\n", i, chi);
- printf("Difference (from MC) is %f for bin 2-100. Limit is %d.\n", fLastChi2MC, fLastChi2MCLimit);
+ if (chi2 < 0.2)
+ fLastChi2MCLimit = i;
- canvas1->cd(2);
+ if (chi < 3)
+ limit = i;
- NormalizeToBinWidth(fCurrentESD);
- gPad->SetLogy();
- fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
- //fCurrentESD->SetLineColor(2);
- fCurrentESD->DrawCopy("HIST");
-
- proj2->SetLineColor(2);
- //proj2->SetMarkerColor(2);
- //proj2->SetMarkerStyle(5);
- proj2->DrawCopy("HIST SAME");
+ }
+ }*/
chi2 = 0;
- for (Int_t i=2; i<=100; ++i)
+ Float_t chi = 0;
+ Int_t limit = 0;
+ for (Int_t i=1; i<=diffMCUnfolded->GetNbinsX(); ++i)
{
- if (fCurrentESD->GetBinContent(i) != 0)
+ if (fMultiplicityESDCorrected[esdCorrId]->GetBinError(i) > 0)
{
- Float_t value = (proj2->GetBinContent(i) - fCurrentESD->GetBinContent(i)) / fCurrentESD->GetBinContent(i);
- chi2 += value * value;
+ Double_t value = diffMCUnfolded->GetBinContent(i) / fMultiplicityESDCorrected[esdCorrId]->GetBinError(i);
+ if (value > 1e8)
+ value = 1e8; //prevent arithmetic exception
+ else if (value < -1e8)
+ value = -1e8;
+ if (i > 1)
+ {
+ chi2 += value * value;
+ chi += TMath::Abs(value);
+ }
+ diffMCUnfolded->SetBinContent(i, value);
}
+ else
+ {
+ //printf("diffMCUnfolded bin %d set to 0\n", i);
+ diffMCUnfolded->SetBinContent(i, 0);
+ }
+ if (chi2 < 1000)
+ fLastChi2MCLimit = i;
+ if (chi < 1000)
+ limit = i;
+ if (i == 150)
+ fLastChi2MC = chi2;
}
- printf("Difference (Residuals) is %f for bin 2-100\n", chi2);
- fLastChi2Residuals = chi2;
+ printf("limits %d %d\n", fLastChi2MCLimit, limit);
+ fLastChi2MCLimit = limit;
+
+ printf("Difference (from MC) is %f for bin 2-150. Limit is %d.\n", fLastChi2MC, fLastChi2MCLimit);
- canvas1->cd(3);
- TH1* clone = dynamic_cast<TH1*> (proj->Clone("clone"));
- clone->Divide(fMultiplicityESDCorrected[esdCorrId]);
- clone->SetTitle("Ratio;Npart;MC/Unfolded Measured");
- clone->GetYaxis()->SetRangeUser(0.8, 1.2);
- clone->GetXaxis()->SetRangeUser(0, 200);
- clone->DrawCopy("HIST");
+ if (!simple)
+ {
+ canvas1->cd(3);
+
+ diffMCUnfolded->SetTitle("#chi^{2};Npart;(MC - Unfolded) / e");
+ //diffMCUnfolded->GetYaxis()->SetRangeUser(-20, 20);
+ diffMCUnfolded->GetXaxis()->SetRangeUser(0, 200);
+ diffMCUnfolded->DrawCopy("HIST");
- /*TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85);
- legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected");
- legend->AddEntry(fMultiplicityMC, "MC");
- legend->AddEntry(fMultiplicityESD, "ESD");
- legend->Draw();*/
+ TH1F* fluctuation = new TH1F("fluctuation", "fluctuation", 20, -5, 5);
+ for (Int_t i=20; i<=diffMCUnfolded->GetNbinsX(); ++i)
+ fluctuation->Fill(diffMCUnfolded->GetBinContent(i));
- canvas1->cd(4);
- residual->GetYaxis()->SetRangeUser(-0.2, 0.2);
- residual->GetXaxis()->SetRangeUser(0, 200);
- residual->DrawCopy();
+ new TCanvas;
+ fluctuation->Draw();
+
+ /*TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85);
+ legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected");
+ legend->AddEntry(fMultiplicityMC, "MC");
+ legend->AddEntry(fMultiplicityESD, "ESD");
+ legend->Draw();*/
+
+ canvas1->cd(4);
+ //residual->GetYaxis()->SetRangeUser(-0.2, 0.2);
+ residual->GetXaxis()->SetRangeUser(0, 200);
+ residual->DrawCopy();
+
+ canvas1->cd(5);
+
+ TH1* ratio = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("ratio");
+ ratio->Divide(mcHist);
+ ratio->SetTitle("Ratio;true multiplicity;Unfolded / MC");
+ ratio->GetYaxis()->SetRangeUser(0.5, 1.5);
+ ratio->GetXaxis()->SetRangeUser(0, 200);
+ ratio->SetStats(kFALSE);
+ ratio->Draw("HIST");
+
+ Double_t ratioChi2 = 0;
+ fLastChi2MCLimit = 0;
+ for (Int_t i=2; i<=150; ++i)
+ {
+ Float_t value = ratio->GetBinContent(i) - 1;
+ if (value > 1e8)
+ value = 1e8; //prevent arithmetic exception
+ else if (value < -1e8)
+ value = -1e8;
+
+ ratioChi2 += value * value;
+
+ if (ratioChi2 < 0.1)
+ fLastChi2MCLimit = i;
+ }
+
+ printf("Sum over (ratio-1)^2 (2..150) is %f\n", ratioChi2);
+ // TODO FAKE
+ fLastChi2MC = ratioChi2;
+ }
canvas1->SaveAs(Form("%s.gif", canvas1->GetName()));
}
Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
- SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
+ SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE);
- // normalize ESD
- fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
+ // smooth efficiency
+ //TH1* tmp = (TH1*) fCurrentEfficiency->Clone("eff_clone");
+ //for (Int_t i=2; i<fCurrentEfficiency->GetNbinsX(); ++i)
+ // fCurrentEfficiency->SetBinContent(i, (tmp->GetBinContent(i-1) + tmp->GetBinContent(i) + tmp->GetBinContent(i+1)) / 3);
+
+ // set efficiency to 1 above 150
+ // FAKE TEST
+ //for (Int_t i=150; i<=fCurrentEfficiency->GetNbinsX(); ++i)
+ // fCurrentEfficiency->SetBinContent(i, 1);
// normalize correction for given nPart
for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
}
}
+ // normalize nTrack
+ /*for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
+ {
+ // with this it is normalized to 1
+ Double_t sum = fCurrentCorrelation->Integral(1, fCurrentCorrelation->GetNbinsX(), j, j);
+
+ for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
+ {
+ if (sum > 0)
+ {
+ fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
+ fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
+ }
+ else
+ {
+ fCurrentCorrelation->SetBinContent(i, j, 0);
+ fCurrentCorrelation->SetBinError(i, j, 0);
+ }
+ }
+ }*/
+
// smooth input spectrum
/*
TH1* esdClone = (TH1*) fCurrentESD->Clone("esdClone");
fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
*/
+ /*new TCanvas;
+ fCurrentEfficiency->Draw();
+
+ new TCanvas;
+ fCurrentCorrelation->DrawCopy("COLZ");
+
+ new TCanvas;
+ ((TH2*) fCurrentCorrelation)->ProjectionX()->DrawCopy();
+
+ new TCanvas;
+ ((TH2*) fCurrentCorrelation)->ProjectionY()->DrawCopy();*/
+
// pick prior distribution
TH1* hPrior = 0;
if (inputDist)
TH1F* convergence = new TH1F("convergence", "convergence", 50, 0.5, 50.5);
+ const Int_t kStartBin = 1;
+
// unfold...
for (Int_t i=0; i<nIterations; i++)
{
for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
{
Float_t norm = 0;
- for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++)
+ for (Int_t t = kStartBin; t<=hResponse->GetNbinsX(); t++)
norm += hResponse->GetBinContent(t,m) * hPrior->GetBinContent(t);
- for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++)
+ for (Int_t t = kStartBin; t<=hResponse->GetNbinsX(); t++)
{
if (norm==0)
hInverseResponseBayes->SetBinContent(t,m,0);
}
}
- for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++)
+ for (Int_t t = kStartBin; t<=hResponse->GetNbinsX(); t++)
{
Float_t value = 0;
for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
}
// this is the last guess, fill before (!) smoothing
- for (Int_t t=1; t<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); t++)
+ for (Int_t t=kStartBin; t<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); t++)
{
+ //as bin error put the difference to the last iteration
+ //fMultiplicityESDCorrected[correlationID]->SetBinError(t, hTemp->GetBinContent(t) - fMultiplicityESDCorrected[correlationID]->GetBinContent(t));
fMultiplicityESDCorrected[correlationID]->SetBinContent(t, hTemp->GetBinContent(t));
- fMultiplicityESDCorrected[correlationID]->SetBinError(t, 0.05 * hTemp->GetBinContent(t)); // TODO
+ fMultiplicityESDCorrected[correlationID]->SetBinError(t, 0.05 * hTemp->GetBinContent(t));
//printf(" bin %d content %f \n", t, fMultiplicityESDCorrected[correlationID]->GetBinContent(t));
}
+ /*new TCanvas;
+ fMultiplicityESDCorrected[correlationID]->DrawCopy();
+ gPad->SetLogy();*/
+
// regularization (simple smoothing)
TH1F* hTrueSmoothed = (TH1F*) hTemp->Clone("truesmoothed");
- for (Int_t t=2; t<hTrueSmoothed->GetNbinsX(); t++)
+ for (Int_t t=kStartBin+2; t<hTrueSmoothed->GetNbinsX(); t++)
{
Float_t average = (hTemp->GetBinContent(t-1)
+ hTemp->GetBinContent(t)
//for (Int_t t=1; t<=hPrior->GetNbinsX(); t++)
// norm = norm + hTrueSmoothed->GetBinContent(t);
- for (Int_t t=1; t<hTrueSmoothed->GetNbinsX(); t++)
+ for (Int_t t=kStartBin; t<=hTrueSmoothed->GetNbinsX(); t++)
{
Float_t newValue = hTrueSmoothed->GetBinContent(t)/norm;
Float_t diff = hPrior->GetBinContent(t) - newValue;
Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
- SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
+ SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE);
// TODO should be taken from correlation map
//TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
- SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx);
+ SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx, kFALSE);
NormalizeToBinWidth((TH2*) fCurrentCorrelation);
}
//____________________________________________________________________
-TH2F* AliMultiplicityCorrection::CalculateMultiplicityESD(TH1* inputMC, Int_t correlationMap, Bool_t normalized)
+TH2F* AliMultiplicityCorrection::CalculateMultiplicityESD(TH1* inputMC, Int_t correlationMap)
{
// runs the distribution given in inputMC through the response matrix identified by
// correlationMap and produces a measured distribution
}
}
- TH1* corr = hist->Project3D("zy");
+ TH2* corr = (TH2*) hist->Project3D("zy");
+ //corr->Rebin2D(2, 1);
corr->Sumw2();
// normalize correction for given nPart
Float_t measured = 0;
Float_t error = 0;
- for (Int_t gen=1; gen<=target->GetNbinsY(); ++gen)
+ for (Int_t gen=1; gen<=corr->GetNbinsX(); ++gen)
{
- Int_t mcGenBin = inputMC->GetXaxis()->FindBin(target->GetYaxis()->GetBinCenter(gen));
-
- Float_t factor = 1;
- if (normalized)
- factor = inputMC->GetBinWidth(mcGenBin);
+ Int_t mcGenBin = inputMC->GetXaxis()->FindBin(corr->GetXaxis()->GetBinCenter(gen));
- measured += inputMC->GetBinContent(mcGenBin) * factor * corr->GetBinContent(gen, meas);
- error += inputMC->GetBinError(mcGenBin) * factor * corr->GetBinContent(gen, meas);
+ measured += inputMC->GetBinContent(mcGenBin) * corr->GetBinContent(gen, meas);
+ error += inputMC->GetBinError(mcGenBin) * corr->GetBinContent(gen, meas);
}
- // bin width of the <measured> axis has to be taken into account
- //measured /= target->GetYaxis()->GetBinWidth(meas);
- //error /= target->GetYaxis()->GetBinWidth(meas);
-
//printf("%f +- %f ; %f +- %f \n", inputMC->GetBinContent(meas), inputMC->GetBinError(meas), measured, error);
- target->SetBinContent(target->GetNbinsX() / 2, meas, measured);
- target->SetBinError(target->GetNbinsX() / 2, meas, error);
+ target->SetBinContent(1 + target->GetNbinsX() / 2, meas, measured);
+ target->SetBinError(1 + target->GetNbinsX() / 2, meas, error);
}
return target;
class TF1;
class TCollection;
-#include <TMatrixF.h>
-#include <TVectorF.h>
+#include <TMatrixD.h>
+#include <TVectorD.h>
class AliMultiplicityCorrection : public TNamed {
public:
enum EventType { kTrVtx = 0, kMB, kINEL };
- enum RegularizationType { kNone = 0, kPol0, kPol1, kCurvature, kEntropy, kTest };
+ enum RegularizationType { kNone = 0, kPol0, kPol1, kEntropy, kCurvature, kTest };
AliMultiplicityCorrection();
AliMultiplicityCorrection(const Char_t* name, const Char_t* title);
Bool_t LoadHistograms(const Char_t* dir);
void SaveHistograms();
void DrawHistograms();
- void DrawComparison(const char* name, Int_t inputRange, Bool_t fullPhaseSpace, Bool_t normalizeESD, TH1* mcHist);
+ void DrawComparison(const char* name, Int_t inputRange, Bool_t fullPhaseSpace, Bool_t normalizeESD, TH1* mcHist, Bool_t simple = kFALSE);
- void ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check = kFALSE, TH1* inputDist = 0);
- void SetRegularizationParameters(RegularizationType type, Float_t weight) { fRegularizationType = type; fRegularizationWeight = weight; };
+ Int_t ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check = kFALSE, TH1* inputDist = 0);
+ void SetRegularizationParameters(RegularizationType type, Float_t weight);
void ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSpace);
- void ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar = 0.07, Int_t nIterations = 30, TH1* inputDist = 0);
+ void ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar = 0.1, Int_t nIterations = 15, TH1* inputDist = 0);
void ApplyGaussianMethod(Int_t inputRange, Bool_t fullPhaseSpace);
void SetMultiplicityESDCorrected(Int_t i, TH1F* hist) { fMultiplicityESDCorrected[i] = hist; }
void SetGenMeasFromFunc(TF1* inputMC, Int_t id);
- TH2F* CalculateMultiplicityESD(TH1* inputMC, Int_t correlationMap, Bool_t normalized = kFALSE);
+ TH2F* CalculateMultiplicityESD(TH1* inputMC, Int_t correlationMap);
static void NormalizeToBinWidth(TH1* hist);
static void NormalizeToBinWidth(TH2* hist);
protected:
enum { kESDHists = 4, kMCHists = 5, kCorrHists = 8 };
- static const Int_t fgMaxParams; // number of fit params
+ static const Int_t fgMaxParams; // bins in unfolded histogram = number of fit params
+ static const Int_t fgMaxInput; // bins in measured histogram
- static Double_t RegularizationPol0(Double_t *params);
- static Double_t RegularizationPol1(Double_t *params);
- static Double_t RegularizationTotalCurvature(Double_t *params);
- static Double_t RegularizationEntropy(Double_t *params);
- static Double_t RegularizationTest(Double_t *params);
+ static Double_t RegularizationPol0(TVectorD& params);
+ static Double_t RegularizationPol1(TVectorD& params);
+ static Double_t RegularizationTotalCurvature(TVectorD& params);
+ static Double_t RegularizationEntropy(TVectorD& params);
+ static Double_t RegularizationTest(TVectorD& params);
static void MinuitFitFunction(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t);
static void MinuitNBD(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3);
- void SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType);
+ void SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t createBigBin);
Float_t BayesCovarianceDerivate(Float_t matrixM[251][251], TH2* hResponse, TH1* fCurrentEfficiency, Int_t k, Int_t i, Int_t r, Int_t u);
static TH1* fCurrentCorrelation; // static variable to be accessed by MINUIT
static TH1* fCurrentEfficiency; // static variable to be accessed by MINUIT
- static TMatrixF* fCorrelationMatrix; // contains fCurrentCorrelation in matrix form
- static TMatrixF* fCorrelationCovarianceMatrix; // contains the errors of fCurrentESD
- static TVectorF* fCurrentESDVector; // contains fCurrentESD
+ static TMatrixD* fCorrelationMatrix; // contains fCurrentCorrelation in matrix form
+ static TMatrixD* fCorrelationCovarianceMatrix; // contains the errors of fCurrentESD
+ static TVectorD* fCurrentESDVector; // contains fCurrentESD
+ static TVectorD* fEntropyAPriori; // a-priori distribution for entropy regularization
static TF1* fNBD; // negative binomial distribution
#include <TH2F.h>
#include <TH3F.h>
#include <TParticle.h>
+#include <TRandom.h>
+#include <TNtuple.h>
#include <AliLog.h>
#include <AliESD.h>
#include "esdTrackCuts/AliESDtrackCuts.h"
#include "AliPWG0Helper.h"
+#include "AliPWG0depHelper.h"
#include "dNdEta/AliMultiplicityCorrection.h"
+#include "AliCorrection.h"
+#include "AliCorrectionMatrix3D.h"
//#define TPCMEASUREMENT
#define ITSMEASUREMENT
AliMultiplicityMCSelector::AliMultiplicityMCSelector() :
AliSelectorRL(),
fMultiplicity(0),
- fEsdTrackCuts(0)
+ fEsdTrackCuts(0),
+ fSystSkipParticles(kFALSE),
+ fSelectProcessType(0),
+ fParticleSpecies(0)
{
//
// Constructor. Initialization of pointers
//
+
+ for (Int_t i = 0; i<4; i++)
+ fParticleCorrection[i] = 0;
}
AliMultiplicityMCSelector::~AliMultiplicityMCSelector()
ReadUserObjects(tree);
fMultiplicity = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+
+ TString option(GetOption());
+ if (option.Contains("skip-particles"))
+ {
+ fSystSkipParticles = kTRUE;
+ AliInfo("WARNING: Systematic study enabled. Particles will be skipped.");
+ }
+
+ if (option.Contains("particle-efficiency"))
+ for (Int_t i = 0; i<4; i++)
+ fParticleCorrection[i] = new AliCorrection(Form("correction_%d", i), Form("correction_%d", i));
+
+ if (option.Contains("only-process-type-nd"))
+ fSelectProcessType = 1;
+
+ if (option.Contains("only-process-type-sd"))
+ fSelectProcessType = 2;
+
+ if (option.Contains("only-process-type-dd"))
+ fSelectProcessType = 3;
+
+ if (fSelectProcessType != 0)
+ AliInfo(Form("WARNING: Systematic study enabled. Only considering process type %d", fSelectProcessType));
+
+ if (option.Contains("particle-species"))
+ fParticleSpecies = new TNtuple("fParticleSpecies", "fParticleSpecies", "Pi_True:K_True:p_True:o_True:Pi_Rec:K_Rec:p_Rec:o_Rec");
}
void AliMultiplicityMCSelector::Init(TTree* tree)
return kFALSE;
}
+ if (fSelectProcessType > 0)
+ {
+ // getting process information; NB: this only works for Pythia
+ Int_t processtype = AliPWG0depHelper::GetPythiaEventProcessType(header);
+
+ Bool_t processEvent = kFALSE;
+
+ // non diffractive
+ if (fSelectProcessType == 1 && processtype !=92 && processtype !=93 && processtype != 94)
+ processEvent = kTRUE;
+
+ // single diffractive
+ if (fSelectProcessType == 2 && (processtype == 92 || processtype == 93))
+ processEvent = kTRUE;
+
+ // double diffractive
+ if (fSelectProcessType == 3 && processtype == 94)
+ processEvent = kTRUE;
+
+ if (!processEvent)
+ {
+ AliDebug(AliLog::kDebug, Form("Skipping this event, because it is not of the requested process type (%d)", processtype));
+ return kTRUE;
+ }
+ }
+
Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD);
+ eventTriggered = kTRUE; // no generated trigger for the simulated events
Bool_t eventVertex = AliPWG0Helper::IsVertexReconstructed(fESD);
// get the ESD vertex
Int_t nMCTracks20 = 0;
Int_t nMCTracksAll = 0;
+ // tracks per particle species, in |eta| < 2 (systematic study)
+ Int_t nMCTracksSpecies[4]; // (pi, K, p, other)
+ for (Int_t i = 0; i<4; ++i)
+ nMCTracksSpecies[i] = 0;
+
for (Int_t iMc = 0; iMc < nPrim; ++iMc)
{
AliDebug(AliLog::kDebug+1, Form("MC Loop: Processing particle %d.", iMc));
continue;
}
- if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
+ Bool_t debug = kFALSE;
+ if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim, debug) == kFALSE)
+ {
+ //printf("%d) DROPPED ", iMC);
+ //particle->Print();
continue;
+ }
+
+ //printf("%d) OK ", iMC);
+ //particle->Print();
//if (particle->Pt() < kPtCut)
// continue;
nMCTracks20++;
nMCTracksAll++;
+
+ Int_t id = -1;
+ switch (TMath::Abs(particle->GetPdgCode()))
+ {
+ case 211: id = 0; break;
+ case 321: id = 1; break;
+ case 2212: id = 2; break;
+ default: id = 3; break;
+ }
+ if (TMath::Abs(particle->Eta()) < 2.0)
+ nMCTracksSpecies[id]++;
+ if (fParticleCorrection[id])
+ fParticleCorrection[id]->GetTrackCorrection()->FillGene(vtxMC[2], particle->Eta(), particle->Pt());
}// end of mc particle
- // FAKE: put back vtxMC[2]
fMultiplicity->FillGenerated(vtxMC[2], eventTriggered, eventVertex, nMCTracks05, nMCTracks10, nMCTracks15, nMCTracks20, nMCTracksAll);
if (!eventTriggered || !eventVertex)
Int_t nESDTracks15 = 0;
Int_t nESDTracks20 = 0;
+ // tracks per particle species, in |eta| < 2 (systematic study)
+ Int_t nESDTracksSpecies[4]; // (pi, K, p, other)
+ for (Int_t i = 0; i<4; ++i)
+ nESDTracksSpecies[i] = 0;
+
#ifdef ITSMEASUREMENT
// get multiplicity from ITS tracklets
for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
if (mult->GetDeltaPhi(i) < -1000)
continue;
+ // systematic study: 10% lower efficiency
+ if (fSystSkipParticles && gRandom->Uniform() < 0.1)
+ continue;
+
Float_t theta = mult->GetTheta(i);
Float_t eta = -TMath::Log(TMath::Tan(theta/2.));
if (TMath::Abs(eta) < 2.0)
nESDTracks20++;
+
+ // TODO define needed, because this only works with new AliRoot
+ Int_t label = mult->GetLabel(i);
+ if (label < 0)
+ continue;
+
+ TParticle* mother = AliPWG0depHelper::FindPrimaryMother(stack, label);
+
+ if (!mother)
+ continue;
+
+ Int_t id = -1;
+ switch (TMath::Abs(mother->GetPdgCode()))
+ {
+ case 211: id = 0; break;
+ case 321: id = 1; break;
+ case 2212: id = 2; break;
+ default: id = 3; break;
+ }
+ if (TMath::Abs(eta) < 2.0)
+ nESDTracksSpecies[id]++;
+ if (fParticleCorrection[id])
+ fParticleCorrection[id]->GetTrackCorrection()->FillMeas(vtxMC[2], mother->Eta(), mother->Pt());
}
#endif
if (TMath::Abs(eta) < 2.0)
nESDTracks20++;
+
+ Int_t label = esdTrack->GetLabel();
+ if (label < 0)
+ continue;
+
+ TParticle* mother = AliPWG0depHelper::FindPrimaryMother(stack, label);
+
+ if (!mother)
+ continue;
+
+ Int_t id = -1;
+ switch (TMath::Abs(mother->GetPdgCode()))
+ {
+ case 211: id = 0; break;
+ case 321: id = 1; break;
+ case 2212: id = 2; break;
+ default: id = 3; break;
+ }
+ if (TMath::Abs(eta) < 2.0)
+ nESDTracksSpecies[id]++;
+ if (fParticleCorrection[id])
+ fParticleCorrection[id]->GetTrackCorrection()->FillMeas(vtxMC[2], mother->Eta(), mother->Pt());
}
#endif
+ if (nMCTracks20 == 0 && nESDTracks20 > 3)
+ printf("WARNING: Event %lld has %d generated and %d reconstructed...\n", entry, nMCTracks05, nESDTracks05);
+
fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks10, nESDTracks15, nESDTracks20);
// TODO should this be vtx[2] or vtxMC[2] ?
fMultiplicity->FillCorrection(vtxMC[2], nMCTracks05, nMCTracks10, nMCTracks15, nMCTracks20, nMCTracksAll, nESDTracks05, nESDTracks10, nESDTracks15, nESDTracks20);
+ if (fParticleSpecies)
+ fParticleSpecies->Fill(nMCTracksSpecies[0], nMCTracksSpecies[1], nMCTracksSpecies[2], nMCTracksSpecies[3], nESDTracksSpecies[0], nESDTracksSpecies[1], nESDTracksSpecies[2], nESDTracksSpecies[3]);
+
return kTRUE;
}
}
fOutput->Add(fMultiplicity);
+ for (Int_t i = 0; i < 4; ++i)
+ fOutput->Add(fParticleCorrection[i]);
+
+ fOutput->Add(fParticleSpecies);
}
void AliMultiplicityMCSelector::Terminate()
AliSelector::Terminate();
fMultiplicity = dynamic_cast<AliMultiplicityCorrection*> (fOutput->FindObject("Multiplicity"));
+ for (Int_t i = 0; i < 4; ++i)
+ fParticleCorrection[i] = dynamic_cast<AliCorrection*> (fOutput->FindObject(Form("correction_%d", i)));
+ fParticleSpecies = dynamic_cast<TNtuple*> (fOutput->FindObject("fParticleSpecies"));
if (!fMultiplicity)
{
TFile* file = TFile::Open("multiplicityMC.root", "RECREATE");
fMultiplicity->SaveHistograms();
+ for (Int_t i = 0; i < 4; ++i)
+ if (fParticleCorrection[i])
+ fParticleCorrection[i]->SaveHistograms();
+ fParticleSpecies->Write();
file->Close();
class AliESDtrackCuts;
class AliMultiplicityCorrection;
+class AliCorrection;
+class TNtuple;
class AliMultiplicityMCSelector : public AliSelectorRL {
public:
void ReadUserObjects(TTree* tree);
AliMultiplicityCorrection* fMultiplicity; // object containing the extracted data
- AliESDtrackCuts* fEsdTrackCuts; // Object containing the parameters of the esd track cuts
+ AliESDtrackCuts* fEsdTrackCuts; // Object containing the parameters of the esd track cuts
+
+ Bool_t fSystSkipParticles; // if true skips particles (systematic study)
+ AliCorrection* fParticleCorrection[4]; // correction from measured to generated particles for trigger, vertex sample in |eta| < 2;
+ // for each of the species: pi, k, p, other; for systematic study of pt cut off
+ Int_t fSelectProcessType; // 0 = all (default), 1 = ND, 2 = SD, 3 = DD (for systematic study)
+ TNtuple *fParticleSpecies; // per event: (pi, k, p, rest (in |eta| < 2)) X (true, recon); (for systematic study)
private:
AliMultiplicityMCSelector(const AliMultiplicityMCSelector&);
// getting process information NB: this only works for Pythia !!!
Int_t processtype = AliPWG0depHelper::GetPythiaEventProcessType(header);
-
+
// can only read pythia headers, either directly or from cocktalil header
AliGenEventHeader* genHeader = (AliGenEventHeader*)(header->GenEventHeader());
Bool_t vertexReconstructed = AliPWG0Helper::IsVertexReconstructed(fESD);
// non diffractive
- if (processtype!=92 && processtype!=93 && processtype!=94) {
+ if (processtype!=92 && processtype!=93 && processtype!=94) {
// NB: passing the wrong process type here (1), since the process type is defined by the index in the array (here non-diffractive)
if (triggerBiasStudy) fdNdEtaCorrectionTriggerBias[0]->FillEvent(vtxMC[2], nGoodTracks, eventTriggered, vertexReconstructed, 1);
if (vertexRecoStudy) fdNdEtaCorrectionVertexReco[0] ->FillEvent(vtxMC[2], nGoodTracks, eventTriggered, vertexReconstructed, 1);
if (vertexRecoStudy) fdNdEtaCorrectionVertexReco[2] ->FillTrackedParticle(vtxMC[2], eta, pt);
}
-
- TParticle* mother = particle;
// find primary particle that created this particle
- while (label >= nPrim)
- {
- //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0));
-
- if (mother->GetMother(0) == -1)
- {
- AliDebug(AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label));
- mother = 0;
- break;
- }
-
- label = mother->GetMother(0);
-
- mother = stack->Particle(label);
- if (!mother)
- {
- AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (find mother loop).", label));
- break;
- }
- }
-
+ TParticle* mother = AliPWG0depHelper::FindPrimaryMother(stack, label);
if (!mother)
continue;
void runMultiplicitySelector(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aMC = kFALSE, Bool_t aDebug = kFALSE, Bool_t aProof = kFALSE, const char* option = "")
{
if (aProof)
- connectProof("jgrosseo@lxb6046");
+ connectProof("proof40@lxb6046");
TString libraries("libEG;libGeom;libESD;libPWG0base");
TString packages("PWG0base");
if (aDebug != kFALSE)
selectorName += "g";
+ //Int_t result = chain->Process(selectorName, option);
Int_t result = executeQuery(chain, &inputList, selectorName, option);
if (result != 0)
mult->LoadHistograms("Multiplicity");
mult->DrawHistograms();
+
+ TH2* hist = (TH2*) gROOT->FindObject("fCorrelation3_zy");
+ canvas = new TCanvas("c1", "c1", 600, 500);
+ hist->SetStats(kFALSE);
+ hist->Draw("COLZ");
+ hist->SetTitle(";true multiplicity in |#eta| < 2;measured multiplicity in |#eta| < 2");
+ hist->GetYaxis()->SetTitleOffset(1.1);
+ gPad->SetRightMargin(0.15);
+ gPad->SetLogz();
+
+ canvas->SaveAs("Plot_Correlation.pdf");
}
void* fit(const char* fileName = "multiplicityMC.root", Int_t hist = 2, Int_t eventType = 0)
//return;
- mult->ApplyBayesianMethod(hist, kFALSE, AliMultiplicityCorrection::kTrVtx);
+ /*mult->ApplyBayesianMethod(hist, kFALSE, AliMultiplicityCorrection::kTrVtx);
mult->DrawComparison("Bayesian", hist, kFALSE, kTRUE, mult->GetMultiplicityMC(hist, AliMultiplicityCorrection::kTrVtx)->ProjectionY());
-
- return;
+ return;*/
TStopwatch timer;
timer.Start();
- mult->SetRegularizationParameters(AliMultiplicityCorrection::kTest, 1);
- mult->ApplyMinuitFit(hist, kFALSE, AliMultiplicityCorrection::kTrVtx);
+ mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 0.2);
+ mult->ApplyMinuitFit(hist, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE);
mult->DrawComparison("MinuitChi2", hist, kFALSE, kFALSE, mult->GetMultiplicityMC(hist, AliMultiplicityCorrection::kTrVtx)->ProjectionY());
timer.Stop();
return mult;
}
-void* fitOther(const char* fileNameMC = "multiplicityMC_3M.root", const char* fileNameESD = "multiplicityMC_3M.root", Int_t histID = 2)
+void fitOther(const char* fileNameMC = "multiplicityMC_3M.root", const char* fileNameESD = "multiplicityMC_3M.root", Bool_t chi2 = kTRUE, Int_t histID = 3, Bool_t fullPhaseSpace = kFALSE)
{
gSystem->Load("libPWG0base");
TFile::Open(fileNameESD);
TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", histID));
- TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", histID));
- //hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityMB%d", histID));
+ TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", ((fullPhaseSpace) ? 4 : histID)));
+ //hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityINEL%d", histID));
mult->SetMultiplicityESD(histID, hist);
- mult->SetMultiplicityVtx(histID, hist2);
- mult->ApplyLaszloMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx);
- return;
-
- mult->SetRegularizationParameters(AliMultiplicityCorrection::kTest, 1.1);
- mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx);
- mult->DrawComparison("MinuitChi2", histID, kFALSE, kTRUE, hist2->ProjectionY());
+ // small hack to get around charge conservation for full phase space ;-)
+ if (fullPhaseSpace)
+ {
+ TH1* corr = mult->GetCorrelation(histID + 4);
- return;
+ for (Int_t i=2; i<=corr->GetNbinsX(); i+=2)
+ for (Int_t j=1; j<=corr->GetNbinsY(); ++j)
+ {
+ corr->SetBinContent(i, j, corr->GetBinContent(i-1, j));
+ corr->SetBinError(i, j, corr->GetBinError(i-1, j));
+ }
+ }
- //mult->ApplyGaussianMethod(histID, kFALSE);
+ /*mult->SetMultiplicityVtx(histID, hist2);
+ mult->ApplyLaszloMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx);
+ return;*/
- for (Float_t f=0.1; f<=0.11; f+=0.05)
+ if (chi2)
+ {
+ mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 1e4);
+ mult->ApplyMinuitFit(histID, fullPhaseSpace, AliMultiplicityCorrection::kTrVtx, kFALSE); //, hist2->ProjectionY("mymchist"));
+ mult->DrawComparison("MinuitChi2", histID, fullPhaseSpace, kTRUE, hist2->ProjectionY("mymchist"));
+ }
+ else
{
- mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, f);
- mult->DrawComparison(Form("Bayesian_%f", f), histID, kFALSE, kTRUE, hist2->ProjectionY());
+ mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
+ mult->DrawComparison("Bayesian", histID, kFALSE, kTRUE, hist2->ProjectionY("mymchist2"));
}
//mult->SetRegularizationParameters(AliMultiplicityCorrection::kEntropy, 1e7);
//mult->ApplyMinuitFit(histID, kFALSE);
//mult->DrawComparison("MinuitChi2", histID, kFALSE, kTRUE, hist2->ProjectionY());
+}
+
+void* fit2Step(const char* fileNameMC = "multiplicityMC_2M.root", const char* fileNameESD = "multiplicityMC_1M_3.root", Int_t histID = 3, Bool_t fullPhaseSpace = kFALSE)
+{
+ gSystem->Load("libPWG0base");
+
+ AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+
+ TFile::Open(fileNameMC);
+ mult->LoadHistograms("Multiplicity");
+
+ TFile::Open(fileNameESD);
+ TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", histID));
+ TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", ((fullPhaseSpace) ? 4 : histID)));
+ //hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityINEL%d", histID));
+
+ mult->SetMultiplicityESD(histID, hist);
+
+ // small hack to get around charge conservation for full phase space ;-)
+ if (fullPhaseSpace)
+ {
+ TH1* corr = mult->GetCorrelation(histID + 4);
+
+ for (Int_t i=2; i<=corr->GetNbinsX(); i+=2)
+ for (Int_t j=1; j<=corr->GetNbinsY(); ++j)
+ {
+ corr->SetBinContent(i, j, corr->GetBinContent(i-1, j));
+ corr->SetBinError(i, j, corr->GetBinError(i-1, j));
+ }
+ }
+
+ mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
+ mult->ApplyMinuitFit(histID, fullPhaseSpace, AliMultiplicityCorrection::kTrVtx, kFALSE);
+ mult->DrawComparison("MinuitChi2", histID, fullPhaseSpace, kTRUE, hist2->ProjectionY("mymchist"));
+
+ TH1* result = (TH1*) mult->GetMultiplicityESDCorrected((fullPhaseSpace) ? 4 : histID))->Clone("firstresult");
+
+ mult->SetRegularizationParameters(AliMultiplicityCorrection::kEntropy, 100000);
+ mult->ApplyMinuitFit(histID, fullPhaseSpace, AliMultiplicityCorrection::kTrVtx, kFALSE, result);
+ mult->DrawComparison("MinuitChi2_Step2", histID, fullPhaseSpace, kTRUE, hist2->ProjectionY("mymchist"));
return mult;
}
return 0;
}
-void EvaluateBayesianMethod(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 2)
+void EvaluateBayesianMethod(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 3)
{
gSystem->mkdir(targetDir);
legend->AddEntry(fitResultsMC, Form("%s MC chi2", GetEventTypeName(type)));
legend->AddEntry(fitResultsRes, Form("%s residual chi2", GetEventTypeName(type)));
- for (Float_t weight = 0; weight < 0.301; weight += 0.02)
+ for (Float_t weight = 0; weight < 1.01; weight += 0.1)
{
Float_t chi2MC = 0;
Float_t residuals = 0;
canvas->SaveAs(Form("%s/%s.C", targetDir, canvas->GetName()));
}
-void EvaluateChi2Method(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 2)
+void EvaluateBayesianMethodIterationsSmoothing(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 3)
+{
+ gSystem->mkdir(targetDir);
+
+ gSystem->Load("libPWG0base");
+
+ AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+ TFile::Open(fileNameMC);
+ mult->LoadHistograms("Multiplicity");
+
+ AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD");
+ TFile::Open(fileNameESD);
+ multESD->LoadHistograms("Multiplicity");
+ mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID));
+
+ Int_t count = 0; // just to order the saved images...
+
+ TFile* graphFile = TFile::Open(Form("%s/EvaluateBayesianMethodIterationsSmoothing.root", targetDir), "RECREATE");
+
+ for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kINEL; ++type)
+ {
+ TString tmp;
+ tmp.Form("EvaluateBayesianMethodIterationsSmoothing_%s", GetEventTypeName(type));
+
+ TCanvas* canvas = new TCanvas(tmp, tmp, 800, 600);
+
+ for (Int_t i = 1; i <= 7; i++)
+ {
+ Int_t iter = i * 20;
+ TGraph* fitResultsMC = new TGraph;
+ fitResultsMC->SetTitle(Form("%d iterations", iter));
+ fitResultsMC->GetXaxis()->SetTitle("smoothing parameter");
+ fitResultsMC->GetYaxis()->SetTitle("P_{1} (2 <= t <= 150)");
+ fitResultsMC->SetName(Form("%s_MC_%d", tmp.Data(), i));
+ TGraph* fitResultsRes = new TGraph;
+ fitResultsRes->SetTitle(Form("%d iterations", iter));
+ fitResultsRes->SetName(Form("%s_Res_%d", tmp.Data(), i));
+ fitResultsRes->GetXaxis()->SetTitle("smoothing parameter");
+ fitResultsRes->GetYaxis()->SetTitle("P_{2}");
+
+ fitResultsMC->SetFillColor(0);
+ fitResultsRes->SetFillColor(0);
+ fitResultsMC->SetMarkerStyle(19+i);
+ fitResultsRes->SetMarkerStyle(19+i);
+ fitResultsRes->SetMarkerColor(kRed);
+ fitResultsRes->SetLineColor(kRed);
+
+ for (Float_t weight = 0.0; weight < 1.01; weight += 0.2)
+ {
+ Float_t chi2MC = 0;
+ Float_t residuals = 0;
+
+ mult->ApplyBayesianMethod(histID, kFALSE, type, weight, iter);
+ mult->DrawComparison(Form("%s/BayesianIterSmooth_%03d_%d_%d_%f", targetDir, count++, type, iter, weight), histID, kFALSE, kTRUE, multESD->GetMultiplicityMC(histID, type)->ProjectionY());
+ mult->GetComparisonResults(&chi2MC, 0, &residuals);
+
+ fitResultsMC->SetPoint(fitResultsMC->GetN(), weight, chi2MC);
+ fitResultsRes->SetPoint(fitResultsRes->GetN(), weight, residuals);
+ }
+
+ fitResultsMC->Write();
+ fitResultsRes->Write();
+ }
+ }
+}
+
+void EvaluateDrawResult(const char* targetDir, Int_t type = 0, Bool_t plotRes = kTRUE)
+{
+ gSystem->Load("libPWG0base");
+
+ TString name;
+ TFile* graphFile = 0;
+ if (type == -1)
+ {
+ name = "EvaluateChi2Method";
+ graphFile = TFile::Open(Form("%s/EvaluateChi2Method.root", targetDir));
+ }
+ else
+ {
+ name.Form("EvaluateBayesianMethodIterationsSmoothing_%s", GetEventTypeName(type));
+ graphFile = TFile::Open(Form("%s/EvaluateBayesianMethodIterationsSmoothing.root", targetDir));
+ }
+
+ TCanvas* canvas = new TCanvas(name, name, 800, 500);
+ if (type == -1)
+ canvas->SetLogx();
+ canvas->SetLogy();
+
+ TLegend* legend = new TLegend(0.6, 0.75, 0.88, 0.88);
+ legend->SetFillColor(0);
+
+ Int_t count = 1;
+
+ Float_t xMin = 1e20;
+ Float_t xMax = 0;
+
+ Float_t yMin = 1e20;
+ Float_t yMax = 0;
+
+ TString xaxis, yaxis;
+
+ while (1)
+ {
+ TGraph* mc = (TGraph*) graphFile->Get(Form("%s_MC_%d", name.Data(), count));
+ TGraph* res = (TGraph*) graphFile->Get(Form("%s_Res_%d", name.Data(), count));
+
+ if (!mc || !res)
+ break;
+
+ xaxis = mc->GetXaxis()->GetTitle();
+ yaxis = mc->GetYaxis()->GetTitle();
+
+ mc->Print();
+ res->Print();
+
+ count++;
+
+ if (plotRes)
+ {
+ xMin = TMath::Min(TMath::Min(xMin, mc->GetXaxis()->GetXmin()), res->GetXaxis()->GetXmin());
+ yMin = TMath::Min(TMath::Min(yMin, mc->GetYaxis()->GetXmin()), res->GetYaxis()->GetXmin());
+
+ xMax = TMath::Max(TMath::Max(xMax, mc->GetXaxis()->GetXmax()), res->GetXaxis()->GetXmax());
+ yMax = TMath::Max(TMath::Max(yMax, mc->GetYaxis()->GetXmax()), res->GetYaxis()->GetXmax());
+ }
+ else
+ {
+ xMin = TMath::Min(xMin, mc->GetXaxis()->GetXmin());
+ yMin = TMath::Min(yMin, mc->GetYaxis()->GetXmin());
+
+ xMax = TMath::Max(xMax, mc->GetXaxis()->GetXmax());
+ yMax = TMath::Max(yMax, mc->GetYaxis()->GetXmax());
+ }
+ }
+
+ if (type >= 0)
+ {
+ xaxis = "smoothing parameter";
+ }
+ else if (type == -1)
+ {
+ xaxis = "weight parameter";
+ //yaxis = "P_{3} (2 <= t <= 150)";
+ }
+ yaxis = "P_{1} (2 <= t <= 150)";
+
+ printf("%f %f %f %f\n", xMin, xMax, yMin, yMax);
+
+ TGraph* dummy = new TGraph;
+ dummy->SetPoint(0, xMin, yMin);
+ dummy->SetPoint(1, xMax, yMax);
+ dummy->SetTitle(Form(";%s;%s", xaxis.Data(), yaxis.Data()));
+
+ dummy->SetMarkerColor(0);
+ dummy->Draw("AP");
+
+ count = 1;
+
+ while (1)
+ {
+ TGraph* mc = (TGraph*) graphFile->Get(Form("%s_MC_%d", name.Data(), count));
+ TGraph* res = (TGraph*) graphFile->Get(Form("%s_Res_%d", name.Data(), count));
+
+ printf("%s_MC_%d %p %p\n", name.Data(), count, mc, res);
+
+ if (!mc || !res)
+ break;
+
+ printf("Loaded %d sucessful.\n", count);
+
+ if (type == -1)
+ {
+ legend->AddEntry(mc, Form("Eq. (%d)", count+9));
+ }
+ else
+ legend->AddEntry(mc);
+
+ mc->Draw("SAME PC");
+
+ if (plotRes)
+ {
+ legend->AddEntry(res);
+ res->Draw("SAME PC");
+ }
+
+ count++;
+ }
+
+ legend->Draw();
+
+ canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName()));
+ canvas->SaveAs(Form("%s/%s.eps", targetDir, canvas->GetName()));
+}
+
+void EvaluateChi2Method(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 3)
{
gSystem->mkdir(targetDir);
mult->SetMultiplicityESD(histID, hist);
- TCanvas* canvas = new TCanvas("EvaluateChi2Method", "EvaluateChi2Method", 800, 600);
- TLegend* legend = new TLegend(0.2, 0.7, 0.58, 0.98);
- legend->SetFillColor(0);
-
- Float_t min = 1e10;
- Float_t max = 0;
-
- TGraph* first = 0;
Int_t count = 0; // just to order the saved images...
- Bool_t firstLoop = kTRUE;
+ TGraph* fitResultsMC = 0;
+ TGraph* fitResultsRes = 0;
+
+ TFile* graphFile = TFile::Open(Form("%s/EvaluateChi2Method.root", targetDir), "RECREATE");
for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kPol0; type <= AliMultiplicityCorrection::kEntropy; ++type)
- //for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kPol1; type <= AliMultiplicityCorrection::kPol1; ++type)
+// for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kPol1; type <= AliMultiplicityCorrection::kPol1; ++type)
{
- TGraph* fitResultsMC = new TGraph;
- fitResultsMC->SetTitle(";Weight Parameter");
- TGraph* fitResultsRes = new TGraph;
- fitResultsRes->SetTitle(";Weight Parameter");
+ fitResultsMC = new TGraph;
+ fitResultsMC->SetTitle(Form("%s", GetRegName(type)));
+ fitResultsMC->GetXaxis()->SetTitle("Weight Parameter");
+ fitResultsMC->SetName(Form("EvaluateChi2Method_MC_%d", type));
+ fitResultsRes = new TGraph;
+ fitResultsRes->SetTitle(Form("%s residual chi2", GetRegName(type)));
+ fitResultsRes->SetName(Form("EvaluateChi2Method_Res_%d", type));
+ fitResultsRes->GetXaxis()->SetTitle("Weight Parameter");
fitResultsMC->SetFillColor(0);
fitResultsRes->SetFillColor(0);
fitResultsRes->SetMarkerColor(kRed);
fitResultsRes->SetLineColor(kRed);
- legend->AddEntry(fitResultsMC, Form("%s MC chi2", GetRegName(type)));
- legend->AddEntry(fitResultsRes, Form("%s residual chi2", GetRegName(type)));
+ for (Int_t i=0; i<5; ++i)
+ {
+ Float_t weight = TMath::Power(TMath::Sqrt(10), i+6);
+ //Float_t weight = TMath::Power(10, i+2);
- if (first == 0)
- first = fitResultsMC;
+ //if (type == AliMultiplicityCorrection::kEntropy) weight = 1e4 * (i+1) * 1.5;
- for (Float_t weight = 1e-4; weight < 2e4; weight *= 10)
- //for (Float_t weight = 0.1; weight < 10; weight *= TMath::Sqrt(TMath::Sqrt(10)))
- {
Float_t chi2MC = 0;
Float_t residuals = 0;
+ Float_t chi2Limit = 0;
+
+ TString runName;
+ runName.Form("MinuitChi2_%02d_%d_%f", count++, type, weight);
mult->SetRegularizationParameters(type, weight);
- mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx);
- mult->DrawComparison(Form("%s/MinuitChi2_%02d_%d_%f", targetDir, count++, type, weight), histID, kFALSE, kTRUE, hist2->ProjectionY());
+ Int_t status = mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx);
+ mult->DrawComparison(Form("%s/%s", targetDir, runName.Data()), histID, kFALSE, kTRUE, hist2->ProjectionY());
+ if (status != 0)
+ {
+ printf("MINUIT did not succeed. Skipping...\n");
+ continue;
+ }
+
mult->GetComparisonResults(&chi2MC, 0, &residuals);
+ TH1* result = mult->GetMultiplicityESDCorrected(histID);
+ result->SetName(runName);
+ result->Write();
fitResultsMC->SetPoint(fitResultsMC->GetN(), weight, chi2MC);
fitResultsRes->SetPoint(fitResultsRes->GetN(), weight, residuals);
-
- min = TMath::Min(min, TMath::Min(chi2MC, residuals));
- max = TMath::Max(max, TMath::Max(chi2MC, residuals));
}
- fitResultsMC->Print();
- fitResultsRes->Print();
-
- canvas->cd();
- fitResultsMC->Draw(Form("%s CP", (firstLoop) ? "A" : "SAME"));
- fitResultsRes->Draw("SAME CP");
-
- firstLoop = kFALSE;
+ graphFile->cd();
+ fitResultsMC->Write();
+ fitResultsRes->Write();
}
- gPad->SetLogx();
- gPad->SetLogy();
- printf("min = %f, max = %f\n", min, max);
- if (min <= 0)
- min = 1e-5;
- first->GetYaxis()->SetRangeUser(min * 0.5, max * 1.5);
-
- legend->Draw();
-
- canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName()));
- canvas->SaveAs(Form("%s/%s.C", targetDir, canvas->GetName()));
+ graphFile->Close();
}
void EvaluateChi2MethodAll()
EvaluateBayesianMethod("multiplicityMC_2M_smoothed.root", "multiplicityMC_3M_NBD.root", "eval-2MS-NBD");
}
-void CompareMethods(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 2)
+void CompareMethods(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 3)
{
gSystem->mkdir(targetDir);
mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID));
- TCanvas* canvas = new TCanvas("CompareMethods", "CompareMethods", 1200, 800);
- canvas->Divide(3, 2);
+ TCanvas* canvas = new TCanvas("CompareMethods", "CompareMethods", 1200, 1200);
+ canvas->Divide(3, 3);
Int_t count = 0;
- for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kINEL; ++type)
+ for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kTrVtx; ++type)
{
- TH1* mc = multESD->GetMultiplicityMC(histID, type)->ProjectionY();
+ TH1* mc = multESD->GetMultiplicityMC(histID, type)->ProjectionY("mymc");
+ mc->Sumw2();
- mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 0.3);
+ mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
mult->ApplyMinuitFit(histID, kFALSE, type);
mult->DrawComparison(Form("%s/CompareMethods_%d_MinuitChi2", targetDir, type), histID, kFALSE, kTRUE, mc);
TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("chi2Result");
mc->GetXaxis()->SetRangeUser(0, 150);
chi2Result->GetXaxis()->SetRangeUser(0, 150);
- // skip errors for now
+/* // skip errors for now
for (Int_t i=1; i<=chi2Result->GetNbinsX(); ++i)
{
chi2Result->SetBinError(i, 0);
bayesResult->SetBinError(i, 0);
- }
+ }*/
canvas->cd(++count);
mc->SetFillColor(kYellow);
gPad->SetLogy();
canvas->cd(count + 3);
- chi2Result->Divide(chi2Result, mc, 1, 1, "B");
- bayesResult->Divide(bayesResult, mc, 1, 1, "B");
+ chi2ResultRatio = (TH1*) chi2Result->Clone("chi2ResultRatio");
+ bayesResultRatio = (TH1*) bayesResult->Clone("bayesResultRatio");
+ chi2ResultRatio->Divide(chi2Result, mc, 1, 1, "");
+ bayesResultRatio->Divide(bayesResult, mc, 1, 1, "");
- // skip errors for now
- for (Int_t i=1; i<=chi2Result->GetNbinsX(); ++i)
- {
- chi2Result->SetBinError(i, 0);
- bayesResult->SetBinError(i, 0);
- }
+ chi2ResultRatio->GetYaxis()->SetRangeUser(0.5, 1.5);
- chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5);
+ chi2ResultRatio->DrawCopy("HIST");
+ bayesResultRatio->DrawCopy("SAME HIST");
- chi2Result->DrawCopy("");
- bayesResult->DrawCopy("SAME");
+ canvas->cd(count + 6);
+ chi2Result->Divide(chi2Result, bayesResult, 1, 1, "");
+ chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5);
+ chi2Result->DrawCopy("HIST");
}
canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName()));
canvas->SaveAs(Form("%s/%s.C", targetDir, canvas->GetName()));
}
-void StatisticsPlot(const char* fileNameMC = "multiplicityMC.root", Int_t histID = 2)
+void StatisticsPlot(const char* fileNameMC = "multiplicityMC_2M.root", Int_t histID = 3)
{
gSystem->Load("libPWG0base");
TFile::Open(fileNameMC);
mult->LoadHistograms("Multiplicity");
- const char* files[] = { "multiplicityMC_0.5M.root", "multiplicityMC_0.75M.root", "multiplicityMC_1M_3.root", "multiplicityMC_1.25M.root", "multiplicityMC_1.5M.root" };
+ const char* files[] = { "multiplicityMC_100k_1.root", "multiplicityMC_200k.root", "multiplicityMC_400k.root", "multiplicityMC_600k.root", "multiplicityMC_800k.root" };
- TGraph* fitResultsChi2 = new TGraph;
- fitResultsChi2->SetTitle(";Nevents;Chi2");
- TGraph* fitResultsBayes = new TGraph;
- fitResultsBayes->SetTitle(";Nevents;Chi2");
- TGraph* fitResultsChi2Limit = new TGraph;
- fitResultsChi2Limit->SetTitle(";Nevents;Multiplicity reach");
- TGraph* fitResultsBayesLimit = new TGraph;
- fitResultsBayesLimit->SetTitle(";Nevents;Multiplicity reach");
+ TGraph* fitResultsChi2 = new TGraph; fitResultsChi2->SetTitle(";Nevents;Chi2");
+ TGraph* fitResultsBayes = new TGraph; fitResultsBayes->SetTitle(";Nevents;Chi2");
+ TGraph* fitResultsChi2Limit = new TGraph; fitResultsChi2Limit->SetTitle(";Nevents;Multiplicity reach");
+ TGraph* fitResultsBayesLimit = new TGraph; fitResultsBayesLimit->SetTitle(";Nevents;Multiplicity reach");
+
+ fitResultsChi2->SetName("fitResultsChi2");
+ fitResultsBayes->SetName("fitResultsBayes");
+ fitResultsChi2Limit->SetName("fitResultsChi2Limit");
+ fitResultsBayesLimit->SetName("fitResultsBayesLimit");
TCanvas* canvas = new TCanvas("StatisticsPlot", "StatisticsPlot", 1200, 600);
canvas->Divide(5, 2);
Float_t min = 1e10;
Float_t max = 0;
+ TFile* file = TFile::Open("StatisticsPlot.root", "RECREATE");
+ file->Close();
+
for (Int_t i=0; i<5; ++i)
{
TFile::Open(files[i]);
mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID));
Int_t nEntries = multESD->GetMultiplicityESD(histID)->GetEntries();
- TH1* mc = multESD->GetMultiplicityMC(histID, AliMultiplicityCorrection::kTrVtx)->ProjectionY();
+ TH1* mc = multESD->GetMultiplicityMC(histID, AliMultiplicityCorrection::kTrVtx)->ProjectionY(Form("mc_%d", i));
- mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 0.3);
+ mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx);
mult->DrawComparison(Form("StatisticsPlot_%d_MinuitChi2", i), histID, kFALSE, kTRUE, mc);
min = TMath::Min(min, chi2MC);
max = TMath::Max(max, chi2MC);
- TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("chi2Result");
+ TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("chi2Result_%d", i));
- mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 0.1);
+ mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 0.07, 10);
mult->DrawComparison(Form("StatisticsPlot_%d_Bayesian", i), histID, kFALSE, kTRUE, mc);
- TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("bayesResult");
+ TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("bayesResult_%d", i));
mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0);
fitResultsBayes->SetPoint(fitResultsBayes->GetN(), nEntries, chi2MC);
fitResultsBayesLimit->SetPoint(fitResultsBayesLimit->GetN(), nEntries, chi2MCLimit);
chi2Result->DrawCopy("");
bayesResult->DrawCopy("SAME");
+
+ TFile* file = TFile::Open("StatisticsPlot.root", "UPDATE");
+ mc->Write();
+ chi2Result->Write();
+ bayesResult->Write();
+ file->Close();
}
canvas->SaveAs(Form("%s.gif", canvas->GetName()));
canvas2->SaveAs(Form("%s.gif", canvas2->GetName()));
canvas2->SaveAs(Form("%s.C", canvas2->GetName()));
+
+ TFile* file = TFile::Open("StatisticsPlot.root", "UPDATE");
+ fitResultsChi2->Write();
+ fitResultsBayes->Write();
+ fitResultsChi2Limit->Write();
+ fitResultsBayesLimit->Write();
+ file->Close();
}
-void StartingConditions(const char* fileNameMC = "multiplicityMC.root", Int_t histID = 2)
+void StartingConditions(const char* fileNameMC = "multiplicityMC_2M.root", Int_t histID = 3)
{
gSystem->Load("libPWG0base");
TFile::Open(fileNameMC);
mult->LoadHistograms("Multiplicity");
- const char* files[] = { "multiplicityMC_0.5M.root", "multiplicityMC_0.75M.root", "multiplicityMC_1M_3.root", "multiplicityMC_1.25M.root", "multiplicityMC_1.5M.root" };
+ const char* files[] = { "multiplicityMC_1M_3.root", "multiplicityMC_100k_1.root", "multiplicityMC_100k_2.root", "multiplicityMC_100k_3.root", "multiplicityMC_100k_4.root" }
// this one we try to unfold
TFile::Open(files[0]);
AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD");
multESD->LoadHistograms("Multiplicity");
mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID));
- TH1* mc = multESD->GetMultiplicityMC(histID, AliMultiplicityCorrection::kTrVtx)->ProjectionY();
+ TH1* mc = multESD->GetMultiplicityMC(histID, AliMultiplicityCorrection::kTrVtx)->ProjectionY("mc");
TGraph* fitResultsChi2 = new TGraph;
fitResultsChi2->SetTitle(";Input Dist ID;Chi2");
TLegend* legend = new TLegend(0.7, 0.7, 1, 1);
+ TFile* file = TFile::Open("StartingConditions.root", "RECREATE");
+ mc->Write();
+ file->Close();
+
for (Int_t i=0; i<8; ++i)
{
if (i == 0)
startCond->SetBinContent(j, 1);
}
- mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 0.3);
+ mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE, startCond);
mult->DrawComparison(Form("StartingConditions_%d_MinuitChi2", i), histID, kFALSE, kTRUE, mc);
min = TMath::Min(min, chi2MC);
max = TMath::Max(max, chi2MC);
- TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("chi2Result");
+ TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("chi2Result_%d", i));
if (!firstChi)
firstChi = (TH1*) chi2Result->Clone("firstChi");
- mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 0.1, 30, startCond);
+ mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, startCond);
mult->DrawComparison(Form("StartingConditions_%d_Bayesian", i), histID, kFALSE, kTRUE, mc);
- TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("bayesResult");
+ TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("bayesResult_%d", i));
if (!firstBayesian)
firstBayesian = (TH1*) bayesResult->Clone("firstBayesian");
fitResultsBayes->SetPoint(fitResultsBayes->GetN(), i, chi2MC);
fitResultsBayesLimit->SetPoint(fitResultsBayesLimit->GetN(), i, chi2MCLimit);
+ TFile* file = TFile::Open("StartingConditions.root", "UPDATE");
+ chi2Result->Write();
+ bayesResult->Write();
+ file->Close();
+
min = TMath::Min(min, chi2MC);
max = TMath::Max(max, chi2MC);
mc->GetXaxis()->SetRangeUser(0, 150);
canvas3->SaveAs(Form("%s.gif", canvas3->GetName()));
}
+void DifferentSamples(const char* fileNameMC = "multiplicityMC_2M.root", Int_t histID = 3)
+{
+ gSystem->Load("libPWG0base");
+
+ AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+
+ TFile::Open(fileNameMC);
+ mult->LoadHistograms("Multiplicity");
+
+ const char* files[] = { "multiplicityMC_100k_1.root", "multiplicityMC_100k_2.root", "multiplicityMC_100k_3.root", "multiplicityMC_100k_4.root", "multiplicityMC_100k_5.root", "multiplicityMC_100k_6.root", "multiplicityMC_100k_7.root", "multiplicityMC_100k_8.root" };
+
+ TGraph* fitResultsChi2 = new TGraph;
+ fitResultsChi2->SetTitle(";Input Dist ID;Chi2");
+ TGraph* fitResultsBayes = new TGraph;
+ fitResultsBayes->SetTitle(";Input Dist ID;Chi2");
+ TGraph* fitResultsChi2Limit = new TGraph;
+ fitResultsChi2Limit->SetTitle(";Input Dist ID;Multiplicity reach");
+ TGraph* fitResultsBayesLimit = new TGraph;
+ fitResultsBayesLimit->SetTitle(";Input Dist ID;Multiplicity reach");
+
+ TCanvas* canvasA = new TCanvas("DifferentSamplesA", "DifferentSamplesA", 1200, 600);
+ canvasA->Divide(4, 2);
+
+ TCanvas* canvasB = new TCanvas("DifferentSamplesB", "DifferentSamplesB", 1200, 600);
+ canvasB->Divide(4, 2);
+
+ TCanvas* canvas4 = new TCanvas("DifferentSamples4", "DifferentSamples4", 1000, 400);
+ canvas4->Divide(2, 1);
+
+ TCanvas* canvas3 = new TCanvas("DifferentSamples3", "DifferentSamples3", 1000, 400);
+ canvas3->Divide(2, 1);
+
+ Float_t min = 1e10;
+ Float_t max = 0;
+
+ TH1* firstChi = 0;
+ TH1* firstBayesian = 0;
+
+ TLegend* legend = new TLegend(0.7, 0.7, 1, 1);
+
+ TFile* file = TFile::Open("DifferentSamples.root", "RECREATE");
+ file->Close();
+
+ for (Int_t i=0; i<8; ++i)
+ {
+ TFile::Open(files[i]);
+ AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD2", "MultiplicityESD2");
+ multESD->LoadHistograms("Multiplicity");
+ mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID));
+ TH1* mc = multESD->GetMultiplicityMC(histID, AliMultiplicityCorrection::kTrVtx)->ProjectionY(Form("mc_%d", i));
+ mc->Sumw2();
+
+ mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
+ mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE);
+ mult->DrawComparison(Form("DifferentSamples_%d_MinuitChi2", i), histID, kFALSE, kTRUE, mc);
+
+ Float_t chi2MC = 0;
+ Int_t chi2MCLimit = 0;
+ mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0);
+ fitResultsChi2->SetPoint(fitResultsChi2->GetN(), i, chi2MC);
+ fitResultsChi2Limit->SetPoint(fitResultsChi2Limit->GetN(), i, chi2MCLimit);
+ min = TMath::Min(min, chi2MC);
+ max = TMath::Max(max, chi2MC);
+
+ TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("chi2Result_%d", i));
+ if (!firstChi)
+ firstChi = (TH1*) chi2Result->Clone("firstChi");
+
+ mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
+ mult->DrawComparison(Form("DifferentSamples_%d_Bayesian", i), histID, kFALSE, kTRUE, mc);
+ TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("bayesResult_%d", i));
+ if (!firstBayesian)
+ firstBayesian = (TH1*) bayesResult->Clone("firstBayesian");
+
+ TFile* file = TFile::Open("DifferentSamples.root", "UPDATE");
+ mc->Write();
+ chi2Result->Write();
+ bayesResult->Write();
+ file->Close();
+
+ mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0);
+ fitResultsBayes->SetPoint(fitResultsBayes->GetN(), i, chi2MC);
+ fitResultsBayesLimit->SetPoint(fitResultsBayesLimit->GetN(), i, chi2MCLimit);
+
+ min = TMath::Min(min, chi2MC);
+ max = TMath::Max(max, chi2MC);
+ mc->GetXaxis()->SetRangeUser(0, 150);
+ chi2Result->GetXaxis()->SetRangeUser(0, 150);
+
+ // skip errors for now
+ for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j)
+ {
+ chi2Result->SetBinError(j, 0);
+ bayesResult->SetBinError(j, 0);
+ }
+
+ canvas4->cd(1);
+ TH1* tmp = (TH1*) chi2Result->Clone("tmp");
+ tmp->SetTitle("Unfolded/MC;Npart;Ratio");
+ tmp->Divide(mc);
+ tmp->GetYaxis()->SetRangeUser(0.5, 1.5);
+ tmp->GetXaxis()->SetRangeUser(0, 200);
+ tmp->SetLineColor(i+1);
+ tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST");
+
+ canvas4->cd(2);
+ tmp = (TH1*) bayesResult->Clone("tmp");
+ tmp->SetTitle("Unfolded/MC;Npart;Ratio");
+ tmp->Divide(mc);
+ tmp->SetLineColor(i+1);
+ tmp->GetYaxis()->SetRangeUser(0.5, 1.5);
+ tmp->GetXaxis()->SetRangeUser(0, 200);
+ tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST");
+
+ canvas3->cd(1);
+ TH1* tmp = (TH1*) chi2Result->Clone("tmp");
+ tmp->SetTitle("Ratio to first result;Npart;Ratio");
+ tmp->Divide(firstChi);
+ tmp->GetYaxis()->SetRangeUser(0.5, 1.5);
+ tmp->GetXaxis()->SetRangeUser(0, 200);
+ tmp->SetLineColor(i+1);
+ legend->AddEntry(tmp, Form("%d", i));
+ tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST");
+
+ canvas3->cd(2);
+ tmp = (TH1*) bayesResult->Clone("tmp");
+ tmp->SetTitle("Ratio to first result;Npart;Ratio");
+ tmp->Divide(firstBayesian);
+ tmp->SetLineColor(i+1);
+ tmp->GetYaxis()->SetRangeUser(0.5, 1.5);
+ tmp->GetXaxis()->SetRangeUser(0, 200);
+ tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST");
+
+ if (i < 4)
+ {
+ canvasA->cd(i+1);
+ }
+ else
+ canvasB->cd(i+1-4);
+
+ mc->SetFillColor(kYellow);
+ mc->DrawCopy();
+ chi2Result->SetLineColor(kRed);
+ chi2Result->DrawCopy("SAME");
+ bayesResult->SetLineColor(kBlue);
+ bayesResult->DrawCopy("SAME");
+ gPad->SetLogy();
+
+ if (i < 4)
+ {
+ canvasA->cd(i+5);
+ }
+ else
+ canvasB->cd(i+5-4);
+
+ chi2Result->Divide(chi2Result, mc, 1, 1, "B");
+ bayesResult->Divide(bayesResult, mc, 1, 1, "B");
+
+ // skip errors for now
+ for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j)
+ {
+ chi2Result->SetBinError(j, 0);
+ bayesResult->SetBinError(j, 0);
+ }
+
+ chi2Result->SetTitle("Ratios;Npart;unfolded measured/MC");
+ chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5);
+
+ chi2Result->DrawCopy("");
+ bayesResult->DrawCopy("SAME");
+ }
+
+ canvas3->cd(1);
+ legend->Draw();
+
+ canvasA->SaveAs(Form("%s.gif", canvasA->GetName()));
+ canvasB->SaveAs(Form("%s.gif", canvasB->GetName()));
+
+ TCanvas* canvas2 = new TCanvas("DifferentSamples2", "DifferentSamples2", 800, 400);
+ canvas2->Divide(2, 1);
+
+ canvas2->cd(1);
+ fitResultsChi2->SetMarkerStyle(20);
+ fitResultsChi2->GetYaxis()->SetRangeUser(0.5 * min, 1.5 * max);
+ fitResultsChi2->Draw("AP");
+
+ fitResultsBayes->SetMarkerStyle(3);
+ fitResultsBayes->SetMarkerColor(2);
+ fitResultsBayes->Draw("P SAME");
+
+ gPad->SetLogy();
+
+ canvas2->cd(2);
+ fitResultsChi2Limit->SetMarkerStyle(20);
+ fitResultsChi2Limit->GetYaxis()->SetRangeUser(0.9 * TMath::Min(fitResultsChi2Limit->GetYaxis()->GetXmin(), fitResultsBayesLimit->GetYaxis()->GetXmin()), 1.1 * TMath::Max(fitResultsChi2Limit->GetYaxis()->GetXmax(), fitResultsBayesLimit->GetYaxis()->GetXmax()));
+ fitResultsChi2Limit->Draw("AP");
+
+ fitResultsBayesLimit->SetMarkerStyle(3);
+ fitResultsBayesLimit->SetMarkerColor(2);
+ fitResultsBayesLimit->Draw("P SAME");
+
+ canvas2->SaveAs(Form("%s.gif", canvas2->GetName()));
+ canvas3->SaveAs(Form("%s.gif", canvas3->GetName()));
+ canvas4->SaveAs(Form("%s.gif", canvas4->GetName()));
+}
+
void Merge(Int_t n, const char** files, const char* output)
{
+ // const char* files[] = { "multiplicityMC_100k_1.root", "multiplicityMC_100k_2.root", "multiplicityMC_100k_3.root", "multiplicityMC_100k_4.root", "multiplicityMC_100k_5.root", "multiplicityMC_100k_6.root", "multiplicityMC_100k_7.root", "multiplicityMC_100k_8.root" };
+
+
gSystem->Load("libPWG0base");
AliMultiplicityCorrection** data = new AliMultiplicityCorrection*[n];
data[0]->Merge(&list);
- data[0]->DrawHistograms();
+ //data[0]->DrawHistograms();
TFile::Open(output, "RECREATE");
data[0]->SaveHistograms();
if (caseNo >= 4)
{
- func = new TF1("nbd", "[0] * TMath::Binomial([2]+TMath::Nint(x)-1, [2]-1) * pow([1] / ([1]+[2]), TMath::Nint(x)) * pow(1 + [1]/[2], -[2])", 0, 50);
+ func = new TF1("nbd", "[0] * TMath::Binomial([2]+TMath::Nint(x)-1, [2]-1) * pow([1] / ([1]+[2]), TMath::Nint(x)) * pow(1 + [1]/[2], -[2])", 0, 250);
func->SetParNames("scaling", "averagen", "k");
}
case 2: func = new TF1("flat", "1000 * 1/(x+1)"); break;
case 3: func = new TF1("flat", "1000 * TMath::Landau(x, 10, 5)"); break;
case 4: func->SetParameters(1e7, 10, 2); break;
- case 5: func->SetParameters(1e7, 20, 3); break;
+ case 5: func->SetParameters(1, 13, 7); break;
case 6: func->SetParameters(1e7, 30, 4); break;
- case 7: func->SetParameters(1e7, 70, 2); break;
+ case 7: func->SetParameters(1e7, 30, 2); break;
case 8: func = new TF1("testlaszlo", "10*1000*x*exp(-0.1*x)"); break;
default: return;
}
- mult->SetGenMeasFromFunc(func, 2);
+ new TCanvas;
+ func->Draw();
+
+ mult->SetGenMeasFromFunc(func, 3);
TFile::Open("out.root", "RECREATE");
mult->SaveHistograms();
//mult->ApplyBayesianMethod(2, kFALSE);
//mult->ApplyMinuitFit(2, kFALSE);
//mult->ApplyGaussianMethod(2, kFALSE);
- mult->ApplyLaszloMethod(2, kFALSE, AliMultiplicityCorrection::kTrVtx);
+ //mult->ApplyLaszloMethod(2, kFALSE, AliMultiplicityCorrection::kTrVtx);
}
void smoothCorrelationMap(const char* fileName = "multiplicityMC.root", Int_t corrMatrix = 2)
}
new TCanvas;
- proj->Draw("COLZ");
+ proj->DrawCopy("COLZ");
TH1* scaling = proj->ProjectionY("scaling", 1, 1);
scaling->Reset();
printf("Fitting %d...\n", i);
TH1* hist = proj->ProjectionY(Form("proj%d", i), i, i, "e");
+
//hist->GetXaxis()->SetRangeUser(0, 50);
//lognormal->SetParameter(0, hist->GetMaximum());
hist->Fit(fitWith, "0 M", "");
TF1* func = hist->GetListOfFunctions()->FindObject(fitWith);
- if (((i-1) % 15 == 0) || ((i % 5 == 0) && i < 30))
+ if (0 && (i % 5 == 0))
{
- new TCanvas;
+ pad = new TCanvas;
hist->Draw();
func->Clone()->Draw("SAME");
- gPad->SetLogy();
+ pad->SetLogy();
}
scaling->Fill(i, func->GetParameter(0));
//TF1* scalingFit = new TF1("mypol0", "[0]");
TF1* scalingFit = over;
- scaling->Fit(scalingFit, "", "", 3, 100);
+ scaling->Fit(scalingFit, "", "", 3, 140);
+ scalingFit->SetRange(0, 200);
+ scalingFit->Draw("SAME");
c1->cd(2);
mean->Draw("P");
//TF1* meanFit = log;
TF1* meanFit = new TF1("mypol1", "[0]+[1]*x");
- mean->Fit(meanFit, "", "", 3, 100);
+ mean->Fit(meanFit, "", "", 3, 140);
+ meanFit->SetRange(0, 200);
+ meanFit->Draw("SAME");
c1->cd(3);
width->Draw("P");
//TF1* widthFit = over;
- TF1* widthFit = new TF1("mypol2", "[0]+[1]*x+[2]*x*x");
- width->Fit(widthFit, "", "", 5, 100);
+ TF1* widthFit = new TF1("mypol", "[0]+[1]*TMath::Sqrt([2]*x)");
+ widthFit->SetParLimits(2, 1e-5, 1e5);
+ width->Fit(widthFit, "", "", 5, 140);
+ widthFit->SetRange(0, 200);
+ widthFit->Draw("SAME");
// build new correction matrix
TH2* new = (TH2*) proj->Clone("new");
for (Int_t j=1; j<=new->GetYaxis()->GetNbins(); j+=1)
{
- if (i < 21)
+ if (i < 11)
{
// leave bins 1..20 untouched
new->SetBinContent(i, j, corr->Integral(1, corr->GetNbinsX(), i, i, j, j));
proj1->Draw();
proj2->Draw("SAME");
}
+
+void GetCrossSections(const char* fileName)
+{
+ gSystem->Load("libPWG0base");
+
+ AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+
+ TFile::Open(fileName);
+ mult->LoadHistograms("Multiplicity");
+
+ TH1* xSection2 = mult->GetCorrelation(3)->Project3D("y")->Clone("xSection2");
+ xSection2->Sumw2();
+ xSection2->Scale(1.0 / xSection2->Integral());
+
+ TH1* xSection15 = mult->GetCorrelation(2)->Project3D("y")->Clone("xSection15");
+ xSection15->Sumw2();
+ xSection15->Scale(1.0 / xSection15->Integral());
+
+ TFile::Open("crosssection.root", "RECREATE");
+ xSection2->Write();
+ xSection15->Write();
+ gFile->Close();
+}
void testAnalysis2(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aMC = kFALSE, Bool_t aDebug = kFALSE, Bool_t aProof = kFALSE, const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction", const char* option = "", const char* proofServer = "jgrosseo@lxb6046")
{
if (aProof)
- {
connectProof(proofServer);
- gProof->AddInput(new TParameter<long>("PROOF_MaxSlavesPerNode", (long)2));
- gProof->AddInput(new TNamed("PROOF_Packetizer", "TAdaptivePacketizer"));
- }
TString libraries("libEG;libGeom;libESD;libPWG0base");
- TString packages("adaptivepacketizer;PWG0base");
+ TString packages("PWG0base");
if (!prepareQuery(libraries, packages, kTRUE))
return;
- //TODO somehow prevent loading several times
gROOT->ProcessLine(".L CreateCuts.C");
gROOT->ProcessLine(".L drawPlots.C");
esdTrackCuts/AliTestESDtrackCutsSelector.h \
TPC/AliROCESDAnalysisSelector.h \
TPC/AliROCRawAnalysisSelector.h \
- TPC/AliROCClusterAnalysisSelector.h
+ TPC/AliROCClusterAnalysisSelector.h \
+ highMultiplicity/AliHighMultiplicitySelector.h
SRCS = $(HDRS:.h=.cxx)
DHDR= PWG0selectorsLinkDef.h
-EINCLUDE=PYTHIA6 EVGEN TPC RAW
+EINCLUDE=PYTHIA6 EVGEN TPC RAW ITS
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff