- Fixed double counting of events by not calliing
authorhdalsgaa <hdalsgaa@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 11 Jan 2011 13:44:16 +0000 (13:44 +0000)
committerhdalsgaa <hdalsgaa@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 11 Jan 2011 13:44:16 +0000 (13:44 +0000)
  AliPhysicsSelection::IsCollisionCandidate in the event inspector
- Fix stupid mistake
- New cuts, etc.
- Better processing of aux data
- Store all (that's _all_) events.

PWG2/FORWARD/analysis2/AddTaskFMD.C
PWG2/FORWARD/analysis2/AliFMDEventInspector.cxx
PWG2/FORWARD/analysis2/AliForwardMCMultiplicityTask.cxx
PWG2/FORWARD/analysis2/AnalyseAOD.C
PWG2/FORWARD/analysis2/Run.sh

index 958e3bb..cc940e1 100644 (file)
@@ -70,10 +70,10 @@ AddTaskFMD(Bool_t mc)
   // Whether to use the merging efficiency correction 
   task->GetCorrections().UseMergingEfficiency(false);
   // Set the number of extra bins (beyond the secondary map border) 
-  task->GetHistCollector().SetNCutBins(1);
+  task->GetHistCollector().SetNCutBins(2);
   // Set the correction cut, that is, when bins in the secondary map 
   // is smaller than this, they are considered empty 
-  task->GetHistCollector().SetCorrectionCut(0.1);
+  task->GetHistCollector().SetCorrectionCut(0.5);
   // Set the overall debug level (1: some output, 3: a lot of output)
   task->SetDebug(0);
   // Set the debug level of a single algorithm 
index d586186..727ec3e 100644 (file)
@@ -367,7 +367,16 @@ AliFMDEventInspector::ReadTriggers(const AliESDEvent* esd, UInt_t& triggers)
     AliWarning("No input handler");
     return kFALSE;
   }
-  
+
+#if 1
+  // Check if this is a collision candidate (INEL)
+  // Note, that we should use the value cached in the input 
+  // handler rather than calling IsCollisionCandiate directly 
+  // on the AliPhysicsSelection obejct.  If we called the latter
+  // then the AliPhysicsSelection object would overcount by a 
+  // factor of 2! :-(
+  Bool_t inel = ih->IsEventSelected();
+#else 
   // Get the physics selection - add that by using the macro 
   // AddTaskPhysicsSelection.C 
   AliPhysicsSelection* ps = 
@@ -379,18 +388,22 @@ AliFMDEventInspector::ReadTriggers(const AliESDEvent* esd, UInt_t& triggers)
   
   // Check if this is a collision candidate (INEL)
   Bool_t inel = ps->IsCollisionCandidate(esd);
+#endif
   if (inel) { 
     triggers |= AliAODForwardMult::kInel;
     fHTriggers->Fill(kInel+0.5);
   }
 
-  // IF this is inel, see if we have a tracklet 
+  // If this is inel, see if we have a tracklet 
   if (inel) { 
     const AliMultiplicity* spdmult = esd->GetMultiplicity();
     if (!spdmult) {
       AliWarning("No SPD multiplicity");
     }
     else { 
+      // Check if we have one or more tracklets 
+      // in the range -1 < eta < 1 to set the INEL>0 
+      // trigger flag. 
       Int_t n = spdmult->GetNumberOfTracklets();
       for (Int_t j = 0; j < n; j++) { 
        if(TMath::Abs(spdmult->GetEta(j)) < 1) { 
index abb0d1a..3a26160 100644 (file)
@@ -321,21 +321,29 @@ AliForwardMCMultiplicityTask::UserExec(Option_t*)
   UShort_t ivz      = 0;
   Double_t vz       = 0;
   UInt_t   found    = fEventInspector.Process(esd, triggers, lowFlux, ivz, vz);
-  if (found & AliFMDEventInspector::kNoEvent)    return;
-  if (found & AliFMDEventInspector::kNoTriggers) return;
+  //Store all events
+  MarkEventForStore();
+  
+  Bool_t isAccepted = true;
+  if (found & AliFMDEventInspector::kNoEvent)    isAccepted = false; // return;
+  if (found & AliFMDEventInspector::kNoTriggers) isAccepted = false; // return;
  
   // Set trigger bits, and mark this event for storage 
   fAODFMD.SetTriggerBits(triggers);
   fMCAODFMD.SetTriggerBits(triggers);
-  MarkEventForStore();
 
-  if (found & AliFMDEventInspector::kNoSPD)     return;
-  if (found & AliFMDEventInspector::kNoFMD)     return;
-  if (found & AliFMDEventInspector::kNoVertex)  return;
-  fAODFMD.SetIpZ(vz);
-  fMCAODFMD.SetIpZ(vz);
+  //All events should be stored - HHD
+  //MarkEventForStore();
 
-  if (found & AliFMDEventInspector::kBadVertex) return;
+  if (found & AliFMDEventInspector::kNoSPD)     isAccepted = false; // return;
+  if (found & AliFMDEventInspector::kNoFMD)     isAccepted = false; // return;
+  if (found & AliFMDEventInspector::kNoVertex)  isAccepted = false; // return;
+
+  if (isAccepted) {
+    fAODFMD.SetIpZ(vz);
+    fMCAODFMD.SetIpZ(vz);
+  }
+  if (found & AliFMDEventInspector::kBadVertex) isAccepted = false; // return;
 
   // We we do not want to use low flux specific code, we disable it here. 
   if (!fEnableLowFlux) lowFlux = false;
@@ -344,7 +352,7 @@ AliForwardMCMultiplicityTask::UserExec(Option_t*)
   AliESDFMD*  esdFMD  = esd->GetFMDData();
   AliMCEvent* mcEvent = MCEvent();
   // Apply the sharing filter (or hit merging or clustering if you like)
-  if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD)) { 
+  if (isAccepted && !fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD)) { 
     AliWarning("Sharing filter failed!");
     return;
   }
@@ -352,6 +360,7 @@ AliForwardMCMultiplicityTask::UserExec(Option_t*)
     AliWarning("MC Sharing filter failed!");
     return;
   }
+  if (!isAccepted) return; // Exit on MC event w/o trigger, vertex, data
   fSharingFilter.CompareResults(fESDFMD, fMCESDFMD);
 
   // Do the energy stuff 
index c05326b..4cb7746 100644 (file)
@@ -71,7 +71,9 @@ public:
   Int_t              fNWithVertex;
   /** Number of events accepted */
   Int_t              fNAccepted;
-
+  /** Number of events accepted */
+  Int_t              fNAll;
+    
   //__________________________________________________________________
   /** 
    * Constructor 
@@ -257,7 +259,12 @@ public:
             fMCSum->GetYaxis()->GetXmin(),
             fMCSum->GetYaxis()->GetXmax());
       }
+      
+      // Add contribution from this event 
+      if (fSumPrimary) fSumPrimary->Add(fPrimary);
+     
+      fNAll++;
+      
       // fAOD->Print();
       // Other trigger/event requirements could be defined 
       if (!fAOD->IsTriggerBits(mask)) continue; 
@@ -277,8 +284,7 @@ public:
       // Add contribution from this event
       if (fMCSum) fMCSum->Add(&(fMCAOD->GetHistogram()));
 
-      // Add contribution from this event 
-      if (fSumPrimary) fSumPrimary->Add(fPrimary);
+      
     }
     printf("\n");
     fVtxEff = Double_t(fNWithVertex)/fNTriggered;
@@ -345,8 +351,9 @@ public:
     if (fSumPrimary) { 
       dndetaTruth = fSumPrimary->ProjectionX("dndetaTruth", -1, -1, "e");
       Info("Finish", "Got dn/deta truth with max: %f, scalling to %d", 
-          dndetaTruth->GetMaximum(), fNTriggered);
-      dndetaTruth->Scale(1./fNTriggered, "width");
+          dndetaTruth->GetMaximum(), fNAll);
+      //dndetaTruth->Scale(1./fNTriggered, "width");
+      dndetaTruth->Scale(1./fNAll, "width");
       dndetaTruth->SetMarkerColor(kGray+3);
       dndetaTruth->SetMarkerStyle(22);
       dndetaTruth->SetMarkerSize(1);
@@ -354,7 +361,7 @@ public:
       Rebin(dndetaTruth, rebin);
     }
 
-    DrawIt(dndeta, dndetaMC, dndetaTruth, mask, energy, doHHD, doComp, rebin);
+    DrawIt(dndeta, dndetaMC, dndetaTruth, mask, energy, doHHD, doComp);
 
     return kTRUE;
   }
@@ -363,60 +370,44 @@ public:
    */
   void DrawIt(TH1* dndeta, TH1* dndetaMC, TH1* dndetaTruth, 
              Int_t mask, Int_t energy, Bool_t doHHD, 
-             Bool_t doComp, Int_t rebin)
+             Bool_t doComp)
   {
     // --- 1st part - prepare data -----------------------------------
-    TH1* sym = Symmetrice(dndeta);
-    // Rebin(sym, rebin);
+    TH1* dndetaSym = Symmetrice(dndeta);
 
     Double_t max = dndeta->GetMaximum();
 
     // Make our histogram stack 
     THStack* stack = new THStack("results", "Results");
 
-    TH1* mc = 0;
-    TH1* mcsym = 0;
-#if 1
-    mc = dndetaTruth;
-#else 
-    // If available, plot the MC truth 
-    if (!gSystem->AccessPathName("AnalysisResults.root")) { 
-      TFile* anares = TFile::Open("AnalysisResults.root", "READ");
-      if (anares) { 
-       TList* list = static_cast<TList*>(anares->Get("Forward"));
-       if (list) { 
-         mc = static_cast<TH1*>(list->FindObject("mcSumEta"));
-         Rebin(mc, rebin);
-       }
-       anares->Close();
-       fOut->cd();
-      }
-    }
-#endif
-    if (mc) {
-      mcsym = Symmetrice(mc);
-      stack->Add(mc);
-      stack->Add(mcsym);
+    TH1* dndetaTruthSym = 0;
+    if (dndetaTruth) {
+      dndetaTruth->SetFillColor(kGray);
+      dndetaTruth->SetFillStyle(3001);
+      dndetaTruthSym = Symmetrice(dndetaTruth);
+      stack->Add(dndetaTruthSym, "e5");
+      stack->Add(dndetaTruth, "e5");
       Info("DrawIt", "Maximum of MC dndeta (truth): %f, was %f", 
-          mc->GetMaximum(),dndeta->GetMaximum());
-      max = TMath::Max(mc->GetMaximum(),max);
+          dndetaTruth->GetMaximum(),dndeta->GetMaximum());
+      max = TMath::Max(dndetaTruth->GetMaximum(),max);
     }
 
     // Get the data from HHD's analysis - if available 
-    TH1* hhd    = 0;
-    TH1* hhdsym = 0;
+    TH1* dndetaHHD    = 0;
+    TH1* dndetaHHDSym = 0;
     Info("DrawIt", "Will %sdraw HHD results", (doHHD ? "" : "not "));
     if (doHHD) {
       TString hhdName(fOut->GetName());
       hhdName.ReplaceAll("dndeta", "hhd");    
-      hhd    = GetHHD(hhdName.Data(), mask & AliAODForwardMult::kNSD);
-      hhdsym = 0;
-      if (hhd) { 
+      dndetaHHD    = GetHHD(hhdName.Data(), mask & AliAODForwardMult::kNSD);
+      dndetaHHDSym = 0;
+      if (dndetaHHD) { 
        // Symmetrice and add to stack 
-       hhdsym = Symmetrice(hhd);
-       stack->Add(hhd);
-       stack->Add(hhdsym);
-       max = TMath::Max(hhd->GetMaximum(),max);
+       dndetaHHD->SetTitle("ALICE Forward (HHD)");
+       dndetaHHDSym = Symmetrice(dndetaHHD);
+       stack->Add(dndetaHHDSym);
+       stack->Add(dndetaHHD);
+       max = TMath::Max(dndetaHHD->GetMaximum(),max);
       }
       else 
        Warning("DrawIt", "No HHD data found");
@@ -427,22 +418,24 @@ public:
     if (dndetaMC) { 
       Info("DrawIt", "Maximum of MC dndeta: %f, was %f", 
           dndetaMC->GetMaximum(),dndeta->GetMaximum());
-      TH1* mcSym = Symmetrice(dndetaMC);
+      TH1* dndetaMCSym = Symmetrice(dndetaMC);
+      stack->Add(dndetaMCSym);
       stack->Add(dndetaMC);
-      stack->Add(mcSym);
       max = TMath::Max(dndetaMC->GetMaximum(),max);
     }
 
     // Add the analysis results to the list 
+    stack->Add(dndetaSym);
     stack->Add(dndeta);
-    stack->Add(sym);
 
     // Get graph of 'other' data - e.g., UA5, CMS, ... - and check if
     // there's any graphs.  Set the pad division based on that.
     Info("DrawIt", "Will %sdraw other results", (doComp ? "" : "not "));
     TMultiGraph* other    = (doComp ? GetData(energy, mask) : 0);
-    THStack*     ratios   = MakeRatios(dndeta, sym, hhd, hhdsym, 
-                                      mc, mcsym, other);
+    THStack*     ratios   = MakeRatios(dndeta,      dndetaSym, 
+                                      dndetaHHD,   dndetaHHDSym, 
+                                      dndetaTruth, dndetaTruthSym, 
+                                      other);
 
     // Check if we have ratios 
 
@@ -569,7 +562,7 @@ public:
 
       // Make a nice band from 0.9 to 1.1 
       TGraphErrors* band = new TGraphErrors(2);
-      band->SetPoint(0, sym->GetXaxis()->GetXmin(), 1);
+      band->SetPoint(0, dndetaSym->GetXaxis()->GetXmin(), 1);
       band->SetPoint(1, dndeta->GetXaxis()->GetXmax(), 1);
       band->SetPointError(0, 0, .1);
       band->SetPointError(1, 0, .1);
@@ -835,6 +828,9 @@ public:
       s->SetBinContent(j, h->GetBinContent(i));
       s->SetBinError(j, h->GetBinError(i));
     }
+    // Fill in overlap bin 
+    s->SetBinContent(l2+1, h->GetBinContent(first));
+    s->SetBinError(l2+1, h->GetBinError(first));
     return s;
   }
   //__________________________________________________________________
index b5a33db..3103fc3 100755 (executable)
@@ -68,7 +68,7 @@ while test $# -gt 0 ; do
        -M|--mc)             mc=`toggle $mc`   ;; 
        -g|--gdb)            gdb=`toggle $gdb`   ;; 
        -H|--hhd)            hhd=`toggle $hhd`   ;; 
-       -O|--other)          other=`toggle $other`   ;; 
+       -O|--other)          comp=`toggle $comp`   ;; 
        -r|--rebin)          rebin=$2         ; shift ;; 
        -v|--vz-min)         vzmin=$2         ; shift ;; 
        -V|--vz-max)         vzmax=$2         ; shift ;;