updated macros for starting correction code and plotting

diff --git a/PWG4/macros/underlyingevent/correct.C b/PWG4/macros/underlyingevent/correct.C
index 48c1f083fe506c2600ee978e9f64dd2e2da7c7e4..095c6e8099155f8829c68b665cc824211c03dc60 100644 (file)
--- a/PWG4/macros/underlyingevent/correct.C
+++ b/PWG4/macros/underlyingevent/correct.C
@@ -127,7 +127,7 @@ void loadlibs()
   gSystem->Load("libPWG4JetTasks");
 }
 
-void CompareBias(const char* mcFile = "PWG4_JetTasksOutput.root", Int_t region)
+void CompareBias(const char* mcFile = "PWG4_JetTasksOutput.root", Int_t region, Int_t ueHist)
 {
   loadlibs();
 
@@ -138,11 +138,22 @@ void CompareBias(const char* mcFile = "PWG4_JetTasksOutput.root", Int_t region)
   h->SetPtRange(0.51, 100);
   h->SetCombineMinMax(kTRUE);
   
-  biasFromData = (TH1*) h->GetNumberDensitypT()->GetBias(AliUEHist::kCFStepBiasStudy, AliUEHist::kCFStepReconstructed, region, "z")->Clone("biasFromData");
-  biasFromData2 = (TH1*) h->GetNumberDensitypT()->GetBias(AliUEHist::kCFStepBiasStudy2, AliUEHist::kCFStepReconstructed, region, "z")->Clone("biasFromData2");
-  //biasFromData = (TH1*) h->GetNumberDensitypT()->GetBias(AliUEHist::kCFStepReconstructed, AliUEHist::kCFStepRealLeading, region, "z")->Clone("biasFromData");
-  biasFromMC   = (TH1*) h->GetNumberDensitypT()->GetBias(AliUEHist::kCFStepReconstructed, AliUEHist::kCFStepTracked, region, "z")->Clone("biasFromMC");
-  //biasFromMC   = (TH1*) h->GetNumberDensitypT()->GetBias(AliUEHist::kCFStepReconstructed, AliUEHist::kCFStepRealLeading, region, "z")->Clone("biasFromMC");
+  const char* axis = "z";
+  Float_t ptLeadMin = 0;
+  Float_t ptLeadMax = -1;
+  
+  if (ueHist == 2)
+  {
+    ptLeadMin = 0.51 + 4;
+    ptLeadMax = 1.99 + 4;
+    axis = "y";
+  }
+  
+  biasFromData = (TH1*) h->GetUEHist(ueHist)->GetBias(AliUEHist::kCFStepBiasStudy, AliUEHist::kCFStepReconstructed, region, axis, ptLeadMin, ptLeadMax)->Clone("biasFromData");
+  biasFromData2 = (TH1*) h->GetUEHist(ueHist)->GetBias(AliUEHist::kCFStepBiasStudy2, AliUEHist::kCFStepReconstructed, region, axis, ptLeadMin, ptLeadMax)->Clone("biasFromData2");
+  //biasFromData = (TH1*) h->GetUEHist(ueHist)->GetBias(AliUEHist::kCFStepReconstructed, AliUEHist::kCFStepRealLeading, region, axis)->Clone("biasFromData");
+  biasFromMC   = (TH1*) h->GetUEHist(ueHist)->GetBias(AliUEHist::kCFStepReconstructed, AliUEHist::kCFStepTracked, region, axis, ptLeadMin, ptLeadMax)->Clone("biasFromMC");
+  //biasFromMC   = (TH1*) h->GetUEHist(ueHist)->GetBias(AliUEHist::kCFStepReconstructed, AliUEHist::kCFStepRealLeading, region, axis)->Clone("biasFromMC");
   
 /*  biasFromData->Draw();
   biasFromMC->SetLineColor(2);
@@ -232,7 +243,7 @@ void SetupRanges(void* obj)
   ((AliUEHistograms*) obj)->SetCombineMinMax(kTRUE);
 }
 
-void DrawRatios(const char* name, void* correctedVoid, void* comparisonVoid, Int_t compareStep = -1, Int_t compareRegion = 2)
+void DrawRatios(const char* name, void* correctedVoid, void* comparisonVoid, Int_t compareStep = -1, Int_t compareRegion = 2, Float_t ptLeadMin = -1, Float_t ptLeadMax = -1)
 {
   AliUEHist* corrected = (AliUEHist*) correctedVoid;
   AliUEHist* comparison = (AliUEHist*) comparisonVoid;
@@ -259,28 +270,70 @@ void DrawRatios(const char* name, void* correctedVoid, void* comparisonVoid, Int
   {
     for (Int_t region=beginRegion; region <= endRegion; region++)
     {
-      TH1* corrHist = corrected->GetUEHist(step, region);
-      TH1* mcHist   = comparison->GetUEHist(step, region);
+      Printf("%f %f", ptLeadMin, ptLeadMax);
+      TH1* corrHist = corrected->GetUEHist(step, region, ptLeadMin, ptLeadMax);
+      TH1* mcHist   = comparison->GetUEHist(step, region, ptLeadMin, ptLeadMax);
       
-      DrawRatio(corrHist, mcHist, Form("%s: step %d %s %s", name, step, corrected->GetStepTitle(step), corrected->GetRegionTitle(region)));
+      DrawRatio(corrHist, mcHist, TString(Form("%s: step %d %s %s", name, step, corrected->GetStepTitle(step), corrected->GetRegionTitle(region))));
     }
   }
 }
 
-void DrawRatios(void* correctedVoid, void* comparisonVoid, Int_t compareStep = -1, Int_t compareRegion = 2)
+void DrawRatios(void* correctedVoid, void* comparisonVoid, Int_t compareStep = -1, Int_t compareRegion = 2, Int_t compareUEHist = 0)
 {
   AliUEHistograms* corrected = (AliUEHistograms*) correctedVoid;
   AliUEHistograms* comparison = (AliUEHistograms*) comparisonVoid;
 
-  DrawRatios("Number density", corrected->GetNumberDensitypT(), comparison->GetNumberDensitypT(), compareStep, compareRegion);
-  if (compareStep != -1)
+  if (compareUEHist == 2)
   {
-    DrawRatios("Pt sum", corrected->GetSumpT(), comparison->GetSumpT(), compareStep, compareRegion);
-    //DrawRatios("Phi correlation", corrected->GetNumberDensityPhi(), comparison->GetNumberDensityPhi(), compareStep, compareRegion);
+    for (Float_t ptLeadMin = 2.51; ptLeadMin < 4; ptLeadMin += 2)
+      DrawRatios(TString(Form("UE %d pT %f", compareUEHist, ptLeadMin)), corrected->GetUEHist(compareUEHist), comparison->GetUEHist(compareUEHist), compareStep, compareRegion, ptLeadMin, ptLeadMin + 1.98);      
+    return;
   }
+
+  if (compareUEHist == -1)
+  {
+    for (Int_t i=0; i<2; i++)
+      DrawRatios(TString(Form("UE %d", i)), corrected->GetUEHist(i), comparison->GetUEHist(i), compareStep, compareRegion);
+  }
+  else
+    DrawRatios(TString(Form("UE %d", compareUEHist)), corrected->GetUEHist(compareUEHist), comparison->GetUEHist(compareUEHist), compareStep, compareRegion);
 }
 
-void correctMC(const char* fileNameCorrections, const char* fileNameESD = 0, Int_t compareStep = -1, Int_t compareRegion = 2)
+void CompareEventsTracks(void* corrVoid, void* mcVoid, Int_t compareStep, Int_t compareRegion, Int_t compareUEHist = 0)
+{
+  AliUEHist* corr = ((AliUEHistograms*) corrVoid)->GetUEHist(compareUEHist);
+  AliUEHist* mc   = ((AliUEHistograms*) mcVoid)->GetUEHist(compareUEHist);
+
+  Float_t ptLeadMin = 0;
+  Float_t ptLeadMax = -1;
+  Int_t axis = 2;
+  
+  if (compareUEHist == 2)
+  {
+    ptLeadMin = 1.01;
+    ptLeadMax = 1.49;
+    axis = 4;
+  }
+
+  TH1* corrHist = corr->GetUEHist(compareStep, compareRegion, ptLeadMin, ptLeadMax);
+  TH1* mcHist   = mc  ->GetUEHist(compareStep, compareRegion, ptLeadMin, ptLeadMax);
+  
+  DrawRatio(corrHist, mcHist, Form("check"));
+  
+  corr->SetBinLimits(corr->GetTrackHist(compareRegion)->GetGrid(compareStep));
+  mc->SetBinLimits(corr->GetTrackHist(compareRegion)->GetGrid(compareStep));
+
+  corrHist =  corr->GetTrackHist(compareRegion)->GetGrid(compareStep)->Project(axis);
+  mcHist   =  mc  ->GetTrackHist(compareRegion)->GetGrid(compareStep)->Project(axis);
+  DrawRatio(corrHist, mcHist, Form("check2"));
+  
+  corrHist =  corr->GetEventHist()->GetGrid(compareStep)->Project(0);
+  mcHist   =  mc  ->GetEventHist()->GetGrid(compareStep)->Project(0);
+  DrawRatio(corrHist, mcHist, Form("check3"));
+}
+
+void correctMC(const char* fileNameCorrections, const char* fileNameESD = 0, Int_t compareStep = -1, Int_t compareRegion = 2, Int_t compareUEHist = 0)
 {
   // corrects the reconstructed step in fileNameESD with fileNameCorr
   // if fileNameESD is 0 data from fileNameCorr is taken
@@ -309,12 +362,14 @@ void correctMC(const char* fileNameCorrections, const char* fileNameESD = 0, Int
   
   esd->Correct(corr);
   
-  DrawRatios(esd, testSample, compareStep, compareRegion);
+  DrawRatios(esd, testSample, compareStep, compareRegion, compareUEHist);
+  
+  //CompareEventsTracks(esd, testSample, compareStep, compareRegion, compareUEHist);
 }
 
 // function to compare only final step for all regions and distributions
 
-void correctData(const char* fileNameCorrections, const char* fileNameESD, Int_t compareStep = -1, Int_t compareRegion = 2)
+void correctData(const char* fileNameCorrections, const char* fileNameESD, Int_t compareStep = -1, Int_t compareRegion = 2, Int_t compareUEHist = 0)
 {
   // corrects fileNameESD with fileNameCorrections and compares the two
   
@@ -337,5 +392,168 @@ void correctData(const char* fileNameCorrections, const char* fileNameESD, Int_t
   esd->Write();
   file3->Close();
   
-  DrawRatios(esd, corr, compareStep, compareRegion);
+  DrawRatios(esd, corr, compareStep, compareRegion, compareUEHist);
 }
+
+void ITSTPCEfficiency(const char* fileNameData, Int_t itsTPC = 0)
+{
+  // its = 0; tpc = 1
+
+  // uncertainty from dN/dpT paper
+  Double_t pTBins[] =  {0.0, 0.1, 0.15,  0.2,  0.25,  0.3,   0.35,  0.4,   0.45,  0.5,   0.6,   0.7,   0.8,   0.9,   1.0,   1.5,   2.0,   2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0, 50.0, 100.0};
+  Float_t effITS[] = {0.,  0.,  0.995, 0.98, 0.986, 0.996, 0.999, 0.997, 0.997, 1.008, 1.005, 1.022, 1.021, 1.021, 1.035, 1.04,  1.03  };  // the last three are the same because i don't have entries
+  Float_t effTPC[] = {0.,  0,   1.042, 1.026,1.021, 1.018, 1.015, 1.015, 1.012, 1.012, 1.007, 1.0075,1.006, 1.006, 1.004, 1.004, 1.009 }; // the last bins put as if they were the same
+
+  TH1F* effHist = new TH1F("effHist", "effHist", 39, pTBins);
+  for (Int_t i=0; i<39; i++)
+  {
+    Int_t bin = i;
+    if (i > 16)
+      bin = 16;
+    effHist->SetBinContent(i+1, (itsTPC == 0) ? effITS[bin] : effTPC[bin]);
+  }
+
+  new TCanvas; effHist->Draw();
+
+  EffectOfModifiedTrackingEfficiency(fileNameData, effHist, (itsTPC == 0) ? "ITS" : "TPC");
+} 
+
+
+void EffectOfModifiedTrackingEfficiency(const char* fileNameData, TH1F* trackingEff, const char* name = "EffectOfModifiedTrackingEfficiency")
+{
+  // trackingEff should contain the change in tracking efficiency, i.e. between before and after in the eta-pT plane
+
+  loadlibs();
+  
+  TFile::Open(fileNameData);
+  list = (TList*) gFile->Get(objectName);
+  AliUEHistograms* corrected = (AliUEHistograms*) list->FindObject("AliUEHistograms");
+  corrected->SetEtaRange(-0.79, 0.79);
+  corrected->SetPtRange(0.51, 100);
+  corrected->SetCombineMinMax(kTRUE);
+  
+  AliUEHist* ueHist = corrected->GetUEHist(0);
+  Int_t region = 2;
+  
+  // histogram before
+  TH1* before = ueHist->GetUEHist(AliUEHist::kCFStepAll, region);
+
+  // copy histogram
+  // the CFStepTriggered step is overwritten here and cannot be used for comparison afterwards anymore
+  ueHist->CorrectTracks(AliUEHist::kCFStepAll, AliUEHist::kCFStepTriggered, (TH1*) 0, 0, -1);
+
+  // reapply tracking efficiency
+  ueHist->CorrectTracks(AliUEHist::kCFStepTriggered, AliUEHist::kCFStepAll, trackingEff, 1);
+
+  // histogram after
+  TH1* after = ueHist->GetUEHist(AliUEHist::kCFStepAll, region);
+  
+  DrawRatio(before, after, name);
+}
+
+void ModifyComposition(const char* fileNameData, const char* fileNameCorrections, Int_t id, Bool_t verbose = kFALSE)
+{
+  loadlibs();
+  
+  TFile::Open(fileNameData);
+  list = (TList*) gFile->Get(objectName);
+  AliUEHistograms* corrected = (AliUEHistograms*) list->FindObject("AliUEHistograms");
+  corrected->SetEtaRange(-0.79, 0.79);
+  corrected->SetPtRange(0.51, 100);
+  corrected->SetCombineMinMax(kTRUE);
+  
+  file2 = TFile::Open(fileNameCorrections);
+  list2 = (TList*) file2->Get(objectName);
+  AliUEHistograms* corrections = (AliUEHistograms*) list2->FindObject("AliUEHistograms");
+  corrections->SetEtaRange(-0.79, 0.79);
+  corrections->SetPtRange(0.51, 100);
+  corrections->SetCombineMinMax(kTRUE);
+  
+  ueHist = (AliUEHist*) corrections->GetUEHist(id);
+  
+  // copy histogram
+  // the CFStepTriggered step is overwritten here and cannot be used for comparison afterwards anymore
+  ueHist->CorrectTracks(AliUEHist::kCFStepAll, AliUEHist::kCFStepTriggered, (TH1*) 0, 0);
+  
+  // histogram before
+  TH1* before[3];
+  for (Int_t region=0; region<3; region++)
+    before[region] = ueHist->GetUEHist(AliUEHist::kCFStepAll, region);
+  
+  defaultEff = ueHist->GetTrackingEfficiency();
+  defaultEffpT = ueHist->GetTrackingEfficiency(1);
+  defaultContainer = ueHist->GetTrackHistEfficiency();
+  
+  c = new TCanvas;
+  defaultEffpT->Draw("");
+  
+  for (Int_t i=0; i<7; i++)
+  {
+    // case 0: // default
+    // case 1: // + 30% kaons
+    // case 2: // - 30% kaons
+    // case 3: // + 30% protons
+    // case 4: // - 30% protons
+    // case 5: // + 30% others
+    // case 6: // - 30% others
+    Int_t correctionIndex = (i+1) / 2 + 1; // bin 1 == protons, bin 2 == kaons, ...
+    Double_t scaleFactor = (i % 2 == 1) ? 1.3 : 0.7;
+    if (i == 0)
+      scaleFactor = 1;
+    
+    newContainer = (AliCFContainer*) defaultContainer->Clone();
+    
+    // modify, change all steps
+    for (Int_t j=0; j<newContainer->GetNStep(); j++)
+    {
+      THnSparse* grid = newContainer->GetGrid(j)->GetGrid();
+      
+      for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
+      {
+        Int_t bins[5];
+        Double_t value = grid->GetBinContent(binIdx, bins);
+        Double_t error = grid->GetBinError(binIdx);
+        
+        if (bins[2] != correctionIndex)
+          continue;
+    
+        value *= scaleFactor;
+        error *= scaleFactor;
+    
+        grid->SetBinContent(bins, value);
+        grid->SetBinError(bins, error);      
+      }
+    }
+    
+    // put in corrections
+    ueHist->SetTrackHistEfficiency(newContainer);
+    
+    // ratio
+    modifiedEff = ueHist->GetTrackingEfficiency();
+    modifiedEff->Divide(modifiedEff, defaultEff);
+    //modifiedEff->Draw("COLZ");
+    
+    c->cd();
+    modifiedEffpT = ueHist->GetTrackingEfficiency(1);
+    modifiedEffpT->SetLineColor(i+1);
+    modifiedEffpT->Draw("SAME");
+    
+    // apply change in tracking efficiency
+    ueHist->CorrectTracks(AliUEHist::kCFStepTriggered, AliUEHist::kCFStepAll, modifiedEff, 0, 1);
+  
+    for (Int_t region=0; region<3; region++)
+    {
+      // histogram after
+      TH1* after = ueHist->GetUEHist(AliUEHist::kCFStepAll, region);
+      
+      if (verbose)
+        DrawRatio(before[region], (TH1*) after->Clone(), Form("Region %d Composition %d", region, i));
+      
+      // ratio is flat, extract deviation
+      after->Divide(before[region]);
+      after->Fit("pol0", "0");
+      Printf("Deviation for region %d case %d is %.2f %%", region, i, 100.0 - 100.0 * after->GetFunction("pol0")->GetParameter(0));
+    }
+    //return;
+  }
+}    

diff --git a/PWG4/macros/underlyingevent/plots.C b/PWG4/macros/underlyingevent/plots.C
new file mode 100644 (file)
index 0000000..e098ba8
--- /dev/null
+++ b/PWG4/macros/underlyingevent/plots.C
@@ -0,0 +1,124 @@
+void loadlibs()
+{
+  gSystem->Load("libANALYSIS");
+  gSystem->Load("libANALYSISalice");
+  gSystem->Load("libCORRFW");
+  gSystem->Load("libJETAN");
+  gSystem->Load("libPWG4JetTasks");
+}
+
+void PlotTrackingEfficiencies()
+{
+  TGrid::Connect("alien://");
+
+  loadlibs();
+
+  ifstream in;
+  in.open("list");
+
+  Int_t i=0;
+
+  legend = new TLegend(0.5, 0.1, 0.9, 0.5);
+  legend->SetNColumns(2);
+
+  TString line;
+  while (in.good())
+  {
+    in >> line;
+
+    if (line.Length() == 0)
+      continue;
+
+    TString fileName;
+    fileName.Form("%s/%s/PWG4_JetTasksOutput.root", "alien:///alice/cern.ch/user/k/kleinb/analysisESD/LHC10d4/output_pwg4train_LHC10d4_101001a", line.Data());
+    Printf("%s", fileName.Data());
+    file = TFile::Open(fileName);
+    if (!file)
+      continue;
+    list = (TList*) file->Get("PWG4_LeadingTrackUE/histosLeadingTrackUE");
+    AliUEHistograms* h = (AliUEHistograms*) list->FindObject("AliUEHistograms");
+    h->SetEtaRange(-0.79, 0.79);
+
+    TH1::AddDirectory(0);
+    eff = (TH1*) h->GetNumberDensitypT()->GetTrackingEfficiency(1)->Clone(Form("trackingefficiency%s", line.Data()));
+
+    delete h;
+    //delete list;
+    file->Close();
+
+    eff->SetLineColor(i+1);
+    eff->GetXaxis()->SetRangeUser(0, 10);
+    eff->GetYaxis()->SetRangeUser(0.6, 0.8);
+    eff->Draw((i == 0) ? "" : "SAME");
+ 
+    legend->AddEntry(eff, line);
+
+    i++;
+    //if (i > 1)
+    //  break;
+  }
+
+  legend->Draw();
+
+  in.close();
+}
+
+void TrackingEfficiencySpecies(const char* fileName = "PWG4_JetTasksOutput.root")
+{
+  loadlibs();
+
+  file = TFile::Open(fileName);
+  list = (TList*) file->Get("PWG4_LeadingTrackUE/histosLeadingTrackUE");
+  AliUEHistograms* h = (AliUEHistograms*) list->FindObject("AliUEHistograms");
+  h->SetEtaRange(-0.79, 0.79);
+
+  eff = (TH2*) h->GetNumberDensitypT()->GetTrackEfficiency(AliUEHist::kCFStepAnaTopology, AliUEHist::kCFStepTrackedOnlyPrim, 1, 2);
+  eff->Draw("colz");
+
+  for (Int_t i=0; i<4; i++)
+  {
+    proj = eff->ProjectionX(Form("p%d", i), i+1, i+1);
+    proj->SetLineColor(i+1);
+    proj->Draw((i==0) ? "" : "SAME");
+  }
+}
+
+void CheckTrackingEfficiency(const char* reference = 0, const char* fileName = "PWG4_JetTasksOutput.root", Bool_t all = kFALSE)
+{
+  loadlibs();
+
+  file = TFile::Open(fileName);
+  list = (TList*) file->Get("PWG4_LeadingTrackUE/histosLeadingTrackUE");
+  tree = (TTree*) list->FindObject("UEAnalysisSettings");
+  
+  if (reference)
+  {
+    file2 = TFile::Open(reference);
+    refHist = (TH1D*) file2->Get("trackingefficiency");
+  }
+  else
+  {
+    AliUEHistograms* h = (AliUEHistograms*) list->FindObject("AliUEHistograms");
+    h->SetEtaRange(-0.79, 0.79);
+    refHist = h->GetNumberDensitypT()->GetTrackingEfficiency(1);
+  }
+
+  refHist->SetLineWidth(3);
+  refHist->GetListOfFunctions()->Clear();
+  refHist->Draw();
+
+  TH1D* hist = 0;
+  
+  tree->SetBranchAddress("fkTrackingEfficiency", &hist);
+  
+  for (Int_t i=0; i<tree->GetEntries(); i++)
+  {
+    tree->GetEntry(i);
+    hist->GetListOfFunctions()->Clear();
+    hist->SetLineColor(i+2);
+    hist->DrawCopy("SAME");
+    if (!all)
+      break;
+  }
+}
+