//------------------------------------------------------------------------- // Implementation of Class AliCollisionNormalization // // This class is used to store the vertex ditributions in the data // and in Monte Carlo, needed to compute the real number of // collisions a given sample is corresponding to. // The strategy matches what described in CERN-THESIS-2009-033 p 119. // // Author: Michele Floris, CERN //------------------------------------------------------------------------- #include "AliCollisionNormalization.h" #include "AliPhysicsSelection.h" #include "AliLog.h" #include "TFile.h" #include "TCanvas.h" #include "AliGenPythiaEventHeader.h" #include "AliGenDPMjetEventHeader.h" #include "AliGenEventHeader.h" #include "AliMCEvent.h" ClassImp(AliCollisionNormalization) const char * AliCollisionNormalization::fProcLabel[] = {"SD","DD","ND", "Unknown"}; AliCollisionNormalization::AliCollisionNormalization() : TObject(), fNbinsVz(0), fMinVz(0), fMaxVz(0), fZRange(9.99), fIsMC(0), fReferenceXS(0), fVerbose(0), fEnergy(900), fHistVzData (0), fHistProcTypes (0), fHistStatBin0 (0), fHistStat (0) { // ctor fNbinsVz = 30; fMinVz = -15; fMaxVz = +15; for(Int_t iproc = 0; iproc < kNProcs; iproc++){ fHistVzMCGen[iproc] = 0; fHistVzMCRec[iproc] = 0; fHistVzMCTrg[iproc] = 0; } BookAllHistos(); } AliCollisionNormalization::AliCollisionNormalization(Int_t nbinz, Float_t minz, Float_t maxz): TObject(), fNbinsVz(0), fMinVz(0), fMaxVz(0), fZRange(9.99), fIsMC(0), fReferenceXS(0), fVerbose(0), fEnergy(900), fHistVzData (0), fHistProcTypes (0), fHistStatBin0 (0), fHistStat (0) { // ctor, allows setting binning fNbinsVz = nbinz; fMinVz = minz; fMaxVz = maxz; for(Int_t iproc = 0; iproc < kNProcs; iproc++){ fHistVzMCGen[iproc] = 0; fHistVzMCRec[iproc] = 0; fHistVzMCTrg[iproc] = 0; } BookAllHistos(); } AliCollisionNormalization::AliCollisionNormalization(const char * dataFile, const char * dataListName, const char * mcFile, const char * mcListName, const char * eventStatFile) : TObject(), fNbinsVz(0), fMinVz(0), fMaxVz(0), fZRange(9.99), fIsMC(0), fReferenceXS(0), fVerbose(0), fEnergy(900), fHistVzData (0), fHistProcTypes (0), fHistStatBin0 (0), fHistStat (0) { // ctor, loads histograms from file for(Int_t iproc = 0; iproc < kNProcs; iproc++){ fHistVzMCGen[iproc] = 0; fHistVzMCRec[iproc] = 0; fHistVzMCTrg[iproc] = 0; } TFile * fdata = new TFile (dataFile); TFile * fmc = new TFile (mcFile ); TFile * fstat = new TFile(eventStatFile); if (!fdata->IsOpen() || !fmc->IsOpen() || !fstat->IsOpen()) { AliFatal("Cannot open input file(s)"); } TList * ldata = (TList*) fdata->Get(dataListName); TList * lmc = (TList*) fmc ->Get(mcListName ); AliCollisionNormalization * cndata = (AliCollisionNormalization*) ldata->FindObject("AliCollisionNormalization"); AliCollisionNormalization * cnmc = (AliCollisionNormalization*) lmc ->FindObject("AliCollisionNormalization"); // Assign or book all histos for(Int_t iproc = 0; iproc < kNProcs; iproc++){ fHistVzMCGen[iproc]= cnmc->GetVzMCGen(iproc); fHistVzMCRec[iproc]= cnmc->GetVzMCRec(iproc); fHistVzMCTrg[iproc]= cnmc->GetVzMCTrg(iproc); } fHistVzData = cndata->GetVzData(); fHistProcTypes = cnmc->GetHistProcTypes(); fHistStatBin0 = (TH1F*) fstat->Get("fHistStatistics_Bin0"); fHistStat = (TH1F*) fstat->Get("fHistStatistics"); } AliCollisionNormalization::~AliCollisionNormalization(){ // dtor for(Int_t iproc = 0; iproc < kNProcs; iproc++){ if(fHistVzMCGen[iproc]) { delete fHistVzMCGen[iproc] ; fHistVzMCGen[iproc] =0;} if(fHistVzMCRec[iproc]) { delete fHistVzMCRec[iproc] ; fHistVzMCRec[iproc] =0;} if(fHistVzMCTrg[iproc]) { delete fHistVzMCTrg[iproc] ; fHistVzMCTrg[iproc] =0;} } if(fHistVzData ) { delete fHistVzData ; fHistVzData =0;} if(fHistStatBin0 ) { delete fHistStatBin0 ; fHistStatBin0 =0;} if(fHistStat ) { delete fHistStat ; fHistStat =0;} if(fHistProcTypes ) { delete fHistProcTypes ; fHistProcTypes =0;} } void AliCollisionNormalization::BookAllHistos(){ // books all histos // Book histos of vz distributions vs multiplicity // if vzOnly == kTRUE, it returns a 1d histo with vz dist only // Do not attach histos to the current directory Bool_t oldStatus = TH1::AddDirectoryStatus(); TH1::AddDirectory(kFALSE); for(Int_t iproc = 0; iproc < kNProcs; iproc++){ fHistVzMCGen [iproc] = (TH2F*) BookVzHisto(TString("fHistVzMCGen")+ fProcLabel[iproc] ,"Vz distribution of generated events vs rec multiplicity "); fHistVzMCRec [iproc] = (TH2F*) BookVzHisto(TString("fHistVzMCRec")+ fProcLabel[iproc] ,"Vz distribution of reconstructed events vs rec multiplicity"); fHistVzMCTrg [iproc] = (TH2F*) BookVzHisto(TString("fHistVzMCTrg")+ fProcLabel[iproc] ,"Vz distribution of triggered events vs rec multiplicity "); } fHistVzData = (TH2F*) BookVzHisto("fHistVzData" ,"Vz distribution of triggered events vs rec multiplicity "); fHistProcTypes = new TH1F ("fHistProcTypes", "Number of events in the different process classes", kNProcs, -0.5 , kNProcs-0.5); fHistProcTypes->GetXaxis()->SetBinLabel(kProcSD+1,"SD"); fHistProcTypes->GetXaxis()->SetBinLabel(kProcND+1,"ND"); fHistProcTypes->GetXaxis()->SetBinLabel(kProcDD+1,"DD"); fHistProcTypes->GetXaxis()->SetBinLabel(kProcUnknown+1,"Unknown"); TH1::AddDirectory(oldStatus); } TH1 * AliCollisionNormalization::BookVzHisto(const char * name , const char * title, Bool_t vzOnly) { // Book histos of vz distributions vs multiplicity // if vzOnly == kTRUE, it returns a 1d histo with vz dist only // Do not attach histos to the current directory Bool_t oldStatus = TH1::AddDirectoryStatus(); TH1::AddDirectory(kFALSE); TH1 * h; Double_t binLimitsVtx[] = {-30,-25,-20,-15,-10,-7,-5.5,-4,-3,-2,-1,0,1,2,3,4,5.5,7,10,15,20,25,30}; if (vzOnly) { h = new TH1F(name,title,22,binLimitsVtx); } else { h = new TH2F(name,title,22,binLimitsVtx,100,Double_t(-0.5),Double_t(99.5)); } h->SetXTitle("V_{z} (cm)"); h->SetYTitle("n_{trk}"); h->Sumw2(); TH1::AddDirectory(oldStatus); return h; } TH2F * AliCollisionNormalization::GetVzMCGen (Int_t procType) { // returns MC gen histo. If proc type < 0 sums over all proc types, reweighting the XS if(procType>=0) return fHistVzMCGen[procType] ; TH2F * sum = (TH2F*) fHistVzMCGen[kProcSD]->Clone(); sum->Reset(); for(Int_t iproc = 0; iproc < kProcUnknown; iproc++){ sum->Add(fHistVzMCGen[iproc],GetProcessWeight(iproc)); } return sum; } TH2F * AliCollisionNormalization::GetVzMCRec (Int_t procType) { // returns MC rec histo. If proc type < 0 sums over all proc types, reweighting the XS if(procType>=0) return fHistVzMCRec[procType] ; TH2F * sum = (TH2F*) fHistVzMCRec[kProcSD]->Clone(); sum->Reset(); for(Int_t iproc = 0; iproc < kProcUnknown; iproc++){ sum->Add(fHistVzMCRec[iproc],GetProcessWeight(iproc)); } return sum; } TH2F * AliCollisionNormalization::GetVzMCTrg (Int_t procType) { // returns MC trg histo. If proc type < 0 sums over all proc types, reweighting the XS if(procType>=0) return fHistVzMCTrg[procType] ; TH2F * sum = (TH2F*) fHistVzMCTrg[kProcSD]->Clone(); sum->Reset(); for(Int_t iproc = 0; iproc < kProcUnknown; iproc++){ sum->Add(fHistVzMCTrg[iproc],GetProcessWeight(iproc)); } return sum; } Double_t AliCollisionNormalization::ComputeNint() { // Compute the number of collisions applying all corrections // TODO: check error propagation TH1 * hVzData = fHistVzData->ProjectionX("_px",2,-1,"E"); // Skip zero bin Int_t allEventsWithVertex = (Int_t) fHistVzData->Integral(0, fHistVzData->GetNbinsX()+1, 2, fHistVzData->GetNbinsY()+1); // include under/overflow!, skip 0 bin // Assign histos reweighted with measured XS TH2F * histVzMCGen = GetVzMCGen(-1); TH2F * histVzMCTrg = GetVzMCTrg(-1); TH2F * histVzMCRec = GetVzMCRec(-1); // Before we start: print number of input events to the physics // selection: this allows to crosscheck that all runs were // successfully processed (useful if running on grid: you may have a // crash without noticing it). AliInfo(Form("Input Events (No cuts: %d, After Phys. Sel.:%d)", Int_t(fHistStat->GetBinContent(1,1)), Int_t(fHistStat->GetBinContent(fHistStat->GetNbinsX(),1)))); // Fixme: will this change with a different trigger scheme? // Get or compute BG. This assumes the CINT1B suite Double_t triggeredEventsWith0MultWithBG = fHistVzData->Integral(0, fHistVzData->GetNbinsX()+1,1, 1); Double_t bg = 0; // This will include beam gas + accidentals if (fHistStatBin0->GetNbinsY() > 4) { // FIXME: we need a better criterion to decide... AliInfo("Using BG computed by Physics Selection"); bg = fHistStatBin0->GetBinContent(fHistStatBin0->GetNbinsX(),AliPhysicsSelection::kStatRowBG); bg += fHistStatBin0->GetBinContent(fHistStatBin0->GetNbinsX(),AliPhysicsSelection::kStatRowAcc); Int_t cint1B = (Int_t) fHistStatBin0->GetBinContent(fHistStatBin0->GetNbinsX(),1); if (cint1B != Int_t(triggeredEventsWith0MultWithBG)) { AliWarning(Form("Events in bin0 from physics selection and local counter not consistent: %d - %d", cint1B, Int_t(triggeredEventsWith0MultWithBG))); } } else { AliInfo("Computing BG using CINT1A/B/C/E, ignoring intensities"); Int_t icol = fHistStatBin0->GetNbinsX(); Int_t cint1B = (Int_t) fHistStatBin0->GetBinContent(icol,1); Int_t cint1A = (Int_t) fHistStatBin0->GetBinContent(icol,2); Int_t cint1C = (Int_t) fHistStatBin0->GetBinContent(icol,3); Int_t cint1E = (Int_t) fHistStatBin0->GetBinContent(icol,4); bg = cint1A + cint1C-2*cint1E ; if (cint1B != triggeredEventsWith0MultWithBG) { AliWarning(Form("Events in bin0 from physics selection and local counter not consistent: %d - %d", cint1B, (Int_t)triggeredEventsWith0MultWithBG)); } } Double_t triggeredEventsWith0Mult = triggeredEventsWith0MultWithBG - bg; if(fVerbose > 0) AliInfo(Form("Measured events with vertex: %d",allEventsWithVertex)); Double_t bin0 = fHistVzData->Integral(0, fHistVzData->GetNbinsX()+1, 1, 1); if(fVerbose > 0) AliInfo(Form("Zero Bin, Meas: %2.2f, BG: %2.2f, Meas - BG: %2.2f", bin0, bg, triggeredEventsWith0Mult)); // This pointers are here so that I use the same names used in Jan Fiete's code (allows for faster/easier comparison) TH2* eTrig = histVzMCTrg; TH2* eTrigVtx = histVzMCRec; TH1* eTrigVtx_projx = eTrigVtx->ProjectionX("eTrigVtx_projx", 2, eTrigVtx->GetNbinsX()+1); // compute trigger and vertex efficiency TH2 * effVtxTrig = (TH2*) histVzMCRec->Clone("effVtxTrig"); effVtxTrig->Reset(); effVtxTrig->Divide(histVzMCRec,histVzMCGen,1,1,"B"); // Apply correction to data to get corrected events TH2 * correctedEvents = (TH2*) fHistVzData->Clone("correctedEvents"); correctedEvents->Divide(effVtxTrig); // TH1 * correctedEvents = fHistVzData->ProjectionX("eTrigVtx_projx", 2, eTrigVtx->GetNbinsX()+1); // loop over vertex bins for (Int_t i = 1; i <= fHistVzData->GetNbinsX(); i++) { Double_t alpha = (Double_t) hVzData->GetBinContent(i) / allEventsWithVertex; Double_t events = alpha * triggeredEventsWith0Mult; if (histVzMCRec->GetBinContent(i,1) == 0) continue; Double_t fZ = eTrigVtx_projx->Integral(0, eTrigVtx_projx->GetNbinsX()+1) / eTrigVtx_projx->GetBinContent(i) * eTrig->GetBinContent(i, 1) / eTrig->Integral(0, eTrig->GetNbinsX()+1, 1, 1); events *= fZ; // multiply with trigger correction Double_t correction = histVzMCGen->GetBinContent(i,1)/histVzMCTrg->GetBinContent(i,1); events *= correction; if (fVerbose > 1) Printf(" Bin %d, alpha is %.2f%%, fZ is %.3f, number of events with 0 mult.: %.2f (MC comparison: %.2f)", i, alpha * 100., fZ, events, histVzMCRec->GetBinContent(i,1)); correctedEvents->SetBinContent(i, 1, events); } // Integrate correctedEvents over full z range Double_t allEvents = correctedEvents->Integral(0, correctedEvents->GetNbinsX()+1,0, correctedEvents->GetNbinsY()+1); // Integrate correctedEvents over needed z range Double_t allEventsZrange = correctedEvents->Integral(correctedEvents->GetXaxis()->FindBin(-fZRange), correctedEvents->GetXaxis()->FindBin(fZRange), 0, correctedEvents->GetNbinsY()+1); if(fVerbose > 1) AliInfo(Form("Results in |Vz| < %3.3f",fZRange)); if(fVerbose > 1) AliInfo(Form(" Events in Bin0: %2.2f, With > 1 track: %2.2f, All corrected: %2.2f", correctedEvents->Integral(correctedEvents->GetXaxis()->FindBin(-fZRange), correctedEvents->GetXaxis()->FindBin(fZRange),1,1), correctedEvents->Integral(correctedEvents->GetXaxis()->FindBin(-fZRange), correctedEvents->GetXaxis()->FindBin(fZRange), 2,correctedEvents->GetNbinsX()+1), allEventsZrange)); if(fVerbose > 1) { Int_t nbin = histVzMCRec->GetNbinsX(); AliInfo(Form("Efficiency in the zero bin: %3.3f", histVzMCRec->Integral(0,nbin+1,1,1)/histVzMCGen->Integral(0,nbin+1,1,1) )); } AliInfo(Form("Number of collisions in full phase space: %2.2f", allEvents)); // effVtxTrig->Draw(); // new TCanvas(); // correctedEvents->Draw(); // FIXME: debug return allEvents; } Long64_t AliCollisionNormalization::Merge(TCollection* list) { // Merge a list of AliCollisionNormalization objects with this // (needed for PROOF). // Returns the number of merged objects (including this). if (!list) return 0; if (list->IsEmpty()) return 1; TIterator* iter = list->MakeIterator(); TObject* obj; // collections of all histograms const Int_t nHists = kNProcs*3+5; TList collections[nHists]; Int_t count = 0; while ((obj = iter->Next())) { AliCollisionNormalization* entry = dynamic_cast (obj); if (entry == 0) continue; Int_t ihist = -1; for(Int_t iproc = 0; iproc < kNProcs; iproc++){ if (entry->fHistVzMCGen[iproc] ) collections[++ihist].Add(entry->fHistVzMCGen[iproc] ); if (entry->fHistVzMCRec[iproc] ) collections[++ihist].Add(entry->fHistVzMCRec[iproc] ); if (entry->fHistVzMCTrg[iproc] ) collections[++ihist].Add(entry->fHistVzMCTrg[iproc] ); } if (entry->fHistVzData ) collections[++ihist].Add(entry->fHistVzData ); if (entry->fHistProcTypes ) collections[++ihist].Add(entry->fHistProcTypes ); if (entry->fHistStatBin0 ) collections[++ihist].Add(entry->fHistStatBin0 ); if (entry->fHistStat ) collections[++ihist].Add(entry->fHistStat ); count++; } Int_t ihist = -1; for(Int_t iproc = 0; iproc < kNProcs; iproc++){ if (fHistVzMCGen[iproc] ) fHistVzMCGen[iproc] ->Merge(&collections[++ihist]); if (fHistVzMCRec[iproc] ) fHistVzMCRec[iproc] ->Merge(&collections[++ihist]); if (fHistVzMCTrg[iproc] ) fHistVzMCTrg[iproc] ->Merge(&collections[++ihist]); } if (fHistVzData ) fHistVzData ->Merge(&collections[++ihist]); if (fHistProcTypes ) fHistProcTypes ->Merge(&collections[++ihist]); if (fHistStatBin0 ) fHistStatBin0 ->Merge(&collections[++ihist]); if (fHistStat ) fHistStat ->Merge(&collections[++ihist]); delete iter; return count+1; } void AliCollisionNormalization::FillVzMCGen(Float_t vz, Int_t ntrk, AliMCEvent * mcEvt) { // Fill MC gen histo and the process types statistics Int_t evtype = GetProcessType(mcEvt); // When I fill gen histos, I also fill statistics of process types (used to detemine ratios afterwards). fHistProcTypes->Fill(evtype); fHistVzMCGen[evtype]->Fill(vz,ntrk); } void AliCollisionNormalization::FillVzMCRec(Float_t vz, Int_t ntrk, AliMCEvent * mcEvt){ // Fill MC rec histo Int_t evtype = GetProcessType(mcEvt); fHistVzMCRec[evtype]->Fill(vz,ntrk); } void AliCollisionNormalization::FillVzMCTrg(Float_t vz, Int_t ntrk, AliMCEvent * mcEvt) { // Fill MC trg histo Int_t evtype = GetProcessType(mcEvt); fHistVzMCTrg[evtype]->Fill(vz,ntrk); } Int_t AliCollisionNormalization::GetProcessType(AliMCEvent * mcEvt) { // Determine if the event was generated with pythia or phojet and return the process type // Check if mcEvt is fine if (!mcEvt) { AliFatal("NULL mc event"); } // Determine if it was a pythia or phojet header, and return the correct process type AliGenPythiaEventHeader * headPy = 0; AliGenDPMjetEventHeader * headPho = 0; AliGenEventHeader * htmp = mcEvt->GenEventHeader(); if(!htmp) { AliFatal("Cannot Get MC Header!!"); } if( TString(htmp->IsA()->GetName()) == "AliGenPythiaEventHeader") { headPy = (AliGenPythiaEventHeader*) htmp; } else if (TString(htmp->IsA()->GetName()) == "AliGenDPMjetEventHeader") { headPho = (AliGenDPMjetEventHeader*) htmp; } else { AliError("Unknown header"); } // Determine process type if(headPy) { if(headPy->ProcessType() == 92 || headPy->ProcessType() == 93) { // single difractive return kProcSD; } else if (headPy->ProcessType() == 94) { // double diffractive return kProcDD; } else if(headPy->ProcessType() != 92 && headPy->ProcessType() != 93 && headPy->ProcessType() != 94) { // non difractive return kProcND; } } else if (headPho) { if(headPho->ProcessType() == 5 || headPho->ProcessType() == 6 ) { // single difractive return kProcSD; } else if (headPho->ProcessType() == 7) { // double diffractive return kProcDD; } else if(headPho->ProcessType() != 5 && headPho->ProcessType() != 6 && headPho->ProcessType() != 7 ) { // non difractive return kProcND; } } // no process type found? AliError(Form("Unknown header: %s", htmp->IsA()->GetName())); return kProcUnknown; } Double_t AliCollisionNormalization::GetProcessWeight(Int_t proc) { // Return a weight to be used when combining the MC histos to // compute efficiency, defined as the ratio XS in generator / XS // measured Float_t ref_SD, ref_DD, ref_ND, error_SD, error_DD, error_ND; GetRelativeFractions(fReferenceXS,ref_SD, ref_DD, ref_ND, error_SD, error_DD, error_ND); static Double_t total = fHistProcTypes->Integral(); if (fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcUnknown)) > 0) { AliError("There were unknown processes!!!"); } static Double_t SD = fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcSD))/total; static Double_t DD = fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcDD))/total; static Double_t ND = fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcND))/total; if (fVerbose > 2) { AliInfo(Form("Total MC evts: %f",total)); AliInfo(Form(" Frac SD %4.4f", SD)); AliInfo(Form(" Frac DD %4.4f", DD)); AliInfo(Form(" Frac ND %4.4f", ND)); } switch(proc) { case kProcSD: return ref_SD/SD; break; case kProcDD: return ref_DD/DD; break; case kProcND: return ref_ND/ND; break; default: AliError("Unknown process"); } return 0; } void AliCollisionNormalization::GetRelativeFractions(Int_t origin, Float_t& ref_SD, Float_t& ref_DD, Float_t& ref_ND, Float_t& error_SD, Float_t& error_DD, Float_t& error_ND) { // Returns fraction of XS (SD, ND, DD) and corresponding error // Stolen from Jan Fiete's drawSystematics macro // origin: // -1 = Pythia (test) // 0 = UA5 // 1 = Data 1.8 TeV // 2 = Tel-Aviv // 3 = Durham // Double_t epsilon = 0.0001; if(TMath::Abs(fEnergy-900) this should be 15.3 // DD: 0.122750 ref_SD = 0.224 * 0.153 / 0.189; error_SD = 0.023 * 0.224 / 0.189; ref_DD = 0.095; error_DD = 0.06; ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0; break; case 11: // UA5 AliInfo("UA5 x-sections hadron level definition for Phojet"); // Fractions in Phojet with UA5 cuts selection for SD // ND: 0.783573 // SD: 0.151601 --> this should be 15.3 // DD: 0.064827 ref_SD = 0.191 * 0.153 / 0.152; error_SD = 0.023 * 0.191 / 0.152; ref_DD = 0.095; error_DD = 0.06; ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0; break; case 2: // tel-aviv model AliInfo("Tel-aviv model x-sections"); ref_SD = 0.171; ref_DD = 0.094; ref_ND = 1 - ref_SD - ref_DD; break; case 3: // durham model AliInfo("Durham model x-sections"); ref_SD = 0.190; ref_DD = 0.125; ref_ND = 1 - ref_SD - ref_DD; break; default: AliFatal(Form("Unknown origin %d, Energy %f", origin, fEnergy)); } } else if(TMath::Abs(fEnergy-1800) this should be 15.9 // DD: 0.127701 ref_SD = 0.217 * 0.159 / 0.167; error_SD = 0.024 * 0.217 / 0.167; ref_DD = 0.075 * 1.43; error_DD = 0.02 * 1.43; ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0; break; case 21: // E710, 1.8 TeV AliInfo("E710 x-sections hadron level definition for Phojet"); // ND: 0.817462 // SD: 0.125506 --> this should be 15.9 // DD: 0.057032 ref_SD = 0.161 * 0.159 / 0.126; error_SD = 0.024 * 0.161 / 0.126; ref_DD = 0.075 * 1.43; error_DD = 0.02 * 1.43; ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0; break; case 1: // data 1.8 TeV AliInfo("??? x-sections"); ref_SD = 0.152; ref_DD = 0.092; ref_ND = 1 - ref_SD - ref_DD; break; default: AliFatal(Form("Unknown origin %d, Energy %f", origin, fEnergy)); } } else { AliFatal(Form("Unknown energy %f", fEnergy)); } }