[u/mrichter/AliRoot.git] / PWG0 / dNdEta / systematicsAnalysis.C

changeProductionCrossSections() {
  //
  // calculates dN/dEta from the standard correction and from
  // corrections that has been evaluated with different relative cross
  // sections (of single diff, double diff and non diff). The ratios
  // (standard to changed x-section) of the different dN/deta
  // distributions are saved to a file.
  //

  gSystem->Load("libPWG0base");

  // folder names in systematics root file
  const Char_t* folderNames[] = {"triggerBiasDD","triggerBiasSD","triggerBiasND","vertexRecoDD","vertexRecoSD","vertexRecoND"};

  const Char_t* changes[]  = {"pythia","ddmore","ddless","sdmore","sdless", "dmore", "dless"};

  Float_t scalesDD[] = {1.0, 1.5, 0.5, 1.0, 1.0, 1.5, 0.5};
  Float_t scalesSD[] = {1.0, 1.0, 1.0, 1.5, 0.5, 1.5, 0.5};

  // cross section from Pythia
  //Float_t sigmaND = 55.2;
  //Float_t sigmaDD = 9.78;
  //Float_t sigmaSD = 14.30;

  // getting data
  TFile* finCorr = TFile::Open("analysis_esd.root");
  dNdEtaAnalysis* analysis = new dNdEtaAnalysis("dndeta", "dndeta");
  analysis->LoadHistograms();

  // getting corrections
  AlidNdEtaCorrection* correctionStandard = new AlidNdEtaCorrection("dndeta_correction","dndeta_correction");

  TFile* finCorr = TFile::Open("correction_map.root");
  correctionStandard->LoadHistograms();

  // getting corrections for pure SD, DD and ND
  TFile* finSyst = TFile::Open("systematics.root");
  AlidNdEtaCorrection* correctionsSyst[7];
  for (Int_t i=0; i<6; i++) {
    correctionsSyst[i] = new AlidNdEtaCorrection(folderNames[i],folderNames[i]);
    correctionsSyst[i]->LoadHistograms();
  }

  TH1F* hRatios[21];
  for (Int_t j=0; j<3; j++) { // j = 0 (change vtx), j = 1 (change trg), j = 2 (change both)    

    for (Int_t i=0; i<7; i++) { // changes in cross section
      
      // temporary correction with changed cross sections
      AlidNdEtaCorrection* correctionTmp = (AlidNdEtaCorrection*)correctionStandard->Clone();;

      correctionTmp->Reset();

      if (j==0 || j==2) {
	//correctionTmp->GetVertexRecoCorrection()->Reset();
	//correctionTmp->GetVertexRecoCorrection()->Add(correctionsSyst[3]->GetVertexRecoCorrection(),scalesDD[i]);
	//correctionTmp->GetVertexRecoCorrection()->Add(correctionsSyst[4]->GetVertexRecoCorrection(),scalesSD[i]);
	//correctionTmp->GetVertexRecoCorrection()->Add(correctionsSyst[5]->GetVertexRecoCorrection(),1);

	correctionTmp->Add(correctionsSyst[3], scalesDD[i]);
	correctionTmp->Add(correctionsSyst[4], scalesSD[i]);
	correctionTmp->Add(correctionsSyst[5], 1);	
      }
      if (j==1 || j==2) {
	//correctionTmp->GetTriggerBiasCorrectionINEL()->Reset();
	//correctionTmp->GetTriggerBiasCorrectionINEL()->Add(correctionsSyst[0]->GetTriggerBiasCorrectionINEL(),scalesDD[i]);
	//correctionTmp->GetTriggerBiasCorrectionINEL()->Add(correctionsSyst[1]->GetTriggerBiasCorrectionINEL(),scalesSD[i]);
	//correctionTmp->GetTriggerBiasCorrectionINEL()->Add(correctionsSyst[2]->GetTriggerBiasCorrectionINEL(),1);

	correctionTmp->Add(correctionsSyst[0], scalesDD[i]);
	correctionTmp->Add(correctionsSyst[1], scalesSD[i]);
	correctionTmp->Add(correctionsSyst[2], 1);
      }
      correctionTmp->Finish();

      dNdEtaAnalysis* analysisTmp = (dNdEtaAnalysis*)analysis->Clone();
      analysisTmp->Finish(correctionTmp, 0.3, AlidNdEtaCorrection::kINEL);

      Int_t counter = i + j*7;
      
      hRatios[counter] = (TH1F*)analysisTmp->GetdNdEtaHistogram(2)->Clone();
      
      TString name("ratio");
      if (j==0) name.Append("_vetexReco_");
      if (j==1) name.Append("_triggerBias_");
      if (j==2) name.Append("_vertexReco_triggerBias_");
      name.Append(changes[i]);
      hRatios[counter]->SetName(name.Data());
      name.Append(Form(" (DD #times %0.1f, SD #times %0.1f)",scalesDD[i],scalesSD[i]));
      hRatios[counter]->SetTitle(name.Data());
      hRatios[counter]->SetYTitle("ratio (Pythia x-sections)/(changed x-sections)");
      
      delete analysisTmp;
      delete correctionTmp;
    }
  }

  for (Int_t i=1; i<21; i++) 
    hRatios[i]->Divide(hRatios[0],hRatios[i],1,1);
  hRatios[0]->Divide(hRatios[0],hRatios[0],1,1);

  TFile* fout = new TFile("systematics_xsections.root","RECREATE");
  
  for (Int_t i=0; i<21; i++) 
    hRatios[i]->Write();
  
  fout->Write();
  fout->Close();
}

changeParticleComposition() {
  //
  // calculates dN/dEta from the standard correction and from
  // corrections that has been evaluated with different relative
  // abundancies of kaons and protons and save the ratios in a file

  gSystem->Load("libPWG0base");

  const Char_t* folderNames[] = {"correction_0","correction_1","correction_2","correction_3"};
  Float_t  scalesPi[] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
  Float_t  scalesK[]  = {0.5, 1.5, 1.0, 1.0, 0.5, 1.5};
  Float_t  scalesP[]  = {1.0, 1.0, 0.5, 1.5, 0.5, 1.5};

  // getting data
  TFile* finCorr = TFile::Open("analysis_esd.root");
  dNdEtaAnalysis* analysis = new dNdEtaAnalysis("dndeta", "dndeta");
  analysis->LoadHistograms();

  // getting corrections for pi, K, p and other particles
  TFile* finSyst = TFile::Open("systematics.root");
  AlidNdEtaCorrection* correctionsSyst[4];
  for (Int_t i=0; i<4; i++) {
    correctionsSyst[i] = new AlidNdEtaCorrection(folderNames[i],folderNames[i]);
    correctionsSyst[i]->LoadHistograms();
  }
  AlidNdEtaCorrection* correctionAll = (AlidNdEtaCorrection*) correctionsSyst[0]->Clone();
  correctionAll->Reset();
  correctionAll->Add(correctionsSyst[0]);
  correctionAll->Add(correctionsSyst[1]);
  correctionAll->Add(correctionsSyst[2]);
  correctionAll->Add(correctionsSyst[3]);
  correctionAll->Finish();
  analysis->Finish(correctionAll, 0.3, AlidNdEtaCorrection::kINEL);

  TH1F* hRatios[6];
  for (Int_t i=0; i<6; i++) { 
    // temporary correction with changed particle composistion
    AlidNdEtaCorrection* correctionTmp = (AlidNdEtaCorrection*)correctionAll->Clone();
    
    correctionTmp->Reset();
    
    correctionTmp->Add(correctionsSyst[0],scalesPi[i]);    
    correctionTmp->Add(correctionsSyst[1],scalesK[i]);    
    correctionTmp->Add(correctionsSyst[2],scalesP[i]);    
    correctionTmp->Add(correctionsSyst[3],1);    
    
    correctionTmp->Finish();

    dNdEtaAnalysis* analysisTmp = (dNdEtaAnalysis*)analysis->Clone();
    analysisTmp->Finish(correctionTmp, 0.3, AlidNdEtaCorrection::kINEL);
    
    hRatios[i] = (TH1F*)analysisTmp->GetdNdEtaHistogram(2)->Clone();
    hRatios[i]->Divide((TH1F*)analysis->GetdNdEtaHistogram(2)->Clone(),hRatios[i],1,1,"B"); 

    TString name(Form("ratio_%d",i));
    hRatios[i]->SetName(name.Data());
    name = TString(Form("#pi #times %0.1f, K #times %0.1f, p #times %0.1f",scalesPi[i],scalesK[i],scalesP[i]));
    hRatios[i]->SetTitle(name.Data());
    hRatios[i]->SetYTitle("ratio (standard/changed compositions)");    
  }

  TFile* fout = new TFile("systematics_composition.root","RECREATE");
  
  for (Int_t i=0; i<6; i++) 
    hRatios[i]->Write();
  
  fout->Write();
  fout->Close();


}
Commit	Line	Data
b541c64c	1	changeProductionCrossSections() {
	2	//
	3	// calculates dN/dEta from the standard correction and from
	4	// corrections that has been evaluated with different relative cross
	5	// sections (of single diff, double diff and non diff). The ratios
	6	// (standard to changed x-section) of the different dN/deta
	7	// distributions are saved to a file.
	8	//
	9
	10	gSystem->Load("libPWG0base");
	11
	12	// folder names in systematics root file
	13	const Char_t* folderNames[] = {"triggerBiasDD","triggerBiasSD","triggerBiasND","vertexRecoDD","vertexRecoSD","vertexRecoND"};
	14
	15	const Char_t* changes[] = {"pythia","ddmore","ddless","sdmore","sdless", "dmore", "dless"};
	16
	17	Float_t scalesDD[] = {1.0, 1.5, 0.5, 1.0, 1.0, 1.5, 0.5};
	18	Float_t scalesSD[] = {1.0, 1.0, 1.0, 1.5, 0.5, 1.5, 0.5};
	19
	20	// cross section from Pythia
	21	//Float_t sigmaND = 55.2;
	22	//Float_t sigmaDD = 9.78;
	23	//Float_t sigmaSD = 14.30;
	24
	25	// getting data
	26	TFile* finCorr = TFile::Open("analysis_esd.root");
	27	dNdEtaAnalysis* analysis = new dNdEtaAnalysis("dndeta", "dndeta");
	28	analysis->LoadHistograms();
	29
	30	// getting corrections
	31	AlidNdEtaCorrection* correctionStandard = new AlidNdEtaCorrection("dndeta_correction","dndeta_correction");
	32
	33	TFile* finCorr = TFile::Open("correction_map.root");
	34	correctionStandard->LoadHistograms();
	35
	36	// getting corrections for pure SD, DD and ND
	37	TFile* finSyst = TFile::Open("systematics.root");
	38	AlidNdEtaCorrection* correctionsSyst[7];
	39	for (Int_t i=0; i<6; i++) {
	40	correctionsSyst[i] = new AlidNdEtaCorrection(folderNames[i],folderNames[i]);
	41	correctionsSyst[i]->LoadHistograms();
	42	}
	43
	44	TH1F* hRatios[21];
	45	for (Int_t j=0; j<3; j++) { // j = 0 (change vtx), j = 1 (change trg), j = 2 (change both)
	46
	47	for (Int_t i=0; i<7; i++) { // changes in cross section
	48
	49	// temporary correction with changed cross sections
	50	AlidNdEtaCorrection* correctionTmp = (AlidNdEtaCorrection*)correctionStandard->Clone();;
	51
	52	correctionTmp->Reset();
	53
	54	if (j==0 \|\| j==2) {
	55	//correctionTmp->GetVertexRecoCorrection()->Reset();
	56	//correctionTmp->GetVertexRecoCorrection()->Add(correctionsSyst[3]->GetVertexRecoCorrection(),scalesDD[i]);
	57	//correctionTmp->GetVertexRecoCorrection()->Add(correctionsSyst[4]->GetVertexRecoCorrection(),scalesSD[i]);
	58	//correctionTmp->GetVertexRecoCorrection()->Add(correctionsSyst[5]->GetVertexRecoCorrection(),1);
	59
	60	correctionTmp->Add(correctionsSyst[3], scalesDD[i]);
	61	correctionTmp->Add(correctionsSyst[4], scalesSD[i]);
	62	correctionTmp->Add(correctionsSyst[5], 1);
	63	}
	64	if (j==1 \|\| j==2) {
65	//correctionTmp->GetTriggerBiasCorrectionINEL()->Reset();
66	//correctionTmp->GetTriggerBiasCorrectionINEL()->Add(correctionsSyst[0]->GetTriggerBiasCorrectionINEL(),scalesDD[i]);
67	//correctionTmp->GetTriggerBiasCorrectionINEL()->Add(correctionsSyst[1]->GetTriggerBiasCorrectionINEL(),scalesSD[i]);
68	//correctionTmp->GetTriggerBiasCorrectionINEL()->Add(correctionsSyst[2]->GetTriggerBiasCorrectionINEL(),1);
69
70	correctionTmp->Add(correctionsSyst[0], scalesDD[i]);
71	correctionTmp->Add(correctionsSyst[1], scalesSD[i]);
72	correctionTmp->Add(correctionsSyst[2], 1);
73	}
74	correctionTmp->Finish();
75
76	dNdEtaAnalysis* analysisTmp = (dNdEtaAnalysis*)analysis->Clone();
77	analysisTmp->Finish(correctionTmp, 0.3, AlidNdEtaCorrection::kINEL);
78
79	Int_t counter = i + j*7;
80
81	hRatios[counter] = (TH1F*)analysisTmp->GetdNdEtaHistogram(2)->Clone();
82
83	TString name("ratio");
84	if (j==0) name.Append("_vetexReco_");
85	if (j==1) name.Append("_triggerBias_");
86	if (j==2) name.Append("_vertexReco_triggerBias_");
87	name.Append(changes[i]);
88	hRatios[counter]->SetName(name.Data());
89	name.Append(Form(" (DD #times %0.1f, SD #times %0.1f)",scalesDD[i],scalesSD[i]));
90	hRatios[counter]->SetTitle(name.Data());
91	hRatios[counter]->SetYTitle("ratio (Pythia x-sections)/(changed x-sections)");
92
93	delete analysisTmp;
94	delete correctionTmp;
95	}
96	}
97
98	for (Int_t i=1; i<21; i++)
99	hRatios[i]->Divide(hRatios[0],hRatios[i],1,1);
100	hRatios[0]->Divide(hRatios[0],hRatios[0],1,1);
101
102	TFile* fout = new TFile("systematics_xsections.root","RECREATE");
103
104	for (Int_t i=0; i<21; i++)
105	hRatios[i]->Write();
106
107	fout->Write();
108	fout->Close();
109	}
110
111	changeParticleComposition() {
112	//
113	// calculates dN/dEta from the standard correction and from
114	// corrections that has been evaluated with different relative
115	// abundancies of kaons and protons and save the ratios in a file
116
117	gSystem->Load("libPWG0base");
118
119	const Char_t* folderNames[] = {"correction_0","correction_1","correction_2","correction_3"};
120	Float_t scalesPi[] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
121	Float_t scalesK[] = {0.5, 1.5, 1.0, 1.0, 0.5, 1.5};
122	Float_t scalesP[] = {1.0, 1.0, 0.5, 1.5, 0.5, 1.5};
123
124	// getting data
125	TFile* finCorr = TFile::Open("analysis_esd.root");
126	dNdEtaAnalysis* analysis = new dNdEtaAnalysis("dndeta", "dndeta");
127	analysis->LoadHistograms();
128
129	// getting corrections for pi, K, p and other particles
130	TFile* finSyst = TFile::Open("systematics.root");
131	AlidNdEtaCorrection* correctionsSyst[4];
132	for (Int_t i=0; i<4; i++) {
133	correctionsSyst[i] = new AlidNdEtaCorrection(folderNames[i],folderNames[i]);
134	correctionsSyst[i]->LoadHistograms();
135	}
136	AlidNdEtaCorrection* correctionAll = (AlidNdEtaCorrection*) correctionsSyst[0]->Clone();
137	correctionAll->Reset();
138	correctionAll->Add(correctionsSyst[0]);
139	correctionAll->Add(correctionsSyst[1]);
140	correctionAll->Add(correctionsSyst[2]);
141	correctionAll->Add(correctionsSyst[3]);
142	correctionAll->Finish();
143	analysis->Finish(correctionAll, 0.3, AlidNdEtaCorrection::kINEL);
144
145	TH1F* hRatios[6];
146	for (Int_t i=0; i<6; i++) {
147	// temporary correction with changed particle composistion
148	AlidNdEtaCorrection* correctionTmp = (AlidNdEtaCorrection*)correctionAll->Clone();
149
150	correctionTmp->Reset();
151
152	correctionTmp->Add(correctionsSyst[0],scalesPi[i]);
153	correctionTmp->Add(correctionsSyst[1],scalesK[i]);
154	correctionTmp->Add(correctionsSyst[2],scalesP[i]);
155	correctionTmp->Add(correctionsSyst[3],1);
156
157	correctionTmp->Finish();
158
159	dNdEtaAnalysis* analysisTmp = (dNdEtaAnalysis*)analysis->Clone();
160	analysisTmp->Finish(correctionTmp, 0.3, AlidNdEtaCorrection::kINEL);
161
162	hRatios[i] = (TH1F*)analysisTmp->GetdNdEtaHistogram(2)->Clone();
163	hRatios[i]->Divide((TH1F*)analysis->GetdNdEtaHistogram(2)->Clone(),hRatios[i],1,1,"B");
164
165	TString name(Form("ratio_%d",i));
166	hRatios[i]->SetName(name.Data());
167	name = TString(Form("#pi #times %0.1f, K #times %0.1f, p #times %0.1f",scalesPi[i],scalesK[i],scalesP[i]));
168	hRatios[i]->SetTitle(name.Data());
169	hRatios[i]->SetYTitle("ratio (standard/changed compositions)");
170	}
171
172	TFile* fout = new TFile("systematics_composition.root","RECREATE");
173
174	for (Int_t i=0; i<6; i++)
175	hRatios[i]->Write();
176
177	fout->Write();
178	fout->Close();
179
180
181	}