ROOT logo
#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 pwglf_forward_scripts_tests Test scripts
 * 
 * 
 *
 * @ingroup pwglf_forward_scripts
 */
/** 
 * 
 * 
 * @param r 
 * @param t 
 * @param oldm 
 * @param newm 
 *
 * @ingroup pwglf_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 pwglf_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 pwglf_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
//
 TestAcc.C:1
 TestAcc.C:2
 TestAcc.C:3
 TestAcc.C:4
 TestAcc.C:5
 TestAcc.C:6
 TestAcc.C:7
 TestAcc.C:8
 TestAcc.C:9
 TestAcc.C:10
 TestAcc.C:11
 TestAcc.C:12
 TestAcc.C:13
 TestAcc.C:14
 TestAcc.C:15
 TestAcc.C:16
 TestAcc.C:17
 TestAcc.C:18
 TestAcc.C:19
 TestAcc.C:20
 TestAcc.C:21
 TestAcc.C:22
 TestAcc.C:23
 TestAcc.C:24
 TestAcc.C:25
 TestAcc.C:26
 TestAcc.C:27
 TestAcc.C:28
 TestAcc.C:29
 TestAcc.C:30
 TestAcc.C:31
 TestAcc.C:32
 TestAcc.C:33
 TestAcc.C:34
 TestAcc.C:35
 TestAcc.C:36
 TestAcc.C:37
 TestAcc.C:38
 TestAcc.C:39
 TestAcc.C:40
 TestAcc.C:41
 TestAcc.C:42
 TestAcc.C:43
 TestAcc.C:44
 TestAcc.C:45
 TestAcc.C:46
 TestAcc.C:47
 TestAcc.C:48
 TestAcc.C:49
 TestAcc.C:50
 TestAcc.C:51
 TestAcc.C:52
 TestAcc.C:53
 TestAcc.C:54
 TestAcc.C:55
 TestAcc.C:56
 TestAcc.C:57
 TestAcc.C:58
 TestAcc.C:59
 TestAcc.C:60
 TestAcc.C:61
 TestAcc.C:62
 TestAcc.C:63
 TestAcc.C:64
 TestAcc.C:65
 TestAcc.C:66
 TestAcc.C:67
 TestAcc.C:68
 TestAcc.C:69
 TestAcc.C:70
 TestAcc.C:71
 TestAcc.C:72
 TestAcc.C:73
 TestAcc.C:74
 TestAcc.C:75
 TestAcc.C:76
 TestAcc.C:77
 TestAcc.C:78
 TestAcc.C:79
 TestAcc.C:80
 TestAcc.C:81
 TestAcc.C:82
 TestAcc.C:83
 TestAcc.C:84
 TestAcc.C:85
 TestAcc.C:86
 TestAcc.C:87
 TestAcc.C:88
 TestAcc.C:89
 TestAcc.C:90
 TestAcc.C:91
 TestAcc.C:92
 TestAcc.C:93
 TestAcc.C:94
 TestAcc.C:95
 TestAcc.C:96
 TestAcc.C:97
 TestAcc.C:98
 TestAcc.C:99
 TestAcc.C:100
 TestAcc.C:101
 TestAcc.C:102
 TestAcc.C:103
 TestAcc.C:104
 TestAcc.C:105
 TestAcc.C:106
 TestAcc.C:107
 TestAcc.C:108
 TestAcc.C:109
 TestAcc.C:110
 TestAcc.C:111
 TestAcc.C:112
 TestAcc.C:113
 TestAcc.C:114
 TestAcc.C:115
 TestAcc.C:116
 TestAcc.C:117
 TestAcc.C:118
 TestAcc.C:119
 TestAcc.C:120
 TestAcc.C:121
 TestAcc.C:122
 TestAcc.C:123
 TestAcc.C:124
 TestAcc.C:125
 TestAcc.C:126
 TestAcc.C:127
 TestAcc.C:128
 TestAcc.C:129
 TestAcc.C:130
 TestAcc.C:131
 TestAcc.C:132
 TestAcc.C:133
 TestAcc.C:134
 TestAcc.C:135
 TestAcc.C:136
 TestAcc.C:137
 TestAcc.C:138
 TestAcc.C:139
 TestAcc.C:140
 TestAcc.C:141
 TestAcc.C:142
 TestAcc.C:143
 TestAcc.C:144
 TestAcc.C:145
 TestAcc.C:146
 TestAcc.C:147
 TestAcc.C:148
 TestAcc.C:149
 TestAcc.C:150
 TestAcc.C:151
 TestAcc.C:152
 TestAcc.C:153
 TestAcc.C:154
 TestAcc.C:155
 TestAcc.C:156
 TestAcc.C:157
 TestAcc.C:158
 TestAcc.C:159
 TestAcc.C:160
 TestAcc.C:161
 TestAcc.C:162
 TestAcc.C:163
 TestAcc.C:164
 TestAcc.C:165
 TestAcc.C:166
 TestAcc.C:167
 TestAcc.C:168
 TestAcc.C:169
 TestAcc.C:170
 TestAcc.C:171
 TestAcc.C:172
 TestAcc.C:173
 TestAcc.C:174
 TestAcc.C:175
 TestAcc.C:176
 TestAcc.C:177
 TestAcc.C:178
 TestAcc.C:179
 TestAcc.C:180
 TestAcc.C:181
 TestAcc.C:182
 TestAcc.C:183
 TestAcc.C:184
 TestAcc.C:185
 TestAcc.C:186
 TestAcc.C:187
 TestAcc.C:188
 TestAcc.C:189
 TestAcc.C:190
 TestAcc.C:191
 TestAcc.C:192
 TestAcc.C:193
 TestAcc.C:194
 TestAcc.C:195
 TestAcc.C:196
 TestAcc.C:197
 TestAcc.C:198
 TestAcc.C:199
 TestAcc.C:200
 TestAcc.C:201
 TestAcc.C:202
 TestAcc.C:203
 TestAcc.C:204
 TestAcc.C:205
 TestAcc.C:206
 TestAcc.C:207
 TestAcc.C:208
 TestAcc.C:209
 TestAcc.C:210
 TestAcc.C:211
 TestAcc.C:212
 TestAcc.C:213
 TestAcc.C:214
 TestAcc.C:215
 TestAcc.C:216
 TestAcc.C:217
 TestAcc.C:218
 TestAcc.C:219
 TestAcc.C:220
 TestAcc.C:221
 TestAcc.C:222
 TestAcc.C:223
 TestAcc.C:224
 TestAcc.C:225
 TestAcc.C:226
 TestAcc.C:227
 TestAcc.C:228
 TestAcc.C:229
 TestAcc.C:230
 TestAcc.C:231
 TestAcc.C:232
 TestAcc.C:233
 TestAcc.C:234
 TestAcc.C:235
 TestAcc.C:236
 TestAcc.C:237
 TestAcc.C:238
 TestAcc.C:239
 TestAcc.C:240
 TestAcc.C:241
 TestAcc.C:242
 TestAcc.C:243
 TestAcc.C:244
 TestAcc.C:245
 TestAcc.C:246
 TestAcc.C:247
 TestAcc.C:248
 TestAcc.C:249
 TestAcc.C:250
 TestAcc.C:251
 TestAcc.C:252
 TestAcc.C:253
 TestAcc.C:254
 TestAcc.C:255
 TestAcc.C:256
 TestAcc.C:257
 TestAcc.C:258
 TestAcc.C:259
 TestAcc.C:260
 TestAcc.C:261
 TestAcc.C:262
 TestAcc.C:263
 TestAcc.C:264
 TestAcc.C:265
 TestAcc.C:266
 TestAcc.C:267
 TestAcc.C:268
 TestAcc.C:269
 TestAcc.C:270
 TestAcc.C:271
 TestAcc.C:272
 TestAcc.C:273
 TestAcc.C:274
 TestAcc.C:275
 TestAcc.C:276