[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / corrs / ExtractAcceptance.C

//_____________________________________________________________________
/** 
 * 
 * 
 * @param d 
 * @param r 
 * @param vz 
 * @param nDead 
 * 
 * @return 
 *
 * @ingroup pwglf_forward_scripts_corr
 */
TH2D* MakeOneRing(UShort_t d, Char_t r, Double_t vz, Int_t& nDead)
{
  AliFMDGeometry*   geom = AliFMDGeometry::Instance();
  AliFMDParameters* pars = AliFMDParameters::Instance();

  UShort_t nS = (r == 'I' || r == 'i' ?  20 :  40);
  UShort_t nT = (r == 'I' || r == 'i' ? 512 : 256);
  
  // Make our two histograms 
  TH2D* hAll = new TH2D("all","All",200,-4,6,nS,0,2*TMath::Pi());
  hAll->SetXTitle("#eta");
  hAll->SetYTitle("#phi");
  hAll->Sumw2();
  hAll->SetDirectory(0);
  TH2D* hOK  = static_cast<TH2D*>(hAll->Clone());
  hOK->SetDirectory(0);
  
  // Loop over all sectors and strips in this ring 
  Int_t nOK  = 0;
  Int_t nAll = 0;
  Int_t nPhi = hAll->GetNbinsY();
  for (UShort_t s = 0; s < nS; s++) { 
    for (UShort_t t = 0; t < nT; t++) { 
      // Get eta,phi by quering the geometry (first for (x,y,z), then
      // correcting for the vertex position, and then calculating 
      // (eta, phi))
      Double_t x, y, z;
      geom->Detector2XYZ(d, r, s, t, x, y, z);
      z -= vz;
      Double_t q, eta, phi, theta;
      AliFMDGeometry::XYZ2REtaPhiTheta(x, y, z, q, eta, phi, theta);
      if (phi < 0) phi += 2*TMath::Pi();

      // Check if this is a dead channel or not 
      Bool_t isDead = pars->IsDead(d, r, s, t);

      // Special check for FMD2i - upper part of sectors 16/17 have 
      // have anomalous gains and/or noise - common sources are 
      // power regulartors for bias currents and the like 
      Int_t VA = t/128;
      if(d==2 && r=='I' && VA>1 && (s==16 || s==17)) isDead =true;

      // Find the eta bin number and corresponding overflow bin
      Int_t etaBin = hAll->GetXaxis()->FindBin(eta);
      Int_t ovrBin = hAll->GetBin(etaBin, nPhi+1); 
      
      // Increment all histogram 
      hAll->Fill(eta, phi);
      hAll->AddBinContent(ovrBin);
      nAll++;

      // If not dead, increment OK histogram 
      if (!isDead) {
	hOK->Fill(eta, phi);
	hOK->AddBinContent(ovrBin);
	nOK++;
      }
      else         nDead++;
    }
  }
  // Divide out the efficiency. 
  // Note, that the overflow bins along eta now contains the ratio 
  // nOK/nAll Strips for a given eta bin. 
  hOK->Divide(hOK,hAll,1,1,"B");

  // Invert overflow bin
  for (Int_t etaBin = 1; etaBin <= hOK->GetNbinsX(); etaBin++) { 
    Double_t ovr    = hOK->GetBinContent(etaBin, nPhi+1);
    Double_t novr   = (ovr < 1e-12 ? 0 : 1./ovr);
    hOK->SetBinContent(etaBin, nPhi+1, novr);
    if (ovr > 0 && ovr != 1)
      Info("", "Setting overflow bin (%3d,%3d) to 1/%f=%f", etaBin, nPhi+1, 
	   ovr, hOK->GetBinContent(etaBin, nPhi+1));
  }
  // Clean up
  delete hAll;

  Info("ExtractAcceptances","Made correction for FMD%d%c at vz=%f - "
       "%d strips out of %d OK (w/overflow)", d, r, vz, nOK, nAll);

  // Return result 
  return hOK;
}

//_____________________________________________________________________
/** 
 * 
 * 
 * @param runNo 
 * @param system 
 * @param energy 
 * @param field 
 * @param nVtxBins 
 * @param vtxLow 
 * @param vtxHigh 
 *
 * @ingroup pwglf_forward_scripts_corr
 */
void ExtractAcceptance(Int_t   runNo=121526, 
		       Int_t   system = 1,
		       Float_t energy = 900,
		       Float_t field  = 5,
		       Int_t   nVtxBins=10, 
		       Float_t vtxLow=-10, 
		       Float_t vtxHigh=10)
{  
  gSystem->Load("libANALYSIS");
  gSystem->Load("libANALYSISalice");
  gSystem->Load("libPWGLFforward2");
  
  // Float_t delta = (vtxHigh - vtxLow) / (Float_t)nVtxBins;
  
  Bool_t kGridOnline = kTRUE; 
  if(!(TGrid::Connect("alien://",0,0,"t")))
    kGridOnline = kFALSE;
  
  // --- Initialisations ------------------------------------------
  //Set up CDB manager
  Info("ExtractAcceptances","Setting up OCDB");
  
  AliCDBManager* cdb = AliCDBManager::Instance();
  if(kGridOnline)
    cdb->SetDefaultStorage("alien://Folder=/alice/data/2012/OCDB");
  else
    cdb->SetDefaultStorage("local://$(ALICE_ROOT)/OCDB");
  cdb->SetRun(runNo);
  
  // Get the geometry 
  Info("ExtractAcceptances","Loading geometry");
  AliGeomManager::LoadGeometry();

  // Get an initialize parameters 
  Info("ExtractAcceptances","Intialising parameters");
  AliFMDParameters* pars = AliFMDParameters::Instance();
  pars->Init();

  // Get an initialise geometry 
  Info("ExtractAcceptances","Initialising geomtry");
  AliFMDGeometry* geom = AliFMDGeometry::Instance();
  geom->Init();
  geom->InitTransformations();

  // --- Output object -----------------------------------------------
  // Make our correction object 
  AliFMDCorrAcceptance* corr = new AliFMDCorrAcceptance();
  corr->SetVertexAxis(nVtxBins, vtxLow, vtxHigh);

  // --- Loop over verticies and rings -------------------------------
  Int_t nDead = 0;
  Float_t dV = (vtxHigh - vtxLow) / nVtxBins;
  for (Double_t v = vtxLow+dV/2; v < vtxHigh; v += dV) { 
    for(UShort_t d = 1; d <= 3;d++) { 
      UShort_t nR = (d == 1 ? 1 : 2);
      for (UShort_t q = 0; q < nR; q++) { 
	Char_t   r  = (q == 0 ? 'I' : 'O');

	// Delegate to other function 
	TH2D* ratio = MakeOneRing(d, r, v, nDead);
	if (!ratio) continue;

	// Set the correction 
	corr->SetCorrection(d, r, v, ratio);
      }
    }
  }

  // Write to a file 
  Info("ExtractAcceptances","Writing to disk");
  AliForwardCorrectionManager& cm = AliForwardCorrectionManager::Instance();
  TString fname = cm.GetFileName(AliForwardCorrectionManager::kAcceptance, 
				 system, energy, field, false);
  TFile* out = TFile::Open(fname.Data(), "RECREATE");
  corr->SetHasOverflow();
  corr->Write(cm.GetObjectName(AliForwardCorrectionManager::kAcceptance));
  out->Write();
  out->Close();

  std::ofstream f("Upload.C");
  if (!f) { 
    Error("ExtractELoss", "Failed to open Upload.C");
    return;
  }
  f << "// Generated by ExtractAcceptance.C\n"
    << "void Upload(const TUrl& url)\n"
    << "{\n"
    << "  if (TString(\"alien\").EqualTo(url.GetProtocol())) {\n"
    << "    if (!TGrid::Connect(\"alien://\")) {\n"
    << "      Error(\"Upload\", \"Failed to connect to AliEn\");\n"
    << "      return;\n"
    << "    }\n"
    << "  }\n\n";
  
  cm.SetPrefix("");
  TString fef(cm.GetFileName(AliForwardCorrectionManager::kAcceptance, 
			      system, energy, field, false));
  TString fep(cm.GetFilePath(AliForwardCorrectionManager::kAcceptance, 
			     system, energy, field, false));
  f << "  TString src  = \"" << fef << "\";\n"
    << "  TString dest = \"" << fep << "\";\n"
    << "  TString out; out.Form(\"%s%s\",url.GetUrl(),dest.Data());\n\n"
    << "  TString dir(gSystem->DirName(out));\n"
    << "  if (gSystem->AccessPathName(dir)) {\n"
    << "    if (gSystem->mkdir(dir, true) < 0) {\n"
    << "      Warning(\"Upload\",\"Failed to make directory %s\","
    << "              dir.Data());\n"
    << "      return;\n"
    << "    }\n"
    << "  }\n"
    << "  if (!TFile::Cp(src,out)) \n"
    << "    Warning(\"Upload\",\"Failed to upload %s -> %s\",\n"
    << "            src.Data(), out.Data());\n"
    << "}\n"
    << "// EOF"
    << std::endl;
  f.close();
  
  Info("ExtracAcceptance", 
       "Run generated Upload.C(DEST) script to copy files in place");
}

//
// EOF
//
Commit	Line	Data
c47f3ae9	1	//_____________________________________________________________________
	2	/**
	3	*
	4	*
	5	* @param d
	6	* @param r
	7	* @param vz
	8	* @param nDead
	9	*
	10	* @return
	11	*
	12	* @ingroup pwglf_forward_scripts_corr
	13	*/
	14	TH2D* MakeOneRing(UShort_t d, Char_t r, Double_t vz, Int_t& nDead)
	15	{
	16	AliFMDGeometry* geom = AliFMDGeometry::Instance();
	17	AliFMDParameters* pars = AliFMDParameters::Instance();
	18
	19	UShort_t nS = (r == 'I' \|\| r == 'i' ? 20 : 40);
	20	UShort_t nT = (r == 'I' \|\| r == 'i' ? 512 : 256);
	21
	22	// Make our two histograms
	23	TH2D* hAll = new TH2D("all","All",200,-4,6,nS,0,2*TMath::Pi());
	24	hAll->SetXTitle("#eta");
	25	hAll->SetYTitle("#phi");
	26	hAll->Sumw2();
	27	hAll->SetDirectory(0);
	28	TH2D* hOK = static_cast<TH2D*>(hAll->Clone());
	29	hOK->SetDirectory(0);
	30
	31	// Loop over all sectors and strips in this ring
	32	Int_t nOK = 0;
	33	Int_t nAll = 0;
	34	Int_t nPhi = hAll->GetNbinsY();
	35	for (UShort_t s = 0; s < nS; s++) {
	36	for (UShort_t t = 0; t < nT; t++) {
	37	// Get eta,phi by quering the geometry (first for (x,y,z), then
	38	// correcting for the vertex position, and then calculating
	39	// (eta, phi))
	40	Double_t x, y, z;
	41	geom->Detector2XYZ(d, r, s, t, x, y, z);
	42	z -= vz;
	43	Double_t q, eta, phi, theta;
	44	AliFMDGeometry::XYZ2REtaPhiTheta(x, y, z, q, eta, phi, theta);
	45	if (phi < 0) phi += 2*TMath::Pi();
	46
	47	// Check if this is a dead channel or not
	48	Bool_t isDead = pars->IsDead(d, r, s, t);
	49
	50	// Special check for FMD2i - upper part of sectors 16/17 have
	51	// have anomalous gains and/or noise - common sources are
	52	// power regulartors for bias currents and the like
	53	Int_t VA = t/128;
	54	if(d==2 && r=='I' && VA>1 && (s==16 \|\| s==17)) isDead =true;
	55
	56	// Find the eta bin number and corresponding overflow bin
	57	Int_t etaBin = hAll->GetXaxis()->FindBin(eta);
	58	Int_t ovrBin = hAll->GetBin(etaBin, nPhi+1);
	59
	60	// Increment all histogram
	61	hAll->Fill(eta, phi);
	62	hAll->AddBinContent(ovrBin);
	63	nAll++;
	64
65	// If not dead, increment OK histogram
66	if (!isDead) {
67	hOK->Fill(eta, phi);
68	hOK->AddBinContent(ovrBin);
69	nOK++;
70	}
71	else nDead++;
72	}
73	}
74	// Divide out the efficiency.
75	// Note, that the overflow bins along eta now contains the ratio
76	// nOK/nAll Strips for a given eta bin.
77	hOK->Divide(hOK,hAll,1,1,"B");
78
79	// Invert overflow bin
80	for (Int_t etaBin = 1; etaBin <= hOK->GetNbinsX(); etaBin++) {
81	Double_t ovr = hOK->GetBinContent(etaBin, nPhi+1);
82	Double_t novr = (ovr < 1e-12 ? 0 : 1./ovr);
83	hOK->SetBinContent(etaBin, nPhi+1, novr);
84	if (ovr > 0 && ovr != 1)
85	Info("", "Setting overflow bin (%3d,%3d) to 1/%f=%f", etaBin, nPhi+1,
86	ovr, hOK->GetBinContent(etaBin, nPhi+1));
87	}
88	// Clean up
89	delete hAll;
90
91	Info("ExtractAcceptances","Made correction for FMD%d%c at vz=%f - "
92	"%d strips out of %d OK (w/overflow)", d, r, vz, nOK, nAll);
93
94	// Return result
95	return hOK;
96	}
97
98	//_____________________________________________________________________
99	/**
100	*
101	*
102	* @param runNo
103	* @param system
104	* @param energy
105	* @param field
106	* @param nVtxBins
107	* @param vtxLow
108	* @param vtxHigh
109	*
110	* @ingroup pwglf_forward_scripts_corr
111	*/
112	void ExtractAcceptance(Int_t runNo=121526,
113	Int_t system = 1,
114	Float_t energy = 900,
115	Float_t field = 5,
116	Int_t nVtxBins=10,
117	Float_t vtxLow=-10,
118	Float_t vtxHigh=10)
119	{
120	gSystem->Load("libANALYSIS");
121	gSystem->Load("libANALYSISalice");
122	gSystem->Load("libPWGLFforward2");
123
124	// Float_t delta = (vtxHigh - vtxLow) / (Float_t)nVtxBins;
125
126	Bool_t kGridOnline = kTRUE;
127	if(!(TGrid::Connect("alien://",0,0,"t")))
128	kGridOnline = kFALSE;
129
130	// --- Initialisations ------------------------------------------
131	//Set up CDB manager
132	Info("ExtractAcceptances","Setting up OCDB");
133
134	AliCDBManager* cdb = AliCDBManager::Instance();
135	if(kGridOnline)
136	cdb->SetDefaultStorage("alien://Folder=/alice/data/2012/OCDB");
137	else
138	cdb->SetDefaultStorage("local://$(ALICE_ROOT)/OCDB");
139	cdb->SetRun(runNo);
140
141	// Get the geometry
142	Info("ExtractAcceptances","Loading geometry");
143	AliGeomManager::LoadGeometry();
144
145	// Get an initialize parameters
146	Info("ExtractAcceptances","Intialising parameters");
147	AliFMDParameters* pars = AliFMDParameters::Instance();
148	pars->Init();
149
150	// Get an initialise geometry
151	Info("ExtractAcceptances","Initialising geomtry");
152	AliFMDGeometry* geom = AliFMDGeometry::Instance();
153	geom->Init();
154	geom->InitTransformations();
155
156	// --- Output object -----------------------------------------------
157	// Make our correction object
158	AliFMDCorrAcceptance* corr = new AliFMDCorrAcceptance();
159	corr->SetVertexAxis(nVtxBins, vtxLow, vtxHigh);
160
161	// --- Loop over verticies and rings -------------------------------
162	Int_t nDead = 0;
163	Float_t dV = (vtxHigh - vtxLow) / nVtxBins;
164	for (Double_t v = vtxLow+dV/2; v < vtxHigh; v += dV) {
165	for(UShort_t d = 1; d <= 3;d++) {
166	UShort_t nR = (d == 1 ? 1 : 2);
167	for (UShort_t q = 0; q < nR; q++) {
168	Char_t r = (q == 0 ? 'I' : 'O');
169
170	// Delegate to other function
171	TH2D* ratio = MakeOneRing(d, r, v, nDead);
172	if (!ratio) continue;
173
174	// Set the correction
175	corr->SetCorrection(d, r, v, ratio);
176	}
177	}
178	}
179
180	// Write to a file
181	Info("ExtractAcceptances","Writing to disk");
182	AliForwardCorrectionManager& cm = AliForwardCorrectionManager::Instance();
183	TString fname = cm.GetFileName(AliForwardCorrectionManager::kAcceptance,
184	system, energy, field, false);
185	TFile* out = TFile::Open(fname.Data(), "RECREATE");
186	corr->SetHasOverflow();
187	corr->Write(cm.GetObjectName(AliForwardCorrectionManager::kAcceptance));
188	out->Write();
189	out->Close();
190
191	std::ofstream f("Upload.C");
192	if (!f) {
193	Error("ExtractELoss", "Failed to open Upload.C");
194	return;
195	}
196	f << "// Generated by ExtractAcceptance.C\n"
197	<< "void Upload(const TUrl& url)\n"
198	<< "{\n"
199	<< " if (TString(\"alien\").EqualTo(url.GetProtocol())) {\n"
200	<< " if (!TGrid::Connect(\"alien://\")) {\n"
201	<< " Error(\"Upload\", \"Failed to connect to AliEn\");\n"
202	<< " return;\n"
203	<< " }\n"
204	<< " }\n\n";
205
206	cm.SetPrefix("");
207	TString fef(cm.GetFileName(AliForwardCorrectionManager::kAcceptance,
208	system, energy, field, false));
209	TString fep(cm.GetFilePath(AliForwardCorrectionManager::kAcceptance,
210	system, energy, field, false));
211	f << " TString src = \"" << fef << "\";\n"
212	<< " TString dest = \"" << fep << "\";\n"
213	<< " TString out; out.Form(\"%s%s\",url.GetUrl(),dest.Data());\n\n"
214	<< " TString dir(gSystem->DirName(out));\n"
215	<< " if (gSystem->AccessPathName(dir)) {\n"
216	<< " if (gSystem->mkdir(dir, true) < 0) {\n"
217	<< " Warning(\"Upload\",\"Failed to make directory %s\","
218	<< " dir.Data());\n"
219	<< " return;\n"
220	<< " }\n"
221	<< " }\n"
222	<< " if (!TFile::Cp(src,out)) \n"
223	<< " Warning(\"Upload\",\"Failed to upload %s -> %s\",\n"
224	<< " src.Data(), out.Data());\n"
225	<< "}\n"
226	<< "// EOF"
227	<< std::endl;
228	f.close();
229
230	Info("ExtracAcceptance",
231	"Run generated Upload.C(DEST) script to copy files in place");
232	}
233
234	//
235	// EOF
236	//
237