Fixes for pA indenfication of events
[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / AliFMDMCEventInspector.cxx
index 3e64b651df92d29339405c068ad6b9381ea1c885..fe5867842e349d98ea71888b8024cb1354210530 100644 (file)
@@ -29,6 +29,8 @@
 #include "AliGenEposEventHeader.h"
 #include "AliGenHerwigEventHeader.h"
 #include "AliGenGeVSimEventHeader.h"
+#include "AliAnalysisManager.h"
+#include "AliMCEventHandler.h"
 #include "AliHeader.h"
 #include <TList.h>
 #include <TH2F.h>
@@ -284,6 +286,7 @@ AliFMDMCEventInspector::ProcessMC(AliMCEvent*       event,
                                  UShort_t&         ivz, 
                                  Double_t&         vz,
                                  Double_t&         b,
+                                 Double_t&         c,
                                  Int_t&            npart, 
                                  Int_t&            nbin,
                                  Double_t&         phiR)
@@ -298,6 +301,7 @@ AliFMDMCEventInspector::ProcessMC(AliMCEvent*       event,
   //                  means outside of the defined vertex range
   //   vz        On return, the z position of the interaction
   //   b         On return, the impact parameter - < 0 if not available
+  //   c         On return, centrality estimate - < 0 if not available 
   //   phiR      On return, the event plane angle - < 0 if not available 
   // 
   // Return:
@@ -336,6 +340,7 @@ AliFMDMCEventInspector::ProcessMC(AliMCEvent*       event,
   npart          = 0;
   nbin           = 0;
   b              = -1;
+  c              = -1;
   if (colGeometry) { 
     b     = colGeometry->ImpactParameter();
     phi   = colGeometry->ReactionPlaneAngle();
@@ -343,6 +348,11 @@ AliFMDMCEventInspector::ProcessMC(AliMCEvent*       event,
             colGeometry->TargetParticipants());
     nbin  = colGeometry->NN();
   }
+  if (fDebug && !colGeometry) { 
+    AliWarningF("Collision header of class %s is not a CollisionHeader", 
+               genHeader->ClassName());
+  }
+    
   if(pythiaHeader) {
     Int_t pythiaType = pythiaHeader->ProcessType();
     // 92 and 93 are SD 
@@ -352,6 +362,15 @@ AliFMDMCEventInspector::ProcessMC(AliMCEvent*       event,
     npart = 2; // Always 2 protons
     nbin  = 1; // Always 1 binary collision 
   }
+  if (b >= 0) { 
+    if      (b <  3.5)  c = 2.5; //0-5%
+    else if (b <  4.95) c = 7.5; //5-10%
+    else if (b <  6.98) c = 15; //10-20%
+    else if (b <  8.55) c = 25; //20-30%
+    else if (b <  9.88) c = 35; //30-40%
+    else if (b < 11.04) c = 45; //40-50%
+    else                c = 55; //50-60%
+  }
   if(dpmHeader) { // Also an AliCollisionGeometry 
     Int_t processType = dpmHeader->ProcessType();
     // 1 & 4 are ND 
@@ -406,12 +425,28 @@ AliFMDMCEventInspector::ProcessMC(AliMCEvent*       event,
   genHeader->PrimaryVertex(vtx);
   vz = vtx[2];
 
+  if (fDebug) { 
+    AliInfoF("vz=%f, phiR=%f, b=%f, npart=%d, nbin=%d", 
+            vz, phiR, b, npart, nbin);
+  }
   fHVertex->Fill(vz);
   fHPhiR->Fill(phiR);
   fHB->Fill(b);
   fHBvsPart->Fill(b, npart);
   fHBvsBin->Fill(b, nbin);
 
+  if(fUseDisplacedVertices) {
+    // Put the vertex at fixed locations 
+    Double_t zvtx  = vz;
+    Double_t ratio = zvtx/37.5;
+    if(ratio > 0) ratio = ratio + 0.5;
+    if(ratio < 0) ratio = ratio - 0.5;
+    Int_t ratioInt = Int_t(ratio);
+    zvtx = 37.5*((Double_t)ratioInt);
+    if(TMath::Abs(zvtx) > 999) 
+      return kBadVertex;
+  }
+
   // Check for the vertex bin 
   ivz = fVtxAxis.FindBin(vz);
   if (ivz <= 0 || ivz > fHEventsTr->GetXaxis()->GetNbins()) { 
@@ -516,43 +551,14 @@ AliFMDMCEventInspector::IsSingleDiffractive(AliStack* stack,
     return false;
   
   // Rapidity shift
-  Double_t m02s = 1 - 2 * p1->Energy() / fEnergy
-  Double_t m12s = 1 - 2 * p2->Energy() / fEnergy;
+  Double_t m02s = (fEnergy > 0 ? 1 - 2 * p1->Energy() / fEnergy : 0)
+  Double_t m12s = (fEnergy > 0 ? 1 - 2 * p2->Energy() / fEnergy : 0);
   
   if (arm == 0 && m02s > xiMin && m02s < xiMax) return true;
   if (arm == 1 && m12s > xiMin && m12s < xiMax) return true;
     
   return false;
 }
-//____________________________________________________________________
-Bool_t
-AliFMDMCEventInspector::ReadCentrality(const AliESDEvent* esd, 
-                                      Double_t& cent, 
-                                      UShort_t& qual) const
-{
-  // 
-  // Read centrality from event 
-  // 
-  // Parameters:
-  //    esd  Event 
-  //    cent On return, the centrality or negative if not found
-  // 
-  // Return:
-  //    False on error, true otherwise 
-  //
-  cent = -1;
-  qual = 0;
-  AliCentrality* centObj = const_cast<AliESDEvent*>(esd)->GetCentrality();
-  if (!centObj) return true;
-
-  qual = centObj->GetQuality();
-  if (qual == 0x8) // Ignore ZDC outliers 
-    cent = centObj->GetCentralityPercentileUnchecked("V0M");  
-  else
-    cent = centObj->GetCentralityPercentile("V0M");        
-
-  return true;
-}
 
 //____________________________________________________________________
 Bool_t