/************************************************************************************************************* This macro analyses the production of SingleMuon and J/Psi as a function of the collision multiplicity in the SPD. It reads and analyses the output of the Analysis Task PWG3/muon/AliAnalysisTaskMuonCollisionMultiplicity. Refer to the corresponding files to know what the output of this task is. Thismacro use a number of input files whose names are hard-coded These are : pTRanges.txt : different ranges in pT in which the study is done. etaRnages.txt : different ranges in eta in which the study is done. *************************************************************************************************************/ // ROOT includes #include #include #include #include #include #include #include #include #include #include // RooFit includes #include #include #include #include #include #include #include #include #include #include // std includes #include // These functions are in charge of doing all the analysis void ProduceTriggerGraph(TFile *inputFile, TFile *outputFile, std::vector multiplicityRanges, Bool_t applyZCut, Bool_t applyPileUpCut); void ProduceSingleMuonRawGraph(TFile *inputFile, TFile *outputFile, std::vector multiplicityRanges, Bool_t applyZCut, Bool_t applyPileUpCut, Bool_t applyMatchTrigCut, Bool_t applyThetaAbsCut, Bool_t applyPDCACut); void AnalysisSingleMuon(TFile *outputFile); void SingleMuonYieldGraph(TFile *outputFile, TGraphErrors *CINT1B); void SingleMuonYieldOverMeanMultGraph(TFile *outputFile, TGraphErrors *CINT1B); void SingleMuonYieldNormalisedGraph(TFile *outputFile, TGraphErrors *CMUS1B); void ProduceFitResults(TFile *inputFile, TFile *outputFile, std::vector multiplicityRanges, Bool_t applyZCut, Bool_t applyPileUpCut, Double_t numberMatchTrig, Bool_t applyThetaAbsCut, Bool_t applyPDCACut); void FitInvariantMass(TList *fitList, TH1D *invariantMass, Bool_t areParametersFixed, TH1D *referenceParameters); void ProduceDimuonGraph(TFile *outputFile, std::vector multiplicityRanges); /////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////// void AnalysisFunctionOfMultiplicity(Bool_t doSingleMuonAnalysis = kTRUE, Bool_t doJPsiAnalysis = kTRUE, TString inputName = "./Input/MUON.Multiplicity.THnSparse.ESDs.root", TString multiplicityFileName = "multiplicityJPsi.txt", TString outputName = "./Output/MultiplicityResults.LHC10e.JPsi.ESDs.root") { // The main function of the macro // Open the input file and create the output file TFile *inputFile = TFile::Open(inputName, "read"); TFile *outputFile = TFile::Open(outputName, "recreate"); // First of all, we need to figure which conditions we are using // For the event : the cut on z_vertex and if we have pile up from the SPD // For the muons : the usual cuts for SMR Bool_t applyZCut = kTRUE; Bool_t applyPileUpCut = kTRUE; Bool_t applyMatchTrigCut = kTRUE; Bool_t applyThetaAbsCut = kTRUE; Bool_t applyPDCACut = kTRUE; Double_t nMatchTrig = 1.0; // Now, we need the multiplicity ranges // We get it from an input file ifstream multiplicityFile(multiplicityFileName); Double_t multiplicityRange = 0.0; std::vector multiplicityRanges; while(multiplicityFile >> multiplicityRange) multiplicityRanges.push_back(multiplicityRange); // With that, we can produce the number of CINT1B in each bin // Specifficaly, we'll save two TGrpahErrors in the output : // - the graph containing the number of minimum bias event in the bin // - the graph containing the mean multiplicity in the bin ProduceTriggerGraph(inputFile, outputFile, multiplicityRanges, applyZCut, applyPileUpCut); // Now, the muon analysis if (doSingleMuonAnalysis) { // First, we produce the single muon raw graph // it is the graph giving the number of muons detected in different the multiplicity bins, and in different eta and pT domains ProduceSingleMuonRawGraph(inputFile, outputFile, multiplicityRanges, applyZCut, applyPileUpCut, applyMatchTrigCut, applyThetaAbsCut, applyPDCACut); // Raw graph are created, it is time to analyse them AnalysisSingleMuon(outputFile); } // Then, the J/Psi analysis if (doJPsiAnalysis) { ProduceFitResults(inputFile, outputFile, multiplicityRanges, applyZCut, applyPileUpCut, nMatchTrig, applyThetaAbsCut, applyPDCACut); ProduceDimuonGraph(outputFile, multiplicityRanges); } inputFile->Close(); outputFile->Close(); cout << "End" << endl; } /////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////// void ProduceTriggerGraph(TFile *inputFile, TFile *outputFile, std::vector multiplicityRanges, Bool_t applyZCut, Bool_t applyPileUpCut) { // Compute the number of CINT1B and CMUS1B in each bin // Along the x axis, the points are placed at the mean multiplicity of the bin, with an error equal to the error on the mean multiplicity // First of all, we need to get the number of CINT1B as a function of the multiplicity in a TH1D // Retrieving the THnSparse in the input THnSparseD *inputCINT1B = (THnSparseD *) ((TList *) (inputFile->Get("Trigger;1"))->FindObject("CINT1B") ); THnSparseD *inputCMUS1B = (THnSparseD *) ((TList *) (inputFile->Get("Trigger;1"))->FindObject("CMUS1B") ); // Reminder of the architecture of this THnSparse // dimension 0 : multiplicity of the event // dimension 1 : is the vertex in the z range (0 for no, 1 for yes)? // dimension 2 : is it an event without pile up (0 for no, 1 for yes)? // Now, applying the cuts on z vertex and pile-up if (applyZCut) { inputCINT1B->GetAxis(1)->SetRangeUser(1.0, 2.0); inputCMUS1B->GetAxis(1)->SetRangeUser(1.0, 2.0); } if (applyPileUpCut) { inputCINT1B->GetAxis(2)->SetRangeUser(1.0, 2.0); inputCMUS1B->GetAxis(2)->SetRangeUser(1.0, 2.0); } // cuts applied, we can project the THnSparse on TH1D, because it is easier to handle TH1D *histoCINT1B = inputCINT1B->Projection(0, "E"); TH1D *histoCMUS1B = inputCMUS1B->Projection(0, "E"); // We can now create and fill the two TGraphErrors TGraphErrors *graphCINT1B = new TGraphErrors(multiplicityRanges.size() - 1); graphCINT1B->SetName("graphCINT1B"); TGraphErrors *graphCMUS1B = new TGraphErrors(multiplicityRanges.size() - 1); graphCMUS1B->SetName("graphCMUS1B"); // Graphs for plotting purposes TGraphAsymmErrors *plotCINT1B = new TGraphAsymmErrors(multiplicityRanges.size() - 1); plotCINT1B->SetName("plotCINT1B"); TGraphAsymmErrors *plotCMUS1B = new TGraphAsymmErrors(multiplicityRanges.size() - 1); plotCMUS1B->SetName("plotCMUS1B"); for (UInt_t iRange = 0; iRange < multiplicityRanges.size() - 1; iRange++) { // CINT1B first Int_t firstBin = histoCINT1B->GetXaxis()->FindBin(multiplicityRanges[iRange]); Int_t lastBin = histoCINT1B->GetXaxis()->FindBin(multiplicityRanges[iRange+1]); histoCINT1B->GetXaxis()->SetRange(firstBin, lastBin); Double_t total = 0.0; Double_t totalErr = 0.0; total = histoCINT1B->IntegralAndError(firstBin, lastBin, totalErr); Double_t mean = histoCINT1B->GetMean(); Double_t meanErr = histoCINT1B->GetMeanError(); graphCINT1B->SetPoint(iRange, mean, total); graphCINT1B->SetPointError(iRange, meanErr, totalErr); plotCINT1B->SetPoint(iRange, mean, total); plotCINT1B->SetPointError(iRange, mean - multiplicityRanges[iRange], multiplicityRanges[iRange+1] - mean, totalErr, totalErr); // CMUS1B second firstBin = histoCMUS1B->GetXaxis()->FindBin(multiplicityRanges[iRange]); lastBin = histoCMUS1B->GetXaxis()->FindBin(multiplicityRanges[iRange+1]); histoCMUS1B->GetXaxis()->SetRange(firstBin, lastBin); total = histoCMUS1B->IntegralAndError(firstBin, lastBin, totalErr); mean = histoCMUS1B->GetMean(); meanErr = histoCMUS1B->GetMeanError(); graphCMUS1B->SetPoint(iRange, mean, total); graphCMUS1B->SetPointError(iRange, meanErr, totalErr); plotCMUS1B->SetPoint(iRange, mean, total); plotCMUS1B->SetPointError(iRange, mean - multiplicityRanges[iRange], multiplicityRanges[iRange+1] - mean, totalErr, totalErr); } // Saving the graphs outputFile->WriteTObject(plotCINT1B); outputFile->WriteTObject(plotCMUS1B); outputFile->WriteTObject(graphCINT1B); outputFile->WriteTObject(graphCMUS1B); delete graphCINT1B; delete histoCINT1B; delete inputCINT1B; delete graphCMUS1B; delete histoCMUS1B; delete inputCMUS1B; } void ProduceSingleMuonRawGraph(TFile *inputFile, TFile *outputFile, std::vector multiplicityRanges, Bool_t applyZCut, Bool_t applyPileUpCut, Bool_t applyMatchTrigCut, Bool_t applyThetaAbsCut, Bool_t applyPDCACut) { // Produce the raw graph for Single Muon that will be used all along the analysis // There are several graph : // - over all the range in both eta and pT // - Single Muons Reference : pT > 1 GeV // - Heavy Flavor Muons : pT > 4 GeV // - bin by bin in pT and over all the range in eta // - bin by bin in eta and with single muon reference // and right now, there is not enough stat to do it bin by bin in both pT and eta // Retrieve the THnSparse and the Minimum bias histo THnSparseD *inputSingleMuon = (THnSparseD *) ((TList *) (inputFile->Get("SingleMuon;1"))->FindObject("muonCMUS1B") ); TGraphErrors *CINT1B = (TGraphErrors *) (outputFile->Get("graphCINT1B;1") ); // Reminder of the architecture of this THnSparse // dimension 0 : multiplicity of the event // dimension 1 : is the vertex in the z range (0 for no, 1 for yes)? // dimension 2 : is it an event without pile up (0 for no, 1 for yes)? // dimension 3 : does the muon match the trigger (0 for no, 1 for yes)? // dimension 4 : theta_abs of the muon // dimension 5 : eta of the muon // dimension 6 : p DCA_x of the muon // dimension 7 : p DCA_y of the muon // dimension 8 : p DCA of the muon // dimension 9 : p of the muon // dimension 10 : pT of the muon // get the pT and eta ranges from files. // Beware, the names of the files are hard-coded, so make sure to have them in your folder. // I might change this in the future, but I felt the main function had enough parameters already. ifstream pTFile ("pTRanges.txt"); Double_t pTRange = 0.0; std::vector pTRanges; while(pTFile >> pTRange) pTRanges.push_back(pTRange); ifstream etaFile ("etaRanges.txt"); Double_t etaRange = 0.0; std::vector etaRanges; while(etaFile >> etaRange) etaRanges.push_back(etaRange); // Apply all the cuts // Beware, theta_abs and pDCA cuts are hard-coded if (applyZCut) inputSingleMuon->GetAxis(1)->SetRangeUser(1.0, 2.0); if (applyPileUpCut) inputSingleMuon->GetAxis(2)->SetRangeUser(1.0, 2.0); if (applyMatchTrigCut) inputSingleMuon->GetAxis(3)->SetRangeUser(1.0, 2.0); if (applyThetaAbsCut) inputSingleMuon->GetAxis(4)->SetRangeUser(2.0, 9.0); if (applyPDCACut) inputSingleMuon->GetAxis(8)->SetRangeUser(0.0, 450.0); // no cut yet, I need to check what are the usual values // First, we get all the raw histos TList *rawSingleMuon = new TList(); rawSingleMuon->SetName("rawSingleMuon"); // All integrated // There are some hard cuts // - pT : 0.5 -> 8 GeV, because we are not confident on what we are seeing outside of these limits // - eta : -4.0, -> -2.5, acceptance of the spectrometer Double_t etaMinAbsolute = -4.0; Double_t etaMaxAbsolute = -2.5; inputSingleMuon->GetAxis(5)->SetRangeUser(etaMinAbsolute, etaMaxAbsolute); Double_t pTMinAbsolute = 0.5; Double_t pTMaxAbsolute = 8.0; inputSingleMuon->GetAxis(10)->SetRangeUser(pTMinAbsolute, pTMaxAbsolute); // Putting this in a TH3D, since it is easier to use and allow for computation of errors on the integral // TH3D is multiplicity, eta and pT TH3D *histoSingleMuon = inputSingleMuon->Projection(0, 5, 10, "E"); // Declare all the histos TGraphAsymmErrors *allSingleMuon = new TGraphAsymmErrors(multiplicityRanges.size() - 1); allSingleMuon->SetName("allSingleMuon"); TGraphAsymmErrors *rawSMR = new TGraphAsymmErrors(multiplicityRanges.size() - 1); // Single Muon Reference (pT > 1.0 GeV) rawSMR->SetName("rawSMR"); TGraphAsymmErrors *rawHFM = new TGraphAsymmErrors(multiplicityRanges.size() - 1); // Heavy Flavor Muon (pT > 4.0 GeV) rawHFM->SetName("rawHFM"); for (UInt_t iRange = 0; iRange < multiplicityRanges.size() - 1; iRange++) { Int_t firstBinX = histoSingleMuon->GetXaxis()->FindBin(multiplicityRanges[iRange]); Int_t lastBinX = histoSingleMuon->GetXaxis()->FindBin(multiplicityRanges[iRange+1]); Int_t firstBinY = histoSingleMuon->GetYaxis()->GetFirst(); Int_t lastBinY = histoSingleMuon->GetYaxis()->GetLast(); Int_t firstBinZ = histoSingleMuon->GetZaxis()->GetFirst(); Int_t lastBinZ = histoSingleMuon->GetZaxis()->GetLast(); Double_t nSingleMuon = 0.0; Double_t errSingleMuon = 0.0; nSingleMuon = histoSingleMuon->IntegralAndError(firstBinX, lastBinX, firstBinY, lastBinY, firstBinZ, lastBinZ, errSingleMuon); allSingleMuon->SetPoint(iRange, CINT1B->GetX()[iRange], nSingleMuon); allSingleMuon->SetPointError(iRange, CINT1B->GetX()[iRange] - multiplicityRanges[iRange], multiplicityRanges[iRange+1] - CINT1B->GetX()[iRange], errSingleMuon, errSingleMuon); // SMR firstBinZ = histoSingleMuon->GetZaxis()->FindBin(1.0); lastBinZ = histoSingleMuon->GetZaxis()->GetLast(); nSingleMuon = histoSingleMuon->IntegralAndError(firstBinX, lastBinX, firstBinY, lastBinY, firstBinZ, lastBinZ, errSingleMuon); rawSMR->SetPoint(iRange, CINT1B->GetX()[iRange], nSingleMuon); rawSMR->SetPointError(iRange, CINT1B->GetX()[iRange] - multiplicityRanges[iRange], multiplicityRanges[iRange+1] - CINT1B->GetX()[iRange], errSingleMuon, errSingleMuon); // HFM firstBinZ = histoSingleMuon->GetZaxis()->FindBin(4.0); lastBinZ = histoSingleMuon->GetZaxis()->GetLast(); nSingleMuon = histoSingleMuon->IntegralAndError(firstBinX, lastBinX, firstBinY, lastBinY, firstBinZ, lastBinZ, errSingleMuon); rawHFM->SetPoint(iRange, CINT1B->GetX()[iRange], nSingleMuon); rawHFM->SetPointError(iRange, CINT1B->GetX()[iRange] - multiplicityRanges[iRange], multiplicityRanges[iRange+1] - CINT1B->GetX()[iRange], errSingleMuon, errSingleMuon); } TList *integratedSingleMuon = new TList(); integratedSingleMuon->SetName("rawSingleMuonIntegrated"); integratedSingleMuon->AddAt(allSingleMuon, 0); integratedSingleMuon->AddAt(rawSMR, 1); integratedSingleMuon->AddAt(rawHFM, 2); rawSingleMuon->AddAt(integratedSingleMuon, 0); //delete allSingleMuon; //delete rawSMR; //delete rawHFM; // Now, for the study in pT TList *rawSingleMuonPt = new TList(); rawSingleMuonPt->SetName("rawSingleMuonPt"); for (UInt_t ipT = 0; ipT < pTRanges.size()-1; ipT++) { TString nameGraph; nameGraph.Form("pTRange_%1f-%1f", pTRanges[ipT], pTRanges[ipT+1]); TGraphAsymmErrors *singleMuonPtRange = new TGraphAsymmErrors(); singleMuonPtRange->SetName(nameGraph); for (UInt_t iRange = 0; iRange < multiplicityRanges.size() - 1; iRange++) { Int_t firstBinX = histoSingleMuon->GetXaxis()->FindBin(multiplicityRanges[iRange]); Int_t lastBinX = histoSingleMuon->GetXaxis()->FindBin(multiplicityRanges[iRange+1]); Int_t firstBinY = histoSingleMuon->GetYaxis()->GetFirst(); Int_t lastBinY = histoSingleMuon->GetYaxis()->GetLast(); Int_t firstBinZ = histoSingleMuon->GetZaxis()->FindBin(pTRanges[ipT]); Int_t lastBinZ = histoSingleMuon->GetZaxis()->FindBin(pTRanges[ipT]+1); Double_t nSingleMuon = 0.0; Double_t errSingleMuon = 0.0; nSingleMuon = histoSingleMuon->IntegralAndError(firstBinX, lastBinX, firstBinY, lastBinY, firstBinZ, lastBinZ, errSingleMuon); singleMuonPtRange->SetPoint(iRange, CINT1B->GetX()[iRange], nSingleMuon); singleMuonPtRange->SetPointError(iRange, CINT1B->GetX()[iRange] - multiplicityRanges[iRange], multiplicityRanges[iRange+1] - CINT1B->GetX()[iRange], errSingleMuon, errSingleMuon); } rawSingleMuonPt->AddAt(singleMuonPtRange, ipT); //delete singleMuonPtRange; } rawSingleMuon->AddAt(rawSingleMuonPt, 1); // At last, the study in eta // Both on SMR and HFM TList *rawSingleMuonEta = new TList(); rawSingleMuonEta->SetName("rawSingleMuonEta"); for (UInt_t iEta = 0; iEta < etaRanges.size()-1; iEta++) { TString nameGraph; nameGraph.Form("SMRetaRange_%1f-%1f", etaRanges[iEta], etaRanges[iEta+1]); TGraphAsymmErrors *SMREtaRange = new TGraphAsymmErrors(); SMREtaRange->SetName(nameGraph); nameGraph.Form("HFMetaRange_%1f-%1f", etaRanges[iEta], etaRanges[iEta+1]); TGraphAsymmErrors *HFMEtaRange = new TGraphAsymmErrors(); HFMEtaRange->SetName(nameGraph); for (UInt_t iRange = 0; iRange < multiplicityRanges.size() - 1; iRange++) { // First, SMR Int_t firstBinX = histoSingleMuon->GetXaxis()->FindBin(multiplicityRanges[iRange]); Int_t lastBinX = histoSingleMuon->GetXaxis()->FindBin(multiplicityRanges[iRange+1]); Int_t firstBinY = histoSingleMuon->GetYaxis()->FindBin(etaRanges[iEta]); Int_t lastBinY = histoSingleMuon->GetYaxis()->FindBin(etaRanges[iEta]+1); Int_t firstBinZ = histoSingleMuon->GetZaxis()->FindBin(1.0); Int_t lastBinZ = histoSingleMuon->GetZaxis()->GetLast(); Double_t nSingleMuon = 0.0; Double_t errSingleMuon = 0.0; nSingleMuon = histoSingleMuon->IntegralAndError(firstBinX, lastBinX, firstBinY, lastBinY, firstBinZ, lastBinZ, errSingleMuon); SMREtaRange->SetPoint(iRange, CINT1B->GetX()[iRange], nSingleMuon); SMREtaRange->SetPointError(iRange, CINT1B->GetX()[iRange] - multiplicityRanges[iRange], multiplicityRanges[iRange+1] - CINT1B->GetX()[iRange], errSingleMuon, errSingleMuon); // Then, HFM firstBinZ = histoSingleMuon->GetZaxis()->FindBin(4.0); lastBinZ = histoSingleMuon->GetZaxis()->GetLast(); nSingleMuon = histoSingleMuon->IntegralAndError(firstBinX, lastBinX, firstBinY, lastBinY, firstBinZ, lastBinZ, errSingleMuon); HFMEtaRange->SetPoint(iRange, CINT1B->GetX()[iRange], nSingleMuon); HFMEtaRange->SetPointError(iRange, CINT1B->GetX()[iRange] - multiplicityRanges[iRange], multiplicityRanges[iRange+1] - CINT1B->GetX()[iRange], errSingleMuon, errSingleMuon); } rawSingleMuonEta->AddAt(SMREtaRange, iEta); rawSingleMuonEta->AddAt(HFMEtaRange, iEta + etaRanges.size()-1); //delete SMREtaRange; //delete HFMEtaRange; } rawSingleMuon->AddAt(rawSingleMuonEta, 2); outputFile->WriteTObject(rawSingleMuon); delete integratedSingleMuon; delete rawSingleMuonPt; delete rawSingleMuonEta; } void AnalysisSingleMuon(TFile *outputFile) { // Analyse the raw data for single muons // First, we retrieve the CINT1B graph, as usual TGraphErrors *CINT1B = (TGraphErrors *) outputFile->Get("graphCINT1B;1"); TGraphErrors *CMUS1B = (TGraphErrors *) outputFile->Get("graphCMUS1B;1"); // Creation of all the yield graphs SingleMuonYieldGraph(outputFile, CINT1B); // Creation of the yield over mean mult graph SingleMuonYieldOverMeanMultGraph(outputFile, CINT1B); // Creation of the yield with the reference bin normalised to 1, for comparison SingleMuonYieldNormalisedGraph(outputFile, CMUS1B); } ////////////////////////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////////////////////// void SingleMuonYieldGraph(TFile *outputFile, TGraphErrors *CINT1B) { // This function use the raw graphs of outputfile and the CINT1B graph to create the yield graph // yield are raw divided by CINT1B TList *yieldSingleMuon = new TList(); yieldSingleMuon->SetName("yieldSingleMuon"); for (Int_t iList = 0; iList < ((TList *) outputFile->Get("rawSingleMuon;1"))->GetEntries(); iList++) { TList *currentList = (TList *) (( TList *) outputFile->Get("rawSingleMuon;1"))->At(iList); TList *newList = new TList(); TString listName = currentList->GetName(); listName.ReplaceAll("raw", 3, "yield", 5); newList->SetName(listName); for (Int_t iGraph = 0; iGraph < currentList->GetEntries(); iGraph++) { TGraphAsymmErrors *rawGraph = (TGraphAsymmErrors *) currentList->At(iGraph); TString yieldName = rawGraph->GetName(); yieldName.ReplaceAll("raw", 3, "yield", 5); TGraphAsymmErrors *yieldGraph = new TGraphAsymmErrors(rawGraph->GetN()); yieldGraph->SetName(yieldName); // fill the yield graph for (Int_t iBin = 0; iBin < yieldGraph->GetN(); iBin++ ) { if (CINT1B->GetY()[iBin] != 0.0) { yieldGraph->SetPoint(iBin, rawGraph->GetX()[iBin], rawGraph->GetY()[iBin]/CINT1B->GetY()[iBin]); if (rawGraph->GetY()[iBin] != 0.0) { Double_t error = yieldGraph->GetY()[iBin] * TMath::Sqrt((rawGraph->GetEYlow()[iBin]/rawGraph->GetY()[iBin])*(rawGraph->GetEYlow()[iBin]/rawGraph->GetY()[iBin]) + (CINT1B->GetEY()[iBin]/CINT1B->GetY()[iBin])*(CINT1B->GetEY()[iBin]/CINT1B->GetY()[iBin])); yieldGraph->SetPointError(iBin, rawGraph->GetEXlow()[iBin], rawGraph->GetEXhigh()[iBin], error, error); } } else yieldGraph->SetPoint(iBin, 0, 0); } newList->AddAt(yieldGraph, iGraph); } yieldSingleMuon->AddAt(newList, iList); } outputFile->WriteTObject(yieldSingleMuon); } ////////////////////////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////////////////////// void SingleMuonYieldOverMeanMultGraph(TFile *outputFile, TGraphErrors *CINT1B) { // This function use the yield graphs of outputfile to create the yield graph divided by the mean multiplicity in each bin // The CINT1B graph is necessary to get the error on the mean multiplicity TList *yieldOverMeanMultSingleMuon = new TList(); yieldOverMeanMultSingleMuon->SetName("yieldOverMeanMultSingleMuon"); for (Int_t iList = 0; iList < ((TList *) outputFile->Get("yieldSingleMuon;1"))->GetEntries(); iList++) { TList *currentList = (TList *) (( TList *) outputFile->Get("yieldSingleMuon;1"))->At(iList); TList *newList = new TList(); TString listName = currentList->GetName(); listName.ReplaceAll("yield", 5, "yieldOverMeanMult", 17); newList->SetName(listName); for (Int_t iGraph = 0; iGraph < currentList->GetEntries(); iGraph++) { TGraphAsymmErrors *yieldGraph = (TGraphAsymmErrors *) currentList->At(iGraph); TString yieldOverMeanMultName = yieldGraph->GetName(); yieldOverMeanMultName.ReplaceAll("yield", 5, "yieldOverMeanMult", 17); TGraphAsymmErrors *yieldOverMeanMultGraph = new TGraphAsymmErrors(yieldGraph->GetN()); yieldOverMeanMultGraph->SetName(yieldOverMeanMultName); // fill the yield graph for (Int_t iBin = 0; iBin < yieldOverMeanMultGraph->GetN(); iBin++ ) { if (CINT1B->GetX()[iBin] != 0.0) { yieldOverMeanMultGraph->SetPoint(iBin, yieldGraph->GetX()[iBin], yieldGraph->GetY()[iBin]/CINT1B->GetX()[iBin]); if (yieldGraph->GetY()[iBin] != 0.0) { Double_t error = yieldOverMeanMultGraph->GetY()[iBin]* TMath::Sqrt((yieldGraph->GetEYlow()[iBin]/yieldGraph->GetY()[iBin])*(yieldGraph->GetEYlow()[iBin]/yieldGraph->GetY()[iBin]) + (CINT1B->GetEX()[iBin]/CINT1B->GetX()[iBin])*(CINT1B->GetEX()[iBin]/CINT1B->GetX()[iBin])); yieldOverMeanMultGraph->SetPointError(iBin, yieldGraph->GetEXlow()[iBin], yieldGraph->GetEXhigh()[iBin], error, error); } } else yieldOverMeanMultGraph->SetPoint(iBin, 0, 0); } newList->AddAt(yieldOverMeanMultGraph, iGraph); } yieldOverMeanMultSingleMuon->AddAt(newList, iList); } outputFile->WriteTObject(yieldOverMeanMultSingleMuon); } ////////////////////////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////////////////////// void SingleMuonYieldNormalisedGraph(TFile *outputFile, TGraphErrors *CMUS1B) { // This function normalise the yield over mean multiplicity so that the reference bin in multiplicity is equal to 1. // The reference bin is defined it as the bin with the maximum number of CMUS1B TList *yieldNormalisedSingleMuon = new TList(); yieldNormalisedSingleMuon->SetName("yieldNormalisedSingleMuon"); // The first step is to find the reference bin in multiplicity // we define it as the bin with the maximum number of CMUS1B Double_t maxCMUS1B = 0.0; Int_t referenceBin = 0; for (Int_t iBin = 0; iBin < CMUS1B->GetN(); iBin++) if (CMUS1B->GetX()[iBin] > 10.0) if (CMUS1B->GetY()[iBin] > maxCMUS1B) { maxCMUS1B = CMUS1B->GetY()[iBin]; referenceBin = iBin; } // Now the loop, as usual // We don't propagate the error on the normalisation factor, since this is artificial for (Int_t iList = 0; iList < ((TList *) outputFile->Get("yieldOverMeanMultSingleMuon;1"))->GetEntries(); iList++) { TList *currentList = (TList *) (( TList *) outputFile->Get("yieldOverMeanMultSingleMuon;1"))->At(iList); TList *newList = new TList(); TString listName = currentList->GetName(); listName.ReplaceAll("yieldOverMeanMult", 17, "yieldNormalised", 15); newList->SetName(listName); for (Int_t iGraph = 0; iGraph < currentList->GetEntries(); iGraph++) { TGraphAsymmErrors *yieldOverMeanMultGraph = (TGraphAsymmErrors *) currentList->At(iGraph); TString yieldNormalisedName = yieldOverMeanMultGraph->GetName(); yieldNormalisedName.ReplaceAll("yieldOverMeanMult", 17, "yieldNormalised", 15); TGraphAsymmErrors *yieldNormalisedGraph = new TGraphAsymmErrors(yieldOverMeanMultGraph->GetN()); yieldNormalisedGraph->SetName(yieldNormalisedName); // fill the yield graph for (Int_t iBin = 0; iBin < yieldNormalisedGraph->GetN(); iBin++ ) { if (yieldOverMeanMultGraph->GetY()[iBin] != 0.0) { yieldNormalisedGraph->SetPoint(iBin, yieldOverMeanMultGraph->GetX()[iBin], yieldOverMeanMultGraph->GetY()[iBin]/yieldOverMeanMultGraph->GetY()[referenceBin]); if (yieldOverMeanMultGraph->GetY()[iBin] != 0.0) yieldNormalisedGraph->SetPointError(iBin, yieldOverMeanMultGraph->GetEXlow()[iBin], yieldOverMeanMultGraph->GetEXhigh()[iBin], yieldOverMeanMultGraph->GetEYlow()[iBin]/yieldOverMeanMultGraph->GetY()[referenceBin], yieldOverMeanMultGraph->GetEYhigh()[iBin]/yieldOverMeanMultGraph->GetY()[referenceBin]); } else yieldNormalisedGraph->SetPoint(iBin, 0, 0); } newList->AddAt(yieldNormalisedGraph, iGraph); } yieldNormalisedSingleMuon->AddAt(newList, iList); } outputFile->WriteTObject(yieldNormalisedSingleMuon); } ///////////////////////////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////////////////////////// void ProduceFitResults(TFile *inputFile, TFile *outputFile, std::vector multiplicityRanges, Bool_t applyZCut, Bool_t applyPileUpCut, Double_t numberMatchTrig, Bool_t applyThetaAbsCut, Bool_t applyPDCACut) { // Produce the raw graph for Dimuons that will be used all along the analysis // There are two graphs : // - J/Psi // - Background in the J/Psi range // Retrieve the THnSparse and the Minimum bias histo THnSparseD *inputDimuon = (THnSparseD *) ((TList *) (inputFile->Get("Dimuon;1"))->FindObject("dimuonCMUS1B") ); // Reminder of the structure of the dimuon histo // dimension 0 : multiplicity of the event // dimension 1 : is the vertex in the z range (0 for no, 1 for yes)? // dimension 2 : is it an event without pile up (0 for no, 1 for yes)? // dimension 3 : does the first muon match the trigger (0 for no, 1 for yes)? // dimension 4 : does the second muon match the trigger (0 for no, 1 for yes)? // dimension 5 : number of muons matching the trigger in the dimuon // dimension 6 : theta_abs of the first muon // dimension 7 : theta_abs of the second muon // dimension 8 : eta of the first muon // dimension 9 : eta of the second muon // dimension 10 : p DCA_x of the first muon // dimension 11 : p DCA_y of the first muon // dimension 12 : p DCA_x of the second muon // dimension 13 : p DCA_y of the second muon // dimension 14 : p of the first muon // dimension 15 : p of the second muon // dimension 16 : p of the dimuon // dimension 17 : pT of the dimuon // dimension 18 : invariant mass of the dimuon // get the pT and eta ranges from files. // Beware, the names of the files are hard-coded, so make sure to have them in your folder. // I might change this in the future, but I felt the main function had enough parameters already. ifstream pTFile ("pTRanges.txt"); Double_t pTRange = 0.0; std::vector pTRanges; while(pTFile >> pTRange) pTRanges.push_back(pTRange); ifstream etaFile ("etaRanges.txt"); Double_t etaRange = 0.0; std::vector etaRanges; while(etaFile >> etaRange) etaRanges.push_back(etaRange); // Apply all the cuts // Beware, theta_abs and pDCA cuts are hard-coded if (applyZCut) inputDimuon->GetAxis(1)->SetRangeUser(1.0, 2.0); if (applyPileUpCut) inputDimuon->GetAxis(2)->SetRangeUser(1.0, 2.0); inputDimuon->GetAxis(5)->SetRangeUser(numberMatchTrig, 3.0); if (applyThetaAbsCut) { inputDimuon->GetAxis(6)->SetRangeUser(2.0, 9.0); inputDimuon->GetAxis(7)->SetRangeUser(2.0, 9.0); } if (applyPDCACut) {// no cut yet, I need to check what are the usual values inputDimuon->GetAxis(12)->SetRangeUser(0.0, 450.0); inputDimuon->GetAxis(15)->SetRangeUser(0.0, 450.0); } // First, we get the invariant mass histos TList *dimuonFitResults = new TList(); dimuonFitResults->SetName("dimuonFitResults"); // There are some hard cuts // - eta : -4.0, -> -2.5, acceptance of the spectrometer Double_t etaMinAbsolute = -4.0; Double_t etaMaxAbsolute = -2.5; inputDimuon->GetAxis(8)->SetRangeUser(etaMinAbsolute, etaMaxAbsolute); inputDimuon->GetAxis(9)->SetRangeUser(etaMinAbsolute, etaMaxAbsolute); // First, we get the invariant mass histo for the whole range in multiplicity inputDimuon->GetAxis(0)->SetRangeUser(multiplicityRanges[0], multiplicityRanges[multiplicityRanges.size()-1]); TH1D *invariantMassIntegrated = inputDimuon->Projection(18, "E"); invariantMassIntegrated->SetName("invariantMassIntegrated"); // Now for the fit on all the multiplicity // This is used to get a value for some parameters TList *fitAll = new TList(); fitAll->SetName("AllMultiplicity"); // Create the container for the reference parameters TH1D *referenceParameters = new TH1D("referenceParameters", "parameters", 8, 1.0, 9.0); referenceParameters->Sumw2(); referenceParameters->GetXaxis()->SetBinLabel(1, "meanJPsi"); referenceParameters->GetXaxis()->SetBinLabel(2, "sigmaJPsi"); referenceParameters->GetXaxis()->SetBinLabel(3, "alphaJPsi"); referenceParameters->GetXaxis()->SetBinLabel(4, "nJPsi"); referenceParameters->GetXaxis()->SetBinLabel(5, "meanPsiPrime"); referenceParameters->GetXaxis()->SetBinLabel(6, "sigmaPsiPrime"); referenceParameters->GetXaxis()->SetBinLabel(7, "expBkg1"); referenceParameters->GetXaxis()->SetBinLabel(8, "expBkg2"); FitInvariantMass(fitAll, invariantMassIntegrated, kFALSE, referenceParameters); dimuonFitResults->AddAt(fitAll, 0); // Now, we make a fit for each range in multiplicity for (UInt_t iMult = 0; iMult < multiplicityRanges.size()-1; iMult++) { TString name; name.Form("Multiplicity_%d-%d", static_cast(multiplicityRanges[iMult]), static_cast(multiplicityRanges[iMult+1])); TList *fitRange = new TList(); fitRange->SetName(name); inputDimuon->GetAxis(0)->SetRangeUser(multiplicityRanges[iMult], multiplicityRanges[iMult+1]); TH1D *invariantMassRange = inputDimuon->Projection(18, "E"); invariantMassRange->SetName("invariantMass"); FitInvariantMass(fitRange, invariantMassRange, kTRUE, referenceParameters); dimuonFitResults->AddAt(fitRange, iMult + 1); } outputFile->WriteTObject(dimuonFitResults); } ///////////////////////////////////////////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////////////////////////////////////////// void FitInvariantMass(TList *fitList, TH1D *invariantMass, Bool_t areParametersFixed, TH1D *referenceParameters) { // Fit the given histo and put the results in outputFile // The fit is : Crystal Ball (J/Psi) + gaussian (Psi') + two exponentials (background) // The parameters histo architecture is the following : // bin 1 : mean J/Psi // bin 2 : sigma J/Psi // bin 3 : alpha J/Psi // bin 4 : n J/Psi // bin 5 : mean Psi' // bin 6 : sigma Psi' // bin 7 : exponential factor background 1 // bin 8 : exponential factor background 2 // The results histo architecture is the following : // bin 1 : J/Psi // bin 2 : Background J/Psi at 2 sigma // bin 3 : signal over background // bin 4 : chi2 // There are some parameters that are not always free : // mean of J/Psi peak // sigma of J/Psi peak // alpha of J/Psi function (Crystall Ball function) // n of J/Psi function // mean of Psi' peak // sigma of Psi' peak // If there are fixed, their values are taken from the parameters histo // First, we have to create two histos, one that will contain the parameters of the fit, and the other the results TH1D *parameters = new TH1D("parameters", "parameters", 8, 1.0, 9.0); parameters->Sumw2(); parameters->GetXaxis()->SetBinLabel(1, "meanJPsi"); parameters->GetXaxis()->SetBinLabel(2, "sigmaJPsi"); parameters->GetXaxis()->SetBinLabel(3, "alphaJPsi"); parameters->GetXaxis()->SetBinLabel(4, "nJPsi"); parameters->GetXaxis()->SetBinLabel(5, "meanPsiPrime"); parameters->GetXaxis()->SetBinLabel(6, "sigmaPsiPrime"); parameters->GetXaxis()->SetBinLabel(7, "expBkg1"); parameters->GetXaxis()->SetBinLabel(8, "expBkg2"); TH1D *results = new TH1D("results", "results", 4, 1.0, 5.0); results->Sumw2(); results->GetXaxis()->SetBinLabel(1, "JPsi"); results->GetXaxis()->SetBinLabel(2, "bkg"); results->GetXaxis()->SetBinLabel(3, "SB"); results->GetXaxis()->SetBinLabel(4, "chi2"); // Declare observable M RooRealVar M("M","Dimuon invariant mass [GeV/c^2]", 2.0, 5.0); // Declare Cristal Ball for J/Psi RooRealVar mean_jpsi("mean_jpsi","mean_jpsi", 3.096, 2.8, 3.2); RooRealVar sigma_jpsi("sigma_jpsi","sigma_jpsi", 0.070, 0, 0.2); RooRealVar alpha_jpsi("alpha","alpha", 5, 0, 10); RooRealVar n_jpsi("n","n",5,0,10); RooCBShape jPsi("jpsi", "crystal ball PDF", M, mean_jpsi, sigma_jpsi, alpha_jpsi, n_jpsi); //Declare gaussian for Psi prime RooRealVar mean_psip("mean_psip", "mean_psip", 3.686, 3.6, 3.75); RooRealVar sigma_psip("sigma_psip","sigma_psip", 0.086); RooGaussian psiPrime("psip","Gaussian", M, mean_psip, sigma_psip); //Declare exponentials for background RooRealVar c1("c1", "c1", -10, -0.1); RooExponential bkg_exp1("bkg_exp1", "background", M, c1); RooRealVar c2("c2", "c2", -5.0, -1); RooExponential bkg_exp2("bkg_exp2", "background", M, c2); // Give a value and fix sigma, alpha and n in the Cristal Ball // The values correspond to a previous fit (ideally, with all the statistic) // I should automatize this one of these day if (areParametersFixed) { Double_t fixMeanJPsi = referenceParameters->GetBinContent(parameters->GetXaxis()->FindBin("meanJPsi")); mean_jpsi.setVal(fixMeanJPsi); mean_jpsi.setError(0.0); mean_jpsi.setRange(fixMeanJPsi, fixMeanJPsi); Double_t fixSigmaJPsi = referenceParameters->GetBinContent(parameters->GetXaxis()->FindBin("sigmaJPsi")); sigma_jpsi.setVal(fixSigmaJPsi); sigma_jpsi.setError(0.0); sigma_jpsi.setRange(fixSigmaJPsi, fixSigmaJPsi); Double_t fixAlphaJPsi = referenceParameters->GetBinContent(parameters->GetXaxis()->FindBin("alphaJPsi")); alpha_jpsi.setVal(fixAlphaJPsi); alpha_jpsi.setError(0.0); alpha_jpsi.setRange(fixAlphaJPsi, fixAlphaJPsi); Double_t fixNJPsi = referenceParameters->GetBinContent(parameters->GetXaxis()->FindBin("nJPsi")); n_jpsi.setVal(fixNJPsi); n_jpsi.setError(0.0); n_jpsi.setRange(fixNJPsi, fixNJPsi); Double_t fixMeanPsiPrime = referenceParameters->GetBinContent(parameters->GetXaxis()->FindBin("meanPsiPrime")); mean_psip.setVal(fixMeanPsiPrime); mean_psip.setError(0.0); mean_psip.setRange(fixMeanPsiPrime, fixMeanPsiPrime); Double_t fixSigmaPsiPrime = referenceParameters->GetBinContent(parameters->GetXaxis()->FindBin("sigmaPsiPrime")); sigma_psip.setVal(fixSigmaPsiPrime); sigma_psip.setError(0.0); sigma_psip.setRange(fixSigmaPsiPrime, fixSigmaPsiPrime); } // Sum the composite signal and background into an extended pdf nsig*sig+nbkg*bkg RooRealVar fitJPsi("fitJPsi", "number of J/Psi events", 1600 ,0.0, 15000); RooRealVar fitPsiPrime("fitPsiPrime", "number of Psi Prime events", 50, 0.0, 200); RooRealVar fitBkg1("fitBkg1", "number of background events", 5000, 0.0, 15000); RooRealVar fitBkg2("fitBkg2", "number of background events", 5000, 0.0, 15000); RooAddPdf *fitFunction = new RooAddPdf("model", "(CB+Gauss+2exp)", RooArgList(bkg_exp1, bkg_exp2, jPsi, psiPrime), RooArgList(fitBkg1, fitBkg2, fitJPsi, fitPsiPrime)); // Define the histo to be fitted RooDataHist *invariantMassFit = new RooDataHist("invariantMassFit", "invariantMassFit", M, invariantMass); // Fit the invariant mass // This is the important part RooFitResult *fitResult = fitFunction->fitTo(*invariantMassFit, RooFit::Save()); //RooFitResult *fitResult = fitFunction->chi2FitTo(*invariantMassFit, RooFit::Save()); RooPlot *plot = M.frame(); TString title = "fittedInvariantMass"; plot->SetTitle(title); plot->SetName(title); invariantMassFit->plotOn(plot, RooFit::Name("invariantMassFit")); fitFunction->plotOn(plot, RooFit::Name("fitFunction"), RooFit::Range(2.0, 5.0)); // Fill the histo with the final values of the parameters parameters->SetBinContent(parameters->GetXaxis()->FindBin("meanJPsi"), mean_jpsi.getVal()); parameters->SetBinError(parameters->GetXaxis()->FindBin("meanJPsi"), mean_jpsi.getError()); parameters->SetBinContent(parameters->GetXaxis()->FindBin("sigmaJPsi"), sigma_jpsi.getVal()); parameters->SetBinError(parameters->GetXaxis()->FindBin("sigmaJPsi"), sigma_jpsi.getError()); parameters->SetBinContent(parameters->GetXaxis()->FindBin("alphaJPsi"), alpha_jpsi.getVal()); parameters->SetBinError(parameters->GetXaxis()->FindBin("alphaJPsi"), alpha_jpsi.getError()); parameters->SetBinContent(parameters->GetXaxis()->FindBin("nJPsi"), n_jpsi.getVal()); parameters->SetBinError(parameters->GetXaxis()->FindBin("nJPsi"), n_jpsi.getError()); parameters->SetBinContent(parameters->GetXaxis()->FindBin("meanPsiPrime"), mean_psip.getVal()); parameters->SetBinError(parameters->GetXaxis()->FindBin("meanPsiPrime"), mean_psip.getError()); parameters->SetBinContent(parameters->GetXaxis()->FindBin("sigmaPsiPrime"), sigma_psip.getVal()); parameters->SetBinError(parameters->GetXaxis()->FindBin("sigmaPsiPrime"), sigma_psip.getError()); // Fill the histo with the reference parameters if (!areParametersFixed) { referenceParameters->SetBinContent(parameters->GetXaxis()->FindBin("meanJPsi"), mean_jpsi.getVal()); referenceParameters->SetBinError(parameters->GetXaxis()->FindBin("meanJPsi"), mean_jpsi.getError()); referenceParameters->SetBinContent(parameters->GetXaxis()->FindBin("sigmaJPsi"), sigma_jpsi.getVal()); referenceParameters->SetBinError(parameters->GetXaxis()->FindBin("sigmaJPsi"), sigma_jpsi.getError()); referenceParameters->SetBinContent(parameters->GetXaxis()->FindBin("alphaJPsi"), alpha_jpsi.getVal()); referenceParameters->SetBinError(parameters->GetXaxis()->FindBin("alphaJPsi"), alpha_jpsi.getError()); referenceParameters->SetBinContent(parameters->GetXaxis()->FindBin("nJPsi"), n_jpsi.getVal()); referenceParameters->SetBinError(parameters->GetXaxis()->FindBin("nJPsi"), n_jpsi.getError()); referenceParameters->SetBinContent(parameters->GetXaxis()->FindBin("meanPsiPrime"), mean_psip.getVal()); referenceParameters->SetBinError(parameters->GetXaxis()->FindBin("meanPsiPrime"), mean_psip.getError()); referenceParameters->SetBinContent(parameters->GetXaxis()->FindBin("sigmaPsiPrime"), sigma_psip.getVal()); referenceParameters->SetBinError(parameters->GetXaxis()->FindBin("sigmaPsiPrime"), sigma_psip.getError()); } parameters->SetBinContent(parameters->GetXaxis()->FindBin("expBkg1"), c1.getVal()); parameters->SetBinError(parameters->GetXaxis()->FindBin("expBkg1"), c1.getError()); parameters->SetBinContent(parameters->GetXaxis()->FindBin("expBkg2"), c2.getVal()); parameters->SetBinError(parameters->GetXaxis()->FindBin("expBkg2"), c2.getError()); // Define the range to compute the results // We choose two sigma // We also need two new variables to compute the background Double_t lowerBoundJPsi = mean_jpsi.getVal() - 2.0*sigma_jpsi.getVal(); Double_t upperBoundJPsi = mean_jpsi.getVal() + 2.0*sigma_jpsi.getVal(); M.setRange("JPsiRange" ,lowerBoundJPsi, upperBoundJPsi); RooAbsReal *fracBkgRange1 = bkg_exp1.createIntegral(M, M, "JPsiRange"); RooAbsReal *fracBkgRange2 = bkg_exp2.createIntegral(M, M, "JPsiRange"); // Get the results, and fill the results histo Double_t nJPsi = fitJPsi.getVal(); Double_t errJPsi = fitJPsi.getError(); results->SetBinContent(results->GetXaxis()->FindBin("JPsi"), nJPsi); results->SetBinError(results->GetXaxis()->FindBin("JPsi"), errJPsi); Double_t nBkg = (fitBkg1.getVal() * fracBkgRange1->getVal() + fitBkg2.getVal() * fracBkgRange2->getVal()); Double_t errBkg = (fitBkg1.getError() * fracBkgRange1->getVal() + fitBkg2.getError() * fracBkgRange2->getVal()); results->SetBinContent(results->GetXaxis()->FindBin("bkg"), nBkg); results->SetBinError(results->GetXaxis()->FindBin("bkg"), errBkg); Double_t SB = 0.0; Double_t errSB = 0.0; if (nJPsi != 0.0 && nBkg != 0.0) { SB = nJPsi/nBkg; errSB = SB * (errJPsi/nJPsi + errBkg/nBkg); } results->SetBinContent(results->GetXaxis()->FindBin("SB"), SB); results->SetBinError(results->GetXaxis()->FindBin("SB"), errSB); Int_t nDF = fitResult->floatParsFinal().getSize(); Double_t chi2 = plot->chiSquare("fitFunction", "invariantMassFit", nDF); results->SetBinContent(results->GetXaxis()->FindBin("chi2"), chi2); results->SetBinError(results->GetXaxis()->FindBin("chi2"), 0.0); fitList->AddAt(invariantMass, 0); fitList->AddAt(plot, 1); fitList->AddAt(parameters, 2); fitList->AddAt(results, 3); } void ProduceDimuonGraph(TFile *outputFile, std::vector multiplicityRanges) { // This function will create all the graphs for dimuons (J/Psi and Background) // - raw // - yield // - yield over mean mult // - yield over mean mult normalised // First, get the CINT1B and CMUS1B graphs TGraphErrors *CINT1B = (TGraphErrors *) outputFile->Get("graphCINT1B"); TGraphErrors *CMUS1B = (TGraphErrors *) outputFile->Get("graphCMUS1B"); TList *JPsiGraph = new TList(); JPsiGraph->SetName("JPsiGraph"); TList *bkgGraph = new TList(); bkgGraph->SetName("bkgGraph"); // Raw Graphs TGraphAsymmErrors *rawJPsiGraph = new TGraphAsymmErrors(multiplicityRanges.size()-1); rawJPsiGraph->SetName("rawJPsiGraph"); TGraphAsymmErrors *rawBkgGraph = new TGraphAsymmErrors(multiplicityRanges.size()-1); rawBkgGraph->SetName("rawBkgGraph"); for (UInt_t iMult = 0; iMult < multiplicityRanges.size()-1; iMult++) { TString name; name.Form("Multiplicity_%d-%d", static_cast(multiplicityRanges[iMult]), static_cast(multiplicityRanges[iMult+1])); TH1D *results = (TH1D *) ((TList *) ((TList *) outputFile->Get("dimuonFitResults;1"))->FindObject(name))->FindObject("results"); // J/Psi Double_t nJPsi = results->GetBinContent(results->GetXaxis()->FindBin("JPsi")); Double_t errJPsi = results->GetBinError(results->GetXaxis()->FindBin("JPsi")); rawJPsiGraph->SetPoint(iMult, CINT1B->GetX()[iMult], nJPsi); rawJPsiGraph->SetPointError(iMult, CINT1B->GetX()[iMult] - multiplicityRanges[iMult], multiplicityRanges[iMult+1] - CINT1B->GetX()[iMult], errJPsi, errJPsi); // Background Double_t nBkg = results->GetBinContent(results->GetXaxis()->FindBin("bkg")); Double_t errBkg = results->GetBinError(results->GetXaxis()->FindBin("bkg")); rawJPsiGraph->SetPoint(iMult, CINT1B->GetX()[iMult], nBkg); rawJPsiGraph->SetPointError(iMult, CINT1B->GetX()[iMult] - multiplicityRanges[iMult], multiplicityRanges[iMult+1] - CINT1B->GetX()[iMult], errBkg, errBkg); } JPsiGraph->AddAt(rawJPsiGraph, 0); bkgGraph->AddAt(rawBkgGraph, 0); // Yield Graphs TGraphAsymmErrors *yieldJPsiGraph = new TGraphAsymmErrors(multiplicityRanges.size()-1); yieldJPsiGraph->SetName("yieldJPsiGraph"); TGraphAsymmErrors *yieldBkgGraph = new TGraphAsymmErrors(multiplicityRanges.size()-1); yieldBkgGraph->SetName("yieldBkgGraph"); for (UInt_t iMult = 0; iMult < multiplicityRanges.size()-1; iMult++) { if (CINT1B->GetY()[iMult] != 0.0) { // J/Psi yieldJPsiGraph->SetPoint(iMult, rawJPsiGraph->GetX()[iMult], rawJPsiGraph->GetY()[iMult]/CINT1B->GetY()[iMult]); if (rawJPsiGraph->GetY()[iMult] != 0.0) { Double_t error = yieldJPsiGraph->GetY()[iMult]* TMath::Sqrt((rawJPsiGraph->GetEYlow()[iMult]/rawJPsiGraph->GetY()[iMult])*(rawJPsiGraph->GetEYlow()[iMult]/rawJPsiGraph->GetY()[iMult]) + (CINT1B->GetEY()[iMult]/CINT1B->GetY()[iMult])*(CINT1B->GetEY()[iMult]/CINT1B->GetY()[iMult])); yieldJPsiGraph->SetPointError(iMult, rawJPsiGraph->GetEXlow()[iMult], rawJPsiGraph->GetEXhigh()[iMult], error, error); } // Background yieldBkgGraph->SetPoint(iMult, rawBkgGraph->GetX()[iMult], rawBkgGraph->GetY()[iMult]/CINT1B->GetY()[iMult]); if (rawBkgGraph->GetY()[iMult] != 0.0) { Double_t error = yieldBkgGraph->GetY()[iMult]* TMath::Sqrt((rawBkgGraph->GetEYlow()[iMult]/rawBkgGraph->GetY()[iMult])*(rawBkgGraph->GetEYlow()[iMult]/rawBkgGraph->GetY()[iMult]) + (CINT1B->GetEY()[iMult]/CINT1B->GetY()[iMult])*(CINT1B->GetEY()[iMult]/CINT1B->GetY()[iMult])); yieldBkgGraph->SetPointError(iMult, rawBkgGraph->GetEXlow()[iMult], rawBkgGraph->GetEXhigh()[iMult], error, error); } } else { yieldJPsiGraph->SetPoint(iMult, rawJPsiGraph->GetX()[iMult], 0); yieldBkgGraph->SetPoint(iMult, rawBkgGraph->GetX()[iMult], 0); } } JPsiGraph->AddAt(yieldJPsiGraph, 1); bkgGraph->AddAt(yieldBkgGraph, 1); // Yield over Mean Mult graph TGraphAsymmErrors *yieldOverMeanMultJPsiGraph = new TGraphAsymmErrors(multiplicityRanges.size()-1); yieldOverMeanMultJPsiGraph->SetName("yieldOverMeanMultJPsiGraph"); TGraphAsymmErrors *yieldOverMeanMultBkgGraph = new TGraphAsymmErrors(multiplicityRanges.size()-1); yieldOverMeanMultBkgGraph->SetName("yieldOverMeanMultBkgGraph"); for (UInt_t iMult = 0; iMult < multiplicityRanges.size()-1; iMult++) { if (CINT1B->GetX()[iMult] != 0.0) { // J/Psi yieldOverMeanMultJPsiGraph->SetPoint(iMult, yieldJPsiGraph->GetX()[iMult], yieldJPsiGraph->GetY()[iMult]/CINT1B->GetX()[iMult]); if (yieldJPsiGraph->GetY()[iMult] != 0.0) { Double_t error = yieldOverMeanMultJPsiGraph->GetY()[iMult]* TMath::Sqrt((yieldJPsiGraph->GetEYlow()[iMult]/yieldJPsiGraph->GetY()[iMult])*(yieldJPsiGraph->GetEYlow()[iMult]/yieldJPsiGraph->GetY()[iMult]) + (CINT1B->GetEX()[iMult]/CINT1B->GetX()[iMult])*(CINT1B->GetEX()[iMult]/CINT1B->GetX()[iMult])); yieldOverMeanMultJPsiGraph->SetPointError(iMult, yieldJPsiGraph->GetEXlow()[iMult], yieldJPsiGraph->GetEXhigh()[iMult], error, error); } // Background yieldOverMeanMultBkgGraph->SetPoint(iMult, yieldBkgGraph->GetX()[iMult], yieldBkgGraph->GetY()[iMult]/CINT1B->GetX()[iMult]); if (yieldBkgGraph->GetY()[iMult] != 0.0) { Double_t error = yieldOverMeanMultBkgGraph->GetY()[iMult]* TMath::Sqrt((yieldBkgGraph->GetEYlow()[iMult]/yieldBkgGraph->GetY()[iMult])*(yieldBkgGraph->GetEYlow()[iMult]/yieldBkgGraph->GetY()[iMult]) + (CINT1B->GetEX()[iMult]/CINT1B->GetX()[iMult])*(CINT1B->GetEX()[iMult]/CINT1B->GetX()[iMult])); yieldOverMeanMultBkgGraph->SetPointError(iMult, yieldBkgGraph->GetEXlow()[iMult], yieldBkgGraph->GetEXhigh()[iMult], error, error); } } else { yieldOverMeanMultJPsiGraph->SetPoint(iMult, yieldJPsiGraph->GetX()[iMult], 0); yieldOverMeanMultBkgGraph->SetPoint(iMult, yieldBkgGraph->GetX()[iMult], 0); } } JPsiGraph->AddAt(yieldOverMeanMultJPsiGraph, 2); bkgGraph->AddAt(yieldOverMeanMultBkgGraph, 2); // Yield over Mean Mult normalised to get the bin with the highest number of CMUS1B equal to 1 Double_t maxCMUS1B = 0.0; Int_t referenceBin = 0; for (Int_t iBin = 0; iBin < CMUS1B->GetN(); iBin++) if (CMUS1B->GetX()[iBin] > 10.0) if (CMUS1B->GetY()[iBin] > maxCMUS1B) { maxCMUS1B = CMUS1B->GetY()[iBin]; referenceBin = iBin; } TGraphAsymmErrors *yieldNormalisedJPsiGraph = new TGraphAsymmErrors(multiplicityRanges.size()-1); yieldNormalisedJPsiGraph->SetName("yieldNormalisedJPsiGraph"); TGraphAsymmErrors *yieldNormalisedBkgGraph = new TGraphAsymmErrors(multiplicityRanges.size()-1); yieldNormalisedBkgGraph->SetName("yieldNormalisedBkgGraph"); for (UInt_t iMult = 0; iMult < multiplicityRanges.size()-1; iMult++) { // JPsi if (yieldOverMeanMultJPsiGraph->GetY()[referenceBin] != 0.0) { yieldNormalisedJPsiGraph->SetPoint(iMult, yieldOverMeanMultJPsiGraph->GetX()[iMult], yieldOverMeanMultJPsiGraph->GetY()[iMult]/yieldOverMeanMultJPsiGraph->GetY()[referenceBin]); yieldNormalisedJPsiGraph->SetPointError(iMult, yieldOverMeanMultJPsiGraph->GetEXlow()[iMult], yieldOverMeanMultJPsiGraph->GetEXhigh()[iMult], yieldOverMeanMultJPsiGraph->GetEYlow()[iMult]/yieldOverMeanMultJPsiGraph->GetY()[referenceBin], yieldOverMeanMultJPsiGraph->GetEYhigh()[iMult]/yieldOverMeanMultJPsiGraph->GetY()[referenceBin]); } // Bkg if (yieldOverMeanMultBkgGraph->GetY()[referenceBin] != 0.0) { yieldNormalisedBkgGraph->SetPoint(iMult, yieldOverMeanMultBkgGraph->GetX()[iMult], yieldOverMeanMultBkgGraph->GetY()[iMult]/yieldOverMeanMultBkgGraph->GetY()[referenceBin]); yieldNormalisedBkgGraph->SetPointError(iMult, yieldOverMeanMultBkgGraph->GetEXlow()[iMult], yieldOverMeanMultBkgGraph->GetEXhigh()[iMult], yieldOverMeanMultBkgGraph->GetEYlow()[iMult]/yieldOverMeanMultBkgGraph->GetY()[referenceBin], yieldOverMeanMultBkgGraph->GetEYhigh()[iMult]/yieldOverMeanMultBkgGraph->GetY()[referenceBin]); } } JPsiGraph->AddAt(yieldNormalisedJPsiGraph, 3); bkgGraph->AddAt(yieldNormalisedBkgGraph, 3); outputFile->WriteTObject(JPsiGraph); outputFile->WriteTObject(bkgGraph); }