PWG2/FORWARD/analysis2/tests/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  * @defgroup pwg2_forward_scripts_tests Test scripts
  12  *
  13  *
  14  *
  15  * @ingroup pwg2_forward_scripts
  16  */
  17 /**
  18  *
  19  *
  20  * @param r
  21  * @param t
  22  * @param oldm
  23  * @param newm
  24  *
  25  * @ingroup pwg2_forward_scripts_tests
  26  */
  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(dx*dx+dy*dy);
  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
 117 /**
 118  *
 119  *
 120  * @param r
 121  * @param dt
 122  * @param offT
 123  * @ingroup pwg2_forward_scripts_tests
 124  */
 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(dx*dx+dy*dy);
 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
 195 /**
 196  *
 197  *
 198  * @ingroup pwg2_forward_scripts_tests
 199  */
 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 }
 273 //
 274 // EOF
 275 //