Some extra tests

diff --git a/FMD/scripts/CompareMatDep.C b/FMD/scripts/CompareMatDep.C
new file mode 100644 (file)
index 0000000..cb2c6d7
--- /dev/null
+++ b/FMD/scripts/CompareMatDep.C
@@ -0,0 +1,145 @@
+Int_t 
+DetIndex(Int_t d, Char_t r)
+{
+  return (d == 1 ? 0 : (d-1) * 2 - 1) + (r == 'I' ? 0 : 1);
+}
+Int_t
+RingColor(Int_t d, Char_t r)
+{
+  return (d == 1 ? kRed : d == 2 ? kGreen : kBlue) + (r == 'I' ? 1 : 3);
+}
+
+void CompareMatDep() {
+  gStyle->SetOptTitle(kFALSE);
+
+  Int_t spread = 10;
+  TMultiGraph* mg = new TMultiGraph;
+
+  for(Int_t d = 1; d<=3;d++) {
+    Int_t nr = (d == 1 ? 1 : 2);
+    for(Int_t q = 0; q < nr; q++) {
+      Char_t        r = (q == 0 ? 'I' : 'O');
+      TGraphErrors* g = new TGraphErrors();
+      g->SetName(Form("fmd%d%c", d, r));
+      g->SetTitle(Form("FMD%d%c (+%d)", d, r, spread*(DetIndex(d,r)+1)));
+      g->SetLineColor(RingColor(d, r));
+      g->SetMarkerColor(RingColor(d, r));
+      g->SetFillColor(RingColor(d, r));
+      g->SetFillStyle(3002);
+      g->SetMarkerStyle(20);
+      g->SetMarkerSize(1.3);
+      mg->Add(g);
+    }
+  }
+  
+  TGraphErrors* g = new TGraphErrors();
+  g->SetName("all");
+  g->SetTitle("All");
+  g->SetLineColor(kBlack);
+  g->SetMarkerColor(kBlack);
+  g->SetFillColor(kGray);
+  g->SetFillStyle(3002);
+  g->SetMarkerStyle(20);
+  g->SetMarkerSize(1.3);
+  mg->Add(g);
+
+  const char* base = "/mnt/samsung/oldhome/ALICE/FMDanalysis/BackgroundCorrection/BackgroundStudy/TestOfMaterialBudget/";
+
+  Int_t  samples[] = { -20, -15, -10, -5, 0, 5, 10 , 15, 20, 999 };
+  Int_t* sample    = samples;
+  Int_t  i         = 0;
+  while (*sample < 999) { 
+    TFile* fin = 0;
+    if(*sample < 0) 
+      fin = TFile::Open(Form("%s/minus%dpercent/fmd_dNdeta.root",
+                            base,-*sample));
+    else if(*sample > 0) 
+      fin = TFile::Open(Form("%s/plus%dpercent/fmd_dNdeta.root",base,*sample));
+    if(*sample == 0) 
+      fin = TFile::Open(Form("%s/refdist/fmd_dNdeta.root", base));
+    if (!fin) { 
+      Error("CompareMatDep", "Failed to open input file for %d", *sample);
+      sample++;
+      i++;
+      continue;
+    }
+    Float_t twa = 0;
+    Float_t tsw = 0;
+
+    for(Int_t d = 1; d<=3;d++) {
+      Int_t nr = (d == 1 ? 1 : 2);
+      for(Int_t q = 0; q < nr; q++) {
+       Char_t        r = (q == 0 ? 'I' : 'O');
+       TH1F* hRingMult = (TH1F*)fin->Get(Form("hRingMult_FMD%d%c",
+                                             d,r));
+       Float_t wav    = 0;
+       Float_t sumofw = 0;
+
+       // Weighted average 
+       for(Int_t n = 1;n<=hRingMult->GetNbinsX();n++) {
+         if(hRingMult->GetBinContent(n) > -1 && 
+            hRingMult->GetBinError(n) > 0 ) {
+           Float_t weight = 1/TMath::Power(hRingMult->GetBinError(n),2);
+           wav            += weight*hRingMult->GetBinContent(n);
+           sumofw         += weight;
+           twa            += weight*hRingMult->GetBinContent(n);
+           tsw             += weight;
+         }
+       }
+       
+       TGraphErrors* g = 
+         static_cast<TGraphErrors*>(mg->GetListOfGraphs()->At(DetIndex(d, r)));
+       g->SetPoint(i, *sample, 100*wav/sumofw+spread*(DetIndex(d,r)+1));
+       g->SetPointError(i, 0, 100*1/TMath::Sqrt(sumofw)); 
+       Info("CompareMatDep", "Setting point %d @ %d of %s to %f +/- %f", 
+            i, *sample, g->GetName(), 100*wav/sumofw, 
+            100*1/TMath::Sqrt(sumofw));
+      } // for q 
+    } // for d
+
+    TGraphErrors* g = 
+      static_cast<TGraphErrors*>(mg->GetListOfGraphs()->At(5));
+    g->SetPoint(i, *sample, 100*twa/tsw);
+    g->SetPointError(i, 0, 100*1/TMath::Sqrt(tsw)); 
+    Info("CompareMatDep", "Setting point %d @ %d of %s to %f +/- %f", 
+        i, *sample, g->GetName(), 100*twa/tsw, 
+        100*1/TMath::Sqrt(tsw));
+
+    sample++;
+    i++;
+  }
+  
+  TCanvas* c  = new TCanvas("c", "C", 800, 800);
+  c->SetTopMargin(0.05);
+  c->SetRightMargin(0.05);
+  c->SetLeftMargin(0.13);
+  c->SetFillColor(0);
+  c->SetBorderSize(0);
+  c->SetBorderMode(0);
+  c->SetGridx();
+  c->SetGridy();
+
+  mg->Draw("a pl");
+  mg->GetHistogram()->SetXTitle("#Delta#rho [%]");
+  mg->GetHistogram()->SetYTitle("Average relative difference [%]");
+  mg->GetHistogram()->GetXaxis()->SetTitleFont(132);
+  mg->GetHistogram()->GetXaxis()->SetLabelFont(132);
+  mg->GetHistogram()->GetXaxis()->SetNdivisions(10);
+  mg->GetHistogram()->GetYaxis()->SetTitleFont(132);
+  mg->GetHistogram()->GetYaxis()->SetLabelFont(132);
+  mg->GetHistogram()->GetYaxis()->SetNdivisions(1010);
+  mg->GetHistogram()->GetYaxis()->SetTitleOffset(1.2);
+  // c->Clear();
+  // mg->Draw("a e3");
+
+  TLegend* leg = gPad->BuildLegend(0.15,0.7,0.4,0.945);
+  leg->SetFillColor(0);
+  leg->SetBorderSize(0);
+  leg->SetTextFont(132);
+  leg->SetFillStyle(0);
+
+  TIter next(mg->GetListOfGraphs());
+  TObject* o = 0;
+
+  c->Print("materialDependence.png");
+}

diff --git a/FMD/scripts/TestAcc.C b/FMD/scripts/TestAcc.C
new file mode 100644 (file)
index 0000000..b012f93
--- /dev/null
+++ b/FMD/scripts/TestAcc.C
@@ -0,0 +1,244 @@
+#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>
+
+//_____________________________________________________________________
+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;
+}
+
+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();
+}
+    
+  
+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();
+}

diff --git a/FMD/scripts/polar.C b/FMD/scripts/polar.C
new file mode 100644 (file)
index 0000000..25458a0
--- /dev/null
+++ b/FMD/scripts/polar.C
@@ -0,0 +1,47 @@
+void
+polar()
+{
+   TCanvas *c1 = new TCanvas("c1","c1",500,500);
+   TGraphPolar * grP1 = new TGraphPolar();
+   grP1->SetTitle("TGraphPolar example");
+
+   for (Int_t i = 0; i < 2; i++) {
+     grP1->SetPoint(i, (i+1) * TMath::Pi() / 4, (i+1) * 0.05);
+     grP1->SetPointError(i, 9*TMath::Pi()/180, 0.0);
+     Double_t rr = grP1->GetY()[i];
+     Double_t tt = grP1->GetX()[i];
+     Double_t x  = rr * TMath::Cos(tt);
+     Double_t y  = rr * TMath::Sin(tt);
+     Printf("(x,y)=(%f,%f)", x, y);
+   }
+   Double_t r = 1;
+   TH2* frame = new TH2F("frame", "Frame", 100, -r,r, 100, -r, r);
+   frame->Draw();
+
+   grP1->SetMarkerStyle(1);
+   grP1->SetMarkerSize(1.);
+   grP1->SetMarkerColor(4);
+   grP1->SetLineColor(4);
+   grP1->SetLineWidth(3);
+   grP1->SetFillColor(kRed+1);
+   grP1->SetFillStyle(3001);
+   grP1->Draw("PNEF same");
+   // grP1->Draw("APNEF");
+
+   // Update, otherwise GetPolargram returns 0
+   c1->Update();
+   TGraphPolargram* gram = grP1->GetPolargram();
+   gram->SetLineWidth(0);
+   // gram->SetLineColor(kWhite);
+   gram->SetNdivPolar(20);
+   gram->SetNdivRadial(10);   
+   gram->SetTextSize(0);
+   gram->SetRadialLabelSize(0);
+   gram->SetPolarLabelSize(0);
+   gram->SetAxisAngle(9*TMath::Pi()/180);
+   gram->SetTwoPi();
+   gram->SetToRadian();
+   c1->SetGridx();
+   c1->SetGridy();
+   c1->Update();
+}