[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / tests / TestAcc.C

#include <TMath.h>
#include <TGraph.h>
#include <TMultiGraph.h>
#include <TMarker.h>
#include <TLine.h>
#include <TArc.h>
#include <TCanvas.h>
#include <TH2F.h>
#include <THStack.h>
/**
 * @defgroup pwg2_forward_scripts_tests Test scripts
 * 
 * 
 *
 * @ingroup pwg2_forward_scripts
 */
/** 
 * 
 * 
 * @param r 
 * @param t 
 * @param oldm 
 * @param newm 
 *
 * @ingroup pwg2_forward_scripts_tests
 */
//_____________________________________________________________________
void AcceptanceCorrection(Char_t r, UShort_t t, Float_t& oldm, Float_t& newm)
{
  // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  
  //AliFMDRing fmdring(ring);
  //fmdring.Init();
  // Float_t   rad       = pars->GetMaxR(r)-pars->GetMinR(r);
  // Float_t   slen      = pars->GetStripLength(r,t);
  // Float_t   sblen     = pars->GetBaseStripLength(r,t);
  const Double_t ic1[] = { 4.9895, 15.3560 };
  const Double_t ic2[] = { 1.8007, 17.2000 };
  const Double_t oc1[] = { 4.2231, 26.6638 };
  const Double_t oc2[] = { 1.8357, 27.9500 };

  Double_t  minR      = (r == 'I' ?  4.5213 : 15.4);
  Double_t  maxR      = (r == 'I' ? 17.2    : 28.0);
  Double_t  rad       = maxR - minR;
  
  Int_t     nStrips   = (r == 'I' ? 512 : 256);
  Int_t     nSec      = (r == 'I' ?  20 :  40);
  Float_t   segment   = rad / nStrips;
  Float_t   basearc   = 2 * TMath::Pi() / (0.5 * nSec);
  Float_t   radius    = minR + t * segment;
  Float_t   baselen   = 0.5 * basearc * radius;

  const Double_t* c1  = (r == 'I' ? ic1 : oc1);
  const Double_t* c2  = (r == 'I' ? ic2 : oc2);

  // If the radius of the strip is smaller than the radius corresponding 
  // to the first corner we have a full strip length 
  Float_t   cr        = TMath::Sqrt(c1[0]*c1[0]+c1[1]*c1[1]);
  if (radius <= cr) newm = 1;
  else {
    // Next, we should find the end-point of the strip - that is, 
    // the coordinates where the line from c1 to c2 intersects a circle 
    // with radius given by the strip. 
    // (See http://mathworld.wolfram.com/Circle-LineIntersection.html)
    Float_t D   = c1[0] * c2[1] - c1[1] * c2[0];
    Float_t dx  = c2[0] - c1[0];
    Float_t dy  = c2[1] - c1[1];
    Float_t dr  = TMath::Sqrt(dx*dx+dy*dy);
    Float_t det = radius * radius * dr * dr - D*D;
    
    if      (det <  0) newm = 1; // No intersection
    else if (det == 0) newm = 1; // Exactly tangent
    else {
      Float_t x   = (+D * dy + dx * TMath::Sqrt(det)) / dr / dr;
      Float_t y   = (-D * dx + dy * TMath::Sqrt(det)) / dr / dr;
      Float_t th  = TMath::ATan2(x, y);

      newm = th / (basearc/2);
    }
  }

  Float_t   slope     = (c1[1] - c2[1]) / (c1[0] - c2[0]);
  Float_t   constant  = (c2[1] * c1[0] - c2[0] * c1[1]) / (c1[0]-c2[0]);
  
  Float_t   d         = (TMath::Power(TMath::Abs(radius*slope),2) + 
			 TMath::Power(radius,2) - TMath::Power(constant,2));
  
  // If below corners return 1
  Float_t arclen = baselen;
  if (d > 0) {
    Float_t   x         = ((-TMath::Sqrt(d) - slope * constant) / 
			   (1+TMath::Power(slope, 2)));
    Float_t   y         = slope*x + constant;
    
    // If x is larger than corner x or y less than corner y, we have full
    // length strip
    if(x < c1[0] && y> c1[1]) {
      //One sector since theta is by definition half-hybrid
      Float_t   theta     = TMath::ATan2(x,y);
      arclen              = radius * theta;
    }
  }
  // Calculate the area of a strip with no cut
  Float_t   basearea1 = 0.5 * baselen * TMath::Power(radius,2);
  Float_t   basearea2 = 0.5 * baselen * TMath::Power((radius-segment),2);
  Float_t   basearea  = basearea1 - basearea2;
  
  // Calculate the area of a strip with cut
  Float_t   area1     = 0.5 * arclen * TMath::Power(radius,2);
  Float_t   area2     = 0.5 * arclen * TMath::Power((radius-segment),2);
  Float_t   area      = area1 - area2;
  
  // Acceptance is ratio 
  oldm = area/basearea;
}

/** 
 * 
 * 
 * @param r 
 * @param dt 
 * @param offT 
 * @ingroup pwg2_forward_scripts_tests
 */
void DrawSolution(Char_t r, UShort_t dt=16, UShort_t offT=128)
{
  TCanvas* c = new TCanvas(Form("c%c", r), r == 'I' ? 
			   "Inner" : "Outer", 800, 500);
  
  const Double_t ic1[] = { 4.9895, 15.3560 };
  const Double_t ic2[] = { 1.8007, 17.2000 };
  const Double_t oc1[] = { 4.2231, 26.6638 };
  const Double_t oc2[] = { 1.8357, 27.9500 };

  const Double_t* c1   = (r == 'I' ? ic1 : oc1);
  const Double_t* c2   = (r == 'I' ? ic2 : oc2);
  Double_t  minR       = (r == 'I' ?  4.5213 : 15.4);
  Double_t  maxR       = (r == 'I' ? 17.2    : 28.0);
  Double_t  rad        = maxR - minR;
  Float_t   cr         = TMath::Sqrt(c1[0]*c1[0]+c1[1]*c1[1]);
  Double_t  theta      = (r == 'I' ? 18 : 9);
  TH2*   frame = new TH2F("frame", "Frame", 100, -10, 10, 100, 0, 30);
  frame->Draw();

  TLine* line = new TLine(c1[0], c1[1], c2[0], c2[1]);
  line->SetLineColor(kBlue+1);
  line->Draw();

  UShort_t nT = (r == 'I' ? 512 : 256);
  for (Int_t t = offT; t < nT; t += dt) {
    Double_t radius = minR + t * rad / nT;

    TArc*    circle = new TArc(0, 0, radius, 90-theta, 90+theta);
    circle->SetFillColor(0);
    circle->SetFillStyle(0);
    circle->Draw();

    // Now find the intersection 
    if (radius <= cr) continue;

    // Next, we should find the end-point of the strip - that is, 
    // the coordinates where the line from c1 to c2 intersects a circle 
    // with radius given by the strip. 
    // (See http://mathworld.wolfram.com/Circle-LineIntersection.html)
    Float_t D   = c1[0] * c2[1] - c1[1] * c2[0];
    Float_t dx  = c2[0] - c1[0];
    Float_t dy  = c2[1] - c1[1];
    Float_t dr  = TMath::Sqrt(dx*dx+dy*dy);
    Float_t det = radius * radius * dr * dr - D*D;
    
    if (det <  0) continue;
    if (det == 0) continue;


    Float_t x   = (+D * dy - dx * TMath::Sqrt(det)) / dr / dr;
    Float_t y   = (-D * dx - dy * TMath::Sqrt(det)) / dr / dr;
    
    TMarker* m = new TMarker(x, y, 20);
    m->SetMarkerColor(kRed+1);
    m->Draw();

    x   = (+D * dy + dx * TMath::Sqrt(det)) / dr / dr;
    y   = (-D * dx + dy * TMath::Sqrt(det)) / dr / dr;

    // Float_t th = TMath::ATan2(x,y);
    // Printf("theta of strip %3d: %f degrees", 180. / TMath::Pi() * th);

    m = new TMarker(x, y, 20);
    m->SetMarkerColor(kGreen+1);
    m->Draw();
  }
  c->cd();
}
    
/** 
 * 
 * 
 * @ingroup pwg2_forward_scripts_tests
 */  
void TestAcc()
{
  TCanvas* c =  new TCanvas("c", "C");
  TGraph*      innerOld = new TGraph(512);
  TGraph*      outerOld = new TGraph(256);
  TGraph*      innerNew = new TGraph(512);
  TGraph*      outerNew = new TGraph(256);
  innerOld->SetName("innerOld");
  innerOld->SetName("Inner type, old");
  innerOld->SetMarkerStyle(21);
  innerOld->SetMarkerColor(kRed+1);
  outerOld->SetName("outerOld");
  outerOld->SetName("Outer type, old");
  outerOld->SetMarkerStyle(21);
  outerOld->SetMarkerColor(kBlue+1);
  innerNew->SetName("innerNew");
  innerNew->SetName("Inner type, new");
  innerNew->SetMarkerStyle(20);
  innerNew->SetMarkerColor(kGreen+1);
  outerNew->SetName("outerNew");
  outerNew->SetName("Outer type, new");
  outerNew->SetMarkerStyle(20);
  outerNew->SetMarkerColor(kMagenta+1);

  TMultiGraph* all      = new TMultiGraph("all", "Acceptances");
  all->Add(innerOld);
  all->Add(outerOld);
  all->Add(innerNew);
  all->Add(outerNew);

  TH2* innerCor = new TH2F("innerCor", "Inner correlation", 100,0,1,100,0,1);
  TH2* outerCor = new TH2F("outerCor", "Outer correlation", 100,0,1,100,0,1);
  innerCor->SetMarkerStyle(20);
  outerCor->SetMarkerStyle(20);
  innerCor->SetMarkerColor(kRed+1);
  outerCor->SetMarkerColor(kBlue+1);
  // innerCor->SetMarkerSize(1.2);
  // outerCor->SetMarkerSize(1.2);
  THStack* stack = new THStack("corr", "Correlations");
  stack->Add(innerCor);
  stack->Add(outerCor);

  for (Int_t i = 0; i < 512; i++) { 
    Float_t oldAcc, newAcc;
    
    AcceptanceCorrection('I', i, oldAcc, newAcc);
    innerOld->SetPoint(i, i, oldAcc);
    innerNew->SetPoint(i, i, newAcc);
    // Printf("Inner strip # %3d:  old=%8.6f, new=%8.6f", i, oldAcc, newAcc);

    innerCor->Fill(oldAcc, newAcc);
    if (i >= 256) continue;
    
    AcceptanceCorrection('O', i, oldAcc, newAcc);
    outerOld->SetPoint(i, i, oldAcc);
    outerNew->SetPoint(i, i, newAcc);
    // Printf("Outer strip # %3d:  old=%8.6f, new=%8.6f", i, oldAcc, newAcc);
    outerCor->Fill(oldAcc, newAcc);
  }

  all->Draw("ap");

  DrawSolution('I',4,256);
  DrawSolution('O',4,128);

  TCanvas* c2 = new TCanvas("cc", "cc");
  stack->Draw("nostack p");
  TLine* l = new TLine(0,0,1,1);
  l->SetLineStyle(2);
  l->Draw();

  c2->cd();
}
//
// EOF
//
Commit	Line	Data
e5f38055	1	#include <TMath.h>
	2	#include <TGraph.h>
	3	#include <TMultiGraph.h>
	4	#include <TMarker.h>
	5	#include <TLine.h>
	6	#include <TArc.h>
	7	#include <TCanvas.h>
	8	#include <TH2F.h>
	9	#include <THStack.h>
56199f2b	10	/**
f0ef1e71	11	* @defgroup pwg2_forward_scripts_tests Test scripts
56199f2b	12	*
	13	*
	14	*
f0ef1e71	15	* @ingroup pwg2_forward_scripts
56199f2b	16	*/
56199f2b	17	/**
	18	*
	19	*
	20	* @param r
	21	* @param t
	22	* @param oldm
	23	* @param newm
	24	*
f0ef1e71	25	* @ingroup pwg2_forward_scripts_tests
56199f2b	26	*/
e5f38055	27	//_____________________________________________________________________
	28	void AcceptanceCorrection(Char_t r, UShort_t t, Float_t& oldm, Float_t& newm)
	29	{
	30	// AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
	31
	32	//AliFMDRing fmdring(ring);
	33	//fmdring.Init();
	34	// Float_t rad = pars->GetMaxR(r)-pars->GetMinR(r);
	35	// Float_t slen = pars->GetStripLength(r,t);
	36	// Float_t sblen = pars->GetBaseStripLength(r,t);
	37	const Double_t ic1[] = { 4.9895, 15.3560 };
	38	const Double_t ic2[] = { 1.8007, 17.2000 };
	39	const Double_t oc1[] = { 4.2231, 26.6638 };
	40	const Double_t oc2[] = { 1.8357, 27.9500 };
	41
	42	Double_t minR = (r == 'I' ? 4.5213 : 15.4);
	43	Double_t maxR = (r == 'I' ? 17.2 : 28.0);
	44	Double_t rad = maxR - minR;
	45
	46	Int_t nStrips = (r == 'I' ? 512 : 256);
	47	Int_t nSec = (r == 'I' ? 20 : 40);
	48	Float_t segment = rad / nStrips;
	49	Float_t basearc = 2 * TMath::Pi() / (0.5 * nSec);
	50	Float_t radius = minR + t * segment;
	51	Float_t baselen = 0.5 * basearc * radius;
	52
	53	const Double_t* c1 = (r == 'I' ? ic1 : oc1);
	54	const Double_t* c2 = (r == 'I' ? ic2 : oc2);
	55
	56	// If the radius of the strip is smaller than the radius corresponding
	57	// to the first corner we have a full strip length
	58	Float_t cr = TMath::Sqrt(c1[0]c1[0]+c1[1]c1[1]);
	59	if (radius <= cr) newm = 1;
	60	else {
	61	// Next, we should find the end-point of the strip - that is,
	62	// the coordinates where the line from c1 to c2 intersects a circle
	63	// with radius given by the strip.
	64	// (See http://mathworld.wolfram.com/Circle-LineIntersection.html)
	65	Float_t D = c1[0] * c2[1] - c1[1] * c2[0];
	66	Float_t dx = c2[0] - c1[0];
	67	Float_t dy = c2[1] - c1[1];
	68	Float_t dr = TMath::Sqrt(dxdx+dydy);
	69	Float_t det = radius * radius * dr * dr - D*D;
	70
	71	if (det < 0) newm = 1; // No intersection
	72	else if (det == 0) newm = 1; // Exactly tangent
	73	else {
	74	Float_t x = (+D * dy + dx * TMath::Sqrt(det)) / dr / dr;
	75	Float_t y = (-D * dx + dy * TMath::Sqrt(det)) / dr / dr;
	76	Float_t th = TMath::ATan2(x, y);
	77
	78	newm = th / (basearc/2);
	79	}
	80	}
	81
	82	Float_t slope = (c1[1] - c2[1]) / (c1[0] - c2[0]);
	83	Float_t constant = (c2[1] * c1[0] - c2[0] * c1[1]) / (c1[0]-c2[0]);
	84
	85	Float_t d = (TMath::Power(TMath::Abs(radius*slope),2) +
	86	TMath::Power(radius,2) - TMath::Power(constant,2));
	87
	88	// If below corners return 1
	89	Float_t arclen = baselen;
	90	if (d > 0) {
91	Float_t x = ((-TMath::Sqrt(d) - slope * constant) /
92	(1+TMath::Power(slope, 2)));
93	Float_t y = slope*x + constant;
94
95	// If x is larger than corner x or y less than corner y, we have full
96	// length strip
97	if(x < c1[0] && y> c1[1]) {
98	//One sector since theta is by definition half-hybrid
99	Float_t theta = TMath::ATan2(x,y);
100	arclen = radius * theta;
101	}
102	}
103	// Calculate the area of a strip with no cut
104	Float_t basearea1 = 0.5 * baselen * TMath::Power(radius,2);
105	Float_t basearea2 = 0.5 * baselen * TMath::Power((radius-segment),2);
106	Float_t basearea = basearea1 - basearea2;
107
108	// Calculate the area of a strip with cut
109	Float_t area1 = 0.5 * arclen * TMath::Power(radius,2);
110	Float_t area2 = 0.5 * arclen * TMath::Power((radius-segment),2);
111	Float_t area = area1 - area2;
112
113	// Acceptance is ratio
114	oldm = area/basearea;
115	}
116
56199f2b	117	/**
	118	*
	119	*
	120	* @param r
	121	* @param dt
	122	* @param offT
f0ef1e71	123	* @ingroup pwg2_forward_scripts_tests
56199f2b	124	*/
e5f38055	125	void DrawSolution(Char_t r, UShort_t dt=16, UShort_t offT=128)
	126	{
	127	TCanvas* c = new TCanvas(Form("c%c", r), r == 'I' ?
	128	"Inner" : "Outer", 800, 500);
	129
	130	const Double_t ic1[] = { 4.9895, 15.3560 };
	131	const Double_t ic2[] = { 1.8007, 17.2000 };
	132	const Double_t oc1[] = { 4.2231, 26.6638 };
	133	const Double_t oc2[] = { 1.8357, 27.9500 };
	134
	135	const Double_t* c1 = (r == 'I' ? ic1 : oc1);
	136	const Double_t* c2 = (r == 'I' ? ic2 : oc2);
	137	Double_t minR = (r == 'I' ? 4.5213 : 15.4);
	138	Double_t maxR = (r == 'I' ? 17.2 : 28.0);
	139	Double_t rad = maxR - minR;
	140	Float_t cr = TMath::Sqrt(c1[0]c1[0]+c1[1]c1[1]);
	141	Double_t theta = (r == 'I' ? 18 : 9);
	142	TH2* frame = new TH2F("frame", "Frame", 100, -10, 10, 100, 0, 30);
	143	frame->Draw();
	144
	145	TLine* line = new TLine(c1[0], c1[1], c2[0], c2[1]);
	146	line->SetLineColor(kBlue+1);
	147	line->Draw();
	148
	149	UShort_t nT = (r == 'I' ? 512 : 256);
	150	for (Int_t t = offT; t < nT; t += dt) {
	151	Double_t radius = minR + t * rad / nT;
	152
	153	TArc* circle = new TArc(0, 0, radius, 90-theta, 90+theta);
	154	circle->SetFillColor(0);
	155	circle->SetFillStyle(0);
	156	circle->Draw();
	157
	158	// Now find the intersection
	159	if (radius <= cr) continue;
	160
	161	// Next, we should find the end-point of the strip - that is,
	162	// the coordinates where the line from c1 to c2 intersects a circle
	163	// with radius given by the strip.
	164	// (See http://mathworld.wolfram.com/Circle-LineIntersection.html)
	165	Float_t D = c1[0] * c2[1] - c1[1] * c2[0];
	166	Float_t dx = c2[0] - c1[0];
	167	Float_t dy = c2[1] - c1[1];
	168	Float_t dr = TMath::Sqrt(dxdx+dydy);
	169	Float_t det = radius * radius * dr * dr - D*D;
	170
	171	if (det < 0) continue;
	172	if (det == 0) continue;
	173
	174
	175	Float_t x = (+D * dy - dx * TMath::Sqrt(det)) / dr / dr;
	176	Float_t y = (-D * dx - dy * TMath::Sqrt(det)) / dr / dr;
	177
	178	TMarker* m = new TMarker(x, y, 20);
	179	m->SetMarkerColor(kRed+1);
	180	m->Draw();
	181
	182	x = (+D * dy + dx * TMath::Sqrt(det)) / dr / dr;
	183	y = (-D * dx + dy * TMath::Sqrt(det)) / dr / dr;
	184
	185	// Float_t th = TMath::ATan2(x,y);
	186	// Printf("theta of strip %3d: %f degrees", 180. / TMath::Pi() * th);
	187
	188	m = new TMarker(x, y, 20);
189	m->SetMarkerColor(kGreen+1);
190	m->Draw();
191	}
192	c->cd();
193	}
194
56199f2b	195	/**
	196	*
	197	*
f0ef1e71	198	* @ingroup pwg2_forward_scripts_tests
56199f2b	199	*/
e5f38055	200	void TestAcc()
	201	{
	202	TCanvas* c = new TCanvas("c", "C");
	203	TGraph* innerOld = new TGraph(512);
	204	TGraph* outerOld = new TGraph(256);
	205	TGraph* innerNew = new TGraph(512);
	206	TGraph* outerNew = new TGraph(256);
	207	innerOld->SetName("innerOld");
	208	innerOld->SetName("Inner type, old");
	209	innerOld->SetMarkerStyle(21);
	210	innerOld->SetMarkerColor(kRed+1);
	211	outerOld->SetName("outerOld");
	212	outerOld->SetName("Outer type, old");
	213	outerOld->SetMarkerStyle(21);
	214	outerOld->SetMarkerColor(kBlue+1);
	215	innerNew->SetName("innerNew");
	216	innerNew->SetName("Inner type, new");
	217	innerNew->SetMarkerStyle(20);
	218	innerNew->SetMarkerColor(kGreen+1);
	219	outerNew->SetName("outerNew");
	220	outerNew->SetName("Outer type, new");
	221	outerNew->SetMarkerStyle(20);
	222	outerNew->SetMarkerColor(kMagenta+1);
	223
	224	TMultiGraph* all = new TMultiGraph("all", "Acceptances");
	225	all->Add(innerOld);
	226	all->Add(outerOld);
	227	all->Add(innerNew);
	228	all->Add(outerNew);
	229
	230	TH2* innerCor = new TH2F("innerCor", "Inner correlation", 100,0,1,100,0,1);
	231	TH2* outerCor = new TH2F("outerCor", "Outer correlation", 100,0,1,100,0,1);
	232	innerCor->SetMarkerStyle(20);
	233	outerCor->SetMarkerStyle(20);
	234	innerCor->SetMarkerColor(kRed+1);
	235	outerCor->SetMarkerColor(kBlue+1);
	236	// innerCor->SetMarkerSize(1.2);
	237	// outerCor->SetMarkerSize(1.2);
	238	THStack* stack = new THStack("corr", "Correlations");
	239	stack->Add(innerCor);
	240	stack->Add(outerCor);
	241
	242	for (Int_t i = 0; i < 512; i++) {
	243	Float_t oldAcc, newAcc;
	244
	245	AcceptanceCorrection('I', i, oldAcc, newAcc);
	246	innerOld->SetPoint(i, i, oldAcc);
	247	innerNew->SetPoint(i, i, newAcc);
	248	// Printf("Inner strip # %3d: old=%8.6f, new=%8.6f", i, oldAcc, newAcc);
	249
	250	innerCor->Fill(oldAcc, newAcc);
	251	if (i >= 256) continue;
	252
	253	AcceptanceCorrection('O', i, oldAcc, newAcc);
	254	outerOld->SetPoint(i, i, oldAcc);
	255	outerNew->SetPoint(i, i, newAcc);
	256	// Printf("Outer strip # %3d: old=%8.6f, new=%8.6f", i, oldAcc, newAcc);
	257	outerCor->Fill(oldAcc, newAcc);
	258	}
	259
	260	all->Draw("ap");
	261
	262	DrawSolution('I',4,256);
	263	DrawSolution('O',4,128);
264
265	TCanvas* c2 = new TCanvas("cc", "cc");
266	stack->Draw("nostack p");
267	TLine* l = new TLine(0,0,1,1);
268	l->SetLineStyle(2);
269	l->Draw();
270
271	c2->cd();
272	}
56199f2b	273	//
	274	// EOF
	275	//