Various upgrades. NSD true trigger for MC, pileup selection for analysis, fixes for...
authorhdalsgaa <hdalsgaa@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 28 Feb 2011 23:56:46 +0000 (23:56 +0000)
committerhdalsgaa <hdalsgaa@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 28 Feb 2011 23:56:46 +0000 (23:56 +0000)
17 files changed:
PWG2/FORWARD/analysis2/AliAODForwardMult.cxx
PWG2/FORWARD/analysis2/AliAODForwardMult.h
PWG2/FORWARD/analysis2/AliBasedNdetaTask.cxx
PWG2/FORWARD/analysis2/AliBasedNdetaTask.h
PWG2/FORWARD/analysis2/AliCentralMultiplicityTask.cxx
PWG2/FORWARD/analysis2/AliCentralMultiplicityTask.h
PWG2/FORWARD/analysis2/AliFMDEnergyFitterTask.cxx
PWG2/FORWARD/analysis2/AliFMDEnergyFitterTask.h
PWG2/FORWARD/analysis2/AliFMDEventInspector.cxx
PWG2/FORWARD/analysis2/AliFMDEventInspector.h
PWG2/FORWARD/analysis2/AliFMDMCCorrector.cxx
PWG2/FORWARD/analysis2/AliFMDMCSharingFilter.cxx
PWG2/FORWARD/analysis2/AliForwardMCMultiplicityTask.cxx
PWG2/FORWARD/analysis2/AliForwardMultiplicityTask.cxx
PWG2/FORWARD/analysis2/AliForwardUtil.cxx
PWG2/FORWARD/analysis2/AliForwarddNdetaTask.cxx
PWG2/FORWARD/analysis2/scripts/RunCopyCentralSecMap.C

index 0b87b4e..5660b8e 100644 (file)
@@ -33,7 +33,8 @@ AliAODForwardMult::AliAODForwardMult()
   : fIsMC(false),
     fHist(),
     fTriggers(0),
-    fIpZ(fgkInvalidIpZ)
+    fIpZ(fgkInvalidIpZ), 
+    fCentrality(-1)
 {
   // 
   // Constructor 
@@ -46,7 +47,8 @@ AliAODForwardMult::AliAODForwardMult(Bool_t isMC)
     fHist("forwardMult", "d^{2}N_{ch}/d#etad#varphi in the forward regions", 
          200, -4, 6, 20, 0, 2*TMath::Pi()),
     fTriggers(0),
-    fIpZ(fgkInvalidIpZ)
+    fIpZ(fgkInvalidIpZ), 
+    fCentrality(-1)
 {
   // 
   // Constructor 
@@ -156,11 +158,14 @@ AliAODForwardMult::Browse(TBrowser* b)
   //   b   Browser to use 
   static TObjString ipz;
   static TObjString trg;
+  static TObjString cnt;
   ipz = Form("ip_z=%fcm", fIpZ);
   trg = GetTriggerString(fTriggers);
+  cnt = Form("%+6.1f%%", fCentrality);
   b->Add(&fHist);
   b->Add(&ipz);
   b->Add(&trg);
+  b->Add(&cnt);
 }
 
 //____________________________________________________________________
@@ -175,13 +180,14 @@ AliAODForwardMult::GetTriggerString(UInt_t mask)
   //   Character string representation of mask 
   static TString trg;
   trg = "";
-  if ((mask & kInel)    != 0x0) trg.Append("INEL ");
-  if ((mask & kInelGt0) != 0x0) trg.Append("INEL>0 ");
-  if ((mask & kNSD)     != 0x0) trg.Append("NSD ");
-  if ((mask & kA)       != 0x0) trg.Append("A ");
-  if ((mask & kB)       != 0x0) trg.Append("B ");
-  if ((mask & kC)       != 0x0) trg.Append("C ");
-  if ((mask & kE)       != 0x0) trg.Append("E ");
+  if ((mask & kInel)        != 0x0) trg.Append("INEL ");
+  if ((mask & kInelGt0)     != 0x0) trg.Append("INEL>0 ");
+  if ((mask & kNSD)         != 0x0) trg.Append("NSD ");
+  if ((mask & kA)           != 0x0) trg.Append("A ");
+  if ((mask & kB)           != 0x0) trg.Append("B ");
+  if ((mask & kC)           != 0x0) trg.Append("C ");
+  if ((mask & kE)           != 0x0) trg.Append("E ");
+  if ((mask & kMCNSD)       != 0x0) trg.Append("MCNSD ");
   return trg.Data();
 }
   
@@ -200,11 +206,12 @@ AliAODForwardMult::Print(Option_t* option) const
   case 1:  str = "pp"; break;
   case 2:  str = "PbPb"; break;
   }
-  std::cout << "Ipz:      " << fIpZ << "cm " << (HasIpZ() ? "" : "in") 
+  std::cout << "Ipz:         " << fIpZ << "cm " << (HasIpZ() ? "" : "in") 
            << "valid\n"
-           << "Triggers: " << GetTriggerString(fTriggers) 
-           << "sNN:      " << GetSNN() << "GeV\n" 
-           << "System:   " << str 
+           << "Triggers:    " << GetTriggerString(fTriggers) 
+           << "sNN:         " << GetSNN() << "GeV\n" 
+           << "System:      " << str 
+           << "Centrality:  " << fCentrality << "%" 
            << std::endl;
 }
 
index e260623..2747dfc 100644 (file)
@@ -109,7 +109,12 @@ public:
     /** C-side trigger */
     kC        = 0x080,  
     /** Empty trigger */
-    kE        = 0x100
+    kE        = 0x100,
+    /** pileup from SPD */
+    kPileUp   = 0x200,    
+    /** true NSD from MC */
+    kMCNSD    = 0x400    
+    
   };
   /** 
    * Default constructor 
@@ -222,6 +227,12 @@ public:
    */
   void SetSystem(UShort_t sys);
   /** 
+   * Set the event centrality 
+   * 
+   * @param c Centrality 
+   */
+  void SetCentrality(Float_t c) { fCentrality = c; }
+  /** 
    * Set the z coordinate of the interaction point
    * 
    * @return Interaction point z coordinate
@@ -260,6 +271,21 @@ public:
    */
   Bool_t InRange(Float_t low, Float_t high) const;
   /** 
+   * Get the event centrality 
+   * 
+   * 
+   * @return 
+   */
+  Float_t GetCentrality() const { return fCentrality; }
+  /** 
+   * Check if we have a valid centrality 
+   * 
+   * 
+   * @return 
+   */
+  Bool_t  HasCentrality() const { return !(fCentrality  < 0); }
+  
+  /** 
    * Get the name of the object 
    * 
    * @return Name of object 
@@ -278,9 +304,10 @@ protected:
   TH2D    fHist;     // Histogram of d^2N_{ch}/(deta dphi) for this event
   UInt_t  fTriggers; // Trigger bit mask 
   Float_t fIpZ;      // Z coordinate of the interaction point
+  Float_t fCentrality; // Event centrality 
 
   static const Float_t fgkInvalidIpZ; // Invalid IpZ value 
-  ClassDef(AliAODForwardMult,1); // AOD forward multiplicity 
+  ClassDef(AliAODForwardMult,2); // AOD forward multiplicity 
 };
 
 //____________________________________________________________________
index f58d810..1091cc2 100644 (file)
@@ -5,6 +5,7 @@
 #include <TH1D.h>
 #include <THStack.h>
 #include <TList.h>
+#include <TParameter.h>
 #include <AliAnalysisManager.h>
 #include <AliAODEvent.h>
 #include <AliAODHandler.h>
@@ -26,7 +27,9 @@ AliBasedNdetaTask::AliBasedNdetaTask()
     fRebin(0),         // Rebinning factor 
     fCutEdges(false), 
     fSymmetrice(true),
-    fCorrEmpty(true)
+    fCorrEmpty(true), 
+    fTriggerEff(1),
+    fShapeCorr(0)
 {}
 
 //____________________________________________________________________
@@ -43,7 +46,9 @@ AliBasedNdetaTask::AliBasedNdetaTask(const char* name)
     fRebin(5),         // Rebinning factor 
     fCutEdges(false), 
     fSymmetrice(true),
-    fCorrEmpty(true)
+    fCorrEmpty(true), 
+    fTriggerEff(1),
+    fShapeCorr(0)
 {
   // Output slot #1 writes into a TH1 container
   DefineOutput(1, TList::Class()); 
@@ -64,7 +69,9 @@ AliBasedNdetaTask::AliBasedNdetaTask(const AliBasedNdetaTask& o)
     fRebin(o.fRebin),          // Int_t - Rebinning factor 
     fCutEdges(o.fCutEdges),    // Bool_t - Whether to cut edges when rebinning
     fSymmetrice(true),
-    fCorrEmpty(true)
+    fCorrEmpty(true), 
+    fTriggerEff(o.fTriggerEff),
+    fShapeCorr(o.fShapeCorr)
 {}
 
 //____________________________________________________________________
@@ -106,6 +113,15 @@ AliBasedNdetaTask::SetTriggerMask(const char* mask)
 
 //________________________________________________________________________
 void 
+AliBasedNdetaTask::SetShapeCorrection(const TH1* c)
+{
+  if (!c) return;
+  fShapeCorr = static_cast<TH1*>(c->Clone());
+  fShapeCorr->SetDirectory(0);
+}
+
+//________________________________________________________________________
+void 
 AliBasedNdetaTask::UserCreateOutputObjects()
 {
   // Create histograms
@@ -184,11 +200,15 @@ AliBasedNdetaTask::CheckEvent(AliAODEvent* aod)
     fTriggers->AddBinContent(kE);
   if (forward->IsTriggerBits(AliAODForwardMult::kInel)) 
     fTriggers->AddBinContent(kMB);
-
+  if (forward->IsTriggerBits(AliAODForwardMult::kMCNSD)) 
+    fTriggers->AddBinContent(kMCNSD);
+  
+  
   // Check if we have an event of interest. 
   if (!forward->IsTriggerBits(fTriggerMask)) return false;
+  
   fTriggers->AddBinContent(kWithTrigger);
-
+  
   // Check that we have a valid vertex
   if (!forward->HasIpZ()) return false;
   fTriggers->AddBinContent(kWithVertex);
@@ -209,9 +229,20 @@ AliBasedNdetaTask::UserExec(Option_t *)
     AliError("Cannot get the AOD event");
     return;
   }  
-
+  
+  TObject* obj = 0;
+  obj = aod->FindListObject("Forward");
+  
+  // AliAODForwardMult* forward = static_cast<AliAODForwardMult*>(obj);
+  // Float_t cent = forward->GetCentrality();
+  
+  // if(cent > -1) {
+  //  if( cent < 40 || cent >60 ) return;
+  //  else std::cout<<"selecting event with cent "<<cent<<std::endl;
+  //  }
+  
   // Get the histogram(s) 
-  TH2D* data   = GetHistogram(aod);
+  TH2D* data   = GetHistogram(aod, false);
   TH2D* dataMC = GetHistogram(aod, true);
 
   // We should have a base object at least 
@@ -220,17 +251,23 @@ AliBasedNdetaTask::UserExec(Option_t *)
     return;
   }
   
-  // Create our sum histograms 
-  if (!fSum)             fSum   = CloneHist(data,  GetName());
-  if (!fSumMC && dataMC) fSumMC = CloneHist(dataMC,Form("%sMC", GetName()));
+  
 
   // Check the event 
   if (!CheckEvent(aod)) return;
-
+  
+  // Create our sum histograms 
+  if (!fSum)  fSum   = CloneHist(data,  GetName()); 
+  else fSum->Add((data));
+  
+  //for(Int_t i = 1; i<=data->GetNbinsX(); i++) {
+  //  std::cout<<fSum->GetBinContent(i,0) <<"   "<<data->GetBinContent(i,0)<<std::endl;
+  // }
+  
+  if (!fSumMC && dataMC) fSumMC = CloneHist(dataMC,Form("%sMC", GetName()));
+  else if (fSumMC) fSumMC->Add((dataMC));
   // Add contribution 
-  fSum->Add((data));
-  if (fSumMC) fSumMC->Add((dataMC));
-
+  
   PostData(1, fSums);
 }
 
@@ -265,9 +302,9 @@ AliBasedNdetaTask::ProjectX(const TH2D* h,
                            bool  error) const
 {
   if (!h) return 0;
-#if USE_ROOT_PROJECT
+  //#if USE_ROOT_PROJECT
   return h->ProjectionX(name, firstbin, lastbin, (error ? "e" : ""));
-#endif
+  //#endif
   
   TAxis* xaxis = h->GetXaxis();
   TAxis* yaxis = h->GetYaxis();
@@ -371,6 +408,32 @@ AliBasedNdetaTask::Terminate(Option_t *)
   Int_t nWithVertex = Int_t(fTriggers->GetBinContent(kWithVertex));
   Int_t nAccepted   = Int_t(fTriggers->GetBinContent(kAccepted));
   Int_t nGood       = nB - nA - nC + 2 * nE;
+  
+  //  Int_t nBg         = nA + nC -nE;
+  //std::cout<<"nBackground "<<nBg<<"   , with vertex in range "<<fNbgValid<<" ,  total "<<fNbgVertex<<"  "<<std::endl<<std::endl;
+  Float_t alpha            = ((Float_t)nAccepted) / ((Float_t)nWithVertex);
+  //Float_t alphaBG  = 1;
+  //if(fNbgVertex > 0) alphaBG= (Float_t)fNbgValid / (Float_t)fNbgVertex;
+  //Float_t nNormalizationJF = nAccepted + alpha*(nMB - nAccepted - nBg);
+  Float_t nNormalizationJF = nAccepted + alpha*(nTriggered - nWithVertex) ; // -alphaBG*(nBg-fNbgVertex);
+  
+  //Float_t nNormalizationBg = fNbgValid + alphaBG*nBg;
+  
+  //std::cout<<nTriggered -nWithVertex<<std::endl;
+  
+  //Double_t vtxEff   = Double_t(nMB) / nTriggered * Double_t(nAccepted) / nGood;
+  Double_t vtxEff   =  (Float_t)nAccepted / nNormalizationJF;
+  vtxEff = vtxEff*fTriggerEff;
+  
+  //if ( fTriggerMask == AliAODForwardMult::kNSD ) vtxEff = 0.97; //From paper...
+  //std::cout<<
+
+  //Double_t vtxEffBg   =  (Float_t)fNbgValid / nNormalizationBg;
+  //Double_t vtxEff   =  1.;//(Float_t)nAccepted / nNormalizationJF;
+  
+  
+  
+  
   if (nTriggered <= 0) { 
     AliError("Number of triggered events <= 0");
     return;
@@ -380,7 +443,8 @@ AliBasedNdetaTask::Terminate(Option_t *)
                    nGood, nB, nA, nC, nE));
     nGood = nMB;
   }
-  Double_t vtxEff   = Double_t(nMB) / nTriggered * Double_t(nAccepted) / nGood;
+  //Double_t vtxEff   =  Double_t(nMB) / nTriggered * Double_t(nAccepted) / nGood;
+  
   Double_t vNorm    = Double_t(nAccepted) / nGood;
   AliInfo(Form("Total of %9d events\n"
               "                   of these %9d are minimum bias\n"
@@ -400,18 +464,56 @@ AliBasedNdetaTask::Terminate(Option_t *)
               nB, nA+nC, nA, nC, nE, nGood, vtxEff, vNorm));
   
   const char* name = GetName();
-
+  
   // Get acceptance normalisation from underflow bins 
   TH1D* norm   = ProjectX(fSum, Form("norm%s",name), 0, 1, fCorrEmpty, false);
+  
+    
+  
+
+  /*
+  std::cout<<norm->GetMaximumBin()<<"    "<< (Float_t)nAccepted / norm->GetBinContent((Float_t)norm->GetMaximumBin()) <<std::endl;
+  Float_t scaleForNorm =  (Float_t)nAccepted / (Float_t)norm->GetBinContent(norm->GetMaximumBin()) ;
+  //Float_t scaleForNorm =  norm->Integral() ;
+  std::cout<<norm->Integral()<<std::endl;
+  
+  norm->Scale(scaleForNorm);
+  //norm->Scale(1., "width");
+  //norm->Scale(norm->GetNbinsX() / (norm->GetXaxis()->GetXmax() - norm->GetXaxis()->GetXmin() ));
+  
+  //norm->Scale(1.,"width");
+  */
   // Project onto eta axis - _ignoring_underflow_bins_!
+  
+  TH2D* tmpNorm = (TH2D*)fSum->Clone("tmpNorm");
+  
+  if(fShapeCorr)
+    fSum->Divide(fShapeCorr);
   TH1D* dndeta = ProjectX(fSum, Form("dndeta%s",name), 1, fSum->GetNbinsY(),
                          fCorrEmpty);
+  
   // Normalize to the acceptance 
-  dndeta->Divide(norm);
+  //dndeta->Divide(norm);
+  
+  for(Int_t i = 1; i<=tmpNorm->GetNbinsX(); i++) {
+    
+    if( tmpNorm->GetBinContent(i,0) > 0 ) {
+      //std::cout<<tmpNorm->GetBinContent(i,0) <<std::endl;
+      dndeta->SetBinContent(i,dndeta->GetBinContent(i) / tmpNorm->GetBinContent(i,0));
+      dndeta->SetBinError(i,dndeta->GetBinError(i) / tmpNorm->GetBinContent(i,0));
+    }
+  }
+  
   // Scale by the vertex efficiency 
-  dndeta->Scale(vNorm, "width");
-  norm->Scale(1. / nAccepted);
-
+  //dndeta->Scale(vNorm, "width");
+    
+  dndeta->Scale(vtxEff, "width");
+  
+  //dndeta->Scale(1./(Float_t)nAccepted);
+  //dndeta->Scale(dndeta->GetNbinsX() / (dndeta->GetXaxis()->GetXmax() - dndeta->GetXaxis()->GetXmin() ));
+  
+  //norm->Scale(1. / nAccepted);
+  //norm->Scale(1. / nNormalizationJF);
   SetHistogramAttributes(dndeta, kRed+1, 20, Form("ALICE %s", name));
   SetHistogramAttributes(norm, kRed+1,20,Form("ALICE %s normalisation", name)); 
 
@@ -455,7 +557,9 @@ AliBasedNdetaTask::Terminate(Option_t *)
   vtxAxis->SetName("vtxAxis");
   vtxAxis->SetTitle(Form("v_{z}#in[%+5.1f,%+5.1f]cm", fVtxMin,fVtxMax));
   fOutput->Add(vtxAxis);
-
+  fOutput->Add(new TParameter<Double_t>("triggerEff", fTriggerEff)); 
+  if (fShapeCorr) fOutput->Add(fShapeCorr);
+  
   PostData(2, fOutput);
 }
 
index 766ce9e..c6ff44b 100644 (file)
@@ -54,6 +54,19 @@ public:
    * @param mask trigger mask 
    */
   void SetTriggerMask(const char* mask);
+  /** 
+   * Trigger efficiency for selected trigger(s)
+   * 
+   * @param e Trigger efficiency 
+   */
+  void SetTriggerEff(Double_t e) { fTriggerEff = e; } 
+  /** 
+   * Set the shape correction (a.k.a., track correction) for selected
+   * trigger(s)
+   * 
+   * @param h Correction
+   */
+  void SetShapeCorrection(const TH1* h);
   /**
    * Destructor
    * 
@@ -171,7 +184,8 @@ protected:
     kMB         = 6, 
     kWithTrigger= 7,
     kWithVertex = 8, 
-    kAccepted   = 9 
+    kAccepted   = 9,
+    kMCNSD      = 10
   };
 
   TH2D*           fSum;          // Sum of histograms 
@@ -189,6 +203,8 @@ protected:
   Bool_t          fCutEdges;     // Whether to cut edges when rebinning
   Bool_t          fSymmetrice;   // Whether to symmetrice data 
   Bool_t          fCorrEmpty;    // Correct for empty bins 
+  Double_t        fTriggerEff;   // Trigger efficiency for selected trigger(s)
+  TH1*            fShapeCorr;    // Shape correction 
 
   ClassDef(AliBasedNdetaTask,1); // Determine multiplicity in base area
 };
index af5a0b5..aa98bf8 100644 (file)
@@ -22,6 +22,7 @@
 #include "AliAnalysisManager.h"
 #include "AliESDEvent.h"
 #include "AliMultiplicity.h"
+#include "AliFMDEventInspector.h"
 #include <TROOT.h>
 #include <TFile.h>
 #include <TError.h>
@@ -35,7 +36,8 @@ AliCentralMultiplicityTask::AliCentralMultiplicityTask(const char* name)
     fList(0),
     fAODCentral(kFALSE),
     fManager(),
-    fUseSecondary(true)
+    fUseSecondary(true),
+    firstEventSeen(false)
 {
   
   DefineOutput(1, TList::Class());
@@ -47,7 +49,8 @@ AliCentralMultiplicityTask::AliCentralMultiplicityTask()
     fList(0),
     fAODCentral(),
     fManager(),
-    fUseSecondary(true)
+    fUseSecondary(true),
+    firstEventSeen(false)
 {
 }
 //____________________________________________________________________
@@ -57,18 +60,20 @@ AliCentralMultiplicityTask::AliCentralMultiplicityTask(const AliCentralMultiplic
     fList(o.fList),
     fAODCentral(o.fAODCentral),
     fManager(o.fManager),
-    fUseSecondary(o.fUseSecondary)
+    fUseSecondary(o.fUseSecondary),
+    firstEventSeen(o.firstEventSeen)
 {
 }
 //____________________________________________________________________
 AliCentralMultiplicityTask&
 AliCentralMultiplicityTask::operator=(const AliCentralMultiplicityTask& o)
 {
-  fData         = o.fData;
-  fList         = o.fList;
-  fAODCentral   = o.fAODCentral;
-  fManager      = o.fManager;
-  fUseSecondary = o.fUseSecondary;
+  fData          = o.fData;
+  fList          = o.fList;
+  fAODCentral    = o.fAODCentral;
+  fManager       = o.fManager;
+  fUseSecondary  = o.fUseSecondary;
+  firstEventSeen = o.firstEventSeen;
   return *this;
 }
 //____________________________________________________________________
@@ -101,6 +106,18 @@ void AliCentralMultiplicityTask::UserExec(Option_t* /*option*/)
   
   AliESDEvent* esd = eventHandler->GetEvent();
   
+  if(!GetManager().IsInit() && !firstEventSeen) {
+    AliFMDEventInspector inspector;
+    inspector.ReadRunDetails(esd);
+    GetManager().Init(inspector.GetCollisionSystem(),
+                     inspector.GetEnergy(),
+                     inspector.GetField());
+    
+    //std::cout<<inspector.GetCollisionSystem()<<"  "<<inspector.GetEnergy()<<"    "<<inspector.GetField()<<std::endl;
+    AliInfo("Manager of corrections in AliCentralMultiplicityTask init");
+    firstEventSeen = kTRUE;
+  }
+  
   //Selecting only events with |valid vertex| < 10 cm
   const AliESDVertex* vertex = esd->GetPrimaryVertexSPD();
   if (!vertex) return;
@@ -190,7 +207,8 @@ AliCentralMultiplicityTask::Manager::Manager() :
   fAcceptance(),
   fSecmap(),
   fAcceptanceName("centralacceptance"),
-  fSecMapName("centralsecmap")
+  fSecMapName("centralsecmap"),
+  fIsInit(kFALSE)
 {
 
 
@@ -202,7 +220,8 @@ AliCentralMultiplicityTask::Manager::Manager(const Manager& o)
    fAcceptance(o.fAcceptance),
    fSecmap(o.fSecmap),
    fAcceptanceName(o.fAcceptanceName),
-   fSecMapName(o.fSecMapName) 
+   fSecMapName(o.fSecMapName),
+   fIsInit(o.fIsInit)
 {}
 //____________________________________________________________________
 AliCentralMultiplicityTask::Manager&
@@ -214,6 +233,7 @@ AliCentralMultiplicityTask::Manager::operator=(const Manager& o)
   fSecmap         = o.fSecmap;
   fAcceptanceName = o.fAcceptanceName;
   fSecMapName     = o.fSecMapName;
+  fIsInit         = o.fIsInit;
   return *this;
 }
 
@@ -284,6 +304,8 @@ AliCentralMultiplicityTask::Manager::Init(UShort_t  sys,
                                          UShort_t  sNN,
                                          Short_t   field) 
 {
+  if(fIsInit) ::Warning("Init","Already initialised - overriding...");
+  
   TFile fsec(GetFullFileName(0,sys,sNN,field));
   fSecmap = 
     dynamic_cast<AliCentralCorrSecondaryMap*>(fsec.Get(fSecMapName.Data()));  
@@ -298,7 +320,14 @@ AliCentralMultiplicityTask::Manager::Init(UShort_t  sys,
     ::Error("Init", "no central Acceptance found!") ;
     return;
   }
+  
+  if(fSecmap && fAcceptance) {
+    fIsInit = kTRUE;
+    ::Info("Init", 
+          "Central Manager initialised for sys %d, energy %d, field %d",sys,sNN,field);
 
+  }
+  
 }
 //
 // EOF
index 2e01b74..0019d2d 100644 (file)
@@ -157,6 +157,14 @@ public:
      * @param field  Magnetic field [kG]
      */
     void Init(UShort_t sys, UShort_t sNN, Short_t field);
+    
+    /** 
+     * Is initialized 
+     * 
+     */
+    Bool_t IsInit() { return fIsInit; }
+    
+    
     /** 
      * Get the acceptance path
      * 
@@ -243,6 +251,7 @@ public:
     AliCentralCorrSecondaryMap* fSecmap;         // Secindary map
     TString                     fAcceptanceName; // Acceptance name
     TString                     fSecMapName;     // Secindary name
+    Bool_t                      fIsInit;         // Are we init
 
     ClassDef(Manager,1); // Manager of data 
   };
@@ -267,7 +276,8 @@ protected:
   TList*                 fList;          //Output List for diagnostics
   AliAODCentralMult      fAODCentral;    // Output object
   Manager                fManager;       //Manager object for corrections
-  Bool_t                 fUseSecondary;   // Whether to secondary map
+  Bool_t                 fUseSecondary;  // Whether to secondary map
+  Bool_t                 firstEventSeen; // Have we seen first event     
   ClassDef(AliCentralMultiplicityTask,1) // Forward multiplicity class
 };
 
index 62f589d..f09dc28 100644 (file)
 #include <TDirectory.h>
 #include <TTree.h>
 #include <TFile.h>
+#include "AliMCEvent.h"
+#include "AliGenHijingEventHeader.h"
+#include "AliHeader.h"
+#include <iostream>
 
 //====================================================================
 AliFMDEnergyFitterTask::AliFMDEnergyFitterTask()
@@ -46,7 +50,9 @@ AliFMDEnergyFitterTask::AliFMDEnergyFitterTask(const char* name)
     fFirstEvent(true),
     fEventInspector("event"),
     fEnergyFitter("energy"),
-    fList(0)
+    fList(0),
+    fbLow(0),
+    fbHigh(100)
 {
   // 
   // Constructor 
@@ -171,6 +177,21 @@ AliFMDEnergyFitterTask::UserExec(Option_t*)
   // static Int_t cnt = 0;
   // cnt++;
   // Get the input data 
+  
+  AliMCEvent* mcevent = MCEvent();
+  if(mcevent) {
+    AliHeader* header            = mcevent->Header();
+    AliGenEventHeader* genHeader = header->GenEventHeader();
+    AliGenHijingEventHeader* hijingHeader = dynamic_cast<AliGenHijingEventHeader*>(header->GenEventHeader());
+    if(hijingHeader) {
+      Float_t b = hijingHeader->ImpactParameter();
+      if(b<fbLow || b>fbHigh) return;
+      else
+       std::cout<<"Selecting event with impact parameter "<<b<<std::endl;
+    }
+    
+  }
+  
   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
   // AliInfo(Form("Event # %6d (esd=%p)", cnt, esd));
   if (!esd) { 
index a5a193a..0324582 100644 (file)
@@ -119,6 +119,10 @@ public:
    * @param dbg Debug level
    */
   void SetDebug(Int_t dbg);
+
+  void SetBLow(Float_t b) {fbLow = b;}
+  void SetBHigh(Float_t b) {fbHigh = b;}
+
 protected: 
   /** 
    * Initialise the sub objects and stuff.  Called on first event 
@@ -130,7 +134,9 @@ protected:
   AliFMDEventInspector fEventInspector; // Algorithm
   AliFMDEnergyFitter   fEnergyFitter;   // Algorithm
   TList*               fList;           // Output list 
-
+  Float_t              fbLow;
+  Float_t              fbHigh;
+  
   ClassDef(AliFMDEnergyFitterTask,1) // Forward multiplicity class
 };
 
index 4a8d3c4..a83788a 100644 (file)
 #include <TROOT.h>
 #include <iostream>
 #include <iomanip>
-
+#include "AliMCEvent.h"
+#include "AliGenPythiaEventHeader.h"
+#include "AliGenDPMjetEventHeader.h"
+#include "AliHeader.h"
+#include "AliMCEventHandler.h"
 //====================================================================
 AliFMDEventInspector::AliFMDEventInspector()
   : TNamed(),
@@ -238,8 +242,8 @@ AliFMDEventInspector::Init(const TAxis& vtxAxis)
   fHTriggers->GetXaxis()->SetBinLabel(kB      +1,"B");
   fHTriggers->GetXaxis()->SetBinLabel(kC      +1,"C");
   fHTriggers->GetXaxis()->SetBinLabel(kE      +1,"E");
-  fHTriggers->GetXaxis()->SetBinLabel(9,         "spare1");
-  fHTriggers->GetXaxis()->SetBinLabel(10,        "spare2");
+  fHTriggers->GetXaxis()->SetBinLabel(kPileUp +1,"Pileup");
+  fHTriggers->GetXaxis()->SetBinLabel(kMCNSD  +1,"nsd");
   fList->Add(fHTriggers);
 
   fHType = new TH1I("type", Form("Event type (cut: SPD mult>%d)", 
@@ -337,10 +341,14 @@ AliFMDEventInspector::Process(const AliESDEvent* event,
 
   fHType->Fill(lowFlux ? 0 : 1);
   
-  cent = -1;
+  cent = -10;
   AliCentrality* centObj = const_cast<AliESDEvent*>(event)->GetCentrality();
-  if (centObj) 
-    cent = centObj->GetCentralityPercentile("VOM");  
+  if (centObj) {
+    AliInfo(Form("Got centrality object %p with quality %d", centObj, centObj->GetQuality()));
+    centObj->Print();
+    cent = centObj->GetCentralityPercentileUnchecked("V0M");  
+  }
+  AliInfo(Form("Centrality is %f", cent));
   fHCent->Fill(cent);
 
   // Check the FMD ESD data 
@@ -372,6 +380,8 @@ AliFMDEventInspector::Process(const AliESDEvent* event,
     ivz = 0;
     return kBadVertex;
   }
+  
+  
   return kOk;
 }
 
@@ -405,7 +415,38 @@ AliFMDEventInspector::ReadTriggers(const AliESDEvent* esd, UInt_t& triggers)
     AliWarning("No input handler");
     return kFALSE;
   }
-
+  
+  AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
+  AliMCEvent* mcEvent = 0;
+  if(mcHandler)
+    mcEvent = mcHandler->MCEvent();
+  
+  if(mcEvent) {
+    //Assign MC only triggers : True NSD etc.
+    AliHeader* header            = mcEvent->Header();
+    AliGenEventHeader* genHeader = header->GenEventHeader();
+    
+    AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
+    AliGenDPMjetEventHeader* dpmHeader = dynamic_cast<AliGenDPMjetEventHeader*>(header->GenEventHeader());
+    Bool_t sd = kFALSE;
+    if(pythiaGenHeader) {
+      Int_t pythiaType = pythiaGenHeader->ProcessType();
+      if(pythiaType==92 || pythiaType==93)
+       sd = kTRUE;
+    }
+    if(dpmHeader) {
+      Int_t processType = dpmHeader->ProcessType();
+      if(processType == 5 || processType == 6)  
+       sd = kTRUE;
+      
+    }
+    if(!sd) {
+      triggers |= AliAODForwardMult::kMCNSD;
+      fHTriggers->Fill(kMCNSD+0.5);
+    }
+    
+  }
+    
   // 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 
@@ -445,7 +486,13 @@ AliFMDEventInspector::ReadTriggers(const AliESDEvent* esd, UInt_t& triggers)
     triggers |= AliAODForwardMult::kNSD;
     fHTriggers->Fill(kNSD+.5);
   }
-
+  //Check pileup
+  Bool_t pileup =  esd->IsPileupFromSPD(3,0.8);
+  if (pileup) {
+    triggers |= AliAODForwardMult::kPileUp;
+    fHTriggers->Fill(kPileUp+.5);
+  }
+    
   // Get trigger stuff 
   TString trigStr = esd->GetFiredTriggerClasses();
   // AliWarning(Form("Fired trigger classes: %s", trigStr.Data()));
@@ -526,6 +573,16 @@ AliFMDEventInspector::ReadVertex(const AliESDEvent* esd, Double_t& vz)
     return kFALSE;
   }
 
+  // HHD test
+  /*
+  if (vertex->IsFromVertexerZ()) { 
+    return kFALSE;
+  }
+  if (TMath::Sqrt(TMath::Power(vertex->GetX(),2) + TMath::Power(vertex->GetY(),2)) > 3 ) { 
+      return kFALSE;
+  }
+  */
+  
   // Get the z coordiante 
   vz = vertex->GetZ();
   return kTRUE;
index d4061a8..4c4b559 100644 (file)
@@ -50,7 +50,11 @@ public:
     /** No vertex found */
     kNoVertex = 0x10, 
     /** Vertex out of range */
-    kBadVertex = 0x20
+    kBadVertex = 0x20,
+    /** Suspected pileup */
+    kPileUp = 0x40,
+    /** Suspected pileup */
+    kMCNSD = 0x80
   };
   /** 
    * Trigger bins 
index c59b2b1..a4520e0 100644 (file)
@@ -25,7 +25,7 @@
 #include <TH2D.h>
 #include <TROOT.h>
 #include <TProfile2D.h>
-
+#include <iostream>
 ClassImp(AliFMDMCCorrector)
 #if 0
 ; // For Emacs
@@ -101,7 +101,7 @@ AliFMDMCCorrector::CorrectMC(AliForwardUtil::Histos& hists,
       }
       if (fUseVertexBias) {
         TH2D*       ef = fcm.GetVertexBias()->GetCorrection(r, uvb);
-        if (!ef) {
+       if (!ef) {
           AliWarning(Form("No event vertex bias correction in vertex bin %d",
                        uvb));
           continue;
index bcd00c5..ddfc934 100644 (file)
@@ -196,7 +196,8 @@ AliFMDMCSharingFilter::FilterMC(const AliESDFMD&  input,
     }
   }
   AliStack* stack = const_cast<AliMCEvent&>(event).Stack();
-  Int_t nTracks   = event.GetNumberOfTracks();
+  Int_t nTracks   = stack->GetNtrack();//event.GetNumberOfTracks();
+  Int_t nPrim     = stack->GetNprimary();//event.GetNumberOfPrimary();
   for (Int_t iTr = 0; iTr < nTracks; iTr++) { 
     AliMCParticle* particle = 
       static_cast<AliMCParticle*>(event.GetTrack(iTr));
@@ -210,7 +211,7 @@ AliFMDMCSharingFilter::FilterMC(const AliESDFMD&  input,
     Bool_t isPrimary = stack->IsPhysicalPrimary(iTr);
 
     // Fill 'dn/deta' histogram 
-    if (isPrimary) { 
+    if (isPrimary && iTr < nPrim) { 
       fSumEta->Fill(particle->Eta());
       primary->Fill(particle->Eta(), particle->Phi());
     }
index d7a9a63..ec948aa 100644 (file)
@@ -272,7 +272,8 @@ AliForwardMCMultiplicityTask::UserExec(Option_t*)
   // Set trigger bits, and mark this event for storage 
   fAODFMD.SetTriggerBits(triggers);
   fMCAODFMD.SetTriggerBits(triggers);
-
+  fAODFMD.SetCentrality(cent);
+  
   //All events should be stored - HHD
   //if (isAccepted) MarkEventForStore();
 
index b6627e9..1c7df91 100644 (file)
@@ -219,7 +219,7 @@ AliForwardMultiplicityTask::UserExec(Option_t*)
   fHistos.Clear();
   fESDFMD.Clear();
   fAODFMD.Clear();
-
+  
   Bool_t   lowFlux  = kFALSE;
   UInt_t   triggers = 0;
   UShort_t ivz      = 0;
@@ -227,7 +227,7 @@ AliForwardMultiplicityTask::UserExec(Option_t*)
   Double_t cent     = -1;
   UInt_t   found    = fEventInspector.Process(esd, triggers, lowFlux, 
                                              ivz, vz, cent);
+  
   if (found & AliFMDEventInspector::kNoEvent)    return;
   if (found & AliFMDEventInspector::kNoTriggers) return;
 
@@ -235,11 +235,15 @@ AliForwardMultiplicityTask::UserExec(Option_t*)
   fAODFMD.SetTriggerBits(triggers);
   fAODFMD.SetSNN(fEventInspector.GetEnergy());
   fAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
+  fAODFMD.SetCentrality(cent);
   MarkEventForStore();
-
+  
   if (found & AliFMDEventInspector::kNoSPD)      return;
   if (found & AliFMDEventInspector::kNoFMD)      return;
   if (found & AliFMDEventInspector::kNoVertex)   return;
+  
+  if (triggers & AliAODForwardMult::kPileUp) return;
+  
   fAODFMD.SetIpZ(vz);
 
   if (found & AliFMDEventInspector::kBadVertex) return;
index b5dd3e5..e1fe70d 100644 (file)
@@ -84,7 +84,7 @@ AliForwardUtil::ParseCenterOfMassEnergy(UShort_t /* sys */, Float_t v)
   // if (sys == AliForwardUtil::kPbPb) energy = energy / 208 * 82;
   if (TMath::Abs(energy - 900.)   < 10)  return 900;
   if (TMath::Abs(energy - 2400.)  < 10)  return 2400;
-  if (TMath::Abs(energy - 2750.)  < 10)  return 2750;
+  if (TMath::Abs(energy - 2750.)  < 20)  return 2750;
   if (TMath::Abs(energy - 5500.)  < 40)  return 5500;
   if (TMath::Abs(energy - 7000.)  < 10)  return 7000;
   if (TMath::Abs(energy - 10000.) < 10)  return 10000;
index 9507db1..7b842a7 100644 (file)
@@ -58,11 +58,29 @@ AliForwarddNdetaTask::GetHistogram(AliAODEvent* aod, Bool_t mc)
     if (oPrimary) {
       TH2D* primary   = static_cast<TH2D*>(oPrimary);
       // Add contribtion from MC 
-      if (primary && !fSumPrimary) fSumPrimary = CloneHist(primary, "truth");
-      if (primary) fSumPrimary->Add(primary);
+      
+      if(fTriggerMask == AliAODForwardMult::kInel) {
+       if (primary && !fSumPrimary) fSumPrimary = CloneHist(primary, "truth");
+       else if (primary) fSumPrimary->Add(primary);
+      }
+      
+      else if(fTriggerMask == AliAODForwardMult::kNSD && forward->IsTriggerBits(AliAODForwardMult::kMCNSD)) {
+       if (primary && !fSumPrimary) fSumPrimary = CloneHist(primary, "truth"); 
+       else if (primary) fSumPrimary->Add(primary);
+      }
+      
+      
+      
+      /*   if(fTriggerMask == AliAODForwardMult::kNSD && forward->IsTriggerBits(AliAODForwardMult::kMCNSD)) {
+       if (primary) fSumPrimary->Add(primary); 
+      
+      }
+      
+      else if (primary) fSumPrimary->Add(primary);*/
+      
     }    
   }
-
+  
   // Here, we get the update 
   if (!fSNNString) { 
     UShort_t sNN = forward->GetSNN();
@@ -106,12 +124,32 @@ AliForwarddNdetaTask::Terminate(Option_t *)
 
   if (fSumPrimary) { 
     Int_t nAll        = Int_t(fTriggers->GetBinContent(kAll));
-    TH1D* dndetaTruth = ProjectX(fSumPrimary,"dndetaTruth",
-                                1,fSumPrimary->GetNbinsY());
-    dndetaTruth->Scale(1./nAll, "width");
-
-    SetHistogramAttributes(dndetaTruth, kGray+3, 22, "Monte-Carlo truth");
-
+    Int_t nNSD        = Int_t(fTriggers->GetBinContent(kMCNSD));
+    
+    //  for(Int_t nn =1; nn<=fSumPrimary->GetNbinsX(); nn++) {
+    //  for(Int_t mm =1; mm<=fSumPrimary->GetNbinsY(); mm++) {
+    // if(fSumPrimary->GetBinContent(nn,mm) > 0) std::cout<<fSumPrimary->GetBinContent(nn,mm)<<" +/-  "<<fSumPrimary->GetBinError(nn,mm)<<std::endl;
+    //  }
+    // }
+    
+    //TH1D* dndetaTruth = ProjectX(fSumPrimary,"dndetaTruth",
+    //                  1,fSumPrimary->GetNbinsY(),false,true);
+    TH1D* dndetaTruth = fSumPrimary->ProjectionX("dndetaTruth",1,fSumPrimary->GetNbinsY(),"e");
+    std::cout<<nAll<<"   "<<nNSD<<std::endl;
+    
+    if(fTriggerMask == AliAODForwardMult::kNSD)
+      dndetaTruth->Scale(1./nNSD, "width");
+    else
+      dndetaTruth->Scale(1./nAll, "width");
+    
+    //for(Int_t nn =1; nn<=dndetaTruth->GetNbinsX(); nn++) {
+    //  if(dndetaTruth->GetBinContent(nn) > 0) std::cout<<dndetaTruth->GetBinContent(nn)<<" +/-  "<<dndetaTruth->GetBinError(nn)<<std::endl;
+    // }
+
+    
+    SetHistogramAttributes(dndetaTruth, kBlue+3, 22, "Monte-Carlo truth");
+    
     fOutput->Add(dndetaTruth);
     fOutput->Add(Rebin(dndetaTruth));
   }
index ce47a0e..83edcfd 100644 (file)
@@ -125,7 +125,7 @@ RunCopyCentralSecMap(UShort_t sys, UShort_t cms, Short_t field, const Char_t* pa
 
   }
   
-  AliCentralMultiplicityTask* task = new AliCentralMultiplicityTask("central");
+  AliCentralMultiplicityTask::Manager* task = new AliCentralMultiplicityTask::Manager();
 
   TFile f(task->GetFullFileName(0.,sys+1,cms,field), "RECREATE");
   m->Write(task->GetSecMapName());