adding nsigma study
authorjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 26 Jul 2006 16:33:58 +0000 (16:33 +0000)
committerjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 26 Jul 2006 16:33:58 +0000 (16:33 +0000)
PWG0/dNdEta/AlidNdEtaSystematicsSelector.cxx
PWG0/dNdEta/AlidNdEtaSystematicsSelector.h
PWG0/dNdEta/drawSystematics.C

index 719a986..afd8666 100644 (file)
@@ -27,6 +27,7 @@ ClassImp(AlidNdEtaSystematicsSelector)
 AlidNdEtaSystematicsSelector::AlidNdEtaSystematicsSelector() :
   AliSelectorRL(),
   fSecondaries(0),
+  fSigmaVertex(0),
   fEsdTrackCuts(0)
 {
   //
@@ -83,7 +84,7 @@ void AlidNdEtaSystematicsSelector::SlaveBegin(TTree* tree)
 
   if (option.Contains("secondaries"))
   {
-    fSecondaries = new TH3F("fSecondaries", "fSecondaries;NacceptedTracks;CratioSecondaries;p_{T}", 201, -0.5, 205.5, 16, 0.45, 2.05, 10, 0, 10);
+    fSecondaries = new TH3F("fSecondaries", "fSecondaries;NacceptedTracks;CratioSecondaries;p_{T}", 2000, -0.5, 205.5, 16, 0.45, 2.05, 10, 0, 10);
   }
 
   if (option.Contains("particle-composition"))
@@ -95,6 +96,12 @@ void AlidNdEtaSystematicsSelector::SlaveBegin(TTree* tree)
       fdNdEtaCorrection[i] = new AlidNdEtaCorrection(name, name);
     }
   }
+
+  if (option.Contains("sigma-vertex"))
+  {
+    fSigmaVertex = new TH1F("fSigmaVertex", "fSigmaVertex;Nsigma2vertex;NacceptedTracks", 10, 0.25, 5.25);
+    printf("WARNING: sigma-vertex analysis enabled. This will produce weird results in the AliESDtrackCuts histograms\n");
+  }
 }
 
 Bool_t AlidNdEtaSystematicsSelector::Process(Long64_t entry)
@@ -148,6 +155,9 @@ Bool_t AlidNdEtaSystematicsSelector::Process(Long64_t entry)
   if (fSecondaries)
     FillSecondaries(list);
 
+  if (fSigmaVertex)
+    FillSigmaVertex();
+
   delete list;
   list = 0;
 
@@ -240,10 +250,10 @@ void AlidNdEtaSystematicsSelector::FillCorrectionMaps(TObjArray* listOfTracks)
     Int_t id = -1;
     switch (TMath::Abs(particle->GetPdgCode()))
     {
-      case 211: id = 0; break;
-      case 321: id = 1; break;
+      case 211:  id = 0; break;
+      case 321:  id = 1; break;
       case 2212: id = 2; break;
-      default: id = 3; break;
+      default:   id = 3; break;
     }
 
     if (vertexReconstructed)
@@ -302,7 +312,11 @@ void AlidNdEtaSystematicsSelector::FillSecondaries(TObjArray* listOfTracks)
       AliDebug(AliLog::kDebug, Form("The ratio between primaries and secondaries is %d:%d = %f", primaries, secondaries, ((secondaries > 0) ? (Double_t) primaries / secondaries : -1)));
 
       for (Double_t factor = 0.5; factor < 2.01; factor += 0.1)
-        fSecondaries->Fill((Double_t) primaries + (Double_t) secondaries * factor, factor, nPrimaries->GetBinCenter(i));
+      {
+        Double_t nTracks = (Double_t) primaries + (Double_t) secondaries * factor;
+        fSecondaries->Fill(nTracks, factor, nPrimaries->GetBinCenter(i));
+        //if (secondaries > 0) printf("We fill: %f %f %f\n", nTracks, factor, nPrimaries->GetBinCenter(i));
+      }
     }
   }
 
@@ -316,6 +330,45 @@ void AlidNdEtaSystematicsSelector::FillSecondaries(TObjArray* listOfTracks)
   iter = 0;
 }
 
+void AlidNdEtaSystematicsSelector::FillSigmaVertex()
+{
+  // fills the fSigmaVertex histogram
+
+  // save the old value
+  Float_t oldSigmaVertex = fEsdTrackCuts->GetMinNsigmaToVertex();
+
+  // set to maximum
+  fEsdTrackCuts->SetMinNsigmaToVertex(5);
+
+  TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
+
+  TIterator* iter = list->MakeIterator();
+  TObject* obj = 0;
+  while ((obj = iter->Next()) != 0)
+  {
+    AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (obj);
+    if (!esdTrack)
+      continue;
+
+    Float_t sigma2Vertex = fEsdTrackCuts->GetSigmaToVertex(esdTrack);
+
+    for (Double_t nSigma = 0.5; nSigma < 5.1; nSigma += 0.5)
+    {
+      if (sigma2Vertex < nSigma)
+        fSigmaVertex->Fill(nSigma);
+    }
+  }
+
+  delete iter;
+  iter = 0;
+
+  delete list;
+  list = 0;
+
+  // set back the old value
+  fEsdTrackCuts->SetMinNsigmaToVertex(oldSigmaVertex);
+}
+
 void AlidNdEtaSystematicsSelector::SlaveTerminate()
 {
   // The SlaveTerminate() function is called after all entries or objects
@@ -337,6 +390,9 @@ void AlidNdEtaSystematicsSelector::SlaveTerminate()
   for (Int_t i=0; i<4; ++i)
     if (fdNdEtaCorrection[i])
       fOutput->Add(fdNdEtaCorrection[i]);
+
+  if (fSigmaVertex)
+    fOutput->Add(fSigmaVertex);
 }
 
 void AlidNdEtaSystematicsSelector::Terminate()
@@ -350,6 +406,7 @@ void AlidNdEtaSystematicsSelector::Terminate()
   fSecondaries = dynamic_cast<TH3F*> (fOutput->FindObject("fSecondaries"));
   for (Int_t i=0; i<4; ++i)
     fdNdEtaCorrection[i] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject(Form("correction_%d", i)));
+  fSigmaVertex = dynamic_cast<TH1F*> (fOutput->FindObject("fSigmaVertex"));
 
   TFile* fout = TFile::Open("systematics.root", "RECREATE");
 
@@ -359,6 +416,9 @@ void AlidNdEtaSystematicsSelector::Terminate()
   if (fSecondaries)
     fSecondaries->Write();
 
+  if (fSigmaVertex)
+    fSigmaVertex->Write();
+
   for (Int_t i=0; i<4; ++i)
     if (fdNdEtaCorrection[i])
       fdNdEtaCorrection[i]->SaveHistograms();
@@ -371,4 +431,10 @@ void AlidNdEtaSystematicsSelector::Terminate()
     new TCanvas;
     fSecondaries->Draw();
   }
+
+  if (fSigmaVertex)
+  {
+    new TCanvas;
+    fSigmaVertex->Draw();
+  }
 }
index 29aef2f..b04d3f3 100644 (file)
@@ -8,6 +8,7 @@
 class AliESDtrackCuts;
 class AlidNdEtaCorrection;
 class TH3F;
+class TH1F;
 
 class AlidNdEtaSystematicsSelector : public AliSelectorRL {
   public:
@@ -25,9 +26,11 @@ class AlidNdEtaSystematicsSelector : public AliSelectorRL {
 
     void FillCorrectionMaps(TObjArray* listOfTracks);
     void FillSecondaries(TObjArray* listOfTracks);
+    void FillSigmaVertex();
 
     TH3F* fSecondaries; // (accepted tracks) vs (tracks from sec)/(n * tracks from sec) vs pT
     AlidNdEtaCorrection* fdNdEtaCorrection[4];      // correction for different particle species: here pi, K, p, others
+    TH1F* fSigmaVertex; // (accepted tracks) vs (n of sigma to vertex cut)
 
     AliESDtrackCuts* fEsdTrackCuts;     // Object containing the parameters of the esd track cuts
 
index e140a12..7d47d7c 100644 (file)
@@ -16,6 +16,8 @@
 
 #endif
 
+extern TPad* gPad;
+
 void SetRanges(TAxis* axis)
 {
   if (strcmp(axis->GetTitle(), "#eta") == 0)
@@ -62,6 +64,9 @@ void Prepare1DPlot(TH1* hist)
 
 void InitPad()
 {
+  if (!gPad)
+    return;
+
   gPad->Range(0, 0, 1, 1);
   gPad->SetLeftMargin(0.15);
   //gPad->SetRightMargin(0.05);
@@ -96,9 +101,9 @@ void Secondaries()
   for (Int_t ptBin=1; ptBin<=secondaries->GetNbinsZ(); ptBin++)
   //for (Int_t ptBin = 1; ptBin<=2; ptBin++)
   {
-    TGraph* graph = new TGraph;
-    graph->Clear();
-    graph->SetTitle(Form("%f < p_{T} < %f", secondaries->GetZaxis()->GetBinLowEdge(ptBin), secondaries->GetZaxis()->GetBinUpEdge(ptBin)));
+    TH1F* hist = new TH1F(Form("secondaries_%d", ptBin), Form("secondaries_%d", ptBin), secondaries->GetNbinsY(), secondaries->GetYaxis()->GetXmin(), secondaries->GetYaxis()->GetXmax());
+
+    hist->SetTitle(Form("%f < p_{T} < %f", secondaries->GetZaxis()->GetBinLowEdge(ptBin), secondaries->GetZaxis()->GetBinUpEdge(ptBin)));
 
     for (Int_t cBin=1; cBin<=secondaries->GetNbinsY(); ++cBin)
     {
@@ -119,58 +124,307 @@ void Secondaries()
       printf("%f %f\n", sum, count);
 
       if (count > 0)
-        graph->SetPoint(graph->GetN(), secondaries->GetYaxis()->GetBinCenter(cBin), sum / count);
+      {
+        hist->SetBinContent(cBin, sum / count);
+        hist->SetBinError(cBin, TMath::Sqrt(sum) / count);
+      }
     }
 
+    hist->Scale(1.0 / hist->GetBinContent(hist->GetXaxis()->FindBin(1)));
+    hist->Add(new TF1("flat", "-1", 0, 2));
+
     new TCanvas;
-    graph->SetMarkerStyle(21);
-    graph->Draw("AP");
-    graph->Print();
+    hist->SetMarkerStyle(21);
+    hist->Draw("");
   }
 }
 
-void Composition()
+void Track2Particle1DComposition(const char* fileName = "correction_map.root", Int_t folderCount, const char** folderNames, Float_t upperPtLimit = 9.9)
 {
   gSystem->Load("libPWG0base");
 
-  AlidNdEtaCorrection* fdNdEtaCorrection[4];
+  TString canvasName;
+  canvasName.Form("Track2Particle1DComposition");
+  TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
+  canvas->Divide(3, 1);
 
-  TFile::Open("systematics.root");
+  TLegend* legend = new TLegend(0.7, 0.7, 0.95, 0.95);
 
-  for (Int_t i=0; i<4; ++i)
+  for (Int_t i=0; i<folderCount; ++i)
   {
-    TString name;
-    name.Form("correction_%d", i);
-    fdNdEtaCorrection[i] = new AlidNdEtaCorrection(name, name);
-    fdNdEtaCorrection[i]->LoadHistograms("systematics.root", name);
+    Track2Particle1DCreatePlots(fileName, folderNames[i], upperPtLimit);
+
+    TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject("gene_nTrackToNPart_x_div_meas_nTrackToNPart_x"));
+    TH1* corrY = dynamic_cast<TH1*> (gROOT->FindObject("gene_nTrackToNPart_y_div_meas_nTrackToNPart_y"));
+    TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject("gene_nTrackToNPart_z_div_meas_nTrackToNPart_z"));
+
+    Prepare1DPlot(corrX);
+    Prepare1DPlot(corrY);
+    Prepare1DPlot(corrZ);
+
+    const char* title = "Track2Particle Correction";
+    corrX->SetTitle(title);
+    corrY->SetTitle(title);
+    corrZ->SetTitle(title);
+
+    corrZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
+
+    corrX->SetLineColor(i+1);
+    corrY->SetLineColor(i+1);
+    corrZ->SetLineColor(i+1);
+
+    canvas->cd(1);
+    InitPad();
+    corrX->DrawCopy(((i>0) ? "SAME" : ""));
+
+    canvas->cd(2);
+    InitPad();
+    corrY->DrawCopy(((i>0) ? "SAME" : ""));
+
+    canvas->cd(3);
+    InitPad();
+    corrZ->DrawCopy(((i>0) ? "SAME" : ""));
+
+    legend->AddEntry(corrZ, folderNames[i]);
   }
 
-  //fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetMeasuredHistogram()->Scale(2);
-  //fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetGeneratedHistogram()->Scale(2);
+  legend->Draw();
+
+  //canvas->SaveAs(Form("Track2Particle1D_%s_%d_%f.gif", fileName, gMax, upperPtLimit));
+  //canvas->SaveAs(Form("Track2Particle1D_%s_%d_%f.eps", fileName, gMax, upperPtLimit));
+}
+
+void DrawDifferentSpecies()
+{
+  gROOT->ProcessLine(".L drawPlots.C");
+
+  const char* folderNames[] = { "correction_0", "correction_1", "correction_2", "correction_3" };
+
+  Track2Particle1DComposition("systematics.root", 4, folderNames);
+}
 
-  AlidNdEtaCorrection* finalCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
+void ScalePtDependent(TH3F* hist, TF1* function)
+{
+  // assumes that pt is the third dimension of hist
+  // scales with function(pt)
 
-  TList* collection = new TList;
+  for (Int_t z=1; z<=hist->GetNbinsZ(); ++z)
+  {
+    Double_t factor = function->Eval(hist->GetZaxis()->GetBinCenter(z));
+    printf("z = %d, pt = %f, scaling with %f\n", z, hist->GetZaxis()->GetBinCenter(z), factor);
+
+    for (Int_t x=1; x<=hist->GetNbinsX(); ++x)
+      for (Int_t y=1; y<=hist->GetNbinsY(); ++y)
+        hist->SetBinContent(x, y, z, hist->GetBinContent(x, y, z) * factor);
+  }
+}
+
+void ChangeComposition(void** correctionPointer, Int_t index)
+{
+  AlidNdEtaCorrection** fdNdEtaCorrection = (AlidNdEtaCorrection**) correctionPointer;
+
+  switch (index)
+  {
+    case 0:
+      break;
+
+    case 1:
+      fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetMeasuredHistogram()->Scale(10);
+      fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetGeneratedHistogram()->Scale(10);
+      break;
+
+    case 2:
+      TF1* ptDependence = new TF1("simple", "x", 0, 100);
+      ScalePtDependent(fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetGeneratedHistogram(), ptDependence);
+      ScalePtDependent(fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetMeasuredHistogram(), ptDependence);
+      break;
+
+  }
+}
+
+void Composition()
+{
+  gSystem->Load("libPWG0base");
+
+  gSystem->Unlink("new_compositions.root");
+
+  for (Int_t comp = 0; comp < 3; ++comp)
+  {
+    AlidNdEtaCorrection* fdNdEtaCorrection[4];
+
+    TFile::Open("systematics.root");
+
+    for (Int_t i=0; i<4; ++i)
+    {
+      TString name;
+      name.Form("correction_%d", i);
+      fdNdEtaCorrection[i] = new AlidNdEtaCorrection(name, name);
+      fdNdEtaCorrection[i]->LoadHistograms("systematics.root", name);
+    }
+
+    ChangeComposition(fdNdEtaCorrection, comp);
+
+    Double_t geneCount[5];
+    Double_t measCount[5];
+    geneCount[4] = 0;
+    measCount[4] = 0;
+
+    for (Int_t i=0; i<4; ++i)
+    {
+      geneCount[i] = fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetGeneratedHistogram()->Integral();
+      measCount[i] = fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetMeasuredHistogram()->Integral();
+
+      geneCount[4] += geneCount[i];
+      measCount[4] += measCount[i];
+
+      printf("Particle %s: %d gene, %d meas\n", ((i == 0) ? "pi" : (i == 1) ? "K" : (i == 2) ? "p" : "others"), (Int_t) geneCount[i], (Int_t) measCount[i]);
+    }
 
-  for (Int_t i=0; i<4; ++i)
-    collection->Add(fdNdEtaCorrection[i]);
+    printf("Generated ratios are:     %f pi, %f K, %f p, %f others\n", geneCount[0] / geneCount[4], geneCount[1] / geneCount[4], geneCount[2] / geneCount[4], geneCount[3] / geneCount[4]);
 
-  finalCorrection->Merge(collection);
+    printf("Reconstructed ratios are: %f pi, %f K, %f p, %f others\n", measCount[0] / measCount[4], measCount[1] / measCount[4], measCount[2] / measCount[4], measCount[3] / measCount[4]);
 
-  delete collection;
+    TList* collection = new TList;
 
-  finalCorrection->Finish();
+    for (Int_t i=0; i<4; ++i)
+      collection->Add(fdNdEtaCorrection[i]);
 
-  TFile* file = TFile::Open("temp.root", "RECREATE");
-  finalCorrection->SaveHistograms();
-  file->Write();
-  file->Close();
+    TString correctionName;
+    correctionName.Form("new_composition_%d", comp);
+
+    AlidNdEtaCorrection* newComposition = new AlidNdEtaCorrection(correctionName, correctionName);
+    newComposition->Merge(collection);
+    newComposition->Finish();
+
+    delete collection;
+
+    TFile* file = TFile::Open("new_compositions.root", "UPDATE");
+    newComposition->SaveHistograms();
+    //file->Write();
+    file->Close();
+  }
 
   gROOT->ProcessLine(".L drawPlots.C");
-  Track2Particle1D("temp.root");
+
+  const char* folderNames[] = { "new_composition_0", "new_composition_1", "new_composition_2" };
+
+  Track2Particle1DComposition("new_compositions.root", 3, folderNames);
+}
+
+Double_t ConvSigma1To2D(Double_t sigma)
+{
+  return TMath::Sqrt( - TMath::Log( 1 - TMath::Erf(sigma / TMath::Sqrt(2)) ) * 2);
+}
+
+Double_t ConvDistance1DTo2D(Double_t distance)
+{
+  return TMath::ErfInverse(1 - TMath::Exp(-distance * distance / 2)) * TMath::Sqrt(2);
+}
+
+Double_t Sigma2VertexCount(TH2F* tracks, Double_t nSigma)
+{
+  Double_t count = 0;
+
+  //nSigma = ConvSigma1To2D(nSigma);
+
+  for (Int_t x=1; x<=tracks->GetNbinsX(); ++x)
+    for (Int_t y=1; y<=tracks->GetNbinsY(); ++y)
+    {
+      Double_t impactX = tracks->GetXaxis()->GetBinCenter(x);
+      Double_t impactY = tracks->GetYaxis()->GetBinCenter(y);
+
+      Float_t d = TMath::Sqrt(impactX*impactX + impactY*impactY);
+
+      d = ConvDistance1DTo2D(d);
+
+      if (d < nSigma)
+        count += tracks->GetBinContent(x, y);
+    }
+
+  return count;
+}
+
+TH2F* Sigma2VertexGaussianTracksHist()
+{
+  TH2F* tracks = new TH2F("Sigma2Vertex_tracks", "Sigma2Vertex_tracks", 200, -5, 5, 200, -5, 5);
+
+  TF2* gaussian2D = new TF2("gaussian2D", "xgausn(0) * ygausn(3)", -5, 5, -5, 5);
+  gaussian2D->SetParameters(1, 0, 1, 1, 0, 1);
+
+  for (Int_t x=1; x<=tracks->GetNbinsX(); ++x)
+    for (Int_t y=1; y<=tracks->GetNbinsY(); ++y)
+      tracks->SetBinContent(x, y, gaussian2D->Eval(tracks->GetXaxis()->GetBinCenter(x), tracks->GetYaxis()->GetBinCenter(y)));
+
+  //normalize
+  tracks->Scale(1.0 / tracks->Integral());
+
+  return tracks;
+}
+
+void Sigma2VertexGaussian()
+{
+  TH2F* tracks = Sigma2VertexGaussianTracksHist();
+
+  TCanvas* canvas = new TCanvas("Sigma2VertexGaussian", "Sigma2VertexGaussian", 1000, 1000);
+  canvas->Divide(2, 2);
+
+  canvas->cd(1);
+  tracks->Draw("COLZ");
+
+  TH1F* ratio = new TH1F("Sigma2Vertex_ratio", "Sigma2Vertex_ratio;n sigma;included", 10, 0.25, 5.25);
+  for (Double_t nSigma = 0.5; nSigma < 5.1; nSigma += 0.5)
+    ratio->Fill(nSigma, Sigma2VertexCount(tracks, nSigma));
+  ratio->SetMarkerStyle(21);
+
+  canvas->cd(2);
+  ratio->Draw("P");
+
+  TH1F* ratio2 = new TH1F("Sigma2Vertex_ratio2", "Sigma2Vertex_ratio2;nSigma;% included 3 sigma / % included n sigma", 10, 0.25, 5.25);
+  Double_t sigma3 = Sigma2VertexCount(tracks, 3);
+  for (Double_t nSigma = 0.5; nSigma < 5.1; nSigma += 0.5)
+    ratio2->Fill(nSigma, sigma3 / Sigma2VertexCount(tracks, nSigma));
+  ratio2->SetMarkerStyle(21);
+
+  canvas->cd(3);
+  ratio2->Draw("P");
+
+  canvas->SaveAs("Sigma2Vertex.eps");
+}
+
+void Sigma2VertexSimulation()
+{
+  TFile* file = TFile::Open("systematics.root");
+
+  TH1F* sigmavertex = dynamic_cast<TH1F*> (file->Get("fSigmaVertex"));
+  if (!sigmavertex)
+  {
+    printf("Could not read histogram\n");
+    return;
+  }
+
+  TH1F* ratio = new TH1F("sigmavertex_ratio", "sigmavertex_ratio;Nsigma;% included 3 sigma / % included n sigma", sigmavertex->GetNbinsX(), sigmavertex->GetXaxis()->GetXmin(), sigmavertex->GetXaxis()->GetXmax());
+
+  for (Int_t i=1; i<=sigmavertex->GetNbinsX(); ++i)
+    ratio->SetBinContent(i, sigmavertex->GetBinContent(sigmavertex->GetXaxis()->FindBin(3)) / sigmavertex->GetBinContent(i));
+
+  TCanvas* canvas = new TCanvas("Sigma2VertexSimulation", "Sigma2VertexSimulation", 1000, 500);
+  canvas->Divide(2, 1);
+
+  canvas->cd(1);
+  sigmavertex->SetMarkerStyle(21);
+  sigmavertex->Draw("P");
+
+  canvas->cd(2);
+  ratio->SetMarkerStyle(21);
+  ratio->Draw("P");
 }
 
+
 void drawSystematics()
 {
-  Composition();
+  //Secondaries();
+  //DrawDifferentSpecies();
+  //Composition();
+
+  Sigma2VertexSimulation();
 }