PWG2/FORWARD/analysis2/scripts/TestAcc.C

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