adding flag for non-field data
authorjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 22 Oct 2009 11:54:31 +0000 (11:54 +0000)
committerjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 22 Oct 2009 11:54:31 +0000 (11:54 +0000)
adding bypass trigger to check cosmics

16 files changed:
PWG0/AliCorrection.cxx
PWG0/AliCorrection.h
PWG0/AliPWG0Helper.cxx
PWG0/AliPWG0Helper.h
PWG0/CreateStandardCuts.C
PWG0/dNdEta/AlidNdEtaCorrection.h
PWG0/dNdEta/AlidNdEtaCorrectionTask.cxx
PWG0/dNdEta/AlidNdEtaCorrectionTask.h
PWG0/dNdEta/AlidNdEtaTask.cxx
PWG0/dNdEta/AlidNdEtaTask.h
PWG0/dNdEta/dNdEtaAnalysis.cxx
PWG0/dNdEta/dNdEtaAnalysis.h
PWG0/dNdEta/drawPlots.C
PWG0/dNdEta/drawSystematics.C
PWG0/dNdEta/drawSystematicsNew.C
PWG0/dNdEta/run.C

index 4b9e9a7..85de079 100644 (file)
@@ -41,13 +41,13 @@ AliCorrection::AliCorrection(const Char_t* name, const Char_t* title, AliPWG0Hel
   Int_t nBinsPt = 0;
 
   // different binnings, better solution could be anticipated...
-  if (analysisMode == AliPWG0Helper::kTPC || analysisMode == AliPWG0Helper::kTPCITS)
+  if ((analysisMode & AliPWG0Helper::kTPC || analysisMode & AliPWG0Helper::kTPCITS) && (analysisMode & AliPWG0Helper::kFieldOn))
   {
     static Float_t binLimitsPtTmp[] = {0.0, 0.05, 0.1, 0.125, 0.15, 0.175, 0.2, 0.225, 0.25, 0.275, 0.3, 0.325, 0.35, 0.375, 0.4, 0.425, 0.45, 0.475, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.5, 2.0, 5.0, 10.0, 100.0};
     binLimitsPt = (Float_t*) binLimitsPtTmp;
     nBinsPt = 28;
   }
-  else if (analysisMode == AliPWG0Helper::kSPD)
+  else if (analysisMode & AliPWG0Helper::kSPD || !(analysisMode & AliPWG0Helper::kFieldOn))
   {
     static Float_t binLimitsPtTmp[] = {-0.5, 100.0};
     binLimitsPt = (Float_t*) binLimitsPtTmp;
@@ -93,7 +93,7 @@ AliCorrection::AliCorrection(const Char_t* name, const Char_t* title, AliPWG0Hel
 
   fEventCorr = new AliCorrectionMatrix2D("EventCorrection", Form("%s EventCorrection", fTitle.Data()), 18, binLimitsVtx, nBinsN, binLimitsN);
 
-  if (analysisMode == AliPWG0Helper::kSPD)
+  if (analysisMode & AliPWG0Helper::kSPD)
   {
     TH3F* dummyBinning = new TH3F("dummyBinning","dummyBinning",18, binLimitsVtx, 60, binLimitsEta, nBinsN2, binLimitsN2);
     fTrackCorr = new AliCorrectionMatrix3D("TrackCorrection", Form("%s TrackCorrection", fTitle.Data()), dummyBinning);
@@ -329,26 +329,38 @@ void AliCorrection::DrawOverview(const char* canvasName)
   if (canvasName)
     canvasNameTmp = canvasName;
 
-  TCanvas* canvas = new TCanvas(canvasNameTmp, canvasNameTmp, 1200, 800);
-  canvas->Divide(3, 2);
+  TCanvas* canvas = new TCanvas(canvasNameTmp, canvasNameTmp, 1200, 1200);
+  canvas->Divide(3, 3);
 
   if (fTrackCorr) {
     canvas->cd(1);
+    fTrackCorr->Get2DCorrectionHistogram("xy", 0, 0)->DrawCopy("COLZ")->GetZaxis()->SetRangeUser(0, 10);
+    
+    canvas->cd(2);
+    fTrackCorr->Get2DCorrectionHistogram("yz", 0, 0)->DrawCopy("COLZ")->GetZaxis()->SetRangeUser(0, 10);
+    
+    canvas->cd(3);
+    fTrackCorr->Get2DCorrectionHistogram("zx", 0, 0)->DrawCopy("COLZ")->GetZaxis()->SetRangeUser(0, 10);
+    
+    canvas->cd(4);
     fTrackCorr->Get1DCorrectionHistogram("x", 0.3, 5, -1, 1)->DrawCopy()->GetYaxis()->SetRangeUser(0, 10);
 
-    canvas->cd(2);
+    canvas->cd(5);
     fTrackCorr->Get1DCorrectionHistogram("y", 0.3, 5, 0, 0)->DrawCopy()->GetYaxis()->SetRangeUser(0, 10);
 
-    canvas->cd(3);
+    canvas->cd(6);
     fTrackCorr->Get1DCorrectionHistogram("z", 0, -1, -1, 1)->DrawCopy()->GetYaxis()->SetRangeUser(0, 10);
   }
 
   if (fEventCorr)
   {
-    canvas->cd(4);
+    canvas->cd(7);
+    fEventCorr->GetCorrectionHistogram()->DrawCopy("COLZ")->GetYaxis()->SetRangeUser(0, 30);
+    
+    canvas->cd(8);
     fEventCorr->Get1DCorrectionHistogram("x")->DrawCopy();
 
-    canvas->cd(5);
+    canvas->cd(9);
     fEventCorr->Get1DCorrectionHistogram("y")->DrawCopy()->GetXaxis()->SetRangeUser(0, 30);
   }
 }
@@ -449,6 +461,6 @@ void AliCorrection::PrintInfo(Float_t ptCut)
 
   Printf("Example centered bin: tracks measured: %.1f tracks generated: %.1f, events measured: %.1f, events generated %.1f", measured->GetBinContent(measured->GetNbinsX() / 2, measured->GetNbinsY() / 2, measured->GetNbinsZ() / 2), generated->GetBinContent(measured->GetNbinsX() / 2, measured->GetNbinsY() / 2, measured->GetNbinsZ() / 2), measuredEvents->GetBinContent(measuredEvents->GetNbinsX() / 2, measuredEvents->GetNbinsY() / 2), generatedEvents->GetBinContent(measuredEvents->GetNbinsX() / 2, measuredEvents->GetNbinsY() / 2));
 
-  PrintStats(10, 1.0, ptCut);
+  PrintStats(10, 0.8, ptCut);
   PrintStats(10, 1.5, ptCut);
 }
index ca038f6..9cf9ed6 100644 (file)
@@ -23,7 +23,7 @@ class AliCorrection : public TNamed
 {
 public:
   AliCorrection();
-  AliCorrection(const Char_t* name, const Char_t* title, AliPWG0Helper::AnalysisMode analysisMode = AliPWG0Helper::kTPC);
+  AliCorrection(const Char_t* name, const Char_t* title, AliPWG0Helper::AnalysisMode analysisMode = (AliPWG0Helper::AnalysisMode) (AliPWG0Helper::kTPC | AliPWG0Helper::kFieldOn));
   AliCorrection(const AliCorrection& c);
 
   virtual ~AliCorrection();
index ef7f859..a2af1ee 100644 (file)
@@ -122,6 +122,11 @@ Bool_t AliPWG0Helper::IsEventTriggered(ULong64_t triggerMask, Trigger trigger)
 
   switch (trigger)
   {
+    case kAcceptAll:
+    {
+      return kTRUE;
+      break;
+    }
     case kMB1:
     {
       if (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)))
@@ -160,11 +165,11 @@ Bool_t AliPWG0Helper::TestVertex(const AliESDVertex* vertex, AnalysisMode analys
     // Checks if a vertex meets the needed quality criteria
 
   Float_t requiredZResolution = -1;
-  if (analysisMode == kSPD || analysisMode == kTPCITS)
+  if (analysisMode & kSPD || analysisMode & kTPCITS)
   {
     requiredZResolution = 0.1;
   }
-  else if (analysisMode == kTPC)
+  else if (analysisMode & kTPC)
     requiredZResolution = 10.;
 
   // check resolution
@@ -188,13 +193,13 @@ const AliESDVertex* AliPWG0Helper::GetVertex(AliESDEvent* aEsd, AnalysisMode ana
   // also the quality criteria that are applied)
 
   const AliESDVertex* vertex = 0;
-  if (analysisMode == kSPD || analysisMode == kTPCITS)
+  if (analysisMode & kSPD || analysisMode & kTPCITS)
   {
     vertex = aEsd->GetPrimaryVertexSPD();
     if (debug)
       Printf("AliPWG0Helper::GetVertex: Returning SPD vertex");
   }
-  else if (analysisMode == kTPC)
+  else if (analysisMode & kTPC)
   {
     if(bRedoTPC){
       if (debug)
@@ -631,20 +636,29 @@ void AliPWG0Helper::PrintConf(AnalysisMode analysisMode, Trigger trigger)
   // Prints the given configuration
   //
 
-  TString str(">>>> Running with ");
+  TString str(">>>> Running with >");
+
+  if (analysisMode & kSPD)
+    str += "SPD-only";
+    
+  if (analysisMode & kTPC)
+     str += "TPC-only";
+    
+  if (analysisMode & kTPCITS)
+     str += "Global tracking";
 
-  switch (analysisMode)
+  if (analysisMode & kFieldOn)
   {
-    case kInvalid: str += "invalid setting"; break;
-    case kSPD : str += "SPD-only"; break;
-    case kTPC : str += "TPC-only"; break;
-    case kTPCITS : str += "Global tracking"; break;
+     str += " (with field)";
   }
-
-  str += " and trigger ";
+  else
+     str += " (WITHOUT field)";
+  
+  str += "< and trigger >";
 
   switch (trigger)
   {
+    case kAcceptAll : str += "ACCEPT ALL (bypass!)"; break;
     case kMB1 : str += "MB1"; break;
     case kMB2 : str += "MB2"; break;
     case kMB3 : str += "MB3"; break;
@@ -655,7 +669,7 @@ void AliPWG0Helper::PrintConf(AnalysisMode analysisMode, Trigger trigger)
     case kOfflineFASTOR : str += "Offline SPD FASTOR"; break;
   }
 
-  str += " <<<<";
+  str += "< <<<<";
 
   Printf("%s", str.Data());
 }
index ff7dde5..b90fca5 100644 (file)
@@ -21,8 +21,8 @@ class TTree;
 class AliPWG0Helper : public TObject
 {
   public:
-    enum Trigger { kMB1 = 0, kMB2, kMB3, kSPDFASTOR, kOfflineMB1, kOfflineMB2, kOfflineMB3, kOfflineFASTOR }; // MB1, MB2, MB3 definition from ALICE-INT-2005-025
-    enum AnalysisMode { kInvalid = -1, kSPD = 0, kTPC, kTPCITS };
+    enum Trigger { kAcceptAll = -1, kMB1 = 0, kMB2, kMB3, kSPDFASTOR, kOfflineMB1, kOfflineMB2, kOfflineMB3, kOfflineFASTOR }; // MB1, MB2, MB3 definition from ALICE-INT-2005-025
+    enum AnalysisMode { kInvalid = -1, kSPD = 0x1, kTPC = 0x2, kTPCITS = 0x4, kFieldOn = 0x8 };
     // in case we want to use bitmaps...
     enum MCProcessType { kInvalidProcess = -1, kND = 0x1, kDD = 0x2, kSD = 0x4 };
 
index fb6c760..1eecb0f 100644 (file)
@@ -2,7 +2,7 @@
 
 // this macro creates the track and event cuts used in this analysis
 
-AliESDtrackCuts* CreateTrackCuts(AliPWG0Helper::AnalysisMode analysisMode, Bool_t fieldOn = kTRUE, Bool_t hists = kTRUE)
+AliESDtrackCuts* CreateTrackCuts(AliPWG0Helper::AnalysisMode analysisMode, Bool_t hists = kTRUE)
 {
   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts");
 
@@ -23,10 +23,10 @@ AliESDtrackCuts* CreateTrackCuts(AliPWG0Helper::AnalysisMode analysisMode, Bool_
   TString tag("Global tracking");
 
   // TPC-only cuts
-  if (analysisMode == AliPWG0Helper::kTPC)
+  if (analysisMode & AliPWG0Helper::kTPC)
   {
-    cov1 = 9;
-    cov2 = 9;
+    cov1 = 1e10;
+    cov2 = 1e10;
     cov3 = 1e10;
     cov4 = 1e10;
     cov5 = 1e10;
@@ -36,9 +36,9 @@ AliESDtrackCuts* CreateTrackCuts(AliPWG0Helper::AnalysisMode analysisMode, Bool_
     
     tag = "TPC-only tracking";
   }
-
+  
   // cuts for data without field
-  if (!fieldOn)
+  if (!(analysisMode & AliPWG0Helper::kFieldOn))
   {
     cov5 = 1e10;
     tag += " without field";
@@ -51,8 +51,8 @@ AliESDtrackCuts* CreateTrackCuts(AliPWG0Helper::AnalysisMode analysisMode, Bool_
   if (sigmaToVertex) {
     esdTrackCuts->SetMaxNsigmaToVertex(nSigma);
   }
-  else{
-    
+  else
+  {
     esdTrackCuts->SetMaxDCAToVertexZ(3.2);
     esdTrackCuts->SetMaxDCAToVertexXY(2.4);
     esdTrackCuts->SetDCAToVertex2D(kTRUE);
@@ -62,9 +62,10 @@ AliESDtrackCuts* CreateTrackCuts(AliPWG0Helper::AnalysisMode analysisMode, Bool_
   esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
 
   esdTrackCuts->SetMinNClustersTPC(50);
-  esdTrackCuts->SetMaxChi2PerClusterTPC(3.5);
+  esdTrackCuts->SetMaxChi2PerClusterTPC(4);
 
   Printf("Created track cuts for: %s", tag.Data());
 
   return esdTrackCuts;
 }
+
index 875a324..dbefa84 100644 (file)
@@ -35,7 +35,7 @@ public:
   };
 
   AlidNdEtaCorrection();
-  AlidNdEtaCorrection(const Char_t* name, const Char_t* title, AliPWG0Helper::AnalysisMode analysis = AliPWG0Helper::kTPC);
+  AlidNdEtaCorrection(const Char_t* name, const Char_t* title, AliPWG0Helper::AnalysisMode analysis = (AliPWG0Helper::AnalysisMode) (AliPWG0Helper::kTPC | AliPWG0Helper::kFieldOn));
 
   virtual Long64_t Merge(TCollection* list);
 
index 377bd8d..4a6bc31 100644 (file)
@@ -37,7 +37,7 @@ AlidNdEtaCorrectionTask::AlidNdEtaCorrectionTask() :
   fESD(0),
   fOutput(0),
   fOption(),
-  fAnalysisMode(AliPWG0Helper::kTPC),
+  fAnalysisMode((AliPWG0Helper::AnalysisMode) (AliPWG0Helper::kTPC | AliPWG0Helper::kFieldOn)),
   fTrigger(AliPWG0Helper::kMB1),
   fFillPhi(kFALSE),
   fDeltaPhiCut(-1),
@@ -86,7 +86,7 @@ AlidNdEtaCorrectionTask::AlidNdEtaCorrectionTask(const char* opt) :
   fESD(0),
   fOutput(0),
   fOption(opt),
-  fAnalysisMode(AliPWG0Helper::kTPC),
+  fAnalysisMode((AliPWG0Helper::AnalysisMode) (AliPWG0Helper::kTPC | AliPWG0Helper::kFieldOn)),
   fTrigger(AliPWG0Helper::kMB1),
   fFillPhi(kFALSE),
   fDeltaPhiCut(0),
@@ -167,10 +167,10 @@ void AlidNdEtaCorrectionTask::ConnectInputData(Option_t *)
     // Enable only the needed branches
     esdH->SetActiveBranches("AliESDHeader Vertex");
 
-    if (fAnalysisMode == AliPWG0Helper::kSPD)
+    if (fAnalysisMode & AliPWG0Helper::kSPD)
       esdH->SetActiveBranches("AliESDHeader Vertex AliMultiplicity");
 
-    if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS) {
+    if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS) {
       esdH->SetActiveBranches("AliESDHeader Vertex Tracks");
     }
   }
@@ -261,13 +261,13 @@ void AlidNdEtaCorrectionTask::CreateOutputObjects()
 
   fVertexCorrelation = new TH2F("fVertexCorrelation", "fVertexCorrelation;MC z-vtx;ESD z-vtx", 120, -30, 30, 120, -30, 30);
   fOutput->Add(fVertexCorrelation);
-  fVertexCorrelationShift = new TH2F("fVertexCorrelationShift", "fVertexCorrelationShift;MC z-vtx;MC z-vtx - ESD z-vtx", 120, -30, 30, 100, -1, 1);
+  fVertexCorrelationShift = new TH3F("fVertexCorrelationShift", "fVertexCorrelationShift;MC z-vtx;MC z-vtx - ESD z-vtx;rec. tracks", 120, -30, 30, 100, -1, 1, 100, -0.5, 99.5);
   fOutput->Add(fVertexCorrelationShift);
   fVertexProfile = new TProfile("fVertexProfile", "fVertexProfile;MC z-vtx;MC z-vtx - ESD z-vtx", 40, -20, 20);
   fOutput->Add(fVertexProfile);
   fVertexShift = new TH1F("fVertexShift", "fVertexShift;(MC z-vtx - ESD z-vtx);Entries", 201, -2, 2);
   fOutput->Add(fVertexShift);
-  fVertexShiftNorm = new TH1F("fVertexShiftNorm", "fVertexShiftNorm;(MC z-vtx - ESD z-vtx) / #sigma_{ESD z-vtx};Entries", 200, -100, 100);
+  fVertexShiftNorm = new TH2F("fVertexShiftNorm", "fVertexShiftNorm;(MC z-vtx - ESD z-vtx) / #sigma_{ESD z-vtx};rec. tracks;Entries", 200, -100, 100, 100, -0.5, 99.5);
   fOutput->Add(fVertexShiftNorm);
 
   fEtaCorrelation = new TH2F("fEtaCorrelation", "fEtaCorrelation;MC #eta;ESD #eta", 120, -3, 3, 120, -3, 3);
@@ -334,9 +334,10 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
     
   if (fStatError > 0)
     Printf("WARNING: Statistical error evaluation active!");
-
+    
   // trigger definition
   Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD, fTrigger);
+  //Printf("Trigger mask: %lld", fESD->GetTriggerMask());
 
   if (!eventTriggered)
     Printf("No trigger");
@@ -384,35 +385,17 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
   genHeader->PrimaryVertex(vtxMC);
 
   // get the ESD vertex
-  const AliESDVertex* vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
   Bool_t eventVertex = kFALSE;
-  if (vtxESD)
+  Double_t vtx[3];
+  const AliESDVertex* vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
+  if (vtxESD && AliPWG0Helper::TestVertex(vtxESD, fAnalysisMode))
   {
-    Double_t vtx[3];
+    eventVertex = kTRUE;
     vtxESD->GetXYZ(vtx);
-
-    Double_t diff = vtxMC[2] - vtx[2];
-    fVertexShift->Fill(diff);
-    if (vtxESD->GetZRes() > 0)
-        fVertexShiftNorm->Fill(diff / vtxESD->GetZRes());
-
-    if (!AliPWG0Helper::TestVertex(vtxESD, fAnalysisMode))
-    {
-        vtxESD = 0;
-    }
-    else
-    {
-      eventVertex = kTRUE;
-
-      if (eventTriggered)
-      {
-        fVertexCorrelation->Fill(vtxMC[2], vtx[2]);
-        fVertexCorrelationShift->Fill(vtxMC[2], vtxMC[2] - vtx[2]);
-        fVertexProfile->Fill(vtxMC[2], vtxMC[2] - vtx[2]);
-      }
-    }
   }
-
+  else
+    vtxESD = 0;
+    
   // fill process type
   Int_t biny = (Int_t) eventTriggered + 2 * (Int_t) eventVertex;
   // INEL
@@ -435,7 +418,7 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
   Float_t* etaArr = 0;
   Float_t* thirdDimArr = 0;
   Float_t* deltaPhiArr = 0;
-  if (fAnalysisMode == AliPWG0Helper::kSPD)
+  if (fAnalysisMode & AliPWG0Helper::kSPD)
   {
     // get tracklets
     const AliMultiplicity* mult = fESD->GetMultiplicity();
@@ -479,18 +462,18 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
       ++inputCount;
     }
   }
-  else if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS)
+  else if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
   {
     if (!fEsdTrackCuts)
     {
       AliDebug(AliLog::kError, "fESDTrackCuts not available");
       return;
     }
-
+    
     if (vtxESD)
     {
       // get multiplicity from ESD tracks
-      TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD, (fAnalysisMode == AliPWG0Helper::kTPC));
+      TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD, (fAnalysisMode & AliPWG0Helper::kTPC));
       Int_t nGoodTracks = list->GetEntries();
   
       Printf("Accepted %d tracks", nGoodTracks);
@@ -511,7 +494,17 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
           continue;
         }
         
-        // TODO fOnlyPrimaries not implemented for TPC
+        //Printf("status is: %u", esdTrack->GetStatus());
+        
+        if (fOnlyPrimaries)
+        {
+          Int_t label = TMath::Abs(esdTrack->GetLabel());
+          if (label == 0)
+            continue;
+          
+          if (stack->IsPhysicalPrimary(label) == kFALSE)
+            continue;
+        }        
   
         etaArr[inputCount] = esdTrack->Eta();
         labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel());
@@ -530,9 +523,9 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
         {
           AliESDtrack* track = 0;
   
-          if (fAnalysisMode == AliPWG0Helper::kTPC)
+          if (fAnalysisMode & AliPWG0Helper::kTPC)
             track = AliESDtrackCuts::GetTPCOnlyTrack(fESD, iTrack);
-          else if (fAnalysisMode == AliPWG0Helper::kTPCITS)
+          else if (fAnalysisMode & AliPWG0Helper::kTPCITS)
             track = fESD->GetTrack(iTrack);
   
           if (!track)
@@ -550,24 +543,24 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
           if (stack->Particle(label)->GetPDG()->Charge() == 0)
             continue;
   
-          if (TMath::Abs(track->Eta()) < 1)
+          if (TMath::Abs(track->Eta()) < 0.8 && track->Pt() > 0.15)
           {
             if (stack->IsPhysicalPrimary(label))
             {
               // primary
               if (fEsdTrackCutsPrim->AcceptTrack(track)) 
               {
-                if (AliESDtrackCuts::GetSigmaToVertex(track) > 900)
-                {
-                  Printf("Track %d has nsigma of %f. Printing track and vertex...", iTrack, AliESDtrackCuts::GetSigmaToVertex(track));
-                  Float_t b[2];
-                  Float_t r[3];
-                  track->GetImpactParameters(b, r);
-                  Printf("Impact parameter %f %f and resolution: %f %f %f", b[0], b[1], r[0], r[1], r[2]);
-                  track->Print("");
-                  if (vtxESD)
-                    vtxESD->Print();
-                }
+//                 if (AliESDtrackCuts::GetSigmaToVertex(track) > 900)
+//                 {
+//                   Printf("Track %d has nsigma of %f. Printing track and vertex...", iTrack, AliESDtrackCuts::GetSigmaToVertex(track));
+//                   Float_t b[2];
+//                   Float_t r[3];
+//                   track->GetImpactParameters(b, r);
+//                   Printf("Impact parameter %f %f and resolution: %f %f %f", b[0], b[1], r[0], r[1], r[2]);
+//                   track->Print("");
+//                   if (vtxESD)
+//                     vtxESD->Print();
+//                 }
               }
             }
             else
@@ -578,7 +571,7 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
           }
   
           // TODO mem leak in the continue statements above
-          if (fAnalysisMode == AliPWG0Helper::kTPC)
+          if (fAnalysisMode & AliPWG0Helper::kTPC)
             delete track;
         }
       }
@@ -615,7 +608,7 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
     Float_t eta = particle->Eta();
     
     Float_t thirdDim = -1;
-    if (fAnalysisMode == AliPWG0Helper::kSPD)
+    if (fAnalysisMode & AliPWG0Helper::kSPD)
     {
       if (fFillPhi)
       {
@@ -743,7 +736,7 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
       // resolutions
       fEtaResolution->Fill(particle->Eta() - etaArr[i]);
 
-      if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS)
+      if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
         if (TMath::Abs(particle->Eta() < 0.9) && particle->Pt() > 0)
           fpTResolution->Fill(particle->Pt(), (particle->Pt() - thirdDimArr[i]) / particle->Pt());
 
@@ -751,12 +744,12 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
       Float_t thirdDim = -1;
 
       Bool_t firstIsPrim = stack->IsPhysicalPrimary(label);
-      // in case of primary the MC values are filled, otherwise (background) the reconstructed values
-      if (label == label2 && firstIsPrim)
+      // in case of same label the MC values are filled, otherwise (background) the reconstructed values
+      if (label == label2)
       {
         eta = particle->Eta();
         
-        if (fAnalysisMode == AliPWG0Helper::kSPD)
+        if (fAnalysisMode & AliPWG0Helper::kSPD)
         {
           if (fFillPhi)
           {
@@ -770,7 +763,7 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
       }
       else
       {
-        if (fAnalysisMode == AliPWG0Helper::kSPD && !fFillPhi)
+        if (fAnalysisMode & AliPWG0Helper::kSPD && !fFillPhi)
         {
           thirdDim = inputCount;
         }
@@ -915,6 +908,18 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
 
   if (eventTriggered && eventVertex)
   {
+    Double_t diff = vtxMC[2] - vtx[2];
+    fVertexShift->Fill(diff);
+    if (vtxESD->GetZRes() > 0)
+        fVertexShiftNorm->Fill(diff / vtxESD->GetZRes(), inputCount);
+
+    fVertexCorrelation->Fill(vtxMC[2], vtx[2]);
+    fVertexCorrelationShift->Fill(vtxMC[2], vtxMC[2] - vtx[2], inputCount);
+    fVertexProfile->Fill(vtxMC[2], vtxMC[2] - vtx[2]);
+  }
+
+  if (eventTriggered && eventVertex)
+  {
     fdNdEtaAnalysisMC->FillEvent(vtxMC[2], inputCount);
     fdNdEtaAnalysisESD->FillEvent(vtxMC[2], inputCount);
   }
@@ -925,7 +930,7 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
   if (fOption.Contains("process-types"))
   {
     // non diffractive
-    if (processType == AliPWG0Helper::kND )
+    if (processType == AliPWG0Helper::kND)
       fdNdEtaCorrectionSpecial[0]->FillEvent(vtxMC[2], inputCount, eventTriggered, eventVertex, processType);
 
     // single diffractive
@@ -1023,25 +1028,34 @@ void AlidNdEtaCorrectionTask::Terminate(Option_t *)
   fVertexCorrelation = dynamic_cast<TH2F*> (fOutput->FindObject("fVertexCorrelation"));
   if (fVertexCorrelation)
     fVertexCorrelation->Write();
-  fVertexCorrelationShift = dynamic_cast<TH2F*> (fOutput->FindObject("fVertexCorrelationShift"));
+  fVertexCorrelationShift = dynamic_cast<TH3F*> (fOutput->FindObject("fVertexCorrelationShift"));
   if (fVertexCorrelationShift)
+  {
+    ((TH2*) fVertexCorrelationShift->Project3D("yx"))->FitSlicesY();
     fVertexCorrelationShift->Write();
+  }
   fVertexProfile = dynamic_cast<TProfile*> (fOutput->FindObject("fVertexProfile"));
   if (fVertexProfile)
     fVertexProfile->Write();
   fVertexShift = dynamic_cast<TH1F*> (fOutput->FindObject("fVertexShift"));
   if (fVertexShift)
     fVertexShift->Write();
-  fVertexShiftNorm = dynamic_cast<TH1F*> (fOutput->FindObject("fVertexShiftNorm"));
+  fVertexShiftNorm = dynamic_cast<TH2F*> (fOutput->FindObject("fVertexShiftNorm"));
   if (fVertexShiftNorm)
+  {
+    fVertexShiftNorm->ProjectionX();
     fVertexShiftNorm->Write();
+  }  
 
   fEtaCorrelation = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaCorrelation"));
   if (fEtaCorrelation)
     fEtaCorrelation->Write();
   fEtaCorrelationShift = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaCorrelationShift"));
   if (fEtaCorrelationShift)
+  {
+    fEtaCorrelationShift->FitSlicesY();
     fEtaCorrelationShift->Write();
+  }
   fEtaProfile = dynamic_cast<TProfile*> (fOutput->FindObject("fEtaProfile"));
   if (fEtaProfile)
     fEtaProfile->Write();
index e8ae725..6021908 100644 (file)
@@ -15,6 +15,7 @@ class TH1F;
 class AliESDEvent;
 class TParticlePDG;
 class TH2F;
+class TH3F;
 class TProfile;
 
 class AlidNdEtaCorrectionTask : public AliAnalysisTask {
@@ -65,10 +66,10 @@ class AlidNdEtaCorrectionTask : public AliAnalysisTask {
     TH1F* fPIDTracks;                            //! pid of reconstructed tracks
 
     TH2F* fVertexCorrelation;                    //! ESD z-vtx vs MC z-vtx
-    TH2F* fVertexCorrelationShift;               //! (MC z-vtx - ESD z-vtx) vs MC z-vtx
+    TH3F* fVertexCorrelationShift;               //! (MC z-vtx - ESD z-vtx) vs MC z-vtx vs n# rec tracks
     TProfile* fVertexProfile;                    //! Profile of MC z-vtx - ESD z-vtx vs. MC z-vtx
     TH1F* fVertexShift;                          //! (MC z-vtx - ESD z-vtx) in +- 10 cm
-    TH1F* fVertexShiftNorm;                      //! (MC z-vtx - ESD z-vtx) / (sigma_ESD-z-vtx) histogrammed
+    TH2F* fVertexShiftNorm;                      //! (MC z-vtx - ESD z-vtx) / (sigma_ESD-z-vtx) vs. no. rec tracks
 
     TH2F* fEtaCorrelation;                       //! ESD eta vs MC eta
     TH2F* fEtaCorrelationShift;                  //! (MC eta - ESD eta) vs MC eta
index 32eceab..b613af8 100644 (file)
@@ -44,7 +44,7 @@ AlidNdEtaTask::AlidNdEtaTask(const char* opt) :
   fESD(0),
   fOutput(0),
   fOption(opt),
-  fAnalysisMode(AliPWG0Helper::kTPC),
+  fAnalysisMode((AliPWG0Helper::AnalysisMode) (AliPWG0Helper::kTPC | AliPWG0Helper::kFieldOn)),
   fTrigger(AliPWG0Helper::kMB1),
   fFillPhi(kFALSE),
   fDeltaPhiCut(-1),
@@ -70,6 +70,7 @@ AlidNdEtaTask::AlidNdEtaTask(const char* opt) :
   fRawPt(0),
   fEtaPhi(0),
   fDeltaPhi(0),
+  fDeltaTheta(0),
   fFiredChips(0),
   fTriggerVsTime(0),
   fStats(0)
@@ -125,16 +126,17 @@ void AlidNdEtaTask::ConnectInputData(Option_t *)
     Printf("ERROR: Could not get ESDInputHandler");
   } else {
     fESD = esdH->GetEvent();
+    
+    TString branches("AliESDHeader Vertex ");
 
+    if (fAnalysisMode & AliPWG0Helper::kSPD || fTrigger == AliPWG0Helper::kOfflineMB1 || fTrigger == AliPWG0Helper::kOfflineMB2 || fTrigger == AliPWG0Helper::kOfflineMB3 || fTrigger == AliPWG0Helper::kOfflineFASTOR)
+      branches += "AliMultiplicity ";
+      
+    if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
+      branches += "Tracks ";
+  
     // Enable only the needed branches
-    esdH->SetActiveBranches("AliESDHeader Vertex");
-
-    if (fAnalysisMode == AliPWG0Helper::kSPD)
-      esdH->SetActiveBranches("AliESDHeader Vertex AliMultiplicity");
-
-    if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS) {
-      esdH->SetActiveBranches("AliESDHeader Vertex Tracks");
-    }
+    esdH->SetActiveBranches(branches);
   }
 
   // disable info messages of AliMCEvent (per event)
@@ -178,7 +180,7 @@ void AlidNdEtaTask::CreateOutputObjects()
   fEvents = new TH1F("dndeta_check_vertex", "dndeta_check_vertex", 800, -40, 40);
   fOutput->Add(fEvents);
 
-  Float_t resMax = (fAnalysisMode == AliPWG0Helper::kSPD) ? 0.2 : 2;
+  Float_t resMax = (fAnalysisMode & AliPWG0Helper::kSPD) ? 0.2 : 2;
   fVertexResolution = new TH1F("dndeta_vertex_resolution_z", "dndeta_vertex_resolution_z", 1000, 0, resMax);
   fOutput->Add(fVertexResolution);
 
@@ -194,15 +196,18 @@ void AlidNdEtaTask::CreateOutputObjects()
   fTriggerVsTime->GetYaxis()->SetTitle("count");
   fOutput->Add(fTriggerVsTime);
 
-  fStats = new TH1F("fStats", "fStats", 2, 0.5, 2.5);
+  fStats = new TH1F("fStats", "fStats", 3, 0.5, 3.5);
   fStats->GetXaxis()->SetBinLabel(1, "vertexer 3d");
   fStats->GetXaxis()->SetBinLabel(2, "vertexer z");
+  fStats->GetXaxis()->SetBinLabel(3, "trigger");
   fOutput->Add(fStats);
 
-  if (fAnalysisMode == AliPWG0Helper::kSPD)
+  if (fAnalysisMode & AliPWG0Helper::kSPD)
   {
     fDeltaPhi = new TH1F("fDeltaPhi", "fDeltaPhi;#Delta #phi;Entries", 500, -0.2, 0.2);
     fOutput->Add(fDeltaPhi);
+    fDeltaTheta = new TH1F("fDeltaTheta", "fDeltaTheta;#Delta #theta;Entries", 500, -0.2, 0.2);
+    fOutput->Add(fDeltaTheta);
     fFiredChips = new TH2F("fFiredChips", "fFiredChips;Chips L1 + L2;tracklets", 1201, -0.5, 1201, 50, -0.5, 49.5);
     fOutput->Add(fFiredChips);
     for (Int_t i=0; i<2; i++)
@@ -212,7 +217,7 @@ void AlidNdEtaTask::CreateOutputObjects()
     }
   }
 
-  if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS)
+  if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
   {
     fRawPt =  new TH1F("fRawPt", "raw pt;p_{T};Count", 2000, 0, 100);
     fOutput->Add(fRawPt);
@@ -286,6 +291,8 @@ void AlidNdEtaTask::Exec(Option_t*)
 
     // trigger definition
     eventTriggered = AliPWG0Helper::IsEventTriggered(fESD, fTrigger);
+    if (eventTriggered)
+      fStats->Fill(3);
 
     // get the ESD vertex
     vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
@@ -368,7 +375,7 @@ void AlidNdEtaTask::Exec(Option_t*)
     Int_t* labelArr = 0;
     Float_t* etaArr = 0;
     Float_t* thirdDimArr = 0;
-    if (fAnalysisMode == AliPWG0Helper::kSPD)
+    if (fAnalysisMode & AliPWG0Helper::kSPD)
     {
       // get tracklets
       const AliMultiplicity* mult = fESD->GetMultiplicity();
@@ -422,6 +429,7 @@ void AlidNdEtaTask::Exec(Option_t*)
         }
 
         fDeltaPhi->Fill(deltaPhi);
+        fDeltaTheta->Fill(mult->GetDeltaTheta(i));
 
         if (fDeltaPhiCut > 0 && TMath::Abs(deltaPhi) > fDeltaPhiCut)
           continue;
@@ -455,7 +463,7 @@ void AlidNdEtaTask::Exec(Option_t*)
       fFiredChips->Fill(firedChips, inputCount);
       Printf("Accepted %d tracklets (%d fired chips)", inputCount, firedChips);
     }
-    else if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS)
+    else if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
     {
       if (!fEsdTrackCuts)
       {
@@ -466,9 +474,9 @@ void AlidNdEtaTask::Exec(Option_t*)
       if (vtxESD)
       {
         // get multiplicity from ESD tracks
-        TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD, (fAnalysisMode == AliPWG0Helper::kTPC));
+        TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD, fAnalysisMode & AliPWG0Helper::kTPC);
         Int_t nGoodTracks = list->GetEntries();
-        Printf("Accepted %d tracks", nGoodTracks);
+        Printf("Accepted %d tracks out of %d total ESD tracks", nGoodTracks, fESD->GetNumberOfTracks());
   
         labelArr = new Int_t[nGoodTracks];
         etaArr = new Float_t[nGoodTracks];
@@ -491,14 +499,24 @@ void AlidNdEtaTask::Exec(Option_t*)
           Float_t eta = esdTrack->Eta();
           Int_t label = TMath::Abs(esdTrack->GetLabel());
           Float_t pT  = esdTrack->Pt();
+          
+          // force pT to fixed value without B field
+          if (!(fAnalysisMode & AliPWG0Helper::kFieldOn))
+            pT = 1;
   
           fPhi->Fill(phi);
           fEtaPhi->Fill(eta, phi);
           if (eventTriggered && vtxESD)
             fRawPt->Fill(pT);
   
-          if (fOnlyPrimaries && label == 0)
-            continue;
+          if (fOnlyPrimaries)
+          {
+            if (label == 0)
+              continue;
+            
+            if (stack->IsPhysicalPrimary(label) == kFALSE)
+              continue;
+          }
   
           if (fUseMCKine)
           {
@@ -584,7 +602,7 @@ void AlidNdEtaTask::Exec(Option_t*)
             }
 
             Float_t thirdDim = -1;
-            if (fAnalysisMode == AliPWG0Helper::kSPD)
+            if (fAnalysisMode & AliPWG0Helper::kSPD)
             {
               if (fFillPhi)
               {
@@ -704,7 +722,7 @@ void AlidNdEtaTask::Exec(Option_t*)
       Float_t eta = particle->Eta();
       Float_t thirdDim = -1;
 
-      if (fAnalysisMode == AliPWG0Helper::kSPD)
+      if (fAnalysisMode & AliPWG0Helper::kSPD)
       {
         if (fFillPhi)
         {
@@ -731,7 +749,7 @@ void AlidNdEtaTask::Exec(Option_t*)
           fdNdEtaAnalysisTrVtx->FillTrack(vtxMC[2], eta, thirdDim);
       }
 
-      if (TMath::Abs(eta) < 1.0)
+      if (TMath::Abs(eta) < 1.0 && particle->Pt() > 0)
         fPartPt->Fill(particle->Pt());
     }
 
@@ -776,6 +794,7 @@ void AlidNdEtaTask::Terminate(Option_t *)
     for (Int_t i=0; i<2; ++i)
       fZPhi[i] = dynamic_cast<TH2F*> (fOutput->FindObject(Form("fZPhi_%d", i)));
     fDeltaPhi = dynamic_cast<TH1F*> (fOutput->FindObject("fDeltaPhi"));
+    fDeltaTheta = dynamic_cast<TH1F*> (fOutput->FindObject("fDeltaTheta"));
     fFiredChips = dynamic_cast<TH2F*> (fOutput->FindObject("fFiredChips"));
     fTriggerVsTime = dynamic_cast<TGraph*> (fOutput->FindObject("fTriggerVsTime"));
     fStats = dynamic_cast<TH1F*> (fOutput->FindObject("fStats"));
@@ -840,7 +859,7 @@ void AlidNdEtaTask::Terminate(Option_t *)
       fdNdEtaAnalysisESD->SaveHistograms();
 
     if (fEsdTrackCuts)
-      fEsdTrackCuts->SaveHistograms("esd_tracks_cuts");
+      fEsdTrackCuts->SaveHistograms("esd_track_cuts");
 
     if (fMult)
       fMult->Write();
@@ -861,6 +880,9 @@ void AlidNdEtaTask::Terminate(Option_t *)
     if (fDeltaPhi)
       fDeltaPhi->Write();
 
+    if (fDeltaTheta)
+      fDeltaTheta->Write();
+    
     if (fPhi)
       fPhi->Write();
 
index b9ee6d0..9d39e93 100644 (file)
@@ -35,7 +35,7 @@ class AlidNdEtaTask : public AliAnalysisTask {
     void SetTrigger(AliPWG0Helper::Trigger trigger) { fTrigger = trigger; }
     void SetFillPhi(Bool_t flag = kTRUE) { fFillPhi = flag; }
     void SetDeltaPhiCut(Float_t cut) { fDeltaPhiCut = cut; }
-
+    
     void SetOption(const char* opt) { fOption = opt; }
 
  protected:
@@ -82,6 +82,7 @@ class AlidNdEtaTask : public AliAnalysisTask {
     TH2F* fEtaPhi;                //! raw eta - phi distribution
     TH2F* fZPhi[2];               //! raw z - phi distribution from tracklets per layer (only SPD)
     TH1F* fDeltaPhi;              //! histogram of delta_phi values for tracklets (only for SPD analysis)
+    TH1F* fDeltaTheta;            //! histogram of delta_theta values for tracklets (only for SPD analysis)
     TH2F* fFiredChips;            //! fired chips l1+l2 vs. number of tracklets (only for SPD analysis)
     TGraph* fTriggerVsTime;       //! trigger as function of event time
     TH1F* fStats;                 //! further statistics : bin 1 = vertexer 3d, bin 2 = vertexer z
index 40d02f3..957f1f2 100644 (file)
@@ -262,13 +262,13 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
     //     alpha = triggered events with vertex at a given vertex position / all triggered events with vertex
     //     triggered events without vertex and 0 multiplicity at a given vertex position = alpha * all triggered events with 0 multiplicity
     //   afterwards we still correct for the trigger efficiency
+    // at the same time calculate expectation from MC (not used, just to check the value)
 
     //TH2* measuredEvents = fData->GetEventCorrection()->GetMeasuredHistogram();
     TH2* correctedEvents = fData->GetEventCorrection()->GetGeneratedHistogram();
 
     TH2* eTrig =    correction->GetVertexRecoCorrection()->GetEventCorrection()->GetGeneratedHistogram();
     TH2* eTrigVtx = correction->GetVertexRecoCorrection()->GetEventCorrection()->GetMeasuredHistogram();
-    //TH1* eTrigVtx_projx = eTrigVtx->ProjectionX("eTrigVtx_projx", 2, rawMeasured->GetNbinsY()+1);
     TH1* eTrigVtx_projx = eTrigVtx->ProjectionX("eTrigVtx_projx", 2, rawMeasured->GetNbinsY()+1);
 
     //new TCanvas; correctedEvents->DrawCopy("TEXT");
@@ -285,6 +285,7 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
     TH1* kineBias = (TH1*) vertexDist->Clone("kineBias");
     kineBias->Reset();
 
+    // loop over vertex bins
     for (Int_t i = 1; i <= rawMeasured->GetNbinsX(); i++)
     {
       Double_t alpha = (Double_t) vertexDist->GetBinContent(i) / allEventsWithVertex;
@@ -302,10 +303,16 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
       // multiply with trigger correction if set above
       events *= fData->GetEventCorrection()->GetCorrectionHistogram()->GetBinContent(i, 1);
 
-      Printf("Bin %d, alpha is %.2f, fZ is %.3f, number of events with 0 mult.: %.2f", i, alpha * 100., fZ, events);
+      // calculate how many events we would have got with a pure MC-based correction
+      //   in the given bin: measured events with vertex (mult > 0) * triggered events with mult 0 (mc) / events with vertex and mult > 0 (mc) * trigger correction for bin 0
+      Double_t mcEvents = vertexDist->GetBinContent(i) * eTrig->GetBinContent(i, 1) / eTrigVtx_projx->GetBinContent(i) * fData->GetEventCorrection()->GetCorrectionHistogram()->GetBinContent(i, 1);
+
+      Printf("Bin %d, alpha is %.2f, fZ is %.3f, number of events with 0 mult.: %.2f (MC comparison: %.2f)", i, alpha * 100., fZ, events, mcEvents);
 
       correctedEvents->SetBinContent(i, 1, events);
     }
+    
+    Printf("In |vtx-z| < 10 cm: %d events have been added", (Int_t) correctedEvents->Integral(vertexDist->FindBin(-9.9), vertexDist->FindBin(9.9), 1, 1));
 
     //new TCanvas; correctedEvents->DrawCopy("TEXT");
     //new TCanvas; kineBias->DrawCopy();
@@ -321,7 +328,7 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
   TH1D* vertexHist = (TH1D*) tmp->ProjectionX("_px", 0, tmp->GetNbinsY() + 1, "e");
 
   // create pt hist
-  if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS)
+  if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
   {
     // reset all ranges
     dataHist->GetXaxis()->SetRange(0, 0);
@@ -362,7 +369,7 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
 
   // integrate over pt (with pt cut) (TPC, TPCITS) or multiplicity (SPD)
   Int_t ptLowBin = 1;
-  if (ptCut > 0 && fAnalysisMode != AliPWG0Helper::kSPD)
+  if (ptCut > 0 && (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS))
     ptLowBin = dataHist->GetZaxis()->FindBin(ptCut);
     
   //new TCanvas; dataHist->DrawCopy();
@@ -409,7 +416,7 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
 
       // adjust acceptance range
       // produce with drawPlots.C: DetermineAcceptance(...)
-      if (fAnalysisMode == AliPWG0Helper::kSPD)
+      if (fAnalysisMode & AliPWG0Helper::kSPD)
       {
         //const Int_t binBeginSPD[30] = { 18, 16, 15, 14, 13, 13, 12, 11, 9, 7, 6, 5, 5, 5, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 1 }; // by eye
         //const Int_t binBeginSPD[30] = { -1, -1, -1, -1, 16, 14, 12, 10, 9, 7, 6, 5, 5, 4, 4, 4, 4, 3, 3, 3, 3, 2, 2, 2, 2, 2, -1, -1, -1, -1}; // limit in correction map is 5
@@ -421,14 +428,14 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
 
         binBegin = binBeginSPD;
       }
-      else if (fAnalysisMode == AliPWG0Helper::kTPC)
+      else if (fAnalysisMode & AliPWG0Helper::kTPC)
       {
         //const Int_t binBeginTPC[30] = {-1, -1, -1, -1, -1, -1, -1, -1, -1, 4, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1}; // limit 5, pt cut off 0.2 mev/c
         const Int_t binBeginTPC[maxBins] = {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 9, 5, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}; // limit 5
 
         binBegin = binBeginTPC;
       }
-      else if (fAnalysisMode == AliPWG0Helper::kTPCITS)
+      else if (fAnalysisMode & AliPWG0Helper::kTPCITS)
       {
         // TODO create map
       }
@@ -501,7 +508,7 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
       Float_t ptCutOffCorrection = 1;
 
       // find pt cut off correction factor
-      if (fAnalysisMode != AliPWG0Helper::kSPD)
+      if ((fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS) && (fAnalysisMode & AliPWG0Helper::kFieldOn))
       {
         if (correction && ptCut > 0)
             ptCutOffCorrection = correction->GetMeasuredFraction(correctionType, ptCut, vtxVsEta->GetYaxis()->GetBinCenter(iEta), vertexBinBegin, vertexBinEnd);
index 8db7e0b..6eea884 100644 (file)
@@ -33,7 +33,7 @@ public:
   enum { kVertexBinning = 1+2 }; // the first is for the whole vertex range, the others divide the vertex range
 
   dNdEtaAnalysis();
-  dNdEtaAnalysis(const Char_t* name, const Char_t* title, AliPWG0Helper::AnalysisMode analysisMode = AliPWG0Helper::kSPD);
+  dNdEtaAnalysis(const Char_t* name, const Char_t* title, AliPWG0Helper::AnalysisMode analysisMode = (AliPWG0Helper::AnalysisMode) (AliPWG0Helper::kTPC | AliPWG0Helper::kFieldOn));
   virtual ~dNdEtaAnalysis();
 
   dNdEtaAnalysis(const dNdEtaAnalysis &c);
index a4a684f..7518edf 100644 (file)
@@ -145,7 +145,8 @@ void PrintInfo(const char* fileName = "correction_map.root", const char* dirName
   for (Int_t i=AlidNdEtaCorrection::kTrack2Particle; i<=AlidNdEtaCorrection::kND; i++)
   {
     Printf("Correction %d", i);
-    dNdEtaCorrection->GetCorrection(i)->PrintInfo(0.3);
+    dNdEtaCorrection->GetCorrection(i)->PrintInfo(0);
+    return;
   }
 }
 
@@ -239,7 +240,7 @@ void CompareVertexDist(Int_t plot = 0, const char* correctionMapFile = "correcti
 
 void CompareTrack2ParticleWithAnalysisData(const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
 {
-  gSystem->Load("libPWG0base");
+  loadlibs();
 
   TFile::Open("analysis_esd.root");
   dNdEtaAnalysis* fdNdEtaAnalysisESD = new dNdEtaAnalysis("dndeta_mbvtx", "dndeta_mbvtx");
@@ -495,8 +496,8 @@ void dNdEta(Bool_t onlyESD = kFALSE, Bool_t save = kTRUE)
   histESDMBVtxNoPt->SetMarkerStyle(22);
   histESDMBTracksNoPt->SetMarkerStyle(23);
   
-  Float_t etaLimit = (fdNdEtaAnalysis->GetAnalysisMode() == AliPWG0Helper::kTPC) ? 0.89 : 1.79;
-  Float_t etaPlotLimit = (fdNdEtaAnalysis->GetAnalysisMode() == AliPWG0Helper::kTPC) ? 1.2 : 2.3;
+  Float_t etaLimit = (fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kTPC) ? 0.89 : 1.79;
+  Float_t etaPlotLimit = (fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kTPC) ? 1.2 : 2.3;
   //Float_t etaLimit = (fdNdEtaAnalysis->GetAnalysisMode() == AliPWG0Helper::kTPC) ? 0.89 : 1.39;
   //Float_t etaPlotLimit = (fdNdEtaAnalysis->GetAnalysisMode() == AliPWG0Helper::kTPC) ? 1.2 : 1.9;
 
@@ -512,6 +513,7 @@ void dNdEta(Bool_t onlyESD = kFALSE, Bool_t save = kTRUE)
 
   Float_t max = TMath::Max(histESDMBVtx->GetMaximum(), histESDMB->GetMaximum());
   max = TMath::Max(max, histESD->GetMaximum());
+  max = TMath::Max(max, histESDMBTracksNoPt->GetMaximum());
 
   TLegend* legend = new TLegend(0.35, 0.05, 0.75, 0.4);
   legend->SetFillColor(0);
@@ -685,7 +687,7 @@ void dNdEta(Bool_t onlyESD = kFALSE, Bool_t save = kTRUE)
   legend7->AddEntry(ratio, "Trigger-bias INEL", "P");
   legend7->AddEntry(ratioTr, "Vertex-reconstruction", "P");
   legend7->AddEntry(ratioNSD, "Trigger-bias NSD", "P");
-  if (fdNdEtaAnalysis->GetAnalysisMode() == AliPWG0Helper::kTPC)
+  if ((fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kFieldOn) && (fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kTPC || fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kTPCITS))
     legend7->AddEntry(ratioTrVtx, "p_{T} cut-off", "P");
   
   TH1* dummy7 = new TH2F("dummy7", ";#eta;Deviation in %", 100, -etaPlotLimit, etaPlotLimit, 100, -5, 7);
@@ -700,7 +702,7 @@ void dNdEta(Bool_t onlyESD = kFALSE, Bool_t save = kTRUE)
   
   ratio->Draw("HIST EP SAME");
   ratioTr->Draw("HIST EP SAME");
-  if (fdNdEtaAnalysis->GetAnalysisMode() == AliPWG0Helper::kTPC)
+  if ((fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kFieldOn) && (fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kTPC || fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kTPCITS))
     ratioTrVtx->Draw("HIST EP SAME");
   ratioTrVtxNoPt->Draw("HIST EP SAME");
   ratioNSD->Draw("HIST EP SAME");
@@ -739,6 +741,7 @@ void dNdEta(Bool_t onlyESD = kFALSE, Bool_t save = kTRUE)
   PrintIntegratedDeviation(histMCnsdNoPt, histESDnsdNoPt, "all events (NSD) (no pt corr)");
   PrintIntegratedDeviation(histMCTrPtCut, histESDMBNoPt, "triggered (no pt corr)");
   PrintIntegratedDeviation(histMCTrVtxPtCut, histESDMBVtxNoPt, "trigger, vertex (no pt corr)");
+  PrintIntegratedDeviation(histESD, histESDNoPt, "pt cut off correction");
 
   TCanvas* canvas3 = new TCanvas("dNdEta", "dNdEta", 600, 600);
   canvas3->Range(0, 0, 1, 1);
@@ -858,7 +861,7 @@ TH1* DrawdNdEtaRatio(TH1* corr, TH1* mc, const char* name, Float_t etaPlotLimit)
   legend->AddEntry(mc, "MC prediction");
   legend->SetTextSize(0.08);
 
-  TH2F* dummy = new TH2F("dummy", "", 100, -etaPlotLimit, etaPlotLimit, 1000, 3.1, corr->GetMaximum() * 1.1);
+  TH2F* dummy = new TH2F("dummy", "", 100, -etaPlotLimit, etaPlotLimit, 1000, 2.7, corr->GetMaximum() * 1.1);
   Prepare1DPlot(dummy);
   dummy->SetStats(kFALSE);
   dummy->SetXTitle("#eta");
@@ -1841,25 +1844,6 @@ void MultiplicityESD()
   fMultiplicityESD->Draw();
 }
 
-void drawPlots(Int_t max)
-{
-  gMax = max;
-
-  ptCutoff();
-  TriggerBias();
-  VtxRecon();
-  Track2Particle2D();
-  Track2Particle3D();
-  ptSpectrum();
-  dNdEta();
-}
-
-void drawPlots()
-{
-  drawPlots(5);
-  drawPlots(2);
-}
-
 void CompareCorrection2Measured(Float_t ptMin = 0.301, const char* dataInput = "analysis_esd_raw.root", const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
 {
   loadlibs();
@@ -1903,6 +1887,21 @@ void CompareCorrection2Measured(Float_t ptMin = 0.301, const char* dataInput = "
   new TCanvas; gROOT->FindObject("mc_zy_div_esd_zy")->Draw("COLZ");
   new TCanvas; gROOT->FindObject("mc_x_div_esd_x")->Draw("COLZ");
   new TCanvas; gROOT->FindObject("mc_y_div_esd_y")->Draw("COLZ");
+
+  TH2* hist3 = (TH2*) dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->GetMeasuredHistogram()->Clone("mc2");
+  hist3->SetTitle("mc2");
+  Printf("mc event contains %f entries", hist3->Integral());
+  Printf("mc event contains %f entries in |vtx-z| < 10", hist3->Integral(hist3->GetXaxis()->FindBin(-9.9), hist3->GetXaxis()->FindBin(9.9), 1, hist3->GetNbinsY()));
+
+  TH2* hist4 = (TH2*) fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Clone("esd2");
+  hist4->SetTitle("esd2");
+  Printf("esd event contains %f entries", hist4->Integral());
+  Printf("esd event contains %f entries in |vtx-z| < 10", hist4->Integral(hist4->GetXaxis()->FindBin(-9.9), hist4->GetXaxis()->FindBin(9.9), 1, hist4->GetNbinsY()));
+  
+  ratio = (TH2*) hist3->Clone("ratio");
+  ratio->Divide(hist4);
+  
+  new TCanvas; ratio->Draw("COLZ");
 }
 
 void CompareCorrection2Generated(Float_t ptMin = 0.301, const char* dataInput = "analysis_mc.root", const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
@@ -2220,10 +2219,44 @@ void DrawTrackletOrigin(const char* fileName = "correction_map.root", Bool_t myF
   legend->Draw();
   
   c->SaveAs("spd_tracklets_efficiency.eps");
+}
+
+void DrawTrackletOrigin_Compare(const char* file1, const char* file2)
+{
+  DrawTrackletOrigin(file1);
+  good1 = (TH1*) gROOT->FindObject("good")->Clone("good1");
+  bad1 = (TH1*) gROOT->FindObject("bad")->Clone("bad1");
+
+  DrawTrackletOrigin(file2);
+  good2 = (TH1*) gROOT->FindObject("good")->Clone("good2");
+  bad2 = (TH1*) gROOT->FindObject("bad")->Clone("bad2");
+     
+  c = new TCanvas("c4", "c4", 600, 600);
+  c->SetRightMargin(0.05);
+  c->SetTopMargin(0.05);
+  c->SetGridx();
+  c->SetGridy();
+  gPad->SetLogy();
+  
+  good1->Draw();
+  bad1->SetLineColor(1);
+  bad1->SetMarkerColor(1);
+  bad1->Draw("SAME");
   
+  Float_t factor = (good1->Integral() + bad1->Integral()) / (good2->Integral() + bad2->Integral());
   
+  good2->Scale(factor);
+  bad2->Scale(factor);
+  
+  good2->SetLineColor(2);
+  bad2->SetMarkerColor(2);
+  
+  good2->Draw("SAME");
+  bad2->Draw("SAME");
+  
+  good1->GetYaxis()->SetRangeUser(1, TMath::Max(good1->GetMaximum(), good2->GetMaximum()) * 1.1);
 }
-
+  
 void Tracklets_Asymmetry()
 {
   TFile::Open("correction_map.root");
@@ -2390,7 +2423,7 @@ void GetAverageCorrectionFactor(Float_t etaRange = 1.5, Float_t vertexRange = 9.
   mcH->Fit("pol0", "", "", -etaRange, etaRange);
 }
 
-void TrackCuts_Comparison(char* histName, Int_t plotWhich = 0, const char* fileName = "correction_map.root", Bool_t mirror = kFALSE)
+void TrackCuts_Comparison_MC(char* histName, Int_t plotWhich = 0, const char* fileName = "correction_map.root", Bool_t mirror = kFALSE)
 {
   // for the nsigmaplot it is needed to run with all cuts except the nsigmatovertex
   //    --> manually disable it in the run.C
@@ -2501,6 +2534,89 @@ void TrackCuts_Comparison(char* histName, Int_t plotWhich = 0, const char* fileN
   //c1->SaveAs(Form("%s.png", histName));
 }
 
+void TrackCuts_Comparison_Data(char* histName, Int_t plotWhich, const char* fileName1, const char* fileName2, Bool_t mirror = kFALSE, const char* label1 = "file1", const char* label2 = "file2")
+{
+  // for the nsigmaplot it is needed to run with all cuts except the nsigmatovertex
+  //    --> manually disable it in the run.C
+  //
+  // plotWhich: 0 = only before
+  //            1 = both
+  //            2 = only after
+  //
+  // mirror: kTRUE --> project negative values on the positive side
+  
+
+  Int_t count = 0;
+  Int_t colors[] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 };
+
+  TLegend* legend = new TLegend(0.5, 0.7, 1, 1);
+  legend->SetFillColor(0);
+  TLegend* legend2 = new TLegend(0.4, 0.6, 1, 1);
+  TLegend* legend3 = new TLegend(0.6, 0.5, 1, 0.7);
+
+  TCanvas* c1 = new TCanvas(histName, histName, 600, 600);
+  //TCanvas* c2 = new TCanvas("c2", "c2", 800, 600);
+  //TCanvas* c3 = new TCanvas("c3", "c3", 800, 600);
+
+  const char* folders2[] = { "before_cuts", "after_cuts" };
+  Bool_t first = kTRUE;
+  for (Int_t j = ((plotWhich < 2) ? 0 : 1); j < ((plotWhich > 0) ? 2 : 1); j++)
+  {
+    const char* folders1[] = { "esd_track_cuts", "esd_track_cuts_primaries", "esd_track_cuts_secondaries" };
+    const char* names[] =    { "all", "primaries", "secondaries" };
+    
+    Float_t normalize[3];
+    
+    for (Int_t i = 0; i < 2; i++)
+    {
+      file = TFile::Open((i == 0) ? fileName1 : fileName2);
+      
+      for (Int_t k = 1; k < 3; k++)
+      {
+        TString folder;
+        folder.Form("%s/%s/%s", folders1[k], folders2[j], histName);
+        Printf("%s", folder.Data());
+        TH1* hist = (TH1*) file->Get(folder);
+        
+        if (mirror)
+        {
+          for (Int_t bin=1; bin<=hist->GetXaxis()->FindBin(-0.0001); bin++)
+          {
+            Int_t newBin = hist->GetXaxis()->FindBin(-hist->GetXaxis()->GetBinCenter(bin));
+            if (bin != newBin)
+            {
+              hist->Fill(-hist->GetXaxis()->GetBinCenter(bin), hist->GetBinContent(bin));
+              hist->SetBinContent(bin, 0);
+            }
+          }
+        }
+      
+        if (i == 0)
+        {
+          normalize[k] = hist->Integral();
+        }
+        else
+          hist->Scale(normalize[k] / hist->Integral());
+        
+        legend->AddEntry(hist, Form("%s %s %s", (i == 0) ? label1 : label2, (k == 1) ? "primaries" : "secondaries", folders2[j]));
+  
+        c1->cd();
+        hist->SetStats(0);
+        hist->SetLineColor(colors[count]);
+        hist->DrawCopy((count == 0) ? "" : "SAME");
+  
+        count++;
+      }
+    }
+
+  }
+
+  //c1->SetLogy();
+  c1->SetGridx();
+  c1->SetGridy();
+  legend->Draw();
+}
+
 void TrackCuts_DCA()
 {
   file = TFile::Open("correction_map.root");
@@ -2723,6 +2839,7 @@ void PrintEventStats(Int_t corrID = 3)
   eta->Scale(1.0 / eta->GetXaxis()->GetBinWidth(1));
 
   Printf("dndeta | eta = 0 is %f", (eta->GetBinContent(eta->FindBin(0.01)) + eta->GetBinContent(eta->FindBin(-0.01))) / 2);
+  Printf("dndeta in |eta| < 0.5 is %f", eta->Integral(eta->FindBin(-0.49), eta->FindBin(0.49)) / (eta->FindBin(0.49) - eta->FindBin(-0.49) + 1));
   Printf("dndeta in |eta| < 1 is %f", eta->Integral(eta->FindBin(-0.99), eta->FindBin(0.99)) / (eta->FindBin(0.99) - eta->FindBin(-0.99) + 1));
   Printf("dndeta in |eta| < 1.5 is %f", eta->Integral(eta->FindBin(-1.49), eta->FindBin(1.49)) / (eta->FindBin(1.49) - eta->FindBin(-1.49) + 1));
 
@@ -2734,11 +2851,19 @@ void PrintEventStats(Int_t corrID = 3)
   
   Printf("+++ TRIGGER EFFICIENCIES +++");
   
-  Printf("INEL = %.1f", 100. * (proj->GetBinContent(1) - stats->GetBinContent(1, 1)) / proj->GetBinContent(1));
-  Printf("NSD  = %.1f", 100. * (proj->GetBinContent(2) - stats->GetBinContent(2, 1)) / proj->GetBinContent(2));
-  Printf("ND  = %.1f",  100. * (proj->GetBinContent(3) - stats->GetBinContent(3, 1)) / proj->GetBinContent(3));
-  Printf("SD  = %.1f",  100. * (proj->GetBinContent(4) - stats->GetBinContent(4, 1)) / proj->GetBinContent(4));
-  Printf("DD  = %.1f",  100. * (proj->GetBinContent(5) - stats->GetBinContent(5, 1)) / proj->GetBinContent(5));
+  Printf("INEL = %.1f", 100. * (proj->GetBinContent(1) - stats->GetBinContent(1, 1) - stats->GetBinContent(1, 3)) / proj->GetBinContent(1));
+  Printf("NSD  = %.1f", 100. * (proj->GetBinContent(2) - stats->GetBinContent(2, 1) - stats->GetBinContent(2, 3)) / proj->GetBinContent(2));
+  Printf("ND  = %.1f",  100. * (proj->GetBinContent(3) - stats->GetBinContent(3, 1) - stats->GetBinContent(3, 3)) / proj->GetBinContent(3));
+  Printf("SD  = %.1f",  100. * (proj->GetBinContent(4) - stats->GetBinContent(4, 1) - stats->GetBinContent(4, 3)) / proj->GetBinContent(4));
+  Printf("DD  = %.1f",  100. * (proj->GetBinContent(5) - stats->GetBinContent(5, 1) - stats->GetBinContent(5, 3)) / proj->GetBinContent(5));
+  
+  Printf("+++ TRIGGER + VERTEX EFFICIENCIES +++");
+  
+  Printf("INEL = %.1f", 100. * stats->GetBinContent(1, 4) / proj->GetBinContent(1));
+  Printf("NSD  = %.1f", 100. * stats->GetBinContent(2, 4) / proj->GetBinContent(2));
+  Printf("ND  = %.1f",  100. * stats->GetBinContent(3, 4) / proj->GetBinContent(3));
+  Printf("SD  = %.1f",  100. * stats->GetBinContent(4, 4) / proj->GetBinContent(4));
+  Printf("DD  = %.1f",  100. * stats->GetBinContent(5, 4) / proj->GetBinContent(5));
   
   for (Int_t i=7; i<=proj->GetNbinsX(); i++)
     if (proj->GetBinContent(i) > 0)
@@ -2833,3 +2958,140 @@ dPhiHist2->Fill(-deltaPhi, hist->GetBinContent(i) / 2);
   tracklets->SetLineColor(4);
   tracklets->Draw("SAME");
 }
+
+void VertexDistributions()
+{
+  loadlibs();
+  
+  TFile::Open("correction_map.root");
+  AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
+  if (!dNdEtaCorrection->LoadHistograms())
+    return;
+  
+  all = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kINEL)->GetEventCorrection()->GetGeneratedHistogram()->ProjectionX("all");
+  trigger = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kINEL)->GetEventCorrection()->GetMeasuredHistogram()->ProjectionX("trigger");
+  vtx = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kVertexReco)->GetEventCorrection()->GetMeasuredHistogram()->ProjectionX("vtx");
+  nottriggered = (TH1*) all->Clone("nottriggered");
+  nottriggered->Add(trigger, -1);
+
+  novertex = (TH1*) trigger->Clone("novertex");
+  novertex->Add(vtx, -1);
+  
+  temphist = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kVertexReco)->GetEventCorrection()->GetMeasuredHistogram();
+  highmult = temphist->ProjectionX("highmult", temphist->GetYaxis()->FindBin(10), temphist->GetNbinsY());
+  //all = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kINEL)->GetEventCorrection()->GetGeneratedHistogram()->ProjectionX("all", temphist->GetYaxis()->FindBin(10), temphist->GetNbinsY());
+  for (Int_t i=1; i<=trigger->GetNbinsX(); i++)
+  {
+    all->SetBinContent(i, all->GetBinContent(i) / all->GetBinWidth(i));
+    trigger->SetBinContent(i, trigger->GetBinContent(i) / trigger->GetBinWidth(i));
+    vtx->SetBinContent(i, vtx->GetBinContent(i) / vtx->GetBinWidth(i));
+    nottriggered->SetBinContent(i, nottriggered->GetBinContent(i) / nottriggered->GetBinWidth(i));
+    novertex->SetBinContent(i, novertex->GetBinContent(i) / novertex->GetBinWidth(i));
+    highmult->SetBinContent(i, highmult->GetBinContent(i) / highmult->GetBinWidth(i));
+  }
+
+  new TCanvas;
+  vtx->SetTitle("");
+  vtx->SetStats(0);
+  vtx->DrawCopy("HIST");
+
+  all->Scale(1.0 / all->Integral());
+  nottriggered->Scale(1.0 / nottriggered->Integral());
+  novertex->Scale(1.0 / novertex->Integral());
+  highmult->Scale(1.0 / highmult->Integral());
+
+  new TCanvas;
+  all->Draw("HIST");
+  novertex->SetLineColor(2);
+  novertex->Draw("HISTSAME");
+  highmult->SetLineColor(3);
+  highmult->Draw("HISTSAME");
+
+  legend = new TLegend(0.5, 0.5, 0.8, 0.8);
+  legend->SetFillColor(0);
+  legend->AddEntry(all, "all");
+  legend->AddEntry(novertex, "no vertex");
+  legend->AddEntry(highmult, "mult > 10");
+  legend->Draw();
+  
+  new TCanvas;
+  trigger->Scale(1.0 / trigger->Integral());
+  vtx->Scale(1.0 / vtx->Integral());
+  
+  trigger->Divide(vtx);
+  
+  trigger->Draw();
+  //vtx->SetLineColor(2);
+  //vtx->Draw("SAME");
+
+  //temphist = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kVertexReco)->GetEventCorrection()->GetMeasuredHistogram();
+  temphist = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kINEL)->GetEventCorrection()->GetGeneratedHistogram();
+  //temphist = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kINEL)->GetEventCorrection()->GetMeasuredHistogram();
+  
+  new TCanvas;
+  legend = new TLegend(0.7, 0.7, 0.99, 0.99);
+  legend->SetFillColor(0);
+  Bool_t first = kTRUE; 
+  for (Int_t i=0; i<20; i+=5)
+  {
+    highmult = temphist->ProjectionX("highmult", i+1, i+1+4);
+    Printf("%f", highmult->Integral());
+    if (highmult->Integral() <= 0)
+      continue;
+  
+    for (Int_t j=1; j<=trigger->GetNbinsX(); j++)
+      highmult->SetBinContent(j, highmult->GetBinContent(j) / highmult->GetBinWidth(j));
+
+    highmult->Scale(1.0 / highmult->Integral());
+    highmult->SetLineColor((i/5)+1);
+    highmult->GetYaxis()->SetRangeUser(0, 0.15);
+    if (first)
+    {
+      highmult->DrawCopy();
+      first = kFALSE;
+    }
+    else
+      highmult->DrawCopy("SAME");
+    legend->AddEntry(highmult->Clone(), Form("%d <= N <= %d", i, i+4));
+  }
+  legend->Draw();
+}
+
+void PlotPt1DCorrection()
+{
+  const char* files[] = { "field.root", "field_onlyprim.root", "nofield.root", "nofield_onlyprim.root" };
+  const char* names[] = { "B: all", "B: primaries", "No B: all", "No B: primaries" };
+  Int_t colors[] = { 1, 2, 3, 4 };
+  
+  loadlibs();
+  
+  dummy = new TH2F("dummy", ";p_{T};correction", 100, 0, 1.4, 100, 0.5, 3);
+  dummy->SetStats(0);
+  //dummy->GetYaxis()->SetTitleOffset(1.3);
+  dummy->Draw();
+  
+  legend = new TLegend(0.48, 0.57, 0.88, 0.88);
+  legend->SetFillColor(0);
+  
+  for (Int_t i=0; i<4; i++)
+  {
+    TFile::Open(files[i]);
+    AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
+    if (!dNdEtaCorrection->LoadHistograms())
+      return;
+      
+    hist = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->Get1DCorrectionHistogram("z", -9.9, 9.9, -0.79, 0.79);
+    hist->SetLineColor(colors[i]);
+    hist->SetLineWidth(2);
+    hist->SetMarkerColor(colors[i]);
+    hist->Draw("SAME");
+    
+    legend->AddEntry(hist, names[i], "L");
+  }
+  
+  legend->Draw();
+}
index 29a114a..7283618 100644 (file)
@@ -798,8 +798,10 @@ void mergeCorrectionsWithDifferentCrosssections(Int_t correctionTarget = 3, Char
 
   const Char_t* changes[]  = { "pythia","ddmore","ddless","sdmore","sdless", "dmore", "dless", "sdlessddmore", "sdmoreddless", "ddmore25","ddless25","sdmore25","sdless25", "dmore25", "dless25", "sdlessddmore25", "sdmoreddless25"};
   Float_t scalesND[] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0 };
-  Float_t scalesDD[] = {1.0, 1.5, 0.5, 1.0, 1.0, 1.5, 0.5, 1.5, 0.5, 1.25, 0.75, 1.0,  1.0,  1.25, 0.75, 1.25, 0.75};
-  Float_t scalesSD[] = {1.0, 1.0, 1.0, 1.5, 0.5, 1.5, 0.5, 0.5, 1.5, 1.0,  1.0,  1.25, 0.75, 1.25, 0.75, 0.75, 1.25};
+  //Float_t scalesDD[] = {1.0, 1.5, 0.5, 1.0, 1.0, 1.5, 0.5, 1.5, 0.5, 1.25, 0.75, 1.0,  1.0,  1.25, 0.75, 1.25, 0.75};
+  //Float_t scalesSD[] = {1.0, 1.0, 1.0, 1.5, 0.5, 1.5, 0.5, 0.5, 1.5, 1.0,  1.0,  1.25, 0.75, 1.25, 0.75, 0.75, 1.25};
+  Float_t scalesDD[] = {1.0, 1.4, 0.6, 1.0, 1.0, 1.4, 0.6, 1.4, 0.6, 1.25, 0.75, 1.0,  1.0,  1.25, 0.75, 1.25, 0.75};
+  Float_t scalesSD[] = {1.0, 1.0, 1.0, 1.4, 0.6, 1.4, 0.6, 0.4, 1.6, 1.0,  1.0,  1.25, 0.75, 1.25, 0.75, 0.75, 1.25};
   Int_t nChanges = 17;
 
   /*
@@ -1545,3 +1547,44 @@ void runStudy(const char* baseCorrectionMapFile = "correction_map.root", const c
 
   legend->Draw();
 }
+
+void ChangePtInCorrection(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction")
+{
+  loadlibs();
+  if (!TFile::Open(fileName))
+    return;
+
+  AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(dirName, dirName);
+  if (!dNdEtaCorrection->LoadHistograms())
+    return;
+
+  dNdEtaCorrection->Finish();
+
+  AliCorrection* corr = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kTrack2Particle);
+  Printf(">>>> Before");
+  corr->PrintInfo(0);
+
+  Float_t factor = 0.5;
+  Float_t ptCutOff = 0.2;
+  
+  TH3* gene = corr->GetTrackCorrection()->GetGeneratedHistogram();
+  TH3* meas = corr->GetTrackCorrection()->GetMeasuredHistogram();
+  
+  for (Int_t z = 1; z <= gene->GetZaxis()->FindBin(ptCutOff - 0.01); z++)
+  {
+    Float_t localFactor = 1 - (ptCutOff - gene->GetZaxis()->GetBinCenter(z)) / ptCutOff * factor;
+    Printf("%f %f", gene->GetZaxis()->GetBinCenter(z), localFactor);
+    for (Int_t x = 1; x <= gene->GetNbinsX(); x++)
+      for (Int_t y = 1; y <= gene->GetNbinsY(); y++)
+      {
+        gene->SetBinContent(x, y, z, gene->GetBinContent(x, y, z) * localFactor);
+        meas->SetBinContent(x, y, z, meas->GetBinContent(x, y, z) * localFactor);
+      }
+  }
+  
+  dNdEtaCorrection->Finish();
+  
+  Printf(">>>> After");
+  corr->PrintInfo(0);
+}
index cf73e98..fbc62e0 100644 (file)
@@ -114,6 +114,59 @@ TPad* DrawChange(Bool_t spd, const char* basename, const char** changes, Int_t n
   return p;
 }
 
+TPad* DrawChangeRatios(Bool_t spd, TFile* file1, TFile* file2, const char* basename, const char** changes, Int_t nChanges, Int_t nDraw, Int_t* colors, const char** names = 0, Float_t scale = 0.10)
+{  
+  Float_t etaMax = 1.05;
+  if (spd)
+    etaMax = 1.79;
+
+  TH1F* hRatios[100];
+  for(Int_t i=0; i<nChanges; i++) {
+    hRatios[i] = (TH1F*)file1->Get(Form("%s%s",basename,changes[i]));
+    hRatios[i]->SetLineWidth(1);
+    hRatios[i]->SetMarkerStyle(22);
+    hRatios[i]->SetMarkerSize(0.8);
+    
+    hRatio2 = (TH1F*)file2->Get(Form("%s%s",basename,changes[i]));
+    hRatios[i]->Divide(hRatio2);
+
+    Float_t average = hRatios[i]->Integral(hRatios[i]->FindBin(-1), hRatios[i]->FindBin(1)) / (hRatios[i]->FindBin(1) - hRatios[i]->FindBin(-1) + 1);
+    Printf("%s: %.2f %%" , hRatios[i]->GetTitle(), (average - 1) * 100);
+  }
+  
+  TPad* p = DrawCanvasAndPad("syst_changeInXsection",700,400);
+  p->SetRightMargin(0.2);
+  p->SetLeftMargin(0.13);
+  
+  TH2F* null = new TH2F("","",100,-etaMax, etaMax,100,1. - scale,1. + scale);
+  null->GetXaxis()->SetTitle("#eta");
+  null->GetYaxis()->SetTitle(hRatios[0]->GetYaxis()->GetTitle());
+  null->Draw();
+  
+  line = new TLine(-etaMax, 1, etaMax, 1);
+  line->Draw();
+
+  TLatex* text[100];
+
+  for(Int_t i=1; i<nDraw; i++) {
+    hRatios[i]->SetLineColor(colors[i]);
+    hRatios[i]->SetMarkerColor(colors[i]);
+    hRatios[i]->Draw("HISTPL SAME");
+    
+    TString str(hRatios[i]->GetTitle());
+    
+    if (names)
+      str = names[i];
+    
+    text[i] = new TLatex(etaMax + 0.03,hRatios[i]->GetBinContent(hRatios[i]->FindBin(etaMax-0.1))-0.002,str.Data());
+    text[i]->SetTextAlign(11);
+    text[i]->SetTextColor(colors[i]);
+    text[i]->Draw();
+  }
+  
+  return p;
+}
+
 void DrawEffectOfChangeInCrossSection(Bool_t spd = kFALSE, const char* fileName = "systematics_vtxtrigger_compositions.root") 
 {
   TFile::Open(fileName);
@@ -127,6 +180,20 @@ void DrawEffectOfChangeInCrossSection(Bool_t spd = kFALSE, const char* fileName
   c->SaveAs("cross_sections.eps");
 }
 
+void DrawEffectOfChangeInCrossSection2Files(Bool_t spd = kFALSE, const char* fileName1 = "systematics_vtxtrigger_compositions.root", const char* fileName2) 
+{
+  file1 = TFile::Open(fileName1);
+  file2 = TFile::Open(fileName2);
+
+  const Char_t* changes[]  = {"pythia","ddmore","ddless","sdmore","sdless", "dmore", "dless", "sdmoreddless", "sdlessddmore", "ddmore25","ddless25","sdmore25","sdless25", "dmore25", "dless25", "sdmoreddless25", "sdlessddmore25" };
+  //const Char_t* changes[]  = { "pythia", "qgsm", "phojet" };
+  //const Int_t nChanges = 3;
+  Int_t colors[] = {1,1,4,1,2,2,4,2,1};
+
+  c = DrawChangeRatios(spd, file1, file2, "ratio_vertexReco_triggerBias_", changes, 17, 9, colors, 0);
+  c->SaveAs("cross_sections.eps");
+}
+
 void DrawEffectOfChangeInComposition(Bool_t spd = kFALSE, const char* fileName = "new_compositions_analysis.root") 
 {
   TFile::Open(fileName);
index fd90f09..f361cfc 100644 (file)
@@ -28,6 +28,10 @@ void run(Int_t runWhat, const Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool
   // aProof option: 0 no proof
   //                1 proof with chain
   //                2 proof with dataset
+  //
+  // option is passed to the task(s)
+  //   option SAVE is removed and results in moving the output files to maps/<ds name>/<trigger>/<det>
+  //
   
   TString taskName;
   if (runWhat == 0 || runWhat == 2)
@@ -49,7 +53,7 @@ void run(Int_t runWhat, const Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool
 
   if (aProof)
   {
-    TProof::Open("alicecaf");
+    TProof::Open("alicecaf"); 
     //gProof->SetParallel(1);
 
     // Enable the needed package
@@ -68,8 +72,8 @@ void run(Int_t runWhat, const Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool
     }
     else
     {
-      gProof->UploadPackage("/afs/cern.ch/alice/caf/sw/ALICE/PARs/v4-16-Release/AF-v4-16");
-      gProof->EnablePackage("AF-v4-16");
+      gProof->UploadPackage("/afs/cern.ch/alice/caf/sw/ALICE/PARs/v4-17-Release/AF-v4-17");
+      gProof->EnablePackage("AF-v4-17");
     }
 
     gProof->UploadPackage("$ALICE_ROOT/PWG0base");
@@ -97,13 +101,13 @@ void run(Int_t runWhat, const Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool
   esdH->SetInactiveBranches("AliESDACORDE FMD ALIESDTZERO ALIESDVZERO ALIESDZDC AliRawDataErrorLogs CaloClusters Cascades EMCALCells EMCALTrigger ESDfriend Kinks Kinks Cascades AliESDTZERO ALIESDACORDE MuonTracks TrdTracks CaloClusters");
   mgr->SetInputEventHandler(esdH);
 
-  AliPWG0Helper::AnalysisMode analysisMode = AliPWG0Helper::kSPD;
+  AliPWG0Helper::AnalysisMode analysisMode = AliPWG0Helper::kTPC | AliPWG0Helper::kFieldOn;
   AliPWG0Helper::Trigger      trigger      = AliPWG0Helper::kMB1;
 
   AliPWG0Helper::PrintConf(analysisMode, trigger);
 
   AliESDtrackCuts* esdTrackCuts = 0;
-  if (analysisMode != AliPWG0Helper::kSPD)
+  if (!(analysisMode & AliPWG0Helper::kSPD))
   {
     // selection of esd tracks
     gROOT->ProcessLine(".L ../CreateStandardCuts.C");
@@ -116,13 +120,22 @@ void run(Int_t runWhat, const Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool
     esdTrackCuts->SetHistogramsOn(kTRUE);
   }
 
-  cInput  = mgr->GetCommonInputContainer();
-
+  cInput = mgr->GetCommonInputContainer();
+  
+  // remove SAVE option if set
+  Bool_t save = kFALSE;
+  TString optStr(option);
+  if (optStr.Contains("SAVE"))
+  {
+    optStr = optStr(0,optStr.Index("SAVE")) + optStr(optStr.Index("SAVE")+4, optStr.Length());
+    save = kTRUE;
+  }
+  
   // Create, add task
   if (runWhat == 0 || runWhat == 2)
   {
     Load("AlidNdEtaTask", aDebug);
-    task = new AlidNdEtaTask(option);
+    task = new AlidNdEtaTask(optStr);
 
     if (mc)
       task->SetReadMC();
@@ -136,7 +149,7 @@ void run(Int_t runWhat, const Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool
     task->SetTrigger(trigger);
     task->SetAnalysisMode(analysisMode);
     task->SetTrackCuts(esdTrackCuts);
-    task->SetDeltaPhiCut(0.05);
+    //task->SetDeltaPhiCut(0.05);
 
     mgr->AddTask(task);
 
@@ -150,7 +163,7 @@ void run(Int_t runWhat, const Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool
   if (runWhat == 1 || runWhat == 2)
   {
     Load("AlidNdEtaCorrectionTask", aDebug);
-    task2 = new AlidNdEtaCorrectionTask(option);
+    task2 = new AlidNdEtaCorrectionTask(optStr);
 
     // syst. error flags
     //task2->SetFillPhi();
@@ -191,6 +204,42 @@ void run(Int_t runWhat, const Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool
     // process dataset
 
     mgr->StartAnalysis("proof", data, nRuns, offset);
+    
+    if (save)
+    {
+      TString path("maps/");
+      path += TString(data).Tokenize("/")->Last()->GetName();
+      
+      switch (trigger)
+      {
+        case AliPWG0Helper::kMB1:
+        case AliPWG0Helper::kOfflineMB1: path += "/mb1"; break;
+        case AliPWG0Helper::kMB2:
+        case AliPWG0Helper::kOfflineMB2: path += "/mb2"; break;
+        case AliPWG0Helper::kMB3:
+        case AliPWG0Helper::kOfflineMB3: path += "/mb3"; break;
+        case AliPWG0Helper::kFASTOR:
+        case AliPWG0Helper::kOfflineFASTOR: path += "/fastor"; break;
+        default: Printf("ERROR: Trigger undefined for path to files"); return;
+      }
+      
+      if (analysisMode & AliPWG0Helper::kSPD)
+        path += "/spd";
+      
+      if (analysisMode & AliPWG0Helper::kTPC)
+        path += "/tpc";
+        
+      gSystem->mkdir(path, kTRUE);
+      if (runWhat == 0 || runWhat == 2)
+        gSystem->Rename("analysis_esd_raw.root", path + "/analysis_esd_raw.root");
+      if (runWhat == 1 || runWhat == 2)
+      {
+        gSystem->Rename("analysis_mc.root", path + "/analysis_mc.root");
+        gSystem->Rename("correction_map.root", path + "/correction_map.root");
+      }
+      
+      Printf(">>>>> Moved files to %s", path.Data());
+    }
   }
   else if (aProof == 3)
   {