#include <TRandom.h>
#include <TProfile.h>
#include <TProfile2D.h>
+#include <AliLog.h>
ClassImp(AliMultiplicityCorrection)
+// Defined where the efficiency drops below 1/3
+// |eta| < 1.4 --> -0.3 ... 0.8
+// |eta| < 1.3 --> -1.9 ... 2.4
+// |eta| < 1.0 --> -5.6 ... 6.1
+//Double_t AliMultiplicityCorrection::fgVtxRangeBegin[kESDHists] = { -10.0, -5.6, -1.9 };
+//Double_t AliMultiplicityCorrection::fgVtxRangeEnd[kESDHists] = { 10.0, 6.1, 2.4 };
+Double_t AliMultiplicityCorrection::fgVtxRangeBegin[kESDHists] = { -10.0, -5.5, -1.9 };
+Double_t AliMultiplicityCorrection::fgVtxRangeEnd[kESDHists] = { 10.0, 5.5, 2.4 };
+
// These are the areas where the quality of the unfolding results are evaluated
// Default defined here, call SetQualityRegions to change them
// unit is in multiplicity (not in bin!)
-
// SPD: peak area - flat area - low stat area
Int_t AliMultiplicityCorrection::fgQualityRegionsB[kQualityRegions] = {1, 20, 70};
Int_t AliMultiplicityCorrection::fgQualityRegionsE[kQualityRegions] = {10, 65, 80};
//
for (Int_t i = 0; i < kESDHists; ++i)
+ {
fMultiplicityESD[i] = 0;
+ fTriggeredEvents[i] = 0;
+ fNoVertexEvents[i] = 0;
+ }
for (Int_t i = 0; i < kMCHists; ++i)
{
return 0;
}
+ Printf("AliMultiplicityCorrection::Open: Reading file %s", fileName);
+
AliMultiplicityCorrection* mult = new AliMultiplicityCorrection(folderName, folderName);
mult->LoadHistograms();
fLastChi2MC(0),
fLastChi2MCLimit(0),
fLastChi2Residuals(0),
- fRatioAverage(0)
+ fRatioAverage(0),
+ fVtxBegin(0),
+ fVtxEnd(0)
{
//
// named constructor
//
-
+
// do not add this hists to the directory
Bool_t oldStatus = TH1::AddDirectoryStatus();
TH1::AddDirectory(kFALSE);
#define NBINNING 201, -0.5, 200.5
- Double_t vtxRange[] = { 15, 6, 2 };
-
for (Int_t i = 0; i < kESDHists; ++i)
- fMultiplicityESD[i] = new TH2F(Form("fMultiplicityESD%d", i), "fMultiplicityESD;vtx-z;Ntracks;Count", 1, -vtxRange[i], vtxRange[i], NBINNING);
+ {
+ fMultiplicityESD[i] = new TH2F(Form("fMultiplicityESD%d", i), "fMultiplicityESD;vtx-z;measured multiplicity;Count", 1, fgVtxRangeBegin[i], fgVtxRangeEnd[i], NBINNING);
+ fTriggeredEvents[i] = new TH1F(Form("fTriggeredEvents%d", i), "fTriggeredEvents;measured multiplicity;Count", NBINNING);
+ fNoVertexEvents[i] = new TH1F(Form("fNoVertexEvents%d", i), "fNoVertexEvents;generated multiplicity;Count", NBINNING);
+ }
for (Int_t i = 0; i < kMCHists; ++i)
{
for (Int_t i = 0; i < kCorrHists; ++i)
{
- fCorrelation[i] = new TH3F(Form("fCorrelation%d", i), "fCorrelation;vtx-z;Npart;Ntracks", 1, -vtxRange[i%3], vtxRange[i%3], NBINNING, NBINNING);
- fMultiplicityESDCorrected[i] = new TH1F(Form("fMultiplicityESDCorrected%d", i), "fMultiplicityESDCorrected;Npart;dN/dN", NBINNING);
+ fCorrelation[i] = new TH3F(Form("fCorrelation%d", i), "fCorrelation;vtx-z;true multiplicity;measured multiplicity", 1, fgVtxRangeBegin[i%3], fgVtxRangeEnd[i%3], NBINNING, NBINNING);
+ fMultiplicityESDCorrected[i] = new TH1F(Form("fMultiplicityESDCorrected%d", i), "fMultiplicityESDCorrected;true multiplicity;Count", NBINNING);
}
TH1::AddDirectory(oldStatus);
AliUnfolding::SetNbins(120, 120);
AliUnfolding::SetSkipBinsBegin(1);
- AliUnfolding::SetNormalizeInput(kTRUE);
+ //AliUnfolding::SetNormalizeInput(kTRUE);
+}
+
+//____________________________________________________________________
+void AliMultiplicityCorrection::Rebin2DY(TH2F*& hist, Int_t nBins, Double_t* newBins) const
+{
+ //
+ // rebins the y axis of a two-dimensional histogram giving variable size binning (missing in ROOT v5/25/02)
+ //
+
+ TH2F* temp = new TH2F(hist->GetName(), Form("%s;%s;%s", hist->GetTitle(), hist->GetXaxis()->GetTitle(), hist->GetYaxis()->GetTitle()), hist->GetNbinsX(), hist->GetXaxis()->GetXmin(), hist->GetXaxis()->GetXmax(), nBins, newBins);
+
+ for (Int_t x=0; x<=hist->GetNbinsX()+1; x++)
+ for (Int_t y=0; y<=hist->GetNbinsY()+1; y++)
+ temp->Fill(hist->GetXaxis()->GetBinCenter(x), hist->GetYaxis()->GetBinCenter(y), hist->GetBinContent(x, y));
+
+ for (Int_t x=0; x<=temp->GetNbinsX()+1; x++)
+ for (Int_t y=0; y<=temp->GetNbinsY()+1; y++)
+ temp->SetBinError(x, y, TMath::Sqrt(temp->GetBinContent(x, y)));
+
+ delete hist;
+ hist = temp;
+}
+
+//____________________________________________________________________
+void AliMultiplicityCorrection::Rebin3DY(TH3F*& hist, Int_t nBins, Double_t* newBins) const
+{
+ //
+ // rebins the y axis of a three-dimensional histogram giving variable size binning (missing in ROOT v5/25/02)
+ // this function is a mess - and it should have been Fons who should have gone through the pain of writing it! (JF)
+ //
+
+ // construct variable size arrays for fixed size binning axes because TH3 lacks some constructors
+ Double_t* xBins = new Double_t[hist->GetNbinsX()+1];
+ Double_t* zBins = new Double_t[hist->GetNbinsZ()+1];
+
+ for (Int_t x=1; x<=hist->GetNbinsX()+1; x++)
+ {
+ xBins[x-1] = hist->GetXaxis()->GetBinLowEdge(x);
+ //Printf("%d %f", x, xBins[x-1]);
+ }
+
+ for (Int_t z=1; z<=hist->GetNbinsZ()+1; z++)
+ {
+ zBins[z-1] = hist->GetZaxis()->GetBinLowEdge(z);
+ //Printf("%d %f", y, yBins[y-1]);
+ }
+
+ TH3F* temp = new TH3F(hist->GetName(), Form("%s;%s;%s;%s", hist->GetTitle(), hist->GetXaxis()->GetTitle(), hist->GetYaxis()->GetTitle(), hist->GetZaxis()->GetTitle()), hist->GetNbinsX(), xBins, nBins, newBins, hist->GetNbinsZ(), zBins);
+
+ for (Int_t x=0; x<=hist->GetNbinsX()+1; x++)
+ for (Int_t y=0; y<=hist->GetNbinsY()+1; y++)
+ for (Int_t z=0; z<=hist->GetNbinsZ()+1; z++)
+ temp->Fill(hist->GetXaxis()->GetBinCenter(x), hist->GetYaxis()->GetBinCenter(y), hist->GetZaxis()->GetBinCenter(z), hist->GetBinContent(x, y, z));
+
+ for (Int_t x=0; x<=temp->GetNbinsX()+1; x++)
+ for (Int_t y=0; y<=temp->GetNbinsY()+1; y++)
+ for (Int_t z=0; z<=hist->GetNbinsZ()+1; z++)
+ temp->SetBinError(x, y, z, TMath::Sqrt(temp->GetBinContent(x, y, z)));
+
+ delete[] xBins;
+ delete[] zBins;
+
+ delete hist;
+ hist = temp;
+}
+
+//____________________________________________________________________
+void AliMultiplicityCorrection::RebinGenerated(Int_t nBins05, Double_t* newBins05, Int_t nBins10, Double_t* newBins10, Int_t nBins13, Double_t* newBins13)
+{
+ //
+ // Rebins the (and only the) generated multiplicity axis
+ //
+
+ Printf("Rebinning generated-multiplicity axis...");
+
+ // do not add this hists to the directory
+ Bool_t oldStatus = TH1::AddDirectoryStatus();
+ TH1::AddDirectory(kFALSE);
+
+ if (kESDHists != 3)
+ AliFatal("This function only works for three ESD hists!");
+
+ for (Int_t i = 0; i < kESDHists; ++i)
+ {
+ Int_t nBins = -1;
+ Double_t* newBins = 0;
+
+ switch (i)
+ {
+ case 0:
+ nBins = nBins05;
+ newBins = newBins05;
+ break;
+ case 1:
+ nBins = nBins10;
+ newBins = newBins10;
+ break;
+ case 2:
+ nBins = nBins13;
+ newBins = newBins13;
+ break;
+ }
+
+ // 1D
+ // TODO mem leak
+ fNoVertexEvents[i] = (TH1F*) fNoVertexEvents[i]->Rebin(nBins, fNoVertexEvents[i]->GetName(), newBins);
+ fMultiplicityESDCorrected[i] = (TH1F*) fMultiplicityESDCorrected[i]->Rebin(nBins, fMultiplicityESDCorrected[i]->GetName(), newBins);
+
+ // 2D
+ Rebin2DY(fMultiplicityVtx[i], nBins, newBins);
+ Rebin2DY(fMultiplicityMB[i], nBins, newBins);
+ Rebin2DY(fMultiplicityINEL[i], nBins, newBins);
+ Rebin2DY(fMultiplicityNSD[i], nBins, newBins);
+
+ // 3D
+ Rebin3DY(fCorrelation[i], nBins, newBins);
+ }
+
+ TH1::AddDirectory(oldStatus);
}
//____________________________________________________________________
if (fMultiplicityESD[i])
delete fMultiplicityESD[i];
fMultiplicityESD[i] = 0;
+
+ if (fTriggeredEvents[i])
+ delete fTriggeredEvents[i];
+ fTriggeredEvents[i]= 0;
+
+ if (fNoVertexEvents[i])
+ delete fNoVertexEvents[i];
+ fNoVertexEvents[i]= 0;
}
for (Int_t i = 0; i < kMCHists; ++i)
TObject* obj;
// collections of all histograms
- TList collections[kESDHists+kMCHists*4+kCorrHists*2];
+ TList collections[3*kESDHists+kMCHists*4+kCorrHists*2];
Int_t count = 0;
while ((obj = iter->Next())) {
continue;
for (Int_t i = 0; i < kESDHists; ++i)
+ {
collections[i].Add(entry->fMultiplicityESD[i]);
+ collections[kESDHists+i].Add(entry->fTriggeredEvents[i]);
+ collections[kESDHists*2+i].Add(entry->fNoVertexEvents[i]);
+ }
for (Int_t i = 0; i < kMCHists; ++i)
{
- collections[kESDHists+i].Add(entry->fMultiplicityVtx[i]);
- collections[kESDHists+kMCHists+i].Add(entry->fMultiplicityMB[i]);
- collections[kESDHists+kMCHists*2+i].Add(entry->fMultiplicityINEL[i]);
- collections[kESDHists+kMCHists*3+i].Add(entry->fMultiplicityNSD[i]);
+ collections[3*kESDHists+i].Add(entry->fMultiplicityVtx[i]);
+ collections[3*kESDHists+kMCHists+i].Add(entry->fMultiplicityMB[i]);
+ collections[3*kESDHists+kMCHists*2+i].Add(entry->fMultiplicityINEL[i]);
+ collections[3*kESDHists+kMCHists*3+i].Add(entry->fMultiplicityNSD[i]);
}
for (Int_t i = 0; i < kCorrHists; ++i)
- collections[kESDHists+kMCHists*4+i].Add(entry->fCorrelation[i]);
+ collections[3*kESDHists+kMCHists*4+i].Add(entry->fCorrelation[i]);
for (Int_t i = 0; i < kCorrHists; ++i)
- collections[kESDHists+kMCHists*4+kCorrHists+i].Add(entry->fMultiplicityESDCorrected[i]);
+ collections[3*kESDHists+kMCHists*4+kCorrHists+i].Add(entry->fMultiplicityESDCorrected[i]);
count++;
}
for (Int_t i = 0; i < kESDHists; ++i)
+ {
fMultiplicityESD[i]->Merge(&collections[i]);
-
+ fTriggeredEvents[i]->Merge(&collections[kESDHists+i]);
+ fNoVertexEvents[i]->Merge(&collections[2*kESDHists+i]);
+ }
+
for (Int_t i = 0; i < kMCHists; ++i)
{
- fMultiplicityVtx[i]->Merge(&collections[kESDHists+i]);
- fMultiplicityMB[i]->Merge(&collections[kESDHists+kMCHists+i]);
- fMultiplicityINEL[i]->Merge(&collections[kESDHists+kMCHists*2+i]);
- fMultiplicityNSD[i]->Merge(&collections[kESDHists+kMCHists*3+i]);
+ fMultiplicityVtx[i]->Merge(&collections[3*kESDHists+i]);
+ fMultiplicityMB[i]->Merge(&collections[3*kESDHists+kMCHists+i]);
+ fMultiplicityINEL[i]->Merge(&collections[3*kESDHists+kMCHists*2+i]);
+ fMultiplicityNSD[i]->Merge(&collections[3*kESDHists+kMCHists*3+i]);
}
for (Int_t i = 0; i < kCorrHists; ++i)
- fCorrelation[i]->Merge(&collections[kESDHists+kMCHists*4+i]);
+ fCorrelation[i]->Merge(&collections[3*kESDHists+kMCHists*4+i]);
for (Int_t i = 0; i < kCorrHists; ++i)
- fMultiplicityESDCorrected[i]->Merge(&collections[kESDHists+kMCHists*4+kCorrHists+i]);
+ fMultiplicityESDCorrected[i]->Merge(&collections[3*kESDHists+kMCHists*4+kCorrHists+i]);
delete iter;
TList oldObjects;
oldObjects.SetOwner(1);
for (Int_t i = 0; i < kESDHists; ++i)
+ {
if (fMultiplicityESD[i])
oldObjects.Add(fMultiplicityESD[i]);
-
+ if (fTriggeredEvents[i])
+ oldObjects.Add(fTriggeredEvents[i]);
+ if (fNoVertexEvents[i])
+ oldObjects.Add(fNoVertexEvents[i]);
+ }
+
for (Int_t i = 0; i < kMCHists; ++i)
{
if (fMultiplicityVtx[i])
for (Int_t i = 0; i < kESDHists; ++i)
{
fMultiplicityESD[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityESD[i]->GetName()));
- if (!fMultiplicityESD[i])
+ fTriggeredEvents[i] = dynamic_cast<TH1F*> (gDirectory->Get(fTriggeredEvents[i]->GetName()));
+ fNoVertexEvents[i] = dynamic_cast<TH1F*> (gDirectory->Get(fNoVertexEvents[i]->GetName()));
+ if (!fMultiplicityESD[i] || !fTriggeredEvents[i] || !fNoVertexEvents[i])
success = kFALSE;
}
gDirectory->cd(dir);
for (Int_t i = 0; i < kESDHists; ++i)
+ {
if (fMultiplicityESD[i])
{
fMultiplicityESD[i]->Write();
fMultiplicityESD[i]->ProjectionY(Form("%s_px", fMultiplicityESD[i]->GetName()), 1, fMultiplicityESD[i]->GetNbinsX())->Write();
}
+ if (fTriggeredEvents[i])
+ fTriggeredEvents[i]->Write();
+ if (fNoVertexEvents[i])
+ fNoVertexEvents[i]->Write();
+ }
for (Int_t i = 0; i < kMCHists; ++i)
{
fMultiplicityESD[2]->Fill(vtx, measured14);
}
+//____________________________________________________________________
+void AliMultiplicityCorrection::FillTriggeredEvent(Int_t measured05, Int_t measured10, Int_t measured14)
+{
+ //
+ // fills raw distribution of triggered events
+ //
+
+ fTriggeredEvents[0]->Fill(measured05);
+ fTriggeredEvents[1]->Fill(measured10);
+ fTriggeredEvents[2]->Fill(measured14);
+}
+
+//____________________________________________________________________
+void AliMultiplicityCorrection::FillNoVertexEvent(Float_t vtx, Bool_t vertexReconstructed, Int_t generated05, Int_t generated10, Int_t generated14, Int_t measured05, Int_t measured10, Int_t measured14)
+{
+ //
+ // fills raw distribution of triggered events
+ //
+
+ if (vtx > fgVtxRangeBegin[0] && vtx < fgVtxRangeEnd[0] && (!vertexReconstructed || measured05 == 0))
+ fNoVertexEvents[0]->Fill(generated05);
+
+ if (vtx > fgVtxRangeBegin[1] && vtx < fgVtxRangeEnd[1] && (!vertexReconstructed || measured10 == 0))
+ fNoVertexEvents[1]->Fill(generated10);
+
+ if (vtx > fgVtxRangeBegin[2] && vtx < fgVtxRangeEnd[2] && (!vertexReconstructed || measured14 == 0))
+ fNoVertexEvents[2]->Fill(generated14);
+}
+
//____________________________________________________________________
void AliMultiplicityCorrection::FillCorrection(Float_t vtx, Int_t generated05, Int_t generated10, Int_t generated14, Int_t generatedAll, Int_t measured05, Int_t measured10, Int_t measured14)
{
fMultiplicityESDCorrected[correlationID]->Sumw2();
// project without under/overflow bins
- fCurrentESD = fMultiplicityESD[inputRange]->ProjectionY("fCurrentESD", 1, fMultiplicityESD[inputRange]->GetXaxis()->GetNbins());
+ Int_t begin = 1;
+ Int_t end = fMultiplicityESD[inputRange]->GetXaxis()->GetNbins();
+ if (fVtxEnd > fVtxBegin)
+ {
+ begin = fVtxBegin;
+ end = fVtxEnd;
+ }
+ fCurrentESD = fMultiplicityESD[inputRange]->ProjectionY("fCurrentESD", begin, end);
fCurrentESD->Sumw2();
// empty under/overflow bins in x, otherwise Project3D takes them into account
}
}
+ if (fVtxEnd > fVtxBegin)
+ hist->GetXaxis()->SetRange(fVtxBegin, fVtxEnd);
+
fCurrentCorrelation = (TH2*) hist->Project3D("zy");
fCurrentCorrelation->Sumw2();
// calculates efficiency for given event type
//
+ TString name1;
+ name1.Form("divisor%d", inputRange);
+
+ TString name2;
+ name2.Form("CurrentEfficiency%d", inputRange);
+
TH1* divisor = 0;
switch (eventType)
{
case kTrVtx : break;
- case kMB: divisor = fMultiplicityMB[inputRange]->ProjectionY("divisor", 1, fMultiplicityMB[inputRange]->GetNbinsX(), "e"); break;
- case kINEL: divisor = fMultiplicityINEL[inputRange]->ProjectionY("divisor", 1, fMultiplicityINEL[inputRange]->GetNbinsX(), "e"); break;
- case kNSD: divisor = fMultiplicityNSD[inputRange]->ProjectionY("divisor", 1, fMultiplicityNSD[inputRange]->GetNbinsX(), "e"); break;
+ case kMB: divisor = fMultiplicityMB[inputRange]->ProjectionY(name1, 1, fMultiplicityMB[inputRange]->GetNbinsX(), "e"); break;
+ case kINEL: divisor = fMultiplicityINEL[inputRange]->ProjectionY(name1, 1, fMultiplicityINEL[inputRange]->GetNbinsX(), "e"); break;
+ case kNSD: divisor = fMultiplicityNSD[inputRange]->ProjectionY(name1, 1, fMultiplicityNSD[inputRange]->GetNbinsX(), "e"); break;
}
- TH1* eff = fMultiplicityVtx[inputRange]->ProjectionY("CurrentEfficiency", 1, fMultiplicityVtx[inputRange]->GetNbinsX(), "e");
+ TH1* eff = fMultiplicityVtx[inputRange]->ProjectionY(name2, 1, fMultiplicityVtx[inputRange]->GetNbinsX(), "e");
if (eventType == kTrVtx)
{
}
//____________________________________________________________________
-TH1* AliMultiplicityCorrection::GetTriggerEfficiency(Int_t inputRange)
+TH1* AliMultiplicityCorrection::GetTriggerEfficiency(Int_t inputRange, EventType eventType)
{
//
// calculates efficiency for given event type
//
- TH1* divisor = fMultiplicityINEL[inputRange]->ProjectionY("divisor", 1, fMultiplicityINEL[inputRange]->GetNbinsX(), "e");
- TH1* eff = fMultiplicityMB[inputRange]->ProjectionY("CurrentEfficiency", 1, fMultiplicityMB[inputRange]->GetNbinsX(), "e");
+ TString name1;
+ name1.Form("divisor%d", inputRange);
+
+ TString name2;
+ name2.Form("CurrentEfficiency%d", inputRange);
+
+ TH1* divisor = 0;
+ switch (eventType)
+ {
+ case kTrVtx : AliFatal("Not supported!"); break;
+ case kMB: divisor = fMultiplicityMB[inputRange]->ProjectionY (name1, 1, fMultiplicityMB[inputRange]->GetNbinsX(), "e"); break;
+ case kINEL: divisor = fMultiplicityINEL[inputRange]->ProjectionY(name1, 1, fMultiplicityINEL[inputRange]->GetNbinsX(), "e"); break;
+ case kNSD: divisor = fMultiplicityNSD[inputRange]->ProjectionY (name1, 1, fMultiplicityNSD[inputRange]->GetNbinsX(), "e"); break;
+ }
+ TH1* eff = fMultiplicityMB[inputRange]->ProjectionY(name2, 1, fMultiplicityMB[inputRange]->GetNbinsX(), "e");
+
eff->Divide(eff, divisor, 1, 1, "B");
return eff;
}
//____________________________________________________________________
-Int_t AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check, TH1* initialConditions)
+Int_t AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Int_t zeroBinEvents, Bool_t check, TH1* initialConditions, Bool_t errorAsBias)
{
//
// correct spectrum using minuit chi2 method
Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
- AliUnfolding::SetCreateOverflowBin(5);
+ //AliUnfolding::SetCreateOverflowBin(5);
AliUnfolding::SetUnfoldingMethod(AliUnfolding::kChi2Minimization);
- SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
+ AliUnfolding::SetMinimumInitialValue(kTRUE, 0.1);
+
+ // use here only vtx efficiency (to MB sample) which is always needed if we use the 0 bin
+ SetupCurrentHists(inputRange, fullPhaseSpace, (eventType == kTrVtx) ? kTrVtx : kMB);
+
+ // TODO set errors on measured with 0.5 * TMath::ChisquareQuantile(0.1, 20000)
+ // see PDG: Statistics / Poission or binomial data / Eq. 32.49a/b in 2004 edition
+
+ Calculate0Bin(inputRange, eventType, zeroBinEvents);
+
+ Int_t resultCode = -1;
+ if (errorAsBias == kFALSE)
+ {
+ resultCode = AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID], check);
+ }
+ else
+ {
+ resultCode = AliUnfolding::UnfoldGetBias(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID]);
+ }
+
+ // HACK store new vertex reco efficiency for bin 0, changing number of events with trigger and vertex in MC map
+ if (zeroBinEvents > 0)
+ {
+ Printf("WARNING: Stored vertex reco efficiency from unfolding for bin 0.");
+ fMultiplicityVtx[inputRange]->SetBinContent(1, 1, fMultiplicityMB[inputRange]->GetBinContent(1, 1) * fCurrentEfficiency->GetBinContent(1));
+ }
+
+ // correct for the trigger bias if requested
+ if (eventType > kMB)
+ {
+ Printf("Applying trigger efficiency");
+ TH1* eff = GetTriggerEfficiency(inputRange, eventType);
+ for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); i++)
+ {
+ fMultiplicityESDCorrected[correlationID]->SetBinContent(i, fMultiplicityESDCorrected[correlationID]->GetBinContent(i) / eff->GetBinContent(i));
+ fMultiplicityESDCorrected[correlationID]->SetBinError(i, fMultiplicityESDCorrected[correlationID]->GetBinError(i) / eff->GetBinContent(i));
+ }
+ }
+
+ return resultCode;
+}
+
+//____________________________________________________________________
+void AliMultiplicityCorrection::Calculate0Bin(Int_t inputRange, EventType eventType, Int_t zeroBinEvents)
+{
+ // fills the 0 bin
- if (!initialConditions)
+ if (eventType == kTrVtx)
+ return;
+
+ Double_t fractionEventsInVertexRange = fMultiplicityESD[inputRange]->Integral(1, fMultiplicityESD[inputRange]->GetXaxis()->GetNbins()) / fMultiplicityESD[inputRange]->Integral(0, fMultiplicityESD[inputRange]->GetXaxis()->GetNbins() + 1);
+
+ // difference of fraction that is inside the considered range between triggered events and events with vertex
+ // Extension to NSD not needed, INEL and NSD vertex distributions are nature-given and unbiased!
+ Double_t differenceVtxDist = (fMultiplicityINEL[inputRange]->Integral(1, fMultiplicityINEL[inputRange]->GetXaxis()->GetNbins()) / fMultiplicityINEL[inputRange]->Integral(0, fMultiplicityINEL[inputRange]->GetXaxis()->GetNbins() + 1)) / (fMultiplicityVtx[inputRange]->Integral(1, fMultiplicityVtx[inputRange]->GetXaxis()->GetNbins()) / fMultiplicityVtx[inputRange]->Integral(0, fMultiplicityVtx[inputRange]->GetXaxis()->GetNbins() + 1));
+
+ Printf("Enabling 0 bin estimate for triggered events without vertex.");
+ Printf(" Events in 0 bin: %d", zeroBinEvents);
+ Printf(" Fraction in range: %.1f%%", fractionEventsInVertexRange * 100);
+ Printf(" Difference Vtx Dist: %f", differenceVtxDist);
+
+ AliUnfolding::SetNotFoundEvents(zeroBinEvents * fractionEventsInVertexRange * differenceVtxDist);
+}
+
+//____________________________________________________________________
+void AliMultiplicityCorrection::FixTriggerEfficiencies(Int_t start)
+{
+ //
+ // sets trigger and vertex efficiencies to 1 for large multiplicities where no event was simulated
+ //
+
+ for (Int_t etaRange = 0; etaRange < kMCHists; etaRange++)
{
- initialConditions = (TH1*) fCurrentESD->Clone("initialConditions");
- initialConditions->Scale(1.0 / initialConditions->Integral());
- if (!check)
+ for (Int_t i = 1; i <= fMultiplicityINEL[etaRange]->GetNbinsY(); i++)
{
- // set minimum value to prevent MINUIT just staying in the small value
- for (Int_t i=1; i<=initialConditions->GetNbinsX(); i++)
- initialConditions->SetBinContent(i, TMath::Max(initialConditions->GetBinContent(i), 1e-3));
+ if (fMultiplicityINEL[etaRange]->GetYaxis()->GetBinCenter(i) < start)
+ continue;
+
+ if (fMultiplicityINEL[etaRange]->GetBinContent(1, i) > 0)
+ continue;
+
+ fMultiplicityINEL[etaRange]->SetBinContent(1, i, 1);
+ fMultiplicityNSD[etaRange]->SetBinContent(1, i, 1);
+ fMultiplicityMB[etaRange]->SetBinContent(1, i, 1);
+ fMultiplicityVtx[etaRange]->SetBinContent(1, i, 1);
}
}
+}
- return AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID], check);
+//____________________________________________________________________
+Float_t AliMultiplicityCorrection::GetFraction0Generated(Int_t inputRange)
+{
+ //
+ // returns the fraction of events that have 0 generated particles in the given range among all events without vertex OR 0 tracklets
+ //
+
+ TH1* multMB = GetNoVertexEvents(inputRange);
+ return multMB->GetBinContent(1) / multMB->Integral();
}
//____________________________________________________________________
void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRange, Bool_t fullPhaseSpace, Bool_t /*normalizeESD*/, TH1* mcHist, Bool_t simple, EventType eventType)
{
// draw comparison plots
-
-
- //mcHist->Rebin(2);
-
+
Int_t esdCorrId = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
-
+
TString tmpStr;
tmpStr.Form("%s_DrawComparison_%d", name, esdCorrId);
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", 1, fMultiplicityESD[inputRange]->GetXaxis()->GetNbins());
- 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);*/
-
- // find bin with <= 50 entries. this is used as limit for the chi2 calculation
- Int_t mcBinLimit = 0;
- for (Int_t i=20; i<mcHist->GetNbinsX(); ++i)
+
+ Int_t begin = 1;
+ Int_t end = fMultiplicityESD[inputRange]->GetXaxis()->GetNbins();
+ if (fVtxEnd > fVtxBegin)
{
- if (mcHist->GetBinContent(i) > 50)
- {
- mcBinLimit = i;
- }
- else
- break;
+ begin = fVtxBegin;
+ end = fVtxEnd;
}
- Printf("AliMultiplicityCorrection::DrawComparison: MC bin limit is %d", mcBinLimit);
-
- // scale to 1
- mcHist->Sumw2();
- if (mcHist->Integral() > 0)
- mcHist->Scale(1.0 / mcHist->Integral());
-
- // calculate residual
+ fCurrentESD = fMultiplicityESD[esdCorrId]->ProjectionY("fCurrentESD", begin, end);
+ fCurrentESD->Sumw2();
- // 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");
-
- // undo trigger, vertex efficiency correction
- fCurrentEfficiency = GetEfficiency(inputRange, eventType);
- tmpESDEfficiencyRecorrected->Multiply(fCurrentEfficiency);
-
- TH2* convoluted = CalculateMultiplicityESD(tmpESDEfficiencyRecorrected, esdCorrId);
- TH1* convolutedProj = convoluted->ProjectionY("convolutedProj", -1, -1, "e");
- if (convolutedProj->Integral() > 0)
+ mcHist->Sumw2();
+ Int_t mcMax = 0;
+ for (Int_t i=5; i<=mcHist->GetNbinsX(); ++i)
{
- convolutedProj->Scale(1.0 / convolutedProj->Integral());
+ if (mcHist->GetBinContent(i) > 0)
+ mcMax = mcHist->GetXaxis()->GetBinCenter(i) + 2;
}
- else
- printf("ERROR: convolutedProj is empty. Something went wrong calculating the convoluted histogram.\n");
-
- 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");
+ if (mcMax == 0)
+ {
+ for (Int_t i=5; i<=fMultiplicityESDCorrected[esdCorrId]->GetNbinsX(); ++i)
+ if (fMultiplicityESDCorrected[esdCorrId]->GetBinContent(i) > 1)
+ mcMax = fMultiplicityESDCorrected[esdCorrId]->GetXaxis()->GetBinCenter(i) + 2;
+ }
+ Printf("AliMultiplicityCorrection::DrawComparison: MC bin limit is %d", mcMax);
+ // calculate residual
+ Float_t tmp;
+ TH1* convolutedProj = (TH1*) GetConvoluted(esdCorrId, eventType)->Clone("convolutedProj");
+ TH1* residual = GetResiduals(esdCorrId, eventType, tmp);
TH1* residualHist = new TH1F("residualHist", "residualHist", 51, -5, 5);
- // find bin limit
- Int_t lastBin = 0;
- for (Int_t i=1; i<=fCurrentESD->GetNbinsX(); ++i)
- {
- if (fCurrentESD->GetBinContent(i) <= 5)
- {
- lastBin = i;
- break;
- }
- }
-
- // TODO fix errors
Float_t chi2 = 0;
- for (Int_t i=1; i<=residual->GetNbinsX(); ++i)
+ for (Int_t i=1; i<=TMath::Min(residual->GetNbinsX(), 75); ++i)
{
- if (fCurrentESD->GetBinError(i) > 0)
- {
- Float_t value = residual->GetBinContent(i) / fCurrentESD->GetBinError(i);
- if (i > 1 && i <= lastBin)
- chi2 += value * value;
- Printf("%d --> %f (%f)", i, value * value, chi2);
- residual->SetBinContent(i, value);
- residualHist->Fill(value);
- }
- else
- {
- //printf("Residual bin %d set to 0\n", i);
- residual->SetBinContent(i, 0);
- }
+ Float_t value = residual->GetBinContent(i);
+ // TODO has to get a parameter (used in Chi2Function, GetResiduals, and here)
+ if (i > 1)
+ chi2 += value * value;
+ Printf("%d --> %f (%f)", i, value * value, chi2);
+ residualHist->Fill(value);
convolutedProj->SetBinError(i, 0);
- residual->SetBinError(i, 0);
}
fLastChi2Residuals = chi2;
- new TCanvas; residualHist->DrawCopy();
-
- //residualHist->Fit("gaus", "N");
- //delete residualHist;
+ //new TCanvas; residualHist->DrawCopy();
- printf("Difference (Residuals) is %f for bin 2-%d\n", fLastChi2Residuals, lastBin);
+ printf("Difference (Residuals) is %f\n", fLastChi2Residuals);
TCanvas* canvas1 = 0;
if (simple)
canvas1->cd(1)->SetLeftMargin(0.12);
canvas1->cd(1)->SetBottomMargin(0.12);
TH1* proj = (TH1*) mcHist->Clone("proj");
+ if (proj->GetEntries() > 0)
+ AliPWG0Helper::NormalizeToBinWidth(proj);
- // normalize without 0 bin
- proj->Scale(1.0 / proj->Integral(2, proj->GetNbinsX()));
- Printf("Normalized without 0 bin!");
- proj->GetXaxis()->SetRangeUser(0, 200);
+ proj->GetXaxis()->SetRangeUser(0, mcMax);
proj->GetYaxis()->SetTitleOffset(1.4);
- //proj->SetLabelSize(0.05, "xy");
- //proj->SetTitleSize(0.05, "xy");
proj->SetTitle(Form(";True multiplicity in |#eta| < %.1f;Entries", (inputRange+1)*0.5));
proj->SetStats(kFALSE);
+ fMultiplicityESDCorrected[esdCorrId]->GetXaxis()->SetRangeUser(0, mcMax);
+ fMultiplicityESDCorrected[esdCorrId]->GetYaxis()->SetRangeUser(0.1, fMultiplicityESDCorrected[esdCorrId]->GetMaximum() * 1.5);
+ fMultiplicityESDCorrected[esdCorrId]->GetYaxis()->SetTitleOffset(1.4);
+ fMultiplicityESDCorrected[esdCorrId]->SetTitle(Form(";True multiplicity in |#eta| < %.1f;Entries", (inputRange+1)*0.5));
+ fMultiplicityESDCorrected[esdCorrId]->SetStats(kFALSE);
+
fMultiplicityESDCorrected[esdCorrId]->SetLineColor(2);
fMultiplicityESDCorrected[esdCorrId]->SetMarkerColor(2);
- //fMultiplicityESDCorrected[esdCorrId]->SetMarkerStyle(5);
- // normalize without 0 bin
- fMultiplicityESDCorrected[esdCorrId]->Scale(1.0 / fMultiplicityESDCorrected[esdCorrId]->Integral(2, fMultiplicityESDCorrected[esdCorrId]->GetNbinsX()));
- Printf("Normalized without 0 bin!");
+ TH1* esdCorrected = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("esdCorrected");
+ AliPWG0Helper::NormalizeToBinWidth(esdCorrected);
+
if (proj->GetEntries() > 0) {
proj->DrawCopy("HIST");
- fMultiplicityESDCorrected[esdCorrId]->DrawCopy("SAME HIST E");
+ esdCorrected->DrawCopy("SAME HIST E");
}
else
- fMultiplicityESDCorrected[esdCorrId]->DrawCopy("HIST E");
-
+ esdCorrected->DrawCopy("HIST E");
+
gPad->SetLogy();
TLegend* legend = new TLegend(0.3, 0.8, 0.93, 0.93);
legend->SetFillColor(0);
legend->SetTextSize(0.04);
legend->Draw();
- // unfortunately does not work. maybe a bug? --> legend->SetTextSizePixels(14);
canvas1->cd(2);
canvas1->cd(2)->SetGridx();
canvas1->cd(2)->SetBottomMargin(0.12);
gPad->SetLogy();
- fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
- //fCurrentESD->SetLineColor(2);
+ fCurrentESD->GetXaxis()->SetRangeUser(0, mcMax);
fCurrentESD->SetTitle(Form(";Measured multiplicity in |#eta| < %.1f;Entries", (inputRange+1)*0.5));
fCurrentESD->SetStats(kFALSE);
fCurrentESD->GetYaxis()->SetTitleOffset(1.4);
- //fCurrentESD->SetLabelSize(0.05, "xy");
- //fCurrentESD->SetTitleSize(0.05, "xy");
fCurrentESD->DrawCopy("HIST E");
convolutedProj->SetLineColor(2);
convolutedProj->SetMarkerColor(2);
convolutedProj->SetMarkerStyle(5);
- //proj2->SetMarkerColor(2);
- //proj2->SetMarkerStyle(5);
convolutedProj->DrawCopy("HIST SAME P");
legend = new TLegend(0.3, 0.8, 0.93, 0.93);
legend->SetTextSize(0.04);
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);
-
- //printf("%d %f\n", i, chi);
-
- if (chi2 < 0.2)
- fLastChi2MCLimit = i;
-
- if (chi < 3)
- limit = i;
-
- }
- }*/
-
- /*chi2 = 0;
- Float_t chi = 0;
- Int_t limit = 0;
- for (Int_t i=1; i<=diffMCUnfolded->GetNbinsX(); ++i)
- {
- if (fMultiplicityESDCorrected[esdCorrId]->GetBinError(i) > 0)
- {
- 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("limits %d %d\n", fLastChi2MCLimit, limit);
- fLastChi2MCLimit = limit;
-
- printf("Difference (from MC) is %f for bin 2-150. Limit is %d.\n", fLastChi2MC, fLastChi2MCLimit);*/
-
if (!simple)
{
- /*canvas1->cd(3);
-
- diffMCUnfolded->SetTitle("#chi^{2};Npart;(MC - Unfolded) / e(unfolded)");
- //diffMCUnfolded->GetYaxis()->SetRangeUser(-20, 20);
- diffMCUnfolded->GetXaxis()->SetRangeUser(0, 200);
- diffMCUnfolded->DrawCopy("HIST");
-
- TH1F* fluctuation = new TH1F("fluctuation", "fluctuation", 20, -5, 5);
- for (Int_t i=20; i<=diffMCUnfolded->GetNbinsX(); ++i)
- fluctuation->Fill(diffMCUnfolded->GetBinContent(i));
-
- //new TCanvas; fluctuation->DrawCopy();
- delete fluctuation;*/
-
- /*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(-5, 5);
- residual->GetXaxis()->SetRangeUser(0, 200);
+ residual->GetXaxis()->SetRangeUser(0, mcMax);
+ residual->SetStats(kFALSE);
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->GetXaxis()->SetRangeUser(0, mcMax);
ratio->SetStats(kFALSE);
ratio->Draw("HIST");
- Double_t ratioChi2 = 0;
- fRatioAverage = 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;
- fRatioAverage += TMath::Abs(value);
-
- if (ratioChi2 < 0.1)
- fLastChi2MCLimit = i;
- }
- fRatioAverage /= 149;
-
- printf("Sum over (ratio-1)^2 (2..150) is %f; average of |ratio-1| is %f\n", ratioChi2, fRatioAverage);
- // TODO FAKE
- fLastChi2MC = ratioChi2;
-
- // FFT of ratio
- canvas1->cd(6);
- const Int_t kFFT = 128;
- Double_t fftReal[kFFT];
- Double_t fftImag[kFFT];
-
- for (Int_t i=0; i<kFFT; ++i)
- {
- fftReal[i] = ratio->GetBinContent(i+1+10);
- // test: ;-)
- //fftReal[i] = cos(TMath::Pi() * 5 * 2 * i / 128);
- fftImag[i] = 0;
- }
-
- FFT(-1, TMath::Nint(TMath::Log(kFFT) / TMath::Log(2)), fftReal, fftImag);
-
- TH1* fftResult = (TH1*) ratio->Clone("fftResult");
- fftResult->SetTitle("FFT;true multiplicity;coeff. (10...137)");
- TH1* fftResult2 = (TH1*) ratio->Clone("fftResult2");
- TH1* fftResult3 = (TH1*) ratio->Clone("fftResult3");
- fftResult->Reset();
- fftResult2->Reset();
- fftResult3->Reset();
- fftResult->GetYaxis()->UnZoom();
- fftResult2->GetYaxis()->UnZoom();
- fftResult3->GetYaxis()->UnZoom();
-
- //Printf("AFTER FFT");
- for (Int_t i=0; i<kFFT; ++i)
- {
- //Printf("%d: %f", i, fftReal[i]);
- fftResult->SetBinContent(i+1, fftReal[i]);
- /*if (i != 0 && TMath::Abs(fftReal[i]) > 0.5)
- {
- Printf("Nulled %d", i);
- fftReal[i] = 0;
- }*/
- }
-
- fftResult->SetLineColor(1);
- fftResult->DrawCopy();
-
-
- //Printf("IMAG");
- for (Int_t i=0; i<kFFT; ++i)
- {
- //Printf("%d: %f", i, fftImag[i]);
- fftResult2->SetBinContent(i+1, fftImag[i]);
-
- fftResult3->SetBinContent(i+1, TMath::Sqrt(fftReal[i] * fftReal[i] + fftImag[i] * fftImag[i]));
- }
-
- fftResult2->SetLineColor(2);
- fftResult2->DrawCopy("SAME");
-
- fftResult3->SetLineColor(4);
- fftResult3->DrawCopy("SAME");
-
- for (Int_t i=1; i<kFFT - 1; ++i)
- {
- if (TMath::Sqrt(fftReal[i] * fftReal[i] + fftImag[i] * fftImag[i]) > 3)
- {
- fftReal[i-1] = 0;
- fftReal[i] = 0;
- fftReal[i+1] = 0;
- fftImag[i-1] = 0;
- fftImag[i] = 0;
- fftImag[i+1] = 0;
- //fftReal[i] = (fftReal[i-1] + fftReal[i+1]) / 2;
- //fftImag[i] = (fftImag[i-1] + fftImag[i+1]) / 2;
- //Printf("Nulled %d to %f %f", i, fftReal[i], fftImag[i]);
- }
- }
-
- FFT(1, TMath::Nint(TMath::Log(kFFT) / TMath::Log(2)), fftReal, fftImag);
-
- TH1* fftResult4 = (TH1*) fftResult3->Clone("fftResult4");
- fftResult4->Reset();
-
- //Printf("AFTER BACK-TRAFO");
- for (Int_t i=0; i<kFFT; ++i)
- {
- //Printf("%d: %f", i, fftReal[i]);
- fftResult4->SetBinContent(i+1+10, fftReal[i]);
- }
-
- canvas1->cd(5);
- fftResult4->SetLineColor(4);
- fftResult4->DrawCopy("SAME");
-
// plot (MC - Unfolded) / error (MC)
canvas1->cd(3);
TH1* diffMCUnfolded2 = dynamic_cast<TH1*> (proj->Clone("diffMCUnfolded2"));
- diffMCUnfolded2->Add(fMultiplicityESDCorrected[esdCorrId], -1);
+ diffMCUnfolded2->Add(esdCorrected, -1);
Int_t ndfQual[kQualityRegions];
for (Int_t region=0; region<kQualityRegions; ++region)
ndfQual[region] = 0;
}
-
Double_t newChi2 = 0;
Double_t newChi2Limit150 = 0;
Int_t ndf = 0;
{
value = diffMCUnfolded2->GetBinContent(i) / proj->GetBinError(i);
newChi2 += value * value;
- if (i > 1 && i <= mcBinLimit)
+ if (i > 1 && i <= mcMax)
newChi2Limit150 += value * value;
++ndf;
Printf("Quality parameter %d (%d <= mult <= %d) is %f with %d df", region, fgQualityRegionsB[region], fgQualityRegionsE[region], fQuality[region], ndfQual[region]);
}
- if (mcBinLimit > 1)
+ if (mcMax > 1)
{
- // TODO another FAKE
- fLastChi2MC = newChi2Limit150 / (mcBinLimit - 1);
- Printf("Chi2 (2..%d) from (MC - Unfolded) / e(MC) is: %.2f ndf is %d --> chi2 / ndf = %.2f", mcBinLimit, newChi2Limit150, mcBinLimit - 1, fLastChi2MC);
+ fLastChi2MC = newChi2Limit150 / (mcMax - 1);
+ Printf("Chi2 (2..%d) from (MC - Unfolded) / e(MC) is: %.2f ndf is %d --> chi2 / ndf = %.2f", mcMax, newChi2Limit150, mcMax - 1, fLastChi2MC);
}
else
fLastChi2MC = -1;
Printf("Chi2 (full range) from (MC - Unfolded) / e(MC) is: %.2f ndf is %d --> chi2 / ndf = %.2f", newChi2, ndf, ((ndf > 0) ? newChi2 / ndf : -1));
- diffMCUnfolded2->SetTitle("#chi^{2};Npart;(MC - Unfolded) / e(MC)");
+ diffMCUnfolded2->SetTitle("#chi^{2};true multiplicity;(MC - Unfolded) / e(MC)");
diffMCUnfolded2->GetYaxis()->SetRangeUser(-5, 5);
- diffMCUnfolded2->GetXaxis()->SetRangeUser(0, 200);
+ diffMCUnfolded2->GetXaxis()->SetRangeUser(0, mcMax);
diffMCUnfolded2->DrawCopy("HIST");
+
+ canvas1->cd(6);
+ // draw penalty factor
+
+ TH1* penalty = AliUnfolding::GetPenaltyPlot(fMultiplicityESDCorrected[esdCorrId]);
+ penalty->SetStats(0);
+ penalty->GetXaxis()->SetRangeUser(0, mcMax);
+ penalty->DrawCopy("HIST");
}
-
- canvas1->SaveAs(Form("%s.gif", canvas1->GetName()));
}
//____________________________________________________________________
}
//____________________________________________________________________
-TH1* AliMultiplicityCorrection::StatisticalUncertainty(AliUnfolding::MethodType methodType, Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t randomizeMeasured, Bool_t randomizeResponse, const TH1* compareTo)
+TH1* AliMultiplicityCorrection::StatisticalUncertainty(AliUnfolding::MethodType methodType, Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Int_t zeroBinEvents, Bool_t randomizeMeasured, Bool_t randomizeResponse, const TH1* compareTo)
{
//
// evaluates the uncertainty that arises from the non-infinite statistics in the response matrix
gRandom->SetSeed(0);
if (methodType == AliUnfolding::kChi2Minimization)
- AliUnfolding::SetCreateOverflowBin(5);
+ {
+ Calculate0Bin(inputRange, eventType, zeroBinEvents);
+ AliUnfolding::SetMinimumInitialValue(kTRUE, 0.1);
+ }
+
AliUnfolding::SetUnfoldingMethod(methodType);
- const Int_t kErrorIterations = 150;
+ const Int_t kErrorIterations = 20;
TH1* maxError = 0;
TH1* firstResult = 0;
Int_t returnCode = AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, measured, 0, result);
if (returnCode != 0)
- return 0;
+ {
+ n--;
+ continue;
+ }
}
// normalize
fMultiplicityESDCorrected[correlationID]->SetBinContent(i, firstResult->GetBinContent(i));
for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
- fMultiplicityESDCorrected[correlationID]->SetBinError(i, maxError->GetBinContent(i) * fMultiplicityESDCorrected[correlationID]->GetBinContent(i));
+ fMultiplicityESDCorrected[correlationID]->SetBinError(i, standardDeviation->GetBinContent(i) * fMultiplicityESDCorrected[correlationID]->GetBinContent(i));
return standardDeviation;
}
//____________________________________________________________________
-void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar, Int_t nIterations, TH1* initialConditions, Bool_t determineError)
+void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar, Int_t nIterations, TH1* initialConditions, Int_t determineError)
{
//
// correct spectrum using bayesian method
//
+ // determineError:
+ // 0 = no errors
+ // 1 = from randomizing
+ // 2 = with UnfoldGetBias
// initialize seed with current time
gRandom->SetSeed(0);
AliUnfolding::SetBayesianParameters(regPar, nIterations);
AliUnfolding::SetUnfoldingMethod(AliUnfolding::kBayesian);
- if (AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID]) != 0)
+
+ if (determineError <= 1)
+ {
+ if (AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID]) != 0)
+ return;
+ }
+ else if (determineError == 2)
+ {
+ AliUnfolding::UnfoldGetBias(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID]);
return;
+ }
- if (!determineError)
+ if (determineError == 0)
{
Printf("AliMultiplicityCorrection::ApplyBayesianMethod: WARNING: No errors calculated.");
return;
const Int_t kErrorIterations = 20;
- printf("Spectrum unfolded. Determining error (%d iterations)...\n", kErrorIterations);
+ Printf("Spectrum unfolded. Determining error (%d iterations)...", kErrorIterations);
TH1* randomized = (TH1*) fCurrentESD->Clone("randomized");
TH1* resultArray[kErrorIterations+1];
TH1* result2 = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone("result2");
result2->Reset();
if (AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, randomized, initialConditions, result2) != 0)
- return;
-
- result2->Scale(1.0 / result2->Integral());
+ {
+ n--;
+ continue;
+ }
resultArray[n+1] = result2;
}
for (Int_t n=0; n<kErrorIterations; ++n)
delete resultArray[n+1];
+ Printf("Comparing bias and error:");
for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
+ {
+ Printf("Bin %d: Content: %f Error: %f Bias: %f", i, fMultiplicityESDCorrected[correlationID]->GetBinContent(i), error->GetBinContent(i) * fMultiplicityESDCorrected[correlationID]->GetBinContent(i), fMultiplicityESDCorrected[correlationID]->GetBinError(i));
fMultiplicityESDCorrected[correlationID]->SetBinError(i, error->GetBinContent(i) * fMultiplicityESDCorrected[correlationID]->GetBinContent(i));
+ }
delete error;
}
DrawComparison("Gaussian", inputRange, fullPhaseSpace, kFALSE, GetMultiplicityMC(mcTarget, kTrVtx)->ProjectionY());
}
+//____________________________________________________________________
+TH1* AliMultiplicityCorrection::GetConvoluted(Int_t i, EventType eventType)
+{
+ // convolutes the corrected histogram i with the response matrix
+
+ TH1* corrected = (TH1*) fMultiplicityESDCorrected[i]->Clone("corrected");
+
+ // undo efficiency correction (hist must not be deleted, is reused)
+ TH1* efficiency = GetEfficiency(i, eventType);
+ //new TCanvas; efficiency->DrawCopy();
+ corrected->Multiply(efficiency);
+
+ TH2* convoluted = CalculateMultiplicityESD(corrected, i);
+ TH1* convolutedProj = convoluted->ProjectionY("GetConvoluted_convolutedProj", 1, convoluted->GetNbinsX());
+
+ delete convoluted;
+ delete corrected;
+
+ return convolutedProj;
+}
+
+//____________________________________________________________________
+TH1* AliMultiplicityCorrection::GetResiduals(Int_t i, EventType eventType, Float_t& residualSum)
+{
+ // creates the residual histogram from the corrected histogram i corresponding to an eventType event sample using the corresponding correlation matrix
+ // residual is : M - UT / eM
+ // residualSum contains the squared sum of the residuals
+
+ TH1* corrected = (TH1*) fMultiplicityESDCorrected[i]->Clone("corrected");
+ TH1* convolutedProj = GetConvoluted(i, eventType);
+
+ Int_t begin = 1;
+ Int_t end = fMultiplicityESD[i]->GetNbinsX();
+ if (fVtxEnd > fVtxBegin)
+ {
+ begin = fVtxBegin;
+ end = fVtxEnd;
+ }
+ TH1* measuredProj = fMultiplicityESD[i]->ProjectionY("measuredProj", begin, end);
+
+ TH1* residuals = (TH1*) measuredProj->Clone("GetResiduals_residuals");
+ residuals->SetTitle(";measured multiplicity;residuals (M-Ut)/e");
+ residuals->Add(convolutedProj, -1);
+
+ residualSum = 0;
+ for (Int_t i=1; i<=residuals->GetNbinsX(); i++)
+ {
+ if (measuredProj->GetBinContent(i) > 0)
+ residuals->SetBinContent(i, residuals->GetBinContent(i) / TMath::Sqrt(measuredProj->GetBinContent(i)));
+ residuals->SetBinError(i, 0);
+
+ if (i > 1)
+ residualSum += residuals->GetBinContent(i) * residuals->GetBinContent(i);
+ }
+
+ delete corrected;
+ delete convolutedProj;
+ delete measuredProj;
+
+ return residuals;
+}
+
//____________________________________________________________________
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
// although it is a TH2F the vertex axis is not used at the moment and all entries are filled in mid-vertex
- // if normalized is set, inputMC is expected to be normalized to the bin width
if (!inputMC)
return 0;
}
}
+ if (fVtxEnd > fVtxBegin)
+ hist->GetXaxis()->SetRange(fVtxBegin, fVtxEnd);
+
TH2* corr = (TH2*) hist->Project3D("zy");
//corr->Rebin2D(2, 1);
corr->Sumw2();
Int_t mcGenBin = inputMC->GetXaxis()->FindBin(corr->GetXaxis()->GetBinCenter(gen));
measured += inputMC->GetBinContent(mcGenBin) * corr->GetBinContent(gen, meas);
+ // TODO fix error
error += inputMC->GetBinError(mcGenBin) * corr->GetBinContent(gen, meas);
}
virtual Long64_t Merge(const TCollection* list);
void FillMeasured(Float_t vtx, Int_t measured05, Int_t measured10, Int_t measured14);
+ void FillTriggeredEvent(Int_t measured05, Int_t measured10, Int_t measured14);
void FillGenerated(Float_t vtx, Bool_t triggered, Bool_t vertex, AliPWG0Helper::MCProcessType processType, Int_t generated05, Int_t generated10, Int_t generated14, Int_t generatedAll);
void FillCorrection(Float_t vtx, Int_t generated05, Int_t generated10, Int_t generated14, Int_t generatedAll, Int_t measured05, Int_t measured10, Int_t measured14);
+ void FillNoVertexEvent(Float_t vtx, Bool_t vertexReconstructed, Int_t generated05, Int_t generated10, Int_t generated14, Int_t measured05, Int_t measured10, Int_t measured14);
Bool_t LoadHistograms(const Char_t* dir = 0);
void SaveHistograms(const char* dir = 0);
void DrawHistograms();
+ void Rebin2DY(TH2F*& hist, Int_t nBins, Double_t* newBins) const;
+ void Rebin3DY(TH3F*& hist, Int_t nBins, Double_t* newBins) const;
+ void RebinGenerated(Int_t nBins05, Double_t* newBins05, Int_t nBins10, Double_t* newBins10, Int_t nBins13, Double_t* newBins13);
void DrawComparison(const char* name, Int_t inputRange, Bool_t fullPhaseSpace, Bool_t normalizeESD, TH1* mcHist, Bool_t simple = kFALSE, EventType eventType = kTrVtx);
- Int_t ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check = kFALSE, TH1* initialConditions = 0);
+ Int_t ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Int_t zeroBinEvents, Bool_t check = kFALSE, TH1* initialConditions = 0, Bool_t errorAsBias = kFALSE);
- void ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar = 1, Int_t nIterations = 100, TH1* initialConditions = 0, Bool_t determineError = kTRUE);
+ void ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar = 1, Int_t nIterations = 100, TH1* initialConditions = 0, Int_t determineError = 1);
static TH1* CalculateStdDev(TH1** results, Int_t max);
- TH1* StatisticalUncertainty(AliUnfolding::MethodType methodType, Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t randomizeMeasured, Bool_t randomizeResponse, const TH1* compareTo = 0);
+ TH1* StatisticalUncertainty(AliUnfolding::MethodType methodType, Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Int_t zeroBinEvents, Bool_t randomizeMeasured, Bool_t randomizeResponse, const TH1* compareTo = 0);
Int_t ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType);
void ApplyGaussianMethod(Int_t inputRange, Bool_t fullPhaseSpace);
void ApplyLaszloMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType);
+ void Calculate0Bin(Int_t inputRange, EventType eventType, Int_t zeroBinEvents);
+ Float_t GetFraction0Generated(Int_t inputRange);
+
TH2F* GetMultiplicityESD(Int_t i) const { return fMultiplicityESD[i]; }
+ TH1F* GetTriggeredEvents(Int_t i) const { return fTriggeredEvents[i]; }
TH2F* GetMultiplicityVtx(Int_t i) const { return fMultiplicityVtx[i]; }
TH2F* GetMultiplicityMB(Int_t i) const { return fMultiplicityMB[i]; }
TH2F* GetMultiplicityINEL(Int_t i) const { return fMultiplicityINEL[i]; }
TH2F* GetMultiplicityNSD(Int_t i) const { return fMultiplicityNSD[i]; }
TH2F* GetMultiplicityMC(Int_t i, EventType eventType) const;
TH3F* GetCorrelation(Int_t i) const { return fCorrelation[i]; }
+ TH1F* GetNoVertexEvents(Int_t i) const { return fNoVertexEvents[i]; }
TH1F* GetMultiplicityESDCorrected(Int_t i) const { return fMultiplicityESDCorrected[i]; }
void SetMultiplicityESD(Int_t i, TH2F* const hist) { fMultiplicityESD[i] = hist; }
+ void SetTriggeredEvents(Int_t i, TH1F* hist) { fTriggeredEvents[i] = hist; }
void SetMultiplicityVtx(Int_t i, TH2F* const hist) { fMultiplicityVtx[i] = hist; }
void SetMultiplicityMB(Int_t i, TH2F* const hist) { fMultiplicityMB[i] = hist; }
void SetMultiplicityINEL(Int_t i, TH2F* const hist) { fMultiplicityINEL[i] = hist; }
void SetMultiplicityNSD(Int_t i, TH2F* const hist) { fMultiplicityNSD[i] = hist; }
void SetMultiplicityMC(Int_t i, EventType eventType, TH2F* const hist);
void SetCorrelation(Int_t i, TH3F* const hist) { fCorrelation[i] = hist; }
+ void SetNoVertexEvents(Int_t i, TH1F* hist) { fNoVertexEvents[i] = hist; }
void SetMultiplicityESDCorrected(Int_t i, TH1F* const hist) { fMultiplicityESDCorrected[i] = hist; }
void SetGenMeasFromFunc(const TF1* inputMC, Int_t id);
TH2F* CalculateMultiplicityESD(TH1* inputMC, Int_t correlationMap);
+
+ void FixTriggerEfficiencies(Int_t start);
+
+ TH1* GetConvoluted(Int_t i, EventType eventType);
+ TH1* GetResiduals(Int_t i, EventType eventType, Float_t& residualSum);
void GetComparisonResults(Float_t* const mc = 0, Int_t* const mcLimit = 0, Float_t* const residuals = 0, Float_t* const ratioAverage = 0) const;
TH1* GetEfficiency(Int_t inputRange, EventType eventType);
- TH1* GetTriggerEfficiency(Int_t inputRange);
+ TH1* GetTriggerEfficiency(Int_t inputRange, EventType eventType);
static void SetQualityRegions(Bool_t SPDStudy);
Float_t GetQuality(Int_t region) const { return fQuality[region]; }
void FFT(Int_t dir, Int_t m, Double_t *x, Double_t *y) const;
+
+ Float_t GetVertexBegin(Int_t i) { return fgVtxRangeBegin[i]; }
+ Float_t GetVertexEnd(Int_t i) { return fgVtxRangeEnd[i]; }
+
+ void SetVertexRange(Int_t begin, Int_t end) { fVtxBegin = begin; fVtxEnd = end; }
protected:
void SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType);
TH1* fCurrentEfficiency; //! current efficiency
TH2F* fMultiplicityESD[kESDHists]; // multiplicity histogram: vtx vs multiplicity; array: |eta| < 0.5, 1.0, 1.4 (0..2)
+ TH1F* fTriggeredEvents[kESDHists]; // (raw) multiplicity distribution of triggered events; array: |eta| < 0.5, 1.0, 1.4 (0..2)
+ TH1F* fNoVertexEvents[kESDHists]; // distribution of true multiplicity just of triggered events without vertex or with 0 tracklets; array: |eta| < 0.5, 1.0, 1.4 (0..2)
TH2F* fMultiplicityVtx[kMCHists]; // multiplicity histogram of events that have a reconstructed vertex : vtx vs multiplicity; array: |eta| < 0.5, 1.0, 1.4, inf (0..3)
TH2F* fMultiplicityMB[kMCHists]; // multiplicity histogram of triggered events : vtx vs multiplicity; array: |eta| < 0.5, 1.0, 1.4, inf (0..3)
TH3F* fCorrelation[kCorrHists]; // vtx vs. (gene multiplicity (trig+vtx)) vs. (meas multiplicity); array: |eta| < 0.5, 1, 1.4, (0..2 and 3..5), the first corrects to the eta range itself, the second to full phase space
TH1F* fMultiplicityESDCorrected[kCorrHists]; // corrected histograms
-
+
Int_t fLastBinLimit; //! last bin limit, determined in SetupCurrentHists()
Float_t fLastChi2MC; //! last Chi2 between MC and unfolded ESD (calculated in DrawComparison)
Int_t fLastChi2MCLimit; //! bin where the last chi2 breached a certain threshold, used to evaluate the multiplicity reach (calc. in DrawComparison)
Float_t fLastChi2Residuals; //! last Chi2 of the ESD and the folded unfolded ESD (calculated in DrawComparison)
Float_t fRatioAverage; //! last average of |ratio-1| where ratio = unfolded / mc (bin 2..150)
+
+ Int_t fVtxBegin; //! vertex range for analysis
+ Int_t fVtxEnd; //! vertex range for analysis
+
+ static Double_t fgVtxRangeBegin[kESDHists]; //! begin of allowed vertex range for this eta bin
+ static Double_t fgVtxRangeEnd[kESDHists]; //! end of allowed vertex range for this eta bin
static Int_t fgQualityRegionsB[kQualityRegions]; //! begin, given in multiplicity units
static Int_t fgQualityRegionsE[kQualityRegions]; //! end
AliMultiplicityCorrection(const AliMultiplicityCorrection&);
AliMultiplicityCorrection& operator=(const AliMultiplicityCorrection&);
- ClassDef(AliMultiplicityCorrection, 5);
+ ClassDef(AliMultiplicityCorrection, 7);
};
#endif
#include "multiplicity/AliMultiplicityCorrection.h"
#include "AliCorrection.h"
#include "AliCorrectionMatrix3D.h"
+#include "AliPhysicsSelection.h"
#include "AliTriggerAnalysis.h"
ClassImp(AliMultiplicityTask)
fAnalysisMode((AliPWG0Helper::AnalysisMode) (AliPWG0Helper::kSPD | AliPWG0Helper::kFieldOn)),
fTrigger(AliTriggerAnalysis::kMB1),
fDeltaPhiCut(-1),
+ fDiffTreatment(AliPWG0Helper::kMCFlags),
fReadMC(kFALSE),
fUseMCVertex(kFALSE),
fMultiplicity(0),
fParticleSpecies(0),
fdNdpT(0),
fPtSpectrum(0),
+ fTemp1(0),
+ fTemp2(0),
fOutput(0)
{
//
for (Int_t i = 0; i<8; i++)
fParticleCorrection[i] = 0;
+
+ for (Int_t i=0; i<3; i++)
+ fEta[i] = 0;
// Define input and output slots here
DefineInput(0, TChain::Class());
DefineOutput(0, TList::Class());
+
+ if (fOption.Contains("only-process-type-nd"))
+ fSelectProcessType = 1;
+
+ if (fOption.Contains("only-process-type-sd"))
+ fSelectProcessType = 2;
+
+ if (fOption.Contains("only-process-type-dd"))
+ fSelectProcessType = 3;
+
+ if (fSelectProcessType != 0)
+ AliInfo(Form("WARNING: Systematic study enabled. Only considering process type %d", fSelectProcessType));
}
AliMultiplicityTask::~AliMultiplicityTask()
fdNdpT->Sumw2();
fOutput->Add(fdNdpT);
- if (fOption.Contains("skip-particles"))
- {
- fSystSkipParticles = kTRUE;
- AliInfo("WARNING: Systematic study enabled. Particles will be skipped.");
- }
-
if (fOption.Contains("particle-efficiency"))
for (Int_t i = 0; i<8; i++)
{
fOutput->Add(fParticleCorrection[i]);
}
- if (fOption.Contains("only-process-type-nd"))
- fSelectProcessType = 1;
-
- if (fOption.Contains("only-process-type-sd"))
- fSelectProcessType = 2;
-
- if (fOption.Contains("only-process-type-dd"))
- fSelectProcessType = 3;
-
- if (fSelectProcessType != 0)
- AliInfo(Form("WARNING: Systematic study enabled. Only considering process type %d", fSelectProcessType));
-
if (fOption.Contains("pt-spectrum-hist"))
{
TFile* file = TFile::Open("ptspectrum_fit.root");
fOutput->Add(fParticleSpecies);
}
+ fTemp1 = new TH2F("fTemp1", "fTemp1", 100, -0.5, 99.5, 100, -0.5, 99.5);
+ fOutput->Add(fTemp1);
+
+ for (Int_t i=0; i<3; i++)
+ {
+ fEta[i] = new TH1F(Form("fEta_%d", i), ";#eta", 100, -2, 2);
+ fOutput->Add(fEta[i]);
+ }
+
// TODO set seed for random generator
}
return;
}
- static AliTriggerAnalysis* triggerAnalysis = new AliTriggerAnalysis;
- Bool_t eventTriggered = triggerAnalysis->IsTriggerFired(fESD, fTrigger);
- //Printf("%lld", fESD->GetTriggerMask());
+ AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
+ if (!inputHandler)
+ {
+ Printf("ERROR: Could not receive input handler");
+ return;
+ }
+
+ Bool_t eventTriggered = inputHandler->IsEventSelected();
+ static AliTriggerAnalysis* triggerAnalysis = 0;
+ if (!triggerAnalysis)
+ {
+ AliPhysicsSelection* physicsSelection = dynamic_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
+ triggerAnalysis = physicsSelection->GetTriggerAnalysis();
+ }
+ if (eventTriggered)
+ eventTriggered = triggerAnalysis->IsTriggerFired(fESD, fTrigger);
+
const AliESDVertex* vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
+ if (vtxESD && !AliPWG0Helper::TestVertex(vtxESD, fAnalysisMode))
+ vtxESD = 0;
+
+ // remove vertices outside +- 15 cm
+ if (vtxESD && TMath::Abs(vtxESD->GetZv()) > 15)
+ vtxESD = 0;
+
Bool_t eventVertex = (vtxESD != 0);
Double_t vtx[3];
Float_t* etaArr = 0;
if (fAnalysisMode & AliPWG0Helper::kSPD)
{
- // get tracklets
- const AliMultiplicity* mult = fESD->GetMultiplicity();
- if (!mult)
- {
- AliDebug(AliLog::kError, "AliMultiplicity not available");
- return;
- }
-
- labelArr = new Int_t[mult->GetNumberOfTracklets()];
- etaArr = new Float_t[mult->GetNumberOfTracklets()];
-
- // get multiplicity from ITS tracklets
- for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
+ if (vtxESD)
{
- //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
-
- Float_t deltaPhi = mult->GetDeltaPhi(i);
-
- if (fDeltaPhiCut > 0 && TMath::Abs(deltaPhi) > fDeltaPhiCut)
- continue;
-
- etaArr[inputCount] = mult->GetEta(i);
- if (mult->GetLabel(i, 0) == mult->GetLabel(i, 1))
+ // get tracklets
+ const AliMultiplicity* mult = fESD->GetMultiplicity();
+ if (!mult)
{
- labelArr[inputCount] = mult->GetLabel(i, 0);
+ AliDebug(AliLog::kError, "AliMultiplicity not available");
+ return;
}
- else
- labelArr[inputCount] = -1;
+
+ labelArr = new Int_t[mult->GetNumberOfTracklets()];
+ etaArr = new Float_t[mult->GetNumberOfTracklets()];
+
+ Bool_t foundInEta10 = kFALSE;
+
+ // get multiplicity from ITS tracklets
+ for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
+ {
+ //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
+
+ Float_t deltaPhi = mult->GetDeltaPhi(i);
+
+ if (fDeltaPhiCut > 0 && TMath::Abs(deltaPhi) > fDeltaPhiCut)
+ continue;
+
+ if (fSystSkipParticles && gRandom->Uniform() < (0.0153))
+ {
+ Printf("Skipped tracklet!");
+ continue;
+ }
+
+ etaArr[inputCount] = mult->GetEta(i);
+ if (mult->GetLabel(i, 0) == mult->GetLabel(i, 1))
+ {
+ labelArr[inputCount] = mult->GetLabel(i, 0);
+ }
+ else
+ labelArr[inputCount] = -1;
+
+ for (Int_t i=0; i<3; i++)
+ {
+ if (vtx[2] > fMultiplicity->GetVertexBegin(i) && vtx[2] < fMultiplicity->GetVertexEnd(i))
+ fEta[i]->Fill(etaArr[inputCount]);
+ }
- ++inputCount;
+ // we have to repeat the trigger here, because the tracklet might have been kicked out fSystSkipParticles
+ if (TMath::Abs(etaArr[inputCount]) < 1)
+ foundInEta10 = kTRUE;
+
+ ++inputCount;
+ }
+
+ if (fSystSkipParticles && (fTrigger & AliTriggerAnalysis::kOneParticle) && !foundInEta10)
+ eventTriggered = kFALSE;
}
}
else if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i));
continue;
}
+
+ if (esdTrack->Pt() < 0.15)
+ continue;
+
+ Float_t d0z0[2],covd0z0[3];
+ esdTrack->GetImpactParameters(d0z0,covd0z0);
+ Float_t sigma= 0.0050+0.0060/TMath::Power(esdTrack->Pt(),0.9);
+ Float_t d0max = 7.*sigma;
+ if (TMath::Abs(d0z0[0]) > d0max)
+ continue;
etaArr[inputCount] = esdTrack->Eta();
labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel());
if (!fReadMC) // Processing of ESD information
{
- if (eventTriggered && eventVertex)
+ Int_t nESDTracks05 = 0;
+ Int_t nESDTracks10 = 0;
+ Int_t nESDTracks14 = 0;
+
+ for (Int_t i=0; i<inputCount; ++i)
{
- Int_t nESDTracks05 = 0;
- Int_t nESDTracks10 = 0;
- Int_t nESDTracks14 = 0;
+ Float_t eta = etaArr[i];
- for (Int_t i=0; i<inputCount; ++i)
- {
- Float_t eta = etaArr[i];
+ if (TMath::Abs(eta) < 0.5)
+ nESDTracks05++;
- if (TMath::Abs(eta) < 0.5)
- nESDTracks05++;
+ if (TMath::Abs(eta) < 1.0)
+ nESDTracks10++;
- if (TMath::Abs(eta) < 1.0)
- nESDTracks10++;
+ if (TMath::Abs(eta) < 1.3)
+ nESDTracks14++;
+ }
+
+ // kick out randomly for combinatorics
+ /*
+ if (gRandom->Uniform() < 1.3e-4 * nESDTracks05 * nESDTracks05)
+ nESDTracks05--;
- if (TMath::Abs(eta) < 1.4)
- nESDTracks14++;
- }
+ if (gRandom->Uniform() < 8.7e-5 * nESDTracks10 * nESDTracks10)
+ nESDTracks10--;
+ if (gRandom->Uniform() < 9.6e-5 * nESDTracks14 * nESDTracks14)
+ nESDTracks14--;
+ */
+
+ //if (nESDTracks05 >= 20 || nESDTracks10 >= 30 || nESDTracks14 >= 32)
+ // Printf("File: %s, IEV: %d, TRG: ---, Orbit: 0x%x, Period: %d, BC: %d; Tracks: %d %d %d", ((TTree*) GetInputData(0))->GetCurrentFile()->GetName(), fESD->GetEventNumberInFile(), fESD->GetOrbitNumber(),fESD->GetPeriodNumber(),fESD->GetBunchCrossNumber(), nESDTracks05, nESDTracks10, nESDTracks14);
+
+ if (eventTriggered)
+ fMultiplicity->FillTriggeredEvent(nESDTracks05, nESDTracks10, nESDTracks14);
+
+ if (eventTriggered && eventVertex)
fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks10, nESDTracks14);
- }
}
else if (fReadMC) // Processing of MC information
{
}
// get process information
- AliPWG0Helper::MCProcessType processType = AliPWG0Helper::GetEventProcessType(header);
+ AliPWG0Helper::MCProcessType processType = AliPWG0Helper::GetEventProcessType(fESD, header, stack, fDiffTreatment);
Bool_t processEvent = kTRUE;
if (fSelectProcessType > 0)
if (TMath::Abs(particle->Eta()) < 1.0)
nMCTracks10 += particleWeight;
- if (TMath::Abs(particle->Eta()) < 1.4)
+ if (TMath::Abs(particle->Eta()) < 1.3)
nMCTracks14 += particleWeight;
nMCTracksAll += particleWeight;
fMultiplicity->FillGenerated(vtxMC[2], eventTriggered, eventVertex, processType, (Int_t) nMCTracks05, (Int_t) nMCTracks10, (Int_t) nMCTracks14, (Int_t) nMCTracksAll);
- if (eventTriggered && eventVertex)
- {
- Int_t nESDTracks05 = 0;
- Int_t nESDTracks10 = 0;
- Int_t nESDTracks14 = 0;
+ // ESD processing
+ Int_t nESDTracks05 = 0;
+ Int_t nESDTracks10 = 0;
+ Int_t nESDTracks14 = 0;
- // tracks per particle species, in |eta| < 2 (systematic study)
- Int_t nESDTracksSpecies[7]; // (pi, K, p, other, nolabel, doublecount_prim, doublecount_all)
- for (Int_t i = 0; i<7; ++i)
- nESDTracksSpecies[i] = 0;
+ // tracks per particle species, in |eta| < 2 (systematic study)
+ Int_t nESDTracksSpecies[7]; // (pi, K, p, other, nolabel, doublecount_prim, doublecount_all)
+ for (Int_t i = 0; i<7; ++i)
+ nESDTracksSpecies[i] = 0;
- Bool_t* foundPrimaries = new Bool_t[nPrim]; // to prevent double counting
- for (Int_t i=0; i<nPrim; i++)
- foundPrimaries[i] = kFALSE;
+ Bool_t* foundPrimaries = new Bool_t[nPrim]; // to prevent double counting
+ for (Int_t i=0; i<nPrim; i++)
+ foundPrimaries[i] = kFALSE;
- Bool_t* foundPrimaries2 = new Bool_t[nPrim]; // to prevent double counting
- for (Int_t i=0; i<nPrim; i++)
- foundPrimaries2[i] = kFALSE;
+ Bool_t* foundPrimaries2 = new Bool_t[nPrim]; // to prevent double counting
+ for (Int_t i=0; i<nPrim; i++)
+ foundPrimaries2[i] = kFALSE;
- Bool_t* foundTracks = new Bool_t[nMCPart]; // to prevent double counting
- for (Int_t i=0; i<nMCPart; i++)
- foundTracks[i] = kFALSE;
+ Bool_t* foundTracks = new Bool_t[nMCPart]; // to prevent double counting
+ for (Int_t i=0; i<nMCPart; i++)
+ foundTracks[i] = kFALSE;
- for (Int_t i=0; i<inputCount; ++i)
- {
- Float_t eta = etaArr[i];
- Int_t label = labelArr[i];
+ for (Int_t i=0; i<inputCount; ++i)
+ {
+ Float_t eta = etaArr[i];
+ Int_t label = labelArr[i];
- Int_t particleWeight = 1;
+ Int_t particleWeight = 1;
- // systematic study: 5% lower efficiency
- if (fSystSkipParticles && (gRandom->Uniform() < 0.05))
- continue;
-
- // in case of systematic study, weight according to the change of the pt spectrum
- if (fPtSpectrum)
- {
- TParticle* mother = 0;
+ // in case of systematic study, weight according to the change of the pt spectrum
+ if (fPtSpectrum)
+ {
+ TParticle* mother = 0;
- // preserve label for later
- Int_t labelCopy = label;
- if (labelCopy >= 0)
- labelCopy = AliPWG0Helper::FindPrimaryMotherLabel(stack, labelCopy);
- if (labelCopy >= 0)
- mother = stack->Particle(labelCopy);
+ // preserve label for later
+ Int_t labelCopy = label;
+ if (labelCopy >= 0)
+ labelCopy = AliPWG0Helper::FindPrimaryMotherLabel(stack, labelCopy);
+ if (labelCopy >= 0)
+ mother = stack->Particle(labelCopy);
- // in case of pt study we do not count particles w/o label, because they cannot be scaled
- if (!mother)
- continue;
+ // in case of pt study we do not count particles w/o label, because they cannot be scaled
+ if (!mother)
+ continue;
- // it cannot be just multiplied because we cannot count "half of a particle"
- // instead a random generator decides if the particle is counted twice (if value > 1)
- // or not (if value < 0)
- Int_t bin = fPtSpectrum->FindBin(mother->Pt());
- if (bin > 0 && bin <= fPtSpectrum->GetNbinsX())
+ // it cannot be just multiplied because we cannot count "half of a particle"
+ // instead a random generator decides if the particle is counted twice (if value > 1)
+ // or not (if value < 0)
+ Int_t bin = fPtSpectrum->FindBin(mother->Pt());
+ if (bin > 0 && bin <= fPtSpectrum->GetNbinsX())
+ {
+ Float_t factor = fPtSpectrum->GetBinContent(bin);
+ if (factor > 0)
{
- Float_t factor = fPtSpectrum->GetBinContent(bin);
- if (factor > 0)
+ Float_t random = gRandom->Uniform();
+ if (factor > 1 && random < factor - 1)
{
- Float_t random = gRandom->Uniform();
- if (factor > 1 && random < factor - 1)
- {
- particleWeight = 2;
- }
- else if (factor < 1 && random < 1 - factor)
- particleWeight = 0;
+ particleWeight = 2;
}
+ else if (factor < 1 && random < 1 - factor)
+ particleWeight = 0;
}
}
+ }
- //Printf("ESD weight is: %d", particleWeight);
+ //Printf("ESD weight is: %d", particleWeight);
- if (TMath::Abs(eta) < 0.5)
- nESDTracks05 += particleWeight;
+ if (TMath::Abs(eta) < 0.5)
+ nESDTracks05 += particleWeight;
- if (TMath::Abs(eta) < 1.0)
- nESDTracks10 += particleWeight;
+ if (TMath::Abs(eta) < 1.0)
+ nESDTracks10 += particleWeight;
- if (TMath::Abs(eta) < 1.4)
- nESDTracks14 += particleWeight;
+ if (TMath::Abs(eta) < 1.3)
+ nESDTracks14 += particleWeight;
- if (fParticleSpecies)
+ if (fParticleSpecies)
+ {
+ Int_t motherLabel = -1;
+ TParticle* mother = 0;
+
+ // find mother
+ if (label >= 0)
+ motherLabel = AliPWG0Helper::FindPrimaryMotherLabel(stack, label);
+ if (motherLabel >= 0)
+ mother = stack->Particle(motherLabel);
+
+ if (!mother)
+ {
+ // count tracks that did not have a label
+ if (TMath::Abs(eta) < etaRange)
+ nESDTracksSpecies[4]++;
+ }
+ else
{
- Int_t motherLabel = -1;
- TParticle* mother = 0;
+ // get particle type (pion, proton, kaon, other) of mother
+ Int_t idMother = -1;
+ switch (TMath::Abs(mother->GetPdgCode()))
+ {
+ case 211: idMother = 0; break;
+ case 321: idMother = 1; break;
+ case 2212: idMother = 2; break;
+ default: idMother = 3; break;
+ }
+
+ // double counting is ok for particle ratio study
+ if (TMath::Abs(eta) < etaRange)
+ nESDTracksSpecies[idMother]++;
- // find mother
- if (label >= 0)
- motherLabel = AliPWG0Helper::FindPrimaryMotherLabel(stack, label);
- if (motherLabel >= 0)
- mother = stack->Particle(motherLabel);
+ // double counting is not ok for efficiency study
- if (!mother)
+ // check if we already counted this particle, this way distinguishes double counted particles (bug/artefact in tracking) or double counted primaries due to secondaries (physics)
+ if (foundTracks[label])
{
- // count tracks that did not have a label
if (TMath::Abs(eta) < etaRange)
- nESDTracksSpecies[4]++;
+ nESDTracksSpecies[6]++;
}
else
{
- // get particle type (pion, proton, kaon, other) of mother
- Int_t idMother = -1;
- switch (TMath::Abs(mother->GetPdgCode()))
- {
- case 211: idMother = 0; break;
- case 321: idMother = 1; break;
- case 2212: idMother = 2; break;
- default: idMother = 3; break;
- }
-
- // double counting is ok for particle ratio study
- if (TMath::Abs(eta) < etaRange)
- nESDTracksSpecies[idMother]++;
+ foundTracks[label] = kTRUE;
- // double counting is not ok for efficiency study
-
- // check if we already counted this particle, this way distinguishes double counted particles (bug/artefact in tracking) or double counted primaries due to secondaries (physics)
- if (foundTracks[label])
+ // particle (primary) already counted?
+ if (foundPrimaries[motherLabel])
{
if (TMath::Abs(eta) < etaRange)
- nESDTracksSpecies[6]++;
+ nESDTracksSpecies[5]++;
}
else
- {
- foundTracks[label] = kTRUE;
-
- // particle (primary) already counted?
- if (foundPrimaries[motherLabel])
- {
- if (TMath::Abs(eta) < etaRange)
- nESDTracksSpecies[5]++;
- }
- else
- foundPrimaries[motherLabel] = kTRUE;
- }
+ foundPrimaries[motherLabel] = kTRUE;
}
}
+ }
- if (fParticleCorrection[0])
+ if (fParticleCorrection[0])
+ {
+ if (label >= 0 && stack->IsPhysicalPrimary(label))
{
- if (label >= 0 && stack->IsPhysicalPrimary(label))
- {
- TParticle* particle = stack->Particle(label);
+ TParticle* particle = stack->Particle(label);
- // get particle type (pion, proton, kaon, other)
- 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;
- }
+ // get particle type (pion, proton, kaon, other)
+ 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;
+ }
- // todo check if values are not completely off??
+ // todo check if values are not completely off??
- // particle (primary) already counted?
- if (!foundPrimaries2[label])
- {
- foundPrimaries2[label] = kTRUE;
- fParticleCorrection[id]->GetTrackCorrection()->FillMeas(vtxMC[2], particle->Eta(), particle->Pt());
- }
+ // particle (primary) already counted?
+ if (!foundPrimaries2[label])
+ {
+ foundPrimaries2[label] = kTRUE;
+ fParticleCorrection[id]->GetTrackCorrection()->FillMeas(vtxMC[2], particle->Eta(), particle->Pt());
}
}
}
-
- if (fParticleCorrection[0])
+ }
+
+ if (fParticleCorrection[0])
+ {
+ // if the particle decays/stops before this radius we do not see it
+ // 8cm larger than SPD layer 2
+ // 123cm TPC radius where a track has about 50 clusters (cut limit)
+ const Float_t endRadius = (fAnalysisMode & AliPWG0Helper::kSPD) ? 8. : 123;
+
+ // loop over all primaries that have not been found
+ for (Int_t i=0; i<nPrim; i++)
{
- // if the particle decays/stops before this radius we do not see it
- // 8cm larger than SPD layer 2
- // 123cm TPC radius where a track has about 50 clusters (cut limit)
- const Float_t endRadius = (fAnalysisMode & AliPWG0Helper::kSPD) ? 8. : 123;
-
- // loop over all primaries that have not been found
- for (Int_t i=0; i<nPrim; i++)
- {
- // already found
- if (foundPrimaries2[i])
- continue;
-
- TParticle* particle = 0;
- TClonesArray* trackrefs = 0;
- mcEvent->GetParticleAndTR(i, particle, trackrefs);
-
- // true primary and charged
- if (!AliPWG0Helper::IsPrimaryCharged(particle, nPrim))
- continue;
-
- //skip particles with larger |eta| than 3, to keep the log clean, is anyway not included in correction map
- if (TMath::Abs(particle->Eta()) > 3)
- continue;
-
- // skipping checking of process type of daughter: Neither kPBrem, kPDeltaRay nor kPCerenkov should appear in the event generation
+ // already found
+ if (foundPrimaries2[i])
+ continue;
- // get particle type (pion, proton, kaon, other)
- Int_t id = -1;
- switch (TMath::Abs(particle->GetPdgCode()))
- {
- case 211: id = 4; break;
- case 321: id = 5; break;
- case 2212: id = 6; break;
- default: id = 7; break;
- }
+ TParticle* particle = 0;
+ TClonesArray* trackrefs = 0;
+ mcEvent->GetParticleAndTR(i, particle, trackrefs);
+
+ // true primary and charged
+ if (!AliPWG0Helper::IsPrimaryCharged(particle, nPrim))
+ continue;
+
+ //skip particles with larger |eta| than 3, to keep the log clean, is anyway not included in correction map
+ if (TMath::Abs(particle->Eta()) > 3)
+ continue;
+
+ // skipping checking of process type of daughter: Neither kPBrem, kPDeltaRay nor kPCerenkov should appear in the event generation
+
+ // get particle type (pion, proton, kaon, other)
+ Int_t id = -1;
+ switch (TMath::Abs(particle->GetPdgCode()))
+ {
+ case 211: id = 4; break;
+ case 321: id = 5; break;
+ case 2212: id = 6; break;
+ default: id = 7; break;
+ }
+
+ if (!fParticleCorrection[id])
+ continue;
- if (!fParticleCorrection[id])
- continue;
-
- // get last track reference
- AliTrackReference* trackref = dynamic_cast<AliTrackReference*> (trackrefs->Last());
+ // get last track reference
+ AliTrackReference* trackref = dynamic_cast<AliTrackReference*> (trackrefs->Last());
+
+ if (!trackref)
+ {
+ Printf("ERROR: Could not get trackref of %d (count %d)", i, trackrefs->GetEntries());
+ particle->Print();
+ continue;
+ }
- if (!trackref)
+ // particle in tracking volume long enough...
+ if (trackref->R() > endRadius)
+ continue;
+
+ if (particle->GetLastDaughter() >= 0)
+ {
+ Int_t uID = stack->Particle(particle->GetLastDaughter())->GetUniqueID();
+ //if (uID != kPBrem && uID != kPDeltaRay && uID < kPCerenkov)
+ if (uID == kPDecay)
{
- Printf("ERROR: Could not get trackref of %d (count %d)", i, trackrefs->GetEntries());
- particle->Print();
- continue;
- }
+ // decayed
- // particle in tracking volume long enough...
- if (trackref->R() > endRadius)
- continue;
-
- if (particle->GetLastDaughter() >= 0)
- {
- Int_t uID = stack->Particle(particle->GetLastDaughter())->GetUniqueID();
- //if (uID != kPBrem && uID != kPDeltaRay && uID < kPCerenkov)
- if (uID == kPDecay)
- {
- // decayed
-
- Printf("Particle %d (%s) decayed at %f, daugher uniqueID: %d:", i, particle->GetName(), trackref->R(), uID);
- particle->Print();
- Printf("Daughers:");
- for (Int_t d = particle->GetFirstDaughter(); d <= particle->GetLastDaughter(); d++)
- stack->Particle(d)->Print();
- Printf("");
-
- fParticleCorrection[id]->GetTrackCorrection()->FillGene(vtxMC[2], particle->Eta(), particle->Pt());
- continue;
- }
- }
-
- if (trackref->DetectorId() == -1)
- {
- // stopped
- Printf("Particle %d stopped at %f:", i, trackref->R());
+ Printf("Particle %d (%s) decayed at %f, daugher uniqueID: %d:", i, particle->GetName(), trackref->R(), uID);
particle->Print();
+ Printf("Daughers:");
+ for (Int_t d = particle->GetFirstDaughter(); d <= particle->GetLastDaughter(); d++)
+ stack->Particle(d)->Print();
Printf("");
- fParticleCorrection[id]->GetTrackCorrection()->FillMeas(vtxMC[2], particle->Eta(), particle->Pt());
+ fParticleCorrection[id]->GetTrackCorrection()->FillGene(vtxMC[2], particle->Eta(), particle->Pt());
continue;
}
-
- Printf("Particle %d simply not tracked", i);
+ }
+
+ if (trackref->DetectorId() == -1)
+ {
+ // stopped
+ Printf("Particle %d stopped at %f:", i, trackref->R());
particle->Print();
Printf("");
+
+ fParticleCorrection[id]->GetTrackCorrection()->FillMeas(vtxMC[2], particle->Eta(), particle->Pt());
+ continue;
}
+
+ Printf("Particle %d simply not tracked", i);
+ particle->Print();
+ Printf("");
}
-
- delete[] foundTracks;
- delete[] foundPrimaries;
- delete[] foundPrimaries2;
+ }
+
+ delete[] foundTracks;
+ delete[] foundPrimaries;
+ delete[] foundPrimaries2;
- if ((Int_t) nMCTracks14 > 10 && nESDTracks14 <= 3)
- {
- TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
- printf("WARNING: Event %lld %s (vtx-z = %f, recon: %f, contrib: %d, res: %f) has %d generated and %d reconstructed...\n", tree->GetReadEntry(), tree->GetCurrentFile()->GetName(), vtxMC[2], vtx[2], vtxESD->GetNContributors(), vtxESD->GetZRes(), nMCTracks14, nESDTracks14);
- }
+// if ((Int_t) nMCTracks14 > 10 && nESDTracks14 <= 3)
+// {
+// TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
+// printf("WARNING: Event %lld %s (vtx-z = %f, recon: %f, contrib: %d, res: %f) has %d generated and %d reconstructed...\n", tree->GetReadEntry(), tree->GetCurrentFile()->GetName(), vtxMC[2], vtx[2], vtxESD->GetNContributors(), vtxESD->GetZRes(), nMCTracks14, nESDTracks14);
+// }
+ if (eventTriggered)
+ {
+ fMultiplicity->FillTriggeredEvent(nESDTracks05, nESDTracks10, nESDTracks14);
+ fMultiplicity->FillNoVertexEvent(vtxMC[2], eventVertex, nMCTracks05, nMCTracks10, nMCTracks14, nESDTracks05, nESDTracks10, nESDTracks14);
+// if (!eventVertex)
+// {
+// if (nESDTracks05 == 0)
+// fTemp1->Fill(nMCTracks05, nESDTracks05);
+// }
+ }
+
+ if (eventTriggered && eventVertex)
+ {
// fill response matrix using vtxMC (best guess)
- fMultiplicity->FillCorrection(vtxMC[2], nMCTracks05, nMCTracks10, nMCTracks14, nMCTracksAll, nESDTracks05, nESDTracks10, nESDTracks14);
+ fMultiplicity->FillCorrection(vtxMC[2], nMCTracks05, nMCTracks10, nMCTracks14, nMCTracksAll, nESDTracks05, nESDTracks10, nESDTracks14);
fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks10, nESDTracks14);
return;
}
- TFile* file = TFile::Open("multiplicity.root", "RECREATE");
+ TString fileName("multiplicity");
+ if (fSelectProcessType == 1)
+ fileName += "ND";
+ if (fSelectProcessType == 2)
+ fileName += "SD";
+ if (fSelectProcessType == 3)
+ fileName += "DD";
+ fileName += ".root";
+ TFile* file = TFile::Open(fileName, "RECREATE");
fMultiplicity->SaveHistograms();
for (Int_t i = 0; i < 8; ++i)
if (fdNdpT)
fdNdpT->Write();
+ fTemp1 = dynamic_cast<TH1*> (fOutput->FindObject("fTemp1"));
+ if (fTemp1)
+ fTemp1->Write();
+
+ for (Int_t i=0; i<3; i++)
+ {
+ fEta[i] = dynamic_cast<TH1*> (fOutput->FindObject(Form("fEta_%d", i)));
+ if (fEta[i])
+ {
+ fEta[i]->Sumw2();
+ Float_t events = fMultiplicity->GetMultiplicityESD(i)->Integral(1, fMultiplicity->GetMultiplicityESD(i)->GetNbinsX());
+ if (events > 0)
+ fEta[i]->Scale(1.0 / events);
+ fEta[i]->Scale(1.0 / fEta[i]->GetXaxis()->GetBinWidth(1));
+ fEta[i]->Write();
+ }
+ }
+
TObjString option(fOption);
option.Write();
file->Close();
- Printf("Written result to multiplicity.root");
+ Printf("Written result to %s", fileName.Data());
}
void SetReadMC(Bool_t flag = kTRUE) { fReadMC = flag; }
void SetUseMCVertex(Bool_t flag = kTRUE) { fUseMCVertex = flag; }
-
+ void SetSkipParticles(Bool_t flag = kTRUE) { fSystSkipParticles = flag; }
+ void SetDiffTreatment(AliPWG0Helper::DiffTreatment diffTreatment) { fDiffTreatment = diffTreatment; }
+
protected:
AliESDEvent *fESD; //! ESD object
AliPWG0Helper::AnalysisMode fAnalysisMode; // detector that is used for analysis
AliTriggerAnalysis::Trigger fTrigger; // trigger that is used
Float_t fDeltaPhiCut; // cut in delta phi (only SPD)
-
+ AliPWG0Helper::DiffTreatment fDiffTreatment; // how to identify SD events (see AliPWG0Helper::GetEventProcessType)
+
Bool_t fReadMC; // if true reads MC data (to build correlation maps)
Bool_t fUseMCVertex; // the MC vtx is used instead of the ESD vertex (for syst. check)
AliMultiplicityCorrection* fMultiplicity; //! object containing the extracted data
AliESDtrackCuts* fEsdTrackCuts; // Object containing the parameters of the esd track cuts
- Bool_t fSystSkipParticles; //! if true skips particles (systematic study)
+ Bool_t fSystSkipParticles; // if true skips particles (systematic study)
AliCorrection* fParticleCorrection[8]; //! correction from measured to generated particles for different particles for trigger, vertex sample in |eta| < 2; switch on with "particle-efficiency"
// for each of the species (0..3): pi, k, p, other; for systematic study of pt cut off
// 4..7 counts for the same species the decayed particles (in generated) and stopped (in measured)
- Int_t fSelectProcessType; //! 0 = all (default), 1 = ND, 2 = SD, 3 = DD (for systematic study)
+ Int_t fSelectProcessType; // 0 = all (default), 1 = ND, 2 = SD, 3 = DD (for systematic study)
TNtuple *fParticleSpecies; //! per event: vtx_mc, (pi, k, p, rest (in |eta| < 1)) X (true, recon) + (nolabel,
// doubleTracks, doublePrimaries) [doubleTracks + doublePrimaries are already part of
// rec. particles!); enable with: particle-species
TH1* fdNdpT; //! true pT spectrum (MC)
TH1D* fPtSpectrum; // function that modifies the pt spectrum (syst. study)
+
+ TH1* fTemp1; //! temp histogram for quick study of variables
+ TH1* fTemp2; //! temp histogram for quick study of variables
+
+ TH1* fEta[3]; //! eta histogram of events in the acceptance region for each of the eta-bins (control histogram)
TList* fOutput; //! list send on output slot 0
// script to correct the multiplicity spectrum + helpers
//
+Bool_t is900GeV = 0;
+Bool_t is2360TeV = 0;
+
+// 900 GeV, MC
+//const Int_t kBinLimits[] = { 42, 57, 60 };
+//const Int_t kTrustLimits[] = { 26, 42, 54 };
+
+// 2.36 TeV
+//const Int_t kBinLimits[] = { 43, 67, 83 };
+//const Int_t kTrustLimits[] = { 33, 57, 73 };
+
+// 7 TeV
+const Int_t kBinLimits[] = { 60, 120, 60 };
+const Int_t kTrustLimits[] = { 26, 70, 54 };
+
void SetTPC()
{
gSystem->Load("libPWG0base");
gSystem->Load("libPWG0base");
}
-void correct(const char* fileNameMC = "multiplicityMC.root", const char* folder = "Multiplicity", const char* fileNameESD = "multiplicityESD.root", Bool_t chi2 = kFALSE, Int_t histID = 1, Bool_t fullPhaseSpace = kFALSE, Float_t beta = 1e5, Int_t eventType = 0 /* AliMultiplicityCorrection::kTrVtx */, const char* targetFile = "unfolded.root", const char* folderESD = "Multiplicity")
+void LoadAndInitialize(void* multVoid, void* esdVoid, void* multTriggerVoid, Int_t histID, Bool_t fullPhaseSpace, Int_t* geneLimits)
{
- loadlibs();
-
- AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(fileNameMC, folder);
- AliMultiplicityCorrection* esd = AliMultiplicityCorrection::Open(fileNameESD, folderESD);
+ AliMultiplicityCorrection* mult = (AliMultiplicityCorrection*) multVoid;
+ AliMultiplicityCorrection* esd = (AliMultiplicityCorrection*) esdVoid;
+ AliMultiplicityCorrection* multTrigger = (AliMultiplicityCorrection*) multTriggerVoid;
- //AliUnfolding::SetNbins(150, 150);
- //AliUnfolding::SetNbins(65, 65);
- //AliUnfolding::SetNbins(35, 35);
- //AliUnfolding::SetDebug(1);
-
- for (Int_t hID = ((histID == -1) ? 0 : histID); hID <= ((histID == -1) ? 2 : histID); hID++)
+ for (Int_t i=0; i<3; i++)
+ geneLimits[i] = kBinLimits[i];
+
+ // REBINNING
+ if (1)
{
- switch (hID)
+ // 900 GeV
+ if (is900GeV)
+ {
+ if (1)
+ {
+ Int_t bins05 = 31;
+ Double_t newBinsEta05[] = { -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
+ 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
+ 20.5, 22.5, 24.5, 26.5, 28.5, 30.5, 34.5, 38.5, 42.5, 50.5,
+ 100.5 };
+ }
+
+ if (0)
+ {
+ Int_t bins05 = 29;
+ Double_t newBinsEta05[] = { -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
+ 10.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
+ 20.5, 22.5, 24.5, 26.5, 28.5, 30.5, 34.5, 38.5, 42.5, 50.5,
+ 100.5 };
+ }
+
+ if (0)
+ {
+ Int_t bins05 = 25;
+ Double_t* newBinsEta05 = new Double_t[bins05+1];
+
+ //newBinsEta05[0] = -0.5;
+ //newBinsEta05[1] = 0.5;
+ //newBinsEta05[2] = 1.5;
+ //newBinsEta05[3] = 2.5;
+
+ for (Int_t i=0; i<bins05+1; i++)
+ newBinsEta05[i] = -0.5 + i*2;
+ //newBinsEta05[i] = 4.5 + (i-4)*2;
+ }
+
+ Int_t bins10 = 54;
+ Double_t newBinsEta10[] = { -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
+ 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
+ 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5,
+ 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5,
+ 40.5, 42.5, 44.5, 46.5, 48.5, 50.5, 52.5, 54.5, 56.5, 58.5,
+ 60.5, 65.5, 70.5, 100.5 };
+
+ geneLimits[0] = bins05;
+ geneLimits[1] = bins10;
+ geneLimits[2] = bins10;
+ }
+
+ // 2.36 TeV
+ if (is2360TeV)
{
- case 0: AliUnfolding::SetNbins(35, 35); break;
- case 1: AliUnfolding::SetNbins(60, 60); break;
- case 2: AliUnfolding::SetNbins(70, 70); beta *= 3; break;
+ Int_t bins05 = 36;
+ Double_t newBinsEta05[] = { -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
+ 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
+ 20.5, 22.5, 24.5, 26.5, 28.5, 30.5, 32.5, 34.5, 36.5, 38.5,
+ 40.5, 45.5, 50.5, 55.5, 60.5, 100.5 };
+ Int_t bins10 = 64;
+ Double_t newBinsEta10[] = { -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
+ 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
+ 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5,
+ 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5,
+ 40.5, 42.5, 44.5, 46.5, 48.5, 50.5, 52.5, 54.5, 56.5, 58.5,
+ 60.5, 62.5, 64.5, 66.5, 68.5, 70.5, 72.5, 74.5, 76.5, 78.5,
+ 80.5, 85.5, 90.5, 100.5 };
+
+ geneLimits[0] = bins05;
+ geneLimits[1] = bins10;
+ geneLimits[2] = bins10;
}
+
+ // 7 TeV
+ if (!is900GeV && !is2360TeV)
+ {
+ Int_t bins05 = 36;
+ Double_t newBinsEta05[] = { -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
+ 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
+ 20.5, 22.5, 24.5, 26.5, 28.5, 30.5, 32.5, 34.5, 36.5, 38.5,
+ 40.5, 45.5, 50.5, 55.5, 60.5, 100.5 };
+
+ Int_t bins10 = 85;
+ Double_t newBinsEta10[] = {
+ -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
+ 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
+ 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5,
+ 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5,
+ 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5,
+ 50.5, 51.5, 52.5, 53.5, 54.5, 55.5, 56.5, 57.5, 58.5, 59.5,
+ 60.5, 62.5, 64.5, 66.5, 68.5, 70.5, 72.5, 74.5, 76.5, 78.5,
+ 80.5, 82.5, 84.5, 86.5, 88.5, 90.5, 92.5, 94.5, 96.5, 98.5,
+ 100.5, 105.5, 110.5, 115.5, 120.5 };
+
+ geneLimits[0] = bins05;
+ geneLimits[1] = bins10;
+ geneLimits[2] = bins10;
+ }
+
+ esd->RebinGenerated(bins05, newBinsEta05, bins10, newBinsEta10, bins10, newBinsEta10);
+ mult->RebinGenerated(bins05, newBinsEta05, bins10, newBinsEta10, bins10, newBinsEta10);
+ multTrigger->RebinGenerated(bins05, newBinsEta05, bins10, newBinsEta10, bins10, newBinsEta10);
+ }
+ for (Int_t hID = ((histID == -1) ? 0 : histID); hID <= ((histID == -1) ? 2 : histID); hID++)
+ {
TH2F* hist = esd->GetMultiplicityESD(hID);
- TH2F* hist2 = esd->GetMultiplicityMC(hID, eventType);
mult->SetMultiplicityESD(hID, hist);
-
+ mult->SetTriggeredEvents(hID, esd->GetTriggeredEvents(hID));
+
+ // insert trigger efficiency in flat response matrix
+ for (Int_t evType = AliMultiplicityCorrection::kTrVtx; evType <= AliMultiplicityCorrection::kNSD; evType++)
+ mult->SetMultiplicityMC(hID, evType, multTrigger->GetMultiplicityMC(hID, evType));
+
+ mult->SetNoVertexEvents(hID, multTrigger->GetNoVertexEvents(hID));
+ mult->FixTriggerEfficiencies(10);
+
// small hack to get around charge conservation for full phase space ;-)
if (fullPhaseSpace)
{
corr->SetBinError(i, j, corr->GetBinError(i-1, j));
}
}
+ }
+}
+
+void correct(const char* fileNameMC = "multiplicityMC.root", const char* folder = "Multiplicity", const char* fileNameESD = "multiplicityESD.root", Bool_t chi2 = 1, Int_t histID = 1, Bool_t fullPhaseSpace = kFALSE, Float_t beta = -25, Int_t eventType = 2 /* AliMultiplicityCorrection::kTrVtx */, const char* targetFile = "unfolded.root", const char* folderESD = "Multiplicity", Bool_t calcBias = kFALSE)
+{
+ loadlibs();
+
+ Int_t geneLimits[] = { 0, 0, 0 };
+ AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(fileNameMC, folder);
+ AliMultiplicityCorrection* esd = AliMultiplicityCorrection::Open(fileNameESD, folderESD);
+ AliMultiplicityCorrection* multTrigger = AliMultiplicityCorrection::Open("multiplicityTrigger.root");
+
+ LoadAndInitialize(mult, esd, multTrigger, histID, fullPhaseSpace, geneLimits);
+
+ //AliUnfolding::SetDebug(1);
+
+ for (Int_t hID = ((histID == -1) ? 0 : histID); hID <= ((histID == -1) ? 2 : histID); hID++)
+ {
+ AliUnfolding::SetNbins(kBinLimits[hID], geneLimits[hID]);
+
+ TH2F* hist2 = esd->GetMultiplicityMC(hID, eventType);
+
if (chi2)
{
- AliUnfolding::SetChi2Regularization(AliUnfolding::kPol0, beta);
+ if (is900GeV)
+ {
+ //Float_t regParamPol0[] = { 2, 2, 10 }; // TPCITS
+ Float_t regParamPol0[] = { 5, 15, 25 }; // SPD
+ Float_t regParamPol1[] = { 0.15, 0.25, 0.5 };
+ }
+ else if (is2360TeV)
+ {
+ Float_t regParamPol0[] = { 10, 12, 40 };
+ Float_t regParamPol1[] = { 0.25, 0.25, 2 };
+ }
+ else
+ {
+ Float_t regParamPol0[] = { 1, 25, 10 };
+ Float_t regParamPol1[] = { 0.15, 0.5, 0.5 };
+ AliUnfolding::SetSkipBinsBegin(3);
+ }
+
+ Int_t reg = AliUnfolding::kPol0;
+ if (beta > 0)
+ reg = AliUnfolding::kPol1;
+
+ Float_t regParam = TMath::Abs(beta);
+ if (histID == -1)
+ {
+ if (beta > 0)
+ regParam = regParamPol1[hID];
+ else
+ regParam = regParamPol0[hID];
+ }
+ AliUnfolding::SetChi2Regularization(reg, regParam);
+
+ //AliUnfolding::SetChi2Regularization(AliUnfolding::kLog, 1000000);
+ //AliUnfolding::SetChi2Regularization(AliUnfolding::kRatio, 10);
+ //TVirtualFitter::SetDefaultFitter("Minuit2");
+
//AliUnfolding::SetCreateOverflowBin(kFALSE);
//mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0); //mult->SetCreateBigBin(kFALSE);
//mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0, 125); mult->SetCreateBigBin(kFALSE);
// AliUnfolding::SetChi2Regularization(AliUnfolding::kEntropy, beta);
//mult->SetRegularizationParameters(AliMultiplicityCorrection::kLog, 1e5);
- //mult->ApplyMinuitFit(histID, fullPhaseSpace, AliMultiplicityCorrection::kTrVtx, kTRUE, mcCompare);
- //mult->SetMultiplicityESDCorrected(histID, (TH1F*) mcCompare);
+ if (0)
+ {
+ // part for checking
+ TH1* mcCompare = hist2->ProjectionY("mcmchist", 1, hist2->GetNbinsX());
+ mcCompare->Sumw2();
+ mult->ApplyMinuitFit(hID, fullPhaseSpace, AliMultiplicityCorrection::kMB, 0, kTRUE, mcCompare);
+ mult->SetMultiplicityESDCorrected(hID, (TH1F*) mcCompare);
+ }
+ else
+ {
+ // Zero Bin
+ Int_t zeroBin = 0;
+ if (is900GeV) // from data
+ {
+ // background subtraction
+ Int_t background = 0;
+
+ //background = 1091 + 4398; // MB1 for run 104824 - 52 (SPD)
+ //background = 1087 + 4398; // MB1 for run 104824 - 52 (TPCITS)
+
+ //background = 417 + 1758; // MB1 for run 104867 - 92 (SPD)
+ //background = 1162+422; // MB1 for run 104892 (TPCITS)
+ //background = 5830 + 1384; // MB1 for run 104824,25,45,52,67,92 (TPC runs) (TPCITS)
+
+ background = 20 + 0; // V0AND for run 104824 - 52
+ //background = 10 + 0; // V0AND for run 104867 - 92
+
+ Printf("NOTE: Subtracting %d background events", background);
+ gSystem->Sleep(1000);
+
+ zeroBin = mult->GetTriggeredEvents(hID)->GetBinContent(1) - background - mult->GetMultiplicityESD(hID)->Integral(0, 2, 1, 1);
+ }
+ else if (is2360TeV)
+ {
+ // from MC
+ Float_t fractionZeroBin = (multTrigger->GetTriggeredEvents(hID)->GetBinContent(1) - multTrigger->GetMultiplicityESD(hID)->Integral(0, 2, 1, 1)) / multTrigger->GetMultiplicityESD(hID)->Integral(1, 1, 1, mult->GetMultiplicityESD(hID)->GetNbinsY());
+ zeroBin = fractionZeroBin * mult->GetMultiplicityESD(hID)->Integral(1, 1, 1, mult->GetMultiplicityESD(hID)->GetNbinsY());
+
+ Printf("Zero bin from MC: Estimating %d events with trigger but without vertex", zeroBin);
+ gSystem->Sleep(1000);
+ }
+ else
+ {
+ AliUnfolding::SetSkip0BinInChi2(kTRUE);
+ }
- mult->ApplyMinuitFit(hID, fullPhaseSpace, eventType, kFALSE, 0, kFALSE); //hist2->ProjectionY("mymchist"));
- //mult->ApplyNBDFit(histID, fullPhaseSpace, eventType);
+ //mult->SetVertexRange(3, 4);
+ mult->ApplyMinuitFit(hID, fullPhaseSpace, eventType, zeroBin, kFALSE, 0, calcBias);
+ }
}
else
{
- mult->ApplyBayesianMethod(hID, fullPhaseSpace, eventType, 1, 10, 0, kTRUE);
+ // HACK
+ //mult->GetMultiplicityESD(hID)->SetBinContent(1, 1, 0);
+ //for (Int_t bin=1; bin<=mult->GetCorrelation(hID)->GetNbinsY(); bin++)
+ // mult->GetCorrelation(hID)->SetBinContent(1, bin, 1, 0);
+ AliUnfolding::SetChi2Regularization(AliUnfolding::kNone, 0);
+ mult->ApplyBayesianMethod(hID, fullPhaseSpace, eventType, 1, beta, 0, kTRUE);
//mult->ApplyBayesianMethod(hID, fullPhaseSpace, eventType, 1, 5, 0, kFALSE);
}
-
- mult->SetMultiplicityMC(hID, eventType, hist2);
}
Printf("Writing result to %s", targetFile);
TFile* file = TFile::Open(targetFile, "RECREATE");
- mult->SaveHistograms();
+ mult->SaveHistograms("Multiplicity");
file->Write();
file->Close();
+ if (histID == -1)
+ return;
+
for (hID = ((histID == -1) ? 0 : histID); hID <= ((histID == -1) ? 2 : histID); hID++)
{
TH2F* hist2 = esd->GetMultiplicityMC(hID, eventType);
TH1* mcCompare = hist2->ProjectionY("mcmchist", 1, hist2->GetNbinsX());
- mult->DrawComparison(Form("%s_%d", (chi2) ? "MinuitChi2" : "Bayesian", hID), hID, fullPhaseSpace, kTRUE, mcCompare);
+ mult->DrawComparison(Form("%s_%d_%f", (chi2) ? "MinuitChi2" : "Bayesian", hID, beta), hID, fullPhaseSpace, kTRUE, mcCompare, kFALSE, eventType);
+
+ Printf("<n> = %f", mult->GetMultiplicityESDCorrected(hID)->GetMean());
}
}
-TH1* GetChi2Bias(Float_t alpha)
+void correctAll(Int_t eventType)
+{
+ const char* suffix = "";
+ switch (eventType)
+ {
+ case 0: suffix = "trvtx"; break;
+ case 1: suffix = "mb"; break;
+ case 2: suffix = "inel"; break;
+ case 3: suffix = "nsd"; break;
+ }
+
+ correct("multiplicityMC.root", "Multiplicity", "multiplicityESD.root", 1, -1, 0, -1, eventType, TString(Form("chi2a_%s.root", suffix)));
+ correct("multiplicityMC.root", "Multiplicity", "multiplicityESD.root", 1, -1, 0, 1, eventType, TString(Form("chi2b_%s.root", suffix)));
+ correct("multiplicityMC.root", "Multiplicity", "multiplicityESD.root", 0, -1, 0, 40, eventType, TString(Form("bayes_%s.root", suffix)));
+
+ if (eventType == 3)
+ drawAll(1);
+ else if (eventType == 2)
+ drawAllINEL();
+}
+
+void drawAll(Bool_t showUA5 = kFALSE)
+{
+ const char* files[] = { "chi2a_nsd.root", "chi2b_nsd.root", "bayes_nsd.root" };
+ drawAll(files, showUA5);
+}
+
+void drawAllINEL()
+{
+ const char* files[] = { "chi2a_inel.root", "chi2b_inel.root", "bayes_inel.root" };
+ drawAll(files);
+}
+
+void drawAllMB()
+{
+ const char* files[] = { "chi2a_mb.root", "chi2b_mb.root", "bayes_mb.root" };
+ drawAll(files);
+}
+
+void drawAll(const char** files, Bool_t showUA5 = kFALSE, Bool_t normalize = kFALSE)
{
loadlibs();
- const char* fileNameMC = "multiplicityMC.root";
- const char* folder = "Multiplicity";
- const char* fileNameESD = "multiplicityESD.root";
- const char* folderESD = "Multiplicity";
+ Int_t colors[] = { 1, 2, 4 };
+
+ c = new TCanvas(Form("c%d", gRandom->Uniform(100)), "c", 1800, 600);
+ c->Divide(3, 1);
+
+ l = new TLegend(0.6, 0.6, 0.99, 0.99);
+ l->SetFillColor(0);
+
+ TH1* hist0 = 0;
+
+ for (Int_t hID=0; hID<3; hID++)
+ {
+ c->cd(hID+1)->SetLogy();
+ gPad->SetGridx();
+ gPad->SetGridy();
+ for (Int_t i=0; i<3; i++)
+ {
+ TFile::Open(files[i]);
+
+ hist = (TH1*) gFile->Get(Form("Multiplicity/fMultiplicityESDCorrected%d", hID));
- AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(fileNameMC, folder);
- AliMultiplicityCorrection* esd = AliMultiplicityCorrection::Open(fileNameESD, folderESD);
+ Float_t average0 = hist->GetMean();
+ hist2 = (TH1*) hist->Clone("temp");
+ hist2->SetBinContent(1, 0);
+ Float_t average1 = hist2->GetMean();
+ Printf("%d: <N> = %.2f <N>/(eta) = %.2f | without 0: <N> = %.2f <N>/(eta) = %.2f", hID, average0, average0 / ((hID+1) - 0.4 * (hID / 2)), average1, average1 / ((hID+1) - 0.4 * (hID / 2)));
+
+ hist->SetLineColor(colors[i]);
+
+ hist->SetStats(0);
+ hist->GetXaxis()->SetRangeUser(0, kBinLimits[hID]);
+
+ Float_t min = 0.1;
+ Float_t max = hist->GetMaximum() * 1.5;
+
+ if (normalize)
+ min = 1e-6;
+
+ hist->GetYaxis()->SetRangeUser(min, max);
+ hist->SetTitle(";unfolded multiplicity;events");
+
+ if (hID == 0)
+ {
+ l->AddEntry(hist, files[i]);
+ }
+
+ if (hist->Integral() <= 0)
+ continue;
+
+ if (normalize)
+ hist->Scale(1.0 / hist->Integral());
+
+ AliPWG0Helper::NormalizeToBinWidth(hist);
+
+ hist->DrawCopy((i == 0) ? "" : "SAME");
+
+ if (!hist0)
+ hist0 = (TH1*) hist->Clone("hist0");
+ }
+
+ if (hist0)
+ {
+ line = new TLine(kTrustLimits[hID], hist0->GetMinimum(), kTrustLimits[hID], hist0->GetMaximum());
+ line->SetLineWidth(3);
+ line->Draw();
+ }
+ }
- //Int_t hID = 0; const Int_t kMaxBins = 35;
- //Int_t hID = 1; const Int_t kMaxBins = 60;
- Int_t hID = 2; const Int_t kMaxBins = 70;
- AliUnfolding::SetNbins(kMaxBins, kMaxBins);
- //AliUnfolding::SetDebug(1);
- //AliUnfolding::SetChi2Regularization(AliUnfolding::kPol0, 50);
- AliUnfolding::SetChi2Regularization(AliUnfolding::kPol0, alpha);
+ c->cd(1);
+ l->Draw();
- TH2F* hist = esd->GetMultiplicityESD(hID);
- mult->SetMultiplicityESD(hID, hist);
+ if (showUA5)
+ {
+ TGraphErrors *gre = new TGraphErrors(23);
+ gre->SetName("Graph");
+ gre->SetTitle("UA5");
+ gre->SetFillColor(1);
+ gre->SetMarkerStyle(24);
+ gre->SetPoint(0,0,0.15);
+ gre->SetPointError(0,0.5,0.01);
+ gre->SetPoint(1,1,0.171);
+ gre->SetPointError(1,0.5,0.007);
+ gre->SetPoint(2,2,0.153);
+ gre->SetPointError(2,0.5,0.007);
+ gre->SetPoint(3,3,0.124);
+ gre->SetPointError(3,0.5,0.006);
+ gre->SetPoint(4,4,0.099);
+ gre->SetPointError(4,0.5,0.005);
+ gre->SetPoint(5,5,0.076);
+ gre->SetPointError(5,0.5,0.005);
+ gre->SetPoint(6,6,0.057);
+ gre->SetPointError(6,0.5,0.004);
+ gre->SetPoint(7,7,0.043);
+ gre->SetPointError(7,0.5,0.004);
+ gre->SetPoint(8,8,0.032);
+ gre->SetPointError(8,0.5,0.003);
+ gre->SetPoint(9,9,0.024);
+ gre->SetPointError(9,0.5,0.003);
+ gre->SetPoint(10,10,0.018);
+ gre->SetPointError(10,0.5,0.002);
+ gre->SetPoint(11,11,0.013);
+ gre->SetPointError(11,0.5,0.002);
+ gre->SetPoint(12,12,0.01);
+ gre->SetPointError(12,0.5,0.002);
+ gre->SetPoint(13,13,0.007);
+ gre->SetPointError(13,0.5,0.001);
+ gre->SetPoint(14,14,0.005);
+ gre->SetPointError(14,0.5,0.001);
+ gre->SetPoint(15,15,0.004);
+ gre->SetPointError(15,0.5,0.001);
+ gre->SetPoint(16,16,0.0027);
+ gre->SetPointError(16,0.5,0.0009);
+ gre->SetPoint(17,17,0.0021);
+ gre->SetPointError(17,0.5,0.0008);
+ gre->SetPoint(18,18,0.0015);
+ gre->SetPointError(18,0.5,0.0007);
+ gre->SetPoint(19,19,0.0011);
+ gre->SetPointError(19,0.5,0.0006);
+ gre->SetPoint(20,20,0.0008);
+ gre->SetPointError(20,0.5,0.0005);
+ gre->SetPoint(21,22,0.00065);
+ gre->SetPointError(21,1,0.0003);
+ gre->SetPoint(22,27.5,0.00013);
+ gre->SetPointError(22,3.5,7.1e-05);
+
+ for (Int_t i=0; i<gre->GetN(); i++)
+ {
+ gre->GetY()[i] *= hist0->Integral();
+ gre->GetEY()[i] *= hist0->Integral();
+ }
+
+ gre->Draw("P");
+
+ ratio = (TGraphErrors*) gre->Clone("ratio");
+
+ for (Int_t i=0; i<gre->GetN(); i++)
+ {
+ //Printf("%f %d", hist0->GetBinContent(hist0->FindBin(ratio->GetX()[i])), hist0->FindBin(ratio->GetX()[i]));
+ Int_t bin = hist0->FindBin(gre->GetX()[i]);
+ if (hist0->GetBinContent(bin) > 0)
+ {
+ ratio->GetEY()[i] = TMath::Sqrt((ratio->GetEY()[i] / ratio->GetY()[i])**2 + (hist0->GetBinError(bin) / hist0->GetBinContent(bin))**2);
+ ratio->GetY()[i] /= hist0->GetBinContent(bin);
+ ratio->GetEY()[i] *= ratio->GetY()[i];
+ }
+ }
+ new TCanvas;
+ gPad->SetGridx();
+ gPad->SetGridy();
+ //ratio->Print();
+ ratio->Draw("AP");
+ ratio->GetYaxis()->SetRangeUser(0.5, 1.5);
+ }
+
+ c->SaveAs("draw_all.png");
+}
+
+TH1* GetChi2Bias(Float_t alpha = -5)
+{
+ loadlibs();
+
+ AliMultiplicityCorrection::EventType eventType = AliMultiplicityCorrection::kMB;
+ Int_t hID = 1;
+ correct("multiplicityMC.root", "Multiplicity", "multiplicityESD.root", 1, hID, 0, alpha, eventType, "nobias.root", "Multiplicity", kFALSE);
+
+ correct("multiplicityMC.root", "Multiplicity", "multiplicityESD.root", 1, hID, 0, alpha, eventType, "withbias.root", "Multiplicity", kTRUE);
+
// without bias calculation
- mult->ApplyMinuitFit(hID, kFALSE, 0, kFALSE);
+ mult = AliMultiplicityCorrection::Open("nobias.root");
baseold = (TH1*) mult->GetMultiplicityESDCorrected(hID)->Clone("baseold");
// with bias calculation
- mult->ApplyMinuitFit(hID, kFALSE, 0, kFALSE, 0, kTRUE);
+ mult = AliMultiplicityCorrection::Open("withbias.root");
base = (TH1*) mult->GetMultiplicityESDCorrected(hID)->Clone("base");
// relative error plots
ratio2 = (TH1*) mult->GetMultiplicityESDCorrected(hID)->Clone("hist1");
ratio2->Reset();
- for (Int_t t = 0; t<kMaxBins; t++)
+ for (Int_t t = 0; t<baseold->GetNbinsX(); t++)
{
Printf("Bin %d; Content: %f; Chi2 Error: %f; Bias: %f; In percent: %.2f %.2f", t+1, base->GetBinContent(t+1), baseold->GetBinError(t+1), base->GetBinError(t+1), (base->GetBinContent(t+1) > 0) ? 100.0 * baseold->GetBinError(t+1) / base->GetBinContent(t+1) : -1, (base->GetBinContent(t+1) > 0) ? 100.0 * base->GetBinError(t+1) / base->GetBinContent(t+1) : -1);
if (base->GetBinContent(t+1) <= 0)
//new TCanvas; base->Draw(); gPad->SetLogy();
- new TCanvas;
+ c = new TCanvas;
ratio1->GetYaxis()->SetRangeUser(0, 1);
- ratio1->GetXaxis()->SetRangeUser(0, kMaxBins);
+ ratio1->GetXaxis()->SetRangeUser(0, 50);
ratio1->Draw();
ratio2->SetLineColor(2);
ratio2->Draw("SAME");
hist2->Draw("SAME");
}
-void DataScan(Bool_t redo = kTRUE)
+void DrawUnfoldingLimit(Int_t hID, Float_t min, Float_t max)
+{
+ line = new TLine(kTrustLimits[hID], min, kTrustLimits[hID], max);
+ line->SetLineWidth(2);
+ line->Draw();
+}
+
+void DataScan(Int_t hID, Bool_t redo = kTRUE)
{
// makes a set of unfolded spectra and compares
// don't forget FindUnfoldedLimit in plots.C
loadlibs();
// files...
- Bool_t mc = kTRUE;
+ Bool_t mc = 1;
const char* fileNameMC = "multiplicityMC.root";
const char* folder = "Multiplicity";
const char* fileNameESD = "multiplicityESD.root";
const char* folderESD = "Multiplicity";
- Int_t hID = 0; const Int_t kMaxBins = 35;
- //Int_t hID = 1; const Int_t kMaxBins = 60;
- //Int_t hID = 2; const Int_t kMaxBins = 70;
+ const Int_t kMaxBins = kBinLimits[hID];
+
AliMultiplicityCorrection::EventType eventType = AliMultiplicityCorrection::kTrVtx;
Bool_t evaluteBias = kFALSE;
// chi2 range
AliUnfolding::RegularizationType regBegin = AliUnfolding::kPol0;
- AliUnfolding::RegularizationType regEnd = AliUnfolding::kPol0;
- Float_t regParamBegin[] = { 0, 1, 10 };
- Float_t regParamEnd[] = { 0, 40, 101 };
- Float_t regParamStep[] = { 0, TMath::Sqrt(10), TMath::Sqrt(10) };
+ AliUnfolding::RegularizationType regEnd = AliUnfolding::kPol1;
+ Float_t regParamBegin[] = { 0, 1, 0.2, 3 };
+ Float_t regParamEnd[] = { 0, 11, 0.5, 31 };
+ Float_t regParamStep[] = { 0, 2, 2, TMath::Sqrt(10) };
// bayesian range
- Int_t iterBegin = 5;
- Int_t iterEnd = 21;
- Int_t iterStep = 5;
+ Int_t iterBegin = 40;
+ Int_t iterEnd = 41;
+ Int_t iterStep = 20;
TList labels;
mult->SetMultiplicityESD(hID, esd->GetMultiplicityESD(hID));
+ // insert trigger efficiency in flat response matrix
+ AliMultiplicityCorrection* multTrigger = AliMultiplicityCorrection::Open("multiplicityTrigger.root");
+ for (Int_t evType = AliMultiplicityCorrection::kTrVtx; evType <= AliMultiplicityCorrection::kNSD; evType++)
+ mult->SetMultiplicityMC(hID, evType, multTrigger->GetMultiplicityMC(hID, evType));
+
+ mult->SetNoVertexEvents(hID, multTrigger->GetNoVertexEvents(hID));
+ mult->FixTriggerEfficiencies(10);
+
+ Float_t fraction0Generated = multTrigger->GetFraction0Generated(hID);
+
if (mc)
mult->SetMultiplicityMC(hID, eventType, esd->GetMultiplicityMC(hID, eventType));
AliUnfolding::SetChi2Regularization(reg, regParam);
- mult->ApplyMinuitFit(hID, kFALSE, eventType, kFALSE, 0, evaluteBias);
+ mult->ApplyMinuitFit(hID, kFALSE, eventType, fraction0Generated, kFALSE, 0, evaluteBias);
file = TFile::Open(Form("datascan_%d.root", count), "RECREATE");
mult->SaveHistograms();
TH1* residualSum = new TH1F("residualSum", ";;residual squared sum", count + 1, 0.5, count + 1.5);
+ Int_t allowedCount = count;
count = 0;
while (1)
{
mult = AliMultiplicityCorrection::Open(Form("datascan_%d.root", count++));
if (!mult)
break;
+ if (count > allowedCount+1)
+ break;
hist = (TH1*) mult->GetMultiplicityESDCorrected(hID);
hist->SetLineColor((count-1) % 4 + 1);
hist->SetLineStyle((count-1) / 4 + 1);
hist->GetXaxis()->SetRangeUser(0, kMaxBins);
+ hist->GetYaxis()->SetRangeUser(0.1, hist->GetMaximum() * 1.5);
hist->SetStats(kFALSE);
hist->SetTitle("");
c->cd(1);
hist->DrawCopy((count == 1) ? "HIST" : "HISTSAME");
+ DrawUnfoldingLimit(hID, 0.1, hist->GetMaximum());
+
if (mc)
{
TH2* mcHist2d = mult->GetMultiplicityMC(hID, eventType);
c->cd(1);
mcHist->SetMarkerStyle(24);
mcHist->Draw("P SAME");
-
+
c->cd(2);
// calculate ratio using only the error on the mc bin
ratio = (TH1*) hist->Clone("ratio");
ratio->GetYaxis()->SetRangeUser(0.5, 1.5);
ratio->GetYaxis()->SetTitle("Ratio unfolded / MC");
ratio->DrawCopy((count == 1) ? "HIST" : "HISTSAME");
+
+ DrawUnfoldingLimit(hID, 0.5, 1.5);
}
c->cd(3);
ratio->GetYaxis()->SetTitle("Ratio unfolded / unfolded reference case");
ratio->DrawCopy((count == 1) ? "" : "SAME");
+ DrawUnfoldingLimit(hID, 0.5, 1.5);
+
c->cd(4);
ratio = (TH1*) hist->Clone("ratio");
for (Int_t bin=1; bin<=ratio->GetNbinsX(); bin++)
ratio->GetYaxis()->SetTitle("Relative error");
ratio->DrawCopy((count == 1) ? "" : "SAME");
+ DrawUnfoldingLimit(hID, 0, 1);
+
c->cd(5);
Float_t sum;
residuals = mult->GetResiduals(hID, AliMultiplicityCorrection::kTrVtx, sum);
name.Form("Multiplicity%d", i);
TFile::Open(files[i]);
+ Printf("Loading from file %s", files[i]);
data[i] = new AliMultiplicityCorrection(name, name);
data[i]->LoadHistograms("Multiplicity");
if (i > 0)
{
//
// builds several response matrices with different particle ratios (systematic study)
+ // particle-species study
//
// WARNING doesn't work uncompiled, see test.C
Int_t secondaries = 0;
Int_t doubleCount = 0;
- for (Int_t num = 0; num < 9; num++)
+ for (Int_t num = 0; num < 7; num++)
{
+ Printf("%d", num);
+
TString name;
name.Form("Multiplicity_%d", num);
AliMultiplicityCorrection* fMultiplicity = new AliMultiplicityCorrection(name, name);
for (Int_t i = 0; i < 4; i++)
ratio[i] = 1;
- switch (num)
- {
- case 1 : ratio[1] = 0.5; break;
- case 2 : ratio[1] = 1.5; break;
- case 3 : ratio[2] = 0.5; break;
- case 4 : ratio[2] = 1.5; break;
- case 5 : ratio[1] = 0.5; ratio[2] = 0.5; break;
- case 6 : ratio[1] = 1.5; ratio[2] = 1.5; break;
- case 7 : ratio[1] = 0.5; ratio[2] = 1.5; break;
- case 8 : ratio[1] = 1.5; ratio[2] = 0.5; break;
- }
-
+ if (num == 1)
+ ratio[1] = 0.7;
+ if (num == 2)
+ ratio[1] = 1.3;
+
+ if (num == 3)
+ ratio[2] = 0.7;
+ if (num == 4)
+ ratio[2] = 1.3;
+
+ if (num == 5)
+ ratio[3] = 0.7;
+ if (num == 6)
+ ratio[3] = 1.3;
+
for (Int_t i=0; i<fParticleSpecies->GetEntries(); i++)
{
fParticleSpecies->GetEvent(i);
+
Float_t* f = fParticleSpecies->GetArgs();
for (Int_t j = 0; j < 4; j++)
{
- gene += ratio[j] * f[j+1];
- meas += ratio[j] * f[j+1+4];
+ if (ratio[j] == 1)
+ {
+ gene += f[j+1];
+ }
+ else
+ {
+ for (Int_t k=0; k<f[j+1]; k++)
+ {
+ if (ratio[j] < 1 && gRandom->Uniform() < ratio[j])
+ {
+ gene += 1;
+ }
+ else if (ratio[j] > 1)
+ {
+ gene += 1;
+ if (gRandom->Uniform() < ratio[j] - 1)
+ gene += 1;
+ }
+ }
+ }
+
+ if (ratio[j] == 1)
+ {
+ meas += f[j+1+4];
+ }
+ else
+ {
+ for (Int_t k=0; k<f[j+1+4]; k++)
+ {
+ if (ratio[j] < 1 && gRandom->Uniform() < ratio[j])
+ {
+ meas += 1;
+ }
+ else if (ratio[j] > 1)
+ {
+ meas += 1;
+ if (gRandom->Uniform() < ratio[j] - 1)
+ meas += 1;
+ }
+ }
+ }
+
tracks += f[j+1+4];
}
//Printf("%.f %.f %.f %.f %.f", f[5], f[6], f[7], f[8], f[9]);
- fMultiplicity->FillCorrection(f[0], gene, gene, gene, gene, 0, meas, meas, meas, meas);
+ fMultiplicity->FillCorrection(f[0], gene, gene, gene, gene, meas, meas, meas);
// HACK all as kND = 1
- fMultiplicity->FillGenerated(f[0], kTRUE, kTRUE, 1, gene, gene, gene, gene, 0);
- fMultiplicity->FillMeasured(f[0], meas, meas, meas, meas);
+ fMultiplicity->FillGenerated(f[0], kTRUE, kTRUE, 1, gene, gene, gene, gene);
+ fMultiplicity->FillMeasured(f[0], meas, meas, meas);
}
//fMultiplicity->DrawHistograms();
{
loadlibs();
- for (Int_t num = 0; num < 9; num++)
+ for (Int_t num = 0; num < 7; num++)
{
TString target;
- target.Form("chi2_species_%d.root", num);
- correct("species.root", Form("Multiplicity_%d", num), "species.root", 1, 1, 0, 1e5, 0, target, "Multiplicity_0");
+ target.Form("chi2a_inel_species_%d.root", num);
+
+ correct("compositions.root", Form("Multiplicity_%d", num), "multiplicityESD.root", 1, -1, 0, -1, 2, target);
}
}
-void MergeModifyCrossSection(const char* output = "multiplicityMC_xsection.root")
+void MergeCrossSection(Int_t xsectionID, const char* output = "multiplicityMC_xsection.root")
{
- const char* files[] = { "multiplicityMC_nd.root", "multiplicityMC_sd.root", "multiplicityMC_dd.root" };
+ const char* files[] = { "multiplicitySD.root", "multiplicityDD.root", "multiplicityND.root" };
loadlibs();
TFile::Open(output, "RECREATE");
gFile->Close();
- for (Int_t num = 0; num < 7; num++)
- {
- AliMultiplicityCorrection* data[3];
- TList list;
-
- Float_t ratio[3];
- switch (num)
- {
- case 0: ratio[0] = 1.0; ratio[1] = 1.0; ratio[2] = 1.0; break;
- case 1: ratio[0] = 1.0; ratio[1] = 1.5; ratio[2] = 1.0; break;
- case 2: ratio[0] = 1.0; ratio[1] = 0.5; ratio[2] = 1.0; break;
- case 3: ratio[0] = 1.0; ratio[1] = 1.0; ratio[2] = 1.5; break;
- case 4: ratio[0] = 1.0; ratio[1] = 1.0; ratio[2] = 0.5; break;
- case 5: ratio[0] = 1.0; ratio[1] = 1.5; ratio[2] = 1.5; break;
- case 6: ratio[0] = 1.0; ratio[1] = 0.5; ratio[2] = 0.5; break;
- default: return;
- }
+ AliMultiplicityCorrection* data[3];
+ TList list;
- for (Int_t i=0; i<3; ++i)
- {
- TString name;
- name.Form("Multiplicity_%d", num);
- if (i > 0)
- name.Form("Multiplicity_%d_%d", num, i);
+ Float_t ref_SD = -1;
+ Float_t ref_DD = -1;
+ Float_t ref_ND = -1;
- TFile::Open(files[i]);
- data[i] = new AliMultiplicityCorrection(name, name);
- data[i]->LoadHistograms("Multiplicity");
+ Float_t error_SD = -1;
+ Float_t error_DD = -1;
+ Float_t error_ND = -1;
- // modify x-section
- for (Int_t j=0; j<AliMultiplicityCorrection::kMCHists; j++)
- {
- data[i]->GetMultiplicityVtx(j)->Scale(ratio[i]);
- data[i]->GetMultiplicityMB(j)->Scale(ratio[i]);
- data[i]->GetMultiplicityINEL(j)->Scale(ratio[i]);
- data[i]->GetMultiplicityNSD(j)->Scale(ratio[i]);
- }
+ gROOT->ProcessLine(gSystem->ExpandPathName(".L $ALICE_ROOT/PWG0/dNdEta/drawSystematics.C"));
+ GetRelativeFractions(xsectionID, ref_SD, ref_DD, ref_ND, error_SD, error_DD, error_ND);
+
+ for (Int_t i=0; i<3; ++i)
+ {
+ TString name;
+ name.Form("Multiplicity");
+ if (i > 0)
+ name.Form("Multiplicity_%d", i);
- for (Int_t j=0; j<AliMultiplicityCorrection::kESDHists; j++)
- data[i]->GetMultiplicityESD(j)->Scale(ratio[i]);
+ TFile::Open(files[i]);
+ data[i] = new AliMultiplicityCorrection(name, name);
+ data[i]->LoadHistograms("Multiplicity");
+ }
+
+ // TODO is the under/overflow bin scaled as well? --> seems like, verify anyway!
- for (Int_t j=0; j<AliMultiplicityCorrection::kCorrHists; j++)
- data[i]->GetCorrelation(j)->Scale(ratio[i]);
+ // calculating relative
+ Float_t sd = data[0]->GetMultiplicityINEL(0)->Integral(0, data[0]->GetMultiplicityINEL(0)->GetNbinsX()+1);
+ Float_t dd = data[1]->GetMultiplicityINEL(0)->Integral(0, data[1]->GetMultiplicityINEL(0)->GetNbinsX()+1);
+ Float_t nd = data[2]->GetMultiplicityINEL(0)->Integral(0, data[2]->GetMultiplicityINEL(0)->GetNbinsX()+1);
+ Float_t total = nd + dd + sd;
+
+ nd /= total;
+ sd /= total;
+ dd /= total;
+
+ Printf("Ratios in the correction map are: ND=%f, DD=%f, SD=%f", nd, dd, sd);
+
+ Float_t ratio[3];
+ ratio[0] = ref_SD / sd;
+ ratio[1] = ref_DD / dd;
+ ratio[2] = ref_ND / nd;
+
+ Printf("SD=%.2f, DD=%.2f, ND=%.2f",ratio[0], ratio[1], ratio[2]);
+
+ for (Int_t i=0; i<3; ++i)
+ {
+ // modify x-section
+ for (Int_t j=0; j<AliMultiplicityCorrection::kMCHists; j++)
+ {
+ data[i]->GetMultiplicityVtx(j)->Scale(ratio[i]);
+ data[i]->GetMultiplicityMB(j)->Scale(ratio[i]);
+ data[i]->GetMultiplicityINEL(j)->Scale(ratio[i]);
+ data[i]->GetMultiplicityNSD(j)->Scale(ratio[i]);
+ }
- if (i > 0)
- list.Add(data[i]);
+ for (Int_t j=0; j<AliMultiplicityCorrection::kESDHists; j++)
+ {
+ data[i]->GetMultiplicityESD(j)->Scale(ratio[i]);
+ data[i]->GetTriggeredEvents(j)->Scale(ratio[i]);
+ data[i]->GetNoVertexEvents(j)->Scale(ratio[i]);
}
+
+ for (Int_t j=0; j<AliMultiplicityCorrection::kCorrHists; j++)
+ data[i]->GetCorrelation(j)->Scale(ratio[i]);
+
+ if (i > 0)
+ list.Add(data[i]);
+ }
- printf("Case %d, %s: Entries in response matrix 3: ND: %.2f SD: %.2f DD: %.2f", num, data[0]->GetName(), data[0]->GetCorrelation(3)->Integral(), data[1]->GetCorrelation(3)->Integral(), data[2]->GetCorrelation(3)->Integral());
+ printf("Case %d, %s: Entries in response matrix: SD: %.2f DD: %.2f ND: %.2f", xsectionID, data[0]->GetName(), data[0]->GetCorrelation(0)->Integral(), data[1]->GetCorrelation(0)->Integral(), data[2]->GetCorrelation(0)->Integral());
- data[0]->Merge(&list);
+ data[0]->Merge(&list);
- Printf(" Total: %.2f", data[0]->GetCorrelation(3)->Integral());
+ Printf(" Total: %.2f", data[0]->GetCorrelation(0)->Integral());
- TFile::Open(output, "UPDATE");
- data[0]->SaveHistograms();
- gFile->Close();
+ TFile::Open(output, "RECREATE");
+ data[0]->SaveHistograms();
+ gFile->Close();
- list.Clear();
+ list.Clear();
- for (Int_t i=0; i<3; ++i)
- delete data[i];
- }
+ for (Int_t i=0; i<3; ++i)
+ delete data[i];
}
void Rebin(const char* fileName = "multiplicityMC_3M.root", Int_t corrMatrix = 3)
canvas->SaveAs(Form("%s.eps", canvas->GetName()));
}
}
+
+void MergeDistributions()
+{
+ loadlibs();
+
+ const char* dir1 = "run000104824-52_pass4";
+ const char* dir2 = "run000104867_90_92_pass4";
+
+ for (Int_t evType = 0; evType < 2; evType++)
+ {
+ Printf("%d", evType);
+
+ const char* evTypeStr = ((evType == 0) ? "inel" : "nsd");
+
+ const char* id = "chi2a";
+ //const char* id = "bayes";
+
+ TString suffix;
+ suffix.Form("/all/spd/%s_", id);
+ if (evType == 1)
+ suffix.Form("/v0and/spd/%s_", id);
+
+ TString file1, file2;
+ file1.Form("%s%s%%s.root", dir1, suffix.Data());
+ file2.Form("%s%s%%s.root", dir2, suffix.Data());
+
+ if (1)
+ {
+ const char* files[] = { Form(file1.Data(), evTypeStr), Form(file2.Data(), evTypeStr) };
+ Merge(2, files, Form("merged/%s_%s.root", id, evTypeStr));
+ }
+ else
+ {
+ AliMultiplicityCorrection* mult1 = AliMultiplicityCorrection::Open(Form(file1.Data(), evTypeStr));
+ AliMultiplicityCorrection* mult2 = AliMultiplicityCorrection::Open(Form(file2.Data(), evTypeStr));
+
+ AliMultiplicityCorrection* target = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+
+ for (Int_t etaRange = 0; etaRange < 3; etaRange++)
+ {
+ hist1 = mult1->GetMultiplicityESDCorrected(etaRange);
+ hist2 = mult2->GetMultiplicityESDCorrected(etaRange);
+ targetHist = target->GetMultiplicityESDCorrected(etaRange);
+
+ //hist1->Scale(1.0 / hist1->Integral());
+ //hist2->Scale(1.0 / hist2->Integral());
+
+ //targetHist->SetBinContent(1, hist1->GetBinContent(1) * 0.5 + hist2->GetBinContent(1) * 0.5);
+ targetHist->SetBinContent(1, hist1->GetBinContent(1) + hist2->GetBinContent(1));
+ for (Int_t i=1; i<=hist1->GetNbinsX(); i++)
+ {
+ if (hist1->GetBinError(i) > 0 && hist2->GetBinError(i) > 0)
+ {
+ Float_t w1 = 1.0 / hist1->GetBinError(i) / hist1->GetBinError(i);
+ Float_t w2 = 1.0 / hist2->GetBinError(i) / hist2->GetBinError(i);
+
+ Float_t average = (hist1->GetBinContent(i) * w1 + hist2->GetBinContent(i) * w2) / (w1 + w2);
+
+ //targetHist->SetBinContent(i, average);
+ //targetHist->SetBinError(i, TMath::Max(mult1->GetBinError(i), mult2->GetBinError(i)));
+
+ targetHist->SetBinContent(i, hist1->GetBinContent(i) + hist2->GetBinContent(i));
+ targetHist->SetBinError(i, hist1->GetBinError(i) + hist2->GetBinError(i));
+ }
+ }
+ }
+
+ file = TFile::Open(Form("merged/%s_%s.root", id, evTypeStr), "RECREATE");
+ target->SaveHistograms();
+ file->Close();
+ }
+
+ const char* files2[] = { Form(file1.Data(), evTypeStr), Form(file2.Data(), evTypeStr), Form("merged/%s_%s.root", id, evTypeStr) };
+ drawAll(files2, (evType == 1), kTRUE);
+ }
+}
+
+void CompareDistributions(Int_t evType)
+{
+ loadlibs();
+ gROOT->ProcessLine(gSystem->ExpandPathName(".L $ALICE_ROOT/PWG0/dNdEta/drawPlots.C"));
+
+ const char* dir1 = "run000104824-52_pass4";
+ const char* dir2 = "run000104867_90_92_pass4";
+
+ const char* evTypeStr = (evType == 0) ? "inel" : "nsd";
+
+ const char* suffix = "/all/spd/chi2a_";
+ if (evType == 1)
+ suffix = "/v0and/spd/chi2a_";
+
+ TString file1, file2;
+ file1.Form("%s%s%s.root", dir1, suffix, evTypeStr);
+ file2.Form("%s%s%s.root", dir2, suffix, evTypeStr);
+
+ for (Int_t hID = 0; hID < 3; hID++)
+ CompareQualityHists(file1, file2, Form("Multiplicity/fMultiplicityESDCorrected%d", hID), 1, 1);
+}
+
+void StatisticalUncertainty(Int_t methodType, Int_t etaRange, Bool_t mc = kFALSE)
+{
+ loadlibs();
+
+ AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open("multiplicityMC.root");
+ AliMultiplicityCorrection* esd = AliMultiplicityCorrection::Open("multiplicityESD.root");
+ AliMultiplicityCorrection* multTrigger = AliMultiplicityCorrection::Open("multiplicityTrigger.root");
+
+ TH1* mcHist = esd->GetMultiplicityMB(etaRange)->ProjectionY("mymc", 1, 1);
+
+ Int_t geneLimits[] = { 0, 0, 0 };
+
+ LoadAndInitialize(mult, esd, multTrigger, etaRange, kFALSE, geneLimits);
+
+ AliUnfolding::SetNbins(kBinLimits[etaRange], geneLimits[etaRange]);
+ AliUnfolding::SetChi2Regularization(AliUnfolding::kPol0, 5);
+ AliUnfolding::SetBayesianParameters(1, 10);
+
+ // background subtraction
+ Int_t background = 0;
+
+ //background = 1091 + 4398; // MB1 for run 104824 - 52
+ background = 417 + 1758; // MB1 for run 104867 - 92
+
+ //background = 20 + 0; // V0AND for run 104824 - 52
+ //background = 10 + 0; // V0AND for run 104867 - 92
+
+ Printf("NOTE: Subtracting %d background events", background);
+
+ Int_t zeroBin = mult->GetTriggeredEvents(etaRange)->GetBinContent(1) - background - mult->GetMultiplicityESD(etaRange)->Integral(0, 2, 1, 1);
+
+ TH1* errorMeasured = (TH1*) mult->StatisticalUncertainty((AliUnfolding::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kMB, zeroBin, kTRUE, kFALSE, ((mc) ? mcHist : 0))->Clone("errorMeasured");
+
+ new TCanvas; errorMeasured->Draw();
+ new TCanvas;
+
+ AliPWG0Helper::NormalizeToBinWidth(mult->GetMultiplicityESDCorrected(etaRange));
+ mult->GetMultiplicityESDCorrected(etaRange)->Draw();
+
+ if (mc)
+ {
+ AliPWG0Helper::NormalizeToBinWidth(mcHist);
+ mcHist->Scale(mult->GetMultiplicityESDCorrected(etaRange)->Integral() / mcHist->Integral());
+ mcHist->SetLineColor(2);
+ mcHist->Draw("SAME");
+ }
+ gPad->SetLogy();
+
+ TH1* errorResponse = (TH1*) mult->StatisticalUncertainty((AliUnfolding::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kMB, zeroBin, kFALSE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorResponse");
+
+ TH1* errorBoth = (TH1*) mult->StatisticalUncertainty((AliUnfolding::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kMB, zeroBin, kTRUE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorBoth");
+
+ if (mc)
+ {
+ TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
+ DrawResultRatio(mcHist, result, "StatisticalUncertainty2.eps");
+ }
+
+ TFile* file = new TFile(Form("StatisticalUncertaintySPD%s.root", (methodType == 0) ? "Chi2" : "Bayesian"), "RECREATE");
+ errorResponse->Write();
+ errorMeasured->Write();
+ errorBoth->Write();
+ file->Close();
+}
+
+void ChangeResponseMatrixEfficiency(Bool_t reduceEff = kTRUE, Float_t factor = 0.01625)
+{
+ loadlibs();
+
+ AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open("multiplicity.root");
+
+ for (Int_t etaR = 0; etaR < 3; etaR++)
+ {
+ TH3* corr = mult->GetCorrelation(etaR);
+ TH3* corrOld = (TH3*) corr->Clone("corrOld");
+
+ corr->Reset();
+
+ Int_t total = 0;
+ Int_t change = 0;
+ for (Int_t x=0; x<=corr->GetNbinsX()+1; x++)
+ {
+ for (Int_t y=0; y<=corr->GetNbinsY()+1; y++)
+ {
+ Printf("%d", y);
+ for (Int_t z=0; z<=corr->GetNbinsZ()+1; z++)
+ {
+ Float_t tracklets = corr->GetZaxis()->GetBinCenter(z);
+
+ for (Int_t n=0; n<corrOld->GetBinContent(x, y, z); n++)
+ {
+ total += tracklets;
+ Float_t trackletsNew = tracklets;
+
+ for (Int_t i=0; i<tracklets; i++)
+ {
+ if (gRandom->Uniform() < factor)
+ {
+ trackletsNew += (reduceEff) ? -1 : 1;
+ change += (reduceEff) ? -1 : 1;
+ }
+ }
+
+ corr->Fill(corr->GetXaxis()->GetBinCenter(x), corr->GetYaxis()->GetBinCenter(y), trackletsNew);
+ }
+ }
+ }
+ }
+
+ Printf("%d: Changed %d out of %d total tracklets", etaR, change, total);
+ }
+
+ file = TFile::Open("out.root", "RECREATE");
+ mult->SaveHistograms();
+ file->Close();
+}
+
{
tmpStr += " (full phase space)";
}
+ else if (etaR == 2)
+ {
+ tmpStr += " in |#eta| < 1.3";
+ }
else
tmpStr += Form(" in |#eta| < %.1f", (etaR+1)* 0.5);
return Form("%s", tmpStr.Data());
delete clone;
}
-void responseMatrixPlot(const char* fileName = 0)
+TH1* responseMatrixPlot(const char* fileName = 0, const char* folder = "Multiplicity")
{
loadlibs();
- AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+ AliMultiplicityCorrection* mult = new AliMultiplicityCorrection(folder, folder);
if (fileName == 0)
fileName = correctionFile;
TFile::Open(fileName);
- mult->LoadHistograms("Multiplicity");
+ mult->LoadHistograms();
// empty under/overflow bins in x, otherwise Project3D takes them into account
- TH1* hist = mult->GetCorrelation(etaRange);
+ TH2* hist = (TH2*) mult->GetCorrelation(etaRange);
for (Int_t y=0; y<=hist->GetYaxis()->GetNbins()+1; ++y)
{
for (Int_t z=0; z<=hist->GetZaxis()->GetNbins()+1; ++z)
hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
}
}
- hist = ((TH3*) mult->GetCorrelation(etaRange))->Project3D("zy");
+
+ hist = (TH2*) (((TH3*) mult->GetCorrelation(etaRange))->Project3D("zy"));
hist->SetStats(kFALSE);
+
+ if (0)
+ {
+ // normalize
+ // normalize correction for given nPart
+ for (Int_t i=1; i<=hist->GetNbinsX(); ++i)
+ {
+ Double_t sum = hist->Integral(i, i, 1, hist->GetNbinsY());
+ if (sum <= 0)
+ continue;
+ for (Int_t j=1; j<=hist->GetNbinsY(); ++j)
+ {
+ // npart sum to 1
+ hist->SetBinContent(i, j, hist->GetBinContent(i, j) / sum);// * correlation->GetXaxis()->GetBinWidth(i));
+ hist->SetBinError(i, j, hist->GetBinError(i, j) / sum);
+ }
+ }
+ }
+
hist->SetTitle(Form(";%s;%s;Entries", GetMultLabel(), GetMultLabel(etaRange, kFALSE)));
hist->GetXaxis()->SetRangeUser(0, longDisplayRange);
hist->GetYaxis()->SetRangeUser(0, longDisplayRange);
hist->Draw("COLZ");
canvas->SaveAs("responsematrix.eps");
+
+ return hist;
}
void multPythiaPhojet()
return canvas;
}
-TCanvas* DrawRatio(TH1* result, Int_t nResultSyst, TH1** resultSyst, TString epsName, Bool_t firstMarker = kFALSE, const char** legendStrings = 0, Bool_t errors = kFALSE)
+void DrawEfficiencyChange()
+{
+ loadlibs();
+
+ const char* fileName = "chi2a_inel.root";
+
+ base = AliMultiplicityCorrection::Open(Form("spd/%s", fileName));
+ low = AliMultiplicityCorrection::Open(Form("spd-loweff/%s", fileName));
+ high = AliMultiplicityCorrection::Open(Form("spd-higheff/%s", fileName));
+
+ const char* legendStrings[] = { "low efficiency", "high efficiency" };
+
+ file = TFile::Open("systunc_detectorefficiency.root", "RECREATE");
+
+ for (Int_t etaR = 1; etaR < 2; etaR++)
+ {
+ base->GetMultiplicityESDCorrected(etaR)->Scale(1.0 / base->GetMultiplicityESDCorrected(etaR)->Integral(2, base->GetMultiplicityESDCorrected(etaR)->GetNbinsX()));
+
+ TH1* hists[2];
+ hists[0] = low->GetMultiplicityESDCorrected(etaR);
+ hists[0]->Scale(1.0 / hists[0]->Integral(2, hists[0]->GetNbinsX()));
+
+ if (high)
+ {
+ hists[1] = high->GetMultiplicityESDCorrected(etaR);
+ hists[1]->Scale(1.0 / hists[1]->Integral(2, hists[1]->GetNbinsX()));
+ }
+
+ largestError = (TH1*) hists[0]->Clone(Form("detectorefficiency_%d", etaR));
+ largestError->Reset();
+
+ DrawRatio(base->GetMultiplicityESDCorrected(etaR), (high) ? 2 : 1, hists, Form("eff_%d.png", etaR), kFALSE, legendStrings, kFALSE, largestError);
+
+ largestError->Write();
+
+ new TCanvas;
+ largestError->DrawCopy();
+ }
+
+ file->Close();
+}
+
+TCanvas* DrawRatio(TH1* result, Int_t nResultSyst, TH1** resultSyst, TString epsName, Bool_t firstMarker = kFALSE, const char** legendStrings = 0, Bool_t errors = kFALSE, TH1* largestErrorLow = 0, TH1* largestErrorHigh = 0)
{
// compares n results with first results. E.g. one gained with the default response, another with a changed one to study
// a systematic effect
// normalize results
- result->Scale(1.0 / result->Integral(2, 200));
+ //result->Scale(1.0 / result->Integral(1, 200));
TCanvas* canvas = new TCanvas(epsName, epsName, 800, 500);
canvas->SetTopMargin(0.05);
TLegend* legend = new TLegend(0.2, 0.7, 0.7, 0.93);
legend->SetFillColor(0);
- legend->SetTextSize(0.04);
+ //legend->SetTextSize(0.04);
if (nResultSyst > 6)
legend->SetNColumns(2);
for (Int_t n=0; n<nResultSyst; ++n)
{
- resultSyst[n]->Scale(1.0 / resultSyst[n]->Integral(2, 200));
+ //resultSyst[n]->Scale(1.0 / resultSyst[n]->Integral(1, 200));
// calculate ratio
TH1* ratio = (TH1*) result->Clone("ratio");
if (legendStrings && legendStrings[n])
legend->AddEntry(ratio, legendStrings[n], "L");
+
+ /*
+ s = new TSpline3(ratio);
+ s->SetNpx(5);
+ s->SetLineColor(kRed);
+ s->Draw("same");
+ */
+
+ if (largestErrorLow && largestErrorHigh)
+ for (Int_t i=1; i<=ratio->GetNbinsX(); ++i)
+ {
+ if (ratio->GetBinContent(i) - 1 > 0)
+ largestErrorHigh->SetBinContent(i, TMath::Max(ratio->GetBinContent(i) - 1, largestErrorHigh->GetBinContent(i)));
+ if (ratio->GetBinContent(i) - 1 < 0)
+ largestErrorLow->SetBinContent(i, TMath::Min(ratio->GetBinContent(i) - 1, largestErrorLow->GetBinContent(i)));
+ }
+
+ if (largestErrorLow && !largestErrorHigh)
+ for (Int_t i=1; i<=ratio->GetNbinsX(); ++i)
+ largestErrorLow->SetBinContent(i, TMath::Max(TMath::Abs(ratio->GetBinContent(i) - 1), largestErrorLow->GetBinContent(i)));
// get average of ratio
Float_t sum = 0;
line->Draw();
canvas->SaveAs(epsName);
- canvas->SaveAs(Form("%s.gif", epsName.Data()));
return canvas;
}
{
loadlibs();
- AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
- TFile::Open(fileName);
- mult->LoadHistograms("Multiplicity");
+ AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(fileName);
TH1* measured = mult->GetMultiplicityESD(etaRange)->ProjectionY("myesd", 1, 1);
- TH1* unfoldedFolded = mult->CalculateMultiplicityESD(mult->GetMultiplicityESDCorrected(etaRange), etaRange)->ProjectionY("myfolded", 1, 1);
+
+ if (0)
+ {
+ // test how far we are from a normal exponential in the unfolded
+ corrected = mult->GetMultiplicityESDCorrected(etaRange);
+ new TCanvas; corrected->DrawCopy(); gPad->SetLogy();
+
+ expo = new TF1("exp", "expo(0)", 0, 50);
+ //expo->SetParameters(10, -0.18);
+ expo->SetParameters(9.96, -0.176);
+ expo->Draw("SAME");
+
+ for (Int_t i=21; i<=50; i++)
+ corrected->SetBinContent(i, expo->Eval(corrected->GetXaxis()->GetBinCenter(i)));
+
+ corrected->SetLineColor(2);
+ corrected->Draw("SAME");
+ }
+
+ TH1* unfoldedFolded = mult->GetConvoluted(etaRange, AliMultiplicityCorrection::kINEL);
+ //mult->CalculateMultiplicityESD(mult->GetMultiplicityESDCorrected(etaRange), etaRange)->ProjectionY("myfolded", 1, 1);
// normalize
- unfoldedFolded->Scale(1.0 / unfoldedFolded->Integral(2, displayRange+1));
- unfoldedFolded->Scale(measured->Integral(2, displayRange+1));
+ //unfoldedFolded->Scale(1.0 / unfoldedFolded->Integral(2, displayRange+1));
+ //unfoldedFolded->Scale(measured->Integral(2, displayRange+1));
TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
canvas->Range(0, 0, 1, 1);
residual->Add(unfoldedFolded, -1);
// projection
- TH1* residualHist = new TH1F("residualHist", ";", 11, -3, 3);
+ TH1* residualHist = new TH1F("residualHist", ";", 16, -3, 3);
residualHist->Sumw2();
Float_t chi2 = 0;
+ Float_t chi2_hump = 0;
+
for (Int_t i=1; i<=displayRange+1; ++i)
{
if (measured->GetBinError(i) > 0)
residual->SetBinError(i, 1);
residualHist->Fill(residual->GetBinContent(i));
+ if (i >= 15 && i <= 23)
+ chi2_hump += residual->GetBinContent(i) * residual->GetBinContent(i);
+
chi2 += residual->GetBinContent(i) * residual->GetBinContent(i);
}
else
}
Printf("chi2 / ndf = %f / %d = %f", chi2, displayRange+1, chi2 / (displayRange+1));
+ Printf("chi2 (hump) / ndf = %f / %d = %f", chi2_hump, 23-15+1, chi2_hump / (23-15+1));
residual->GetYaxis()->SetTitle("Residuals: (1/e) (M - R #otimes U)");
residual->GetYaxis()->SetRangeUser(-4.5, 4.5);
// SPD TPC
//const char* fileName[] = { "multiplicityMC_400k_syst.root", "multiplicityMC_TPC_4kfiles_syst.root" };
- const char* fileName[] = { "LHC09b9_0.9TeV_0T/mb1/spd/multiplicity.root", "LHC09b8_0.9TeV_0.5T/mb1/tpc/multiplicity.root" };
+ //const char* fileName[] = { "LHC09b9_0.9TeV_0T/mb1/spd/multiplicity.root", "LHC09b8_0.9TeV_0.5T/mb1/tpc/multiplicity.root" };
//const char* fileName[] = { "spd/multiplicity.root", "tpc/multiplicity.root" };
- //const char* fileName[] = { "multiplicity.root", "multiplicity.root" };
+ const char* fileName[] = { "multiplicityparticle-efficiency.root", "multiplicity.root" };
Float_t etaRangeArr[] = {0.49, 0.9};
const char* titles[] = { "SPD Tracklets", "TPC Tracks" };
TLegend* legends[2];
- for (Int_t loop=0; loop<2; ++loop)
+ for (Int_t loop=0; loop<1; ++loop)
{
Printf("%s", fileName[loop]);
TH1* genePt = gene->Project3D(Form("z_%d", i));
TH1* measPt = meas->Project3D(Form("z_%d", i));
- if (i == 2)
- {
- Printf("WARNING: Rebinning for protons!");
- genePt->Rebin(2);
- measPt->Rebin(2);
- }
+// if (i == 2)
+// {
+// Printf("WARNING: Rebinning for protons!");
+// genePt->Rebin(2);
+// measPt->Rebin(2);
+// }
genePt->Sumw2();
measPt->Sumw2();
legend->Draw();
canvas2->SaveAs(Form("%s.eps", canvas2->GetName()));
}
+
+ return;
canvas->cd();
canvas->SaveAs(Form("%s.eps", canvas->GetName()));
TH1* mc = 0;
// loop over cases (normal, enhanced/reduced ratios)
- Int_t nMax = 9;
+ Int_t nMax = 7;
for (Int_t i = 0; i<nMax; ++i)
{
- AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(Form("chi2_species_%d.root", i), Form("Multiplicity_%d", i));
- if (i == 0)
- mc = mult->GetMultiplicityVtx(etaRange)->ProjectionY("mymchist", 1, 1);
+ AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(Form("chi2a_inel_species_%d.root", i), "Multiplicity");
results[i] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i));
}
- DrawResultRatio(mc, results[0], "ParticleSpeciesComparison1_1.eps");
-
for (Int_t i=1; i<=results[0]->GetNbinsX(); i++)
- {
results[0]->SetBinError(i, 0);
- mc->SetBinError(i, 0);
- }
- const char* legendStrings[] = { "K #times 0.5", "K #times 1.5", "p #times 0.5", "p #times 1.5", "K #times 0.5, p #times 0.5", "K #times 1.5, p #times 1.5", "K #times 0.5, p #times 1.5", "K #times 1.5, p #times 0.5" };
+ const char* legendStrings[] = { "K #times 0.7", "K #times 1.3", "p #times 0.7", "p #times 1.3", "others #times 0.7", "others #times 1.3", };
DrawRatio(results[0], nMax-1, results+1, "ParticleSpeciesComparison1_2.eps", kFALSE, legendStrings);
-
- //not valid: draw chi2 uncertainty on top!
- /*TFile::Open("bayesianUncertainty_400k_100k_syst.root");
- TH1* errorHist = (TH1*) gFile->Get("errorBoth");
- errorHist->SetLineColor(1);
- errorHist->SetLineWidth(2);
- TH1* errorHist2 = (TH1*) errorHist->Clone("errorHist2");
- for (Int_t i=1; i<=errorHist->GetNbinsX(); i++)
- {
- errorHist->SetBinContent(i, errorHist->GetBinContent(i) + 1);
- errorHist2->SetBinContent(i, 1 - errorHist2->GetBinContent(i));
- }
-
- errorHist->DrawCopy("SAME");
- errorHist2->DrawCopy("SAME");*/
-
- //canvas->SaveAs(canvas->GetName());
-
- DrawRatio(mc, nMax, results, "ParticleSpeciesComparison1_3.eps", kTRUE, 0);
-
- //errorHist->DrawCopy("SAME");
- //errorHist2->DrawCopy("SAME");
-
- //canvas2->SaveAs(canvas2->GetName());
}
/*void ParticleSpeciesComparison2()
return corr;
}
+void CompareVertexRecoEfficiencies()
+{
+ loadlibs();
+
+ const char* file1 = "LHC10a12_run10482X";
+ const char* file2 = "LHC10a14_run10482X_Phojet";
+
+ const char* suffix = "/all/spd/multiplicityMC_xsection.root";
+
+ hist1 = AliMultiplicityCorrection::Open(Form("%s%s", file1, suffix))->GetEfficiency(etaRange, AliMultiplicityCorrection::kMB);
+ hist2 = AliMultiplicityCorrection::Open(Form("%s%s", file2, suffix))->GetEfficiency(etaRange, AliMultiplicityCorrection::kMB);
+
+ ratio = (TH1*) hist1->Clone("ratio");
+ ratio->Divide(hist2);
+
+ new TCanvas;
+ hist1->Draw();
+ hist2->SetLineColor(2);
+ hist2->Draw("SAME");
+
+ new TCanvas;
+ ratio->Draw();
+}
+
void TriggerVertexCorrection()
{
//
TH1* corrINEL = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kINEL));
TH1* corrNSD = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kNSD));
TH1* corrMB = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kMB));
+
+ TH1* triggerEff = Invert(mult->GetTriggerEfficiency(etaRange, AliMultiplicityCorrection::kNSD));
TCanvas* canvas = new TCanvas("TriggerVertexCorrection", "TriggerVertexCorrection", 800, 500);
gPad->SetGridx();
corrNSD->SetMarkerColor(4);
corrNSD->Draw("SAME PE");
- Printf(" MB INEL NSD");
- Printf("bin 0: %f %f %f", corrMB->GetBinContent(1), corrINEL->GetBinContent(1), corrNSD->GetBinContent(1));
- Printf("bin 1: %f %f %f", corrMB->GetBinContent(2), corrINEL->GetBinContent(2), corrNSD->GetBinContent(2));
+ triggerEff->SetLineColor(6);
+ triggerEff->SetMarkerStyle(25);
+ triggerEff->SetMarkerColor(6);
+ //triggerEff->Draw("SAME PE");
+
+ Printf(" MB INEL NSD TRIGINEL");
+ Printf("bin 0: %.2f %.2f %.2f %.2f", corrMB->GetBinContent(1), corrINEL->GetBinContent(1), corrNSD->GetBinContent(1), triggerEff->GetBinContent(1));
+ Printf("bin 1: %.2f %.2f %.2f %.2f", corrMB->GetBinContent(2), corrINEL->GetBinContent(2), corrNSD->GetBinContent(2), triggerEff->GetBinContent(2));
TLegend* legend = new TLegend(0.3, 0.6, 0.85, 0.85);
legend->SetFillColor(0);
canvas->SaveAs(Form("%s.eps", canvas->GetName()));
}
-void StatisticalUncertainty(Int_t methodType, Bool_t mc = kFALSE)
+void DrawStatisticalUncertainty(const char* fileName = "StatisticalUncertainty.root")
{
- loadlibs();
-
- TFile::Open(correctionFile);
- AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
- mult->LoadHistograms("Multiplicity");
-
- TFile::Open(measuredFile);
- AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
- mult2->LoadHistograms("Multiplicity");
-
- mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
-
- TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, 1);
-
- AliUnfolding::SetNbins(70, 70);
- //AliUnfolding::SetNbins(35, 35);
- AliUnfolding::SetChi2Regularization(AliUnfolding::kPol0, 1e3);
- AliUnfolding::SetBayesianParameters(1, 10);
-
- TH1* errorMeasured = (TH1*) mult->StatisticalUncertainty((AliUnfolding::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kFALSE, ((mc) ? mcHist : 0))->Clone("errorMeasured");
-
- new TCanvas; errorMeasured->Draw();
- new TCanvas;
-
- mult->GetMultiplicityESDCorrected(etaRange)->Scale(1.0 / mult->GetMultiplicityESDCorrected(etaRange)->Integral());
- mult->GetMultiplicityESDCorrected(etaRange)->Draw();
- mcHist->Scale(1.0 / mcHist->Integral());
- mcHist->SetLineColor(2);
- mcHist->Draw("SAME");
- gPad->SetLogy();
-
- return;
-
- TH1* errorResponse = (TH1*) mult->StatisticalUncertainty((AliUnfolding::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorResponse");
-
- TH1* errorBoth = (TH1*) mult->StatisticalUncertainty((AliUnfolding::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorBoth");
-
- if (!mc)
- {
- TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
- DrawResultRatio(mcHist, result, "StatisticalUncertainty2.eps");
- }
-
- TFile* file = new TFile(Form("StatisticalUncertaintySPD%s.root", (methodType == 0) ? "Chi2" : "Bayesian"), "RECREATE");
- errorResponse->Write();
- errorMeasured->Write();
- errorBoth->Write();
- file->Close();
-}
-
-void DrawStatisticalUncertainty()
-{
- TFile::Open("StatisticalUncertainty.root");
+ TFile::Open(fileName);
errorResponse = (TH1*) gFile->Get("errorResponse");
errorMeasured = (TH1*) gFile->Get("errorMeasured");
errorResponse->SetLineColor(1);
errorResponse->GetXaxis()->SetRangeUser(0, longDisplayRange);
- errorResponse->GetYaxis()->SetRangeUser(0, 0.3);
+ errorResponse->GetYaxis()->SetRangeUser(0, 1);
errorResponse->SetStats(kFALSE);
errorResponse->SetTitle(";true multiplicity;Uncertainty");
errorBoth->SetLineColor(4);
errorBoth->Draw("SAME");
+ errorResponse->Scale(1.0 / sqrt(2));
+ errorMeasured->Scale(1.0 / sqrt(2));
+ errorBoth->Scale(1.0 / sqrt(2));
+
Printf("Average errorResponse: %f", errorResponse->Integral(2, displayRange) / (displayRange - 1));
Printf("Average errorMeasured: %f", errorMeasured->Integral(2, displayRange) / (displayRange - 1));
Printf("Average errorBoth: %f", errorBoth->Integral(2, displayRange) / (displayRange - 1));
new TCanvas; esd_proj->DrawCopy();
TH1* percentage = (TH1*) (esd_proj->Clone("percentage"));
+ percentage->SetTitle("percentage");
percentage->Reset();
for (Int_t i=1; i<=esd_proj->GetNbinsX(); i++)
if (esd_proj->GetBinContent(i) > 0)
esd_proj->SetBinContent(i, 1);
+ Int_t limit = -1;
for (Int_t i=1; i<=percentage->GetNbinsX(); i++)
{
TH1* binResponse = corr->ProjectionY("proj", i, i);
binResponse->Scale(1.0 / binResponse->Integral());
binResponse->Multiply(esd_proj);
//new TCanvas; binResponse->Draw();
- percentage->SetBinContent(i, binResponse->Integral());
+ Float_t value = binResponse->Integral();
+ percentage->SetBinContent(i, value);
+ if (limit == -1 && value < 0.9)
+ limit = percentage->GetXaxis()->GetBinCenter(i-1);
//return;
}
+ Printf("Limit is %d", limit);
+
new TCanvas; percentage->Draw();
+
new TCanvas;
mc = esd->GetMultiplicityVtx(etaRange)->ProjectionY("mc", 1, esd->GetMultiplicityVtx(etaRange)->GetNbinsX());
mc->SetLineColor(2);
mc->Draw("");
}
+void FindUnfoldedLimitAll()
+{
+ loadlibs();
+ for (etaRange = 0; etaRange <= 2; etaRange++)
+ FindUnfoldedLimit();
+}
+
void CompareUnfoldedWithMC()
{
loadlibs();
mcHist->Draw("SAME");
gPad->SetLogy();
}
+
+void PaperNumbers()
+{
+ const char* label[] = { "SD", "DD", "ND" };
+ const char* files[] = { "multiplicity_syst_sd.root", "multiplicity_syst_dd.root", "multiplicity_syst_nd.root" };
+
+ loadlibs();
+
+ Printf("vertex reco");
+
+ Float_t oneParticle[3];
+ for (Int_t i=0; i<3; i++)
+ {
+ AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(files[i]);
+
+ eff = mult->GetEfficiency(2, AliMultiplicityCorrection::kMB);
+ //eff = mult->GetTriggerEfficiency(2);
+
+ oneParticle[i] = 100.0 * eff->GetBinContent(2);
+
+ Printf("%s: %.2f", label[i], oneParticle[i]);
+ }
+
+ Float_t ua5_SD = 0.153;
+ Float_t ua5_DD = 0.080;
+ Float_t ua5_ND = 0.767;
+
+ Float_t vtxINELUA5 = ua5_SD * oneParticle[0] + ua5_DD * oneParticle[1] + ua5_ND * oneParticle[2];
+ Float_t vtxNSDUA5 = (ua5_DD * oneParticle[1] + ua5_ND * oneParticle[2]) / (ua5_DD + ua5_ND);
+
+ Printf("INEL (UA5) = %.1f", vtxINELUA5);
+ Printf("NSD (UA5) = %.1f", vtxNSDUA5);
+
+ Printf("total for >= 2 charged tracks");
+
+ AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open("multiplicity.root");
+
+ vtx = mult->GetMultiplicityVtx(2)->ProjectionY("vtx", 1, mult->GetMultiplicityVtx(2)->GetNbinsX(), "e");
+ all = mult->GetMultiplicityINEL(2)->ProjectionY("all", 1, mult->GetMultiplicityINEL(2)->GetNbinsX(), "e");
+
+ Int_t begin = vtx->GetXaxis()->FindBin(2);
+ Int_t end = vtx->GetNbinsX();
+
+ Float_t above2 = vtx->Integral(begin, end) / all->Integral(begin, end);
+
+ Printf("%.1f", 100.0 * above2);
+}
+
+void DrawRawEta()
+{
+ TFile::Open(measuredFile);
+
+ Int_t colors[] = {1, 2, 4, 6, 7, 8, 9, 10};
+
+ for (Int_t i=0; i<3; i++)
+ {
+ hist = dynamic_cast<TH1*> (gFile->Get(Form("fEta_%d", i)));
+
+ hist->GetYaxis()->SetRangeUser(0, 8);
+ hist->SetLineColor(colors[i]);
+ hist->SetStats(kFALSE);
+ hist->Draw((i == 0) ? "HIST" : "HIST SAME");
+ }
+
+ gPad->SetGridx();
+ gPad->SetGridy();
+}
+
+void CrossSectionUncertainties(Int_t xsectionID, Int_t eventType)
+{
+ const char* files[] = { "multiplicitySD.root", "multiplicityDD.root", "multiplicityND.root" };
+
+ loadlibs();
+
+ AliMultiplicityCorrection* data[3];
+
+ Float_t ref_SD = -1;
+ Float_t ref_DD = -1;
+ Float_t ref_ND = -1;
+
+ Float_t error_SD = -1;
+ Float_t error_DD = -1;
+ Float_t error_ND = -1;
+
+ gROOT->ProcessLine(gSystem->ExpandPathName(".L $ALICE_ROOT/PWG0/dNdEta/drawSystematics.C"));
+ GetRelativeFractions(xsectionID, ref_SD, ref_DD, ref_ND, error_SD, error_DD, error_ND);
+
+ TH1* unc[9*3];
+ TH1* unc2[9*3];
+ const char* legendStrings[8];
+
+ for (Int_t i=0; i<9; i++)
+ {
+ Float_t factorSD = 0;
+ Float_t factorDD = 0;
+
+ if (i > 0 && i < 3)
+ factorSD = (i % 2 == 0) ? 1 : -1;
+ else if (i >= 3 && i < 5)
+ factorDD = (i % 2 == 0) ? 1 : -1;
+ else if (i >= 5 && i < 9)
+ {
+ factorSD = ((i % 2 == 0) ? 1.0 : -1.0) / TMath::Sqrt(2);
+ if (i == 5 || i == 6)
+ factorDD = factorSD;
+ else
+ factorDD = -factorSD;
+ }
+
+ Float_t scalesSD = ref_SD + factorSD * error_SD;
+ Float_t scalesDD = ref_DD + factorDD * error_DD;
+ Float_t scalesND = 1.0 - scalesDD[i] - scalesSD[i];
+
+ str = new TString;
+ str->Form("Case %d: SD: %.2f DD: %.2f ND: %.2f", i, scalesSD, scalesDD, scalesND);
+ Printf("%s", str->Data());
+ if (i > 0)
+ legendStrings[i-1] = str->Data();
+
+ for (Int_t j=0; j<3; ++j)
+ {
+ TString name;
+ name.Form("Multiplicity_%d", j);
+
+ TFile::Open(files[j]);
+ data[j] = new AliMultiplicityCorrection(name, name);
+ data[j]->LoadHistograms("Multiplicity");
+ }
+
+ // calculating relative
+ Float_t sd = data[0]->GetMultiplicityINEL(0)->Integral(0, data[0]->GetMultiplicityINEL(0)->GetNbinsX()+1);
+ Float_t dd = data[1]->GetMultiplicityINEL(0)->Integral(0, data[1]->GetMultiplicityINEL(0)->GetNbinsX()+1);
+ Float_t nd = data[2]->GetMultiplicityINEL(0)->Integral(0, data[2]->GetMultiplicityINEL(0)->GetNbinsX()+1);
+ Float_t total = nd + dd + sd;
+
+ nd /= total;
+ sd /= total;
+ dd /= total;
+
+ Printf("Ratios in the correction map are: ND=%f, DD=%f, SD=%f", nd, dd, sd);
+
+ Float_t ratio[3];
+ ratio[0] = scalesSD / sd;
+ ratio[1] = scalesDD / dd;
+ ratio[2] = scalesND / nd;
+
+ Printf("SD=%.2f, DD=%.2f, ND=%.2f", ratio[0], ratio[1], ratio[2]);
+
+ TList list;
+ for (Int_t k=0; k<3; ++k)
+ {
+ // modify x-section
+ for (Int_t j=0; j<AliMultiplicityCorrection::kMCHists; j++)
+ {
+ data[k]->GetMultiplicityVtx(j)->Scale(ratio[k]);
+ data[k]->GetMultiplicityMB(j)->Scale(ratio[k]);
+ data[k]->GetMultiplicityINEL(j)->Scale(ratio[k]);
+ data[k]->GetMultiplicityNSD(j)->Scale(ratio[k]);
+ }
+
+ if (k > 0)
+ list.Add(data[k]);
+ }
+
+ printf("Case %d, %s: Entries in response matrix: SD: %.2f DD: %.2f ND: %.2f", xsectionID, data[0]->GetName(), data[0]->GetCorrelation(0)->Integral(), data[1]->GetCorrelation(0)->Integral(), data[2]->GetCorrelation(0)->Integral());
+
+ data[0]->Merge(&list);
+ data[0]->FixTriggerEfficiencies(20);
+
+ Printf(" Total: %.2f", data[0]->GetCorrelation(0)->Integral());
+
+ for (Int_t etaR = 0; etaR < 3; etaR++)
+ {
+ unc[i+(etaR*9)] = (TH1*) data[0]->GetTriggerEfficiency(etaR, eventType)->Clone(Form("eff_%d_%d", etaR, i));
+ unc2[i+(etaR*9)] = (TH1*) data[0]->GetEfficiency(etaR, eventType)->Clone(Form("eff_%d_%d", etaR, i));
+ }
+ }
+
+ file = TFile::Open(Form("systunc_fractions_%s.root", (eventType == 2) ? "inel" : "nsd"), "RECREATE");
+
+ for (Int_t etaR = 0; etaR < 3; etaR++)
+ {
+ largestError1 = (TH1*) unc[0]->Clone(Form("fractions_trig_%d", etaR));
+ largestError1->Reset();
+ largestError2 = (TH1*) unc[0]->Clone(Form("fractions_trigvtx_%d", etaR));
+ largestError2->Reset();
+
+ etaRange = etaR; // correct labels in DrawRatio
+ DrawRatio(unc[etaR*9], 8, unc+1+etaR*9, Form("xsections_trig%d.png", etaR), kFALSE, legendStrings, kFALSE, largestError1);
+ DrawRatio(unc2[etaR*9], 8, unc2+1+etaR*9, Form("xsections_trigvtx%d.png", etaR), kFALSE, legendStrings, kFALSE, largestError2);
+
+ largestError2->SetBinContent(1, largestError1->GetBinContent(1));
+
+ largestError2->GetYaxis()->UnZoom();
+ largestError2->Write();
+
+ new TCanvas;
+ largestError2->DrawCopy();
+ }
+ file->Close();
+}
+
+void PythiaPhojetUncertainties()
+{
+ PythiaPhojetUncertainties(900, 2);
+ PythiaPhojetUncertainties(900, 3);
+ PythiaPhojetUncertainties(2360, 2);
+ PythiaPhojetUncertainties(2360, 3);
+}
+
+void PythiaPhojetUncertainties(Int_t energy, Int_t eventType)
+{
+ loadlibs();
+
+ AliMultiplicityCorrection* mult[2];
+
+ const char* trigger = "all";
+ if (eventType == 3)
+ trigger = "v0and";
+
+ if (energy == 2360)
+ trigger = "spdgfobits";
+
+ if (energy == 7000)
+ trigger = "all-onepart";
+
+ if (energy == 900)
+ TFile::Open(Form("LHC10a12_run10482X/%s/spd/multiplicityMC_xsection.root", trigger));
+ else if (energy == 2360)
+ TFile::Open(Form("LHC10a10_run105054_7/%s/spd/multiplicityMC_xsection.root", trigger));
+ else
+ TFile::Open(Form("LHC10b2_run114931/%s/spd/multiplicity.root", trigger));
+ mult[0] = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+ mult[0]->LoadHistograms("Multiplicity");
+
+ if (energy == 900)
+ TFile::Open(Form("LHC10a14_run10482X_Phojet/%s/spd/multiplicityMC_xsection.root", trigger));
+ else if (energy == 2360)
+ TFile::Open(Form("LHC10a15_run10505X_Phojet/%s/spd/multiplicityMC_xsection.root", trigger));
+ else
+ TFile::Open(Form("LHC10b1_run114931/%s/spd/multiplicity.root", trigger));
+ mult[1] = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+ mult[1]->LoadHistograms("Multiplicity");
+
+ TH1* unc[6];
+ TH1* unc2[6];
+
+ for (Int_t i=0; i<2; i++)
+ {
+ mult[i]->FixTriggerEfficiencies(20);
+ for (Int_t etaR = 0; etaR < 3; etaR++)
+ {
+ unc[i+(etaR*2)] = (TH1*) mult[i]->GetTriggerEfficiency(etaR, eventType)->Clone(Form("eff_%d_%d", etaR, i));
+ unc2[i+(etaR*2)] = (TH1*) mult[i]->GetEfficiency(etaR, eventType)->Clone(Form("eff_%d_%d", etaR, i));
+ }
+ }
+
+ file = TFile::Open(Form("systunc_pythiaphojet_%d_%s.root", energy, (eventType == 2) ? "inel" : "nsd"), "RECREATE");
+
+ for (Int_t etaR = 0; etaR < 3; etaR++)
+ {
+ largestError1Low = (TH1*) unc[0]->Clone(Form("pythiaphojet_trig_%d_low", etaR));
+ largestError1Low->Reset();
+ largestError1High = (TH1*) largestError1Low->Clone(Form("pythiaphojet_trig_%d_high", etaR));
+ largestError2Low = (TH1*) largestError1Low->Clone(Form("pythiaphojet_trigvtx_%d_low", etaR));
+ largestError2High = (TH1*) largestError1Low->Clone(Form("pythiaphojet_trigvtx_%d_high", etaR));
+
+ DrawRatio(unc[etaR*2], 1, unc+1+etaR*2, Form("pythiaphojet_trig%d.png", etaR), kFALSE, 0, kFALSE, largestError1Low, largestError1High);
+ DrawRatio(unc2[etaR*2], 1, unc2+1+etaR*2, Form("pythiaphojet_trigvtx%d.png", etaR), kFALSE, 0, kFALSE, largestError2Low, largestError2High);
+
+ largestError2Low->SetBinContent(1, largestError1Low->GetBinContent(1));
+ largestError2High->SetBinContent(1, largestError1High->GetBinContent(1));
+
+ largestError2Low->GetYaxis()->UnZoom();
+ largestError2High->GetYaxis()->UnZoom();
+ largestError2Low->Write();
+ largestError2High->Write();
+
+ new TCanvas;
+ largestError2Low->DrawCopy();
+
+ new TCanvas;
+ largestError2High->DrawCopy();
+ }
+ file->Close();
+}
+
+void CombinatoricsAsFunctionOfN(Int_t etaR)
+{
+ loadlibs();
+
+ mult = AliMultiplicityCorrection::Open("multiplicityESD.root");
+
+ //Float_t integratedCombinatorics[] = { 0.00062, 0.00074, 0.001 }; // 900 GeV
+ Float_t integratedCombinatorics[] = { 7.5e-4, 1.1e-3, 1.2e-3 }; // 2.36 TeV
+
+ multHist = mult->GetMultiplicityESD(etaR)->ProjectionY();
+ multHist->Scale(1.0 / multHist->Integral());
+
+ Float_t integral = 0;
+ Float_t tracks = 0;
+ for (Int_t i=1; i<= multHist->GetNbinsX(); i++)
+ {
+ Float_t x = multHist->GetXaxis()->GetBinCenter(i);
+ integral += x*x * multHist->GetBinContent(i+1);
+ tracks += x * multHist->GetBinContent(i+1);
+ Printf("%f %f %f", x, multHist->GetBinContent(i+1), integral);
+ }
+
+ Printf("Integral is %f; %f", integral, integral / tracks);
+ Printf("alpha is %e", integratedCombinatorics[etaR] / (integral / tracks));
+}
+
+void TryCombinatoricsCorrection(Float_t etaR, Float_t alpha)
+{
+ loadlibs();
+
+ mult = AliMultiplicityCorrection::Open("multiplicity.root");
+
+ multHist = mult->GetMultiplicityESD(etaR)->ProjectionY();
+
+ target = (TH1*) multHist->Clone("target");
+ target->Reset();
+
+ Int_t totalTracks1 = 0;
+ Int_t totalTracks2 = 0;
+
+ for (Int_t i=1; i<=multHist->GetNbinsX(); i++)
+ {
+ for (Int_t j=0; j<multHist->GetBinContent(i); j++)
+ {
+ Int_t origTracks = (Int_t) multHist->GetXaxis()->GetBinCenter(i);
+ tracks = origTracks;
+ //tracks -= gRandom->Poisson(1.3e-4 * origTracks * origTracks);
+ if (gRandom->Uniform() < alpha * origTracks * origTracks)
+ tracks--;
+ target->Fill(tracks);
+
+ totalTracks1 += origTracks;
+ totalTracks2 += tracks;
+ }
+ }
+
+ new TCanvas; multHist->Draw(); target->Draw("SAME"); target->SetLineColor(2);
+
+ Printf("%d %d %f %f", totalTracks1, totalTracks2, 1. - (Float_t) totalTracks2 / totalTracks1, 1. - target->GetMean() / multHist->GetMean());
+}
+
+void ApplyCombinatoricsCorrection()
+{
+ loadlibs();
+
+ mult = AliMultiplicityCorrection::Open("multiplicity.root");
+
+ fractionMC = new TF1("fractionMC", "[0]+[1]*x", 0, 200);
+ fractionMC->SetParameters(0.0004832, 5.82e-5); // mariella's presentation 16.04.10
+
+ fractionData = new TF1("fractionData", "[0]+[1]*x", 0, 200);
+ fractionData->SetParameters(0.00052, 0.0001241); // mariella's presentation 16.04.10
+
+ Int_t etaR = 1;
+ esd = mult->GetMultiplicityESD(etaR);
+
+ for (Int_t x=1; x<=esd->GetNbinsX(); x++)
+ {
+ for (Int_t y=2; y<=esd->GetNbinsY(); y++)
+ {
+ printf("Before: %f %f ", esd->GetBinContent(x, y-1), esd->GetBinContent(x, y));
+
+ Float_t n = esd->GetYaxis()->GetBinCenter(y);
+ Int_t addComb = (fractionData->Eval(n) - fractionMC->Eval(n)) * n * esd->GetBinContent(x, y);
+
+ // shift addComb down one bin
+ esd->SetBinContent(x, y-1, esd->GetBinContent(x, y-1) + addComb);
+ esd->SetBinContent(x, y, esd->GetBinContent(x, y) - addComb);
+
+ // recalc errors
+ esd->SetBinError(x, y-1, TMath::Sqrt(esd->GetBinContent(x, y-1)));
+ esd->SetBinError(x, y, TMath::Sqrt(esd->GetBinContent(x, y)));
+
+ Printf("comb: %d; after: %f +- %f %f +- %f", addComb, esd->GetBinContent(x, y-1), esd->GetBinError(x, y-1), esd->GetBinContent(x, y), esd->GetBinError(x, y));
+ }
+ }
+
+ TFile::Open("out.root", "RECREATE");
+ mult->SaveHistograms();
+}
-void run(Char_t* data, Long64_t nRuns = -1, Long64_t offset = 0, Bool_t aDebug = kFALSE, Int_t aProof = 0, Bool_t mc = kTRUE, const char* option = "")
+void run(Char_t* data, Long64_t nRuns = -1, Long64_t offset = 0, Bool_t aDebug = kFALSE, Int_t aProof = 0, Int_t requiredData = 1, const char* option = "")
{
// aProof option: 0 no proof
// 1 proof with chain
// 2 proof with dataset
+ //
+ // requiredData option: 0 = only ESD
+ // 1 = ESD+MC
+ // 2 = RAW (ESD+check on event type)
+ //
if (nRuns < 0)
nRuns = 1234567890;
{
gEnv->SetValue("XSec.GSI.DelegProxy", "2");
TProof::Open("alicecaf");
- //gProof->SetParallel(1);
-
+
// Enable the needed package
if (1)
{
// Create the analysis manager
mgr = new AliAnalysisManager;
- AliPWG0Helper::AnalysisMode analysisMode = AliPWG0Helper::kTPC | AliPWG0Helper::kFieldOn;
- AliTriggerAnalysis::Trigger trigger = AliTriggerAnalysis::kMB1;
+ // Add ESD handler
+ AliESDInputHandler* esdH = new AliESDInputHandler;
+ esdH->SetInactiveBranches("AliESDACORDE FMD ALIESDTZERO ALIESDZDC AliRawDataErrorLogs CaloClusters Cascades EMCALCells EMCALTrigger ESDfriend Kinks AliESDTZERO ALIESDACORDE MuonTracks TrdTracks");
+ mgr->SetInputEventHandler(esdH);
+
+ // physics selection
+ gROOT->ProcessLine(".L $ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
+ physicsSelectionTask = AddTaskPhysicsSelection((requiredData == 2) ? kFALSE : kTRUE);
+
+ // FO efficiency (for MC)
+ if (1 && requiredData != 2)
+ {
+ //const char* fastORFile = "../dNdEta/spdFOEff_run104824_52.root";
+ //const char* fastORFile = "../dNdEta/spdFOEff_run104867_92.root";
+ //const char* fastORFile = "../dNdEta/spdFOEff_run105054_7.root";
+ const char* fastORFile = "../dNdEta/spdFOEff_run114931.root";
+
+ Printf("NOTE: Simulating FAST-OR efficiency on the analysis level using file %s", fastORFile);
+ TFile::Open(fastORFile);
+
+ spdFOEff = (TH1F*) gFile->Get("spdFOEff");
+ physicsSelectionTask->GetPhysicsSelection()->Initialize(114931);
+ physicsSelectionTask->GetPhysicsSelection()->GetTriggerAnalysis()->SetSPDGFOEfficiency(spdFOEff);
+ }
+
+ AliPWG0Helper::AnalysisMode analysisMode = AliPWG0Helper::kSPD | AliPWG0Helper::kFieldOn;
+ //AliPWG0Helper::AnalysisMode analysisMode = AliPWG0Helper::kTPCITS | AliPWG0Helper::kFieldOn;
+
+ //AliTriggerAnalysis::Trigger trigger = AliTriggerAnalysis::kAcceptAll | AliTriggerAnalysis::kOfflineFlag;
+ AliTriggerAnalysis::Trigger trigger = AliTriggerAnalysis::kAcceptAll | AliTriggerAnalysis::kOfflineFlag | AliTriggerAnalysis::kOneParticle;
+
+ //AliTriggerAnalysis::Trigger trigger = AliTriggerAnalysis::kMB1Prime | AliTriggerAnalysis::kOfflineFlag;
+ //AliTriggerAnalysis::Trigger trigger = AliTriggerAnalysis::kSPDGFOBits | AliTriggerAnalysis::kOfflineFlag;
+ //AliTriggerAnalysis::Trigger trigger = AliTriggerAnalysis::kV0AND | AliTriggerAnalysis::kOfflineFlag;
- AliPWG0Helper::PrintConf(analysisMode, trigger);
+ AliPWG0Helper::DiffTreatment diffTreatment = AliPWG0Helper::kMCFlags;
+ //AliPWG0Helper::DiffTreatment diffTreatment = AliPWG0Helper::kE710Cuts;
+
+ AliPWG0Helper::PrintConf(analysisMode, trigger, diffTreatment);
TString taskName("AliMultiplicityTask.cxx+");
if (aDebug)
} else
gROOT->Macro(taskName);
- task = new AliMultiplicityTask(option);
+ // 0 bin calculation
+ if (0)
+ {
+ }
+
+ // V0 syst. study
+ if (0)
+ {
+ Printf("NOTE: Systematic study for VZERO enabled!");
+ //physicsSelectionTask->GetPhysicsSelection()->Initialize(104867);
+ for (Int_t i=0; i<1; i++)
+ {
+ // for MC and data
+ physicsSelectionTask->GetPhysicsSelection()->GetTriggerAnalysis(i)->SetV0HwPars(15, 61.5, 86.5);
+ physicsSelectionTask->GetPhysicsSelection()->GetTriggerAnalysis(i)->SetV0AdcThr(15);
+ // only for MC
+ //physicsSelectionTask->GetPhysicsSelection()->GetTriggerAnalysis(i)->SetV0HwPars(0, 0, 125);
+ //physicsSelectionTask->GetPhysicsSelection()->GetTriggerAnalysis(i)->SetV0AdcThr(0);
+ }
+ }
+
+ TString optionStr(option);
+
+ // remove SAVE option if set
+ Bool_t save = kFALSE;
+ TString optionStr(option);
+ if (optionStr.Contains("SAVE"))
+ {
+ optionStr = optionStr(0,optionStr.Index("SAVE")) + optionStr(optionStr.Index("SAVE")+4, optionStr.Length());
+ save = kTRUE;
+ }
+
+ task = new AliMultiplicityTask(optionStr);
if (!(analysisMode & AliPWG0Helper::kSPD))
{
task->SetTrackCuts(esdTrackCuts);
}
- else
- task->SetDeltaPhiCut(0.05);
+ //else
+ // task->SetDeltaPhiCut(0.05);
task->SetAnalysisMode(analysisMode);
task->SetTrigger(trigger);
+ task->SetDiffTreatment(diffTreatment);
- if (mc)
+ if (requiredData == 1)
task->SetReadMC();
//task->SetUseMCVertex();
+
+ if (requiredData != 2)
+ task->SetSkipParticles();
mgr->AddTask(task);
- TString optionStr(option);
-
- if (mc) {
+ if (requiredData == 1) {
// Enable MC event handler
AliMCEventHandler* handler = new AliMCEventHandler;
if (!optionStr.Contains("particle-efficiency"))
task->SetPtSpectrum((TH1D*) hist->Clone("pt-spectrum"));
}
- // Add ESD handler
- AliESDInputHandler* esdH = new AliESDInputHandler;
- esdH->SetInactiveBranches("AliESDACORDE FMD ALIESDTZERO ALIESDZDC AliRawDataErrorLogs CaloClusters Cascades EMCALCells EMCALTrigger ESDfriend Kinks AliESDTZERO ALIESDACORDE MuonTracks TrdTracks");
- mgr->SetInputEventHandler(esdH);
-
// Attach input
cInput = mgr->GetCommonInputContainer();
mgr->ConnectInput(task, 0, cInput);
// Attach output
cOutput = mgr->CreateContainer("cOutput", TList::Class(), AliAnalysisManager::kOutputContainer);
- //cOutput->SetDataOwned(kTRUE);
mgr->ConnectOutput(task, 0, cOutput);
// Enable debug printouts
// process dataset
mgr->StartAnalysis("proof", data, nRuns, offset);
+
+ if (save)
+ {
+ TString path("maps/");
+ path += TString(data).Tokenize("/")->Last()->GetName();
+
+ UInt_t triggerNoFlags = (UInt_t) trigger % (UInt_t) AliTriggerAnalysis::kStartOfFlags;
+ switch (triggerNoFlags)
+ {
+ case AliTriggerAnalysis::kAcceptAll: path += "/all"; break;
+ case AliTriggerAnalysis::kMB1: path += "/mb1"; break;
+ case AliTriggerAnalysis::kMB2: path += "/mb2"; break;
+ case AliTriggerAnalysis::kMB3: path += "/mb3"; break;
+ case AliTriggerAnalysis::kSPDGFO: path += "/spdgfo"; break;
+ case AliTriggerAnalysis::kSPDGFOBits: path += "/spdgfobits"; break;
+ case AliTriggerAnalysis::kV0AND: path += "/v0and"; break;
+ case AliTriggerAnalysis::kNSD1: path += "/nsd1"; break;
+ case AliTriggerAnalysis::kMB1Prime: path += "/mb1prime"; break;
+ default: Printf("ERROR: Trigger undefined for path to files"); return;
+ }
+
+ if (trigger & AliTriggerAnalysis::kOneParticle)
+ path += "-onepart";
+
+ if (analysisMode & AliPWG0Helper::kSPD)
+ path += "/spd";
+
+ if (analysisMode & AliPWG0Helper::kTPC)
+ path += "/tpc";
+
+ if (analysisMode & AliPWG0Helper::kTPCITS)
+ path += "/tpcits";
+
+ gSystem->mkdir(path, kTRUE);
+
+ TString fileName("multiplicity");
+ if (optionStr.Contains("only-process-type-nd"))
+ fileName += "ND";
+ if (optionStr.Contains("only-process-type-sd"))
+ fileName += "SD";
+ if (optionStr.Contains("only-process-type-dd"))
+ fileName += "DD";
+ fileName += ".root";
+
+ gSystem->Rename(fileName, path + "/" + fileName);
+ gSystem->Rename("event_stat.root", path + "/event_stat.root");
+
+ Printf(">>>>> Moved files to %s", path.Data());
+ }
}
else if (aProof == 3)
{
gROOT->ProcessLine(".L CreateChainFromDataSet.C");
ds = gProof->GetDataSet(data)->GetStagedSubset();
- chain = CreateChainFromDataSet(ds);
+ chain = CreateChainFromDataSet(ds, "esdTree", nRuns);
mgr->StartAnalysis("local", chain, nRuns, offset);
}
else
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff