Fixed up default value to be those in stock ForwardAODConfig.C
authorcholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 3 Mar 2011 12:13:45 +0000 (12:13 +0000)
committercholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 3 Mar 2011 12:13:45 +0000 (12:13 +0000)
Added script ForwardAODConfig.C to configure task - user can
override this by copying and editing in working directory

Added calculation of N_ch from Poisson method - needs a bit of
work.

Fixed up some of the flow code.   Marked obvious memory leaks and
such.

18 files changed:
PWG2/FORWARD/analysis2/AddTaskFMD.C
PWG2/FORWARD/analysis2/AliCentralMultiplicityTask.cxx
PWG2/FORWARD/analysis2/AliFMDCorrector.cxx
PWG2/FORWARD/analysis2/AliFMDDensityCalculator.cxx
PWG2/FORWARD/analysis2/AliFMDDensityCalculator.h
PWG2/FORWARD/analysis2/AliFMDEnergyFitter.cxx
PWG2/FORWARD/analysis2/AliFMDEventInspector.cxx
PWG2/FORWARD/analysis2/AliFMDHistCollector.h
PWG2/FORWARD/analysis2/AliForwardFlowBase.cxx
PWG2/FORWARD/analysis2/AliForwardFlowBase.h
PWG2/FORWARD/analysis2/AliForwardFlowTaskQC.cxx
PWG2/FORWARD/analysis2/AliForwardFlowTaskQC.h
PWG2/FORWARD/analysis2/AliForwardMultiplicityBase.cxx
PWG2/FORWARD/analysis2/AliForwardMultiplicityTask.cxx
PWG2/FORWARD/analysis2/ForwardAODConfig.C [new file with mode: 0644]
PWG2/FORWARD/analysis2/MakeAOD.C
PWG2/FORWARD/analysis2/MakeFlow.C [new file with mode: 0644]
PWG2/FORWARD/analysis2/RunManagerFlow.C [deleted file]

index eab2982a0fa76cf71d32d9eaf0cf4673d7fcf19a..12ade395a885e4cb9a0e2baa890229f80b080285 100644 (file)
@@ -33,71 +33,22 @@ AddTaskFMD(Bool_t mc, UShort_t sys=0, UShort_t sNN=0, Short_t field)
   // --- Do a local initialisation with assumed values ---------------
   if (sys > 0 && sNN > 0) 
     AliForwardCorrectionManager::Instance().Init(sys,sNN,field);
-  
-  // Whether to enable low flux specific code 
-  task->SetEnableLowFlux(kFALSE);
-  // Set the number of SPD tracklets for which we consider the event a
-  // low flux event
-  task->GetEventInspector().SetLowFluxCut(1000); 
-  // Set the maximum error on v_z [cm]
-  task->GetEventInspector().SetMaxVzErr(0.2);
-  // Set the eta axis to use - note, this overrides whatever is used
-  // by the rest of the algorithms - but only for the energy fitter
-  // algorithm. 
-  task->GetEnergyFitter().SetEtaAxis(200, -4, 6);
-  // Set maximum energy loss to consider 
-  task->GetEnergyFitter().SetMaxE(10); 
-  // Set number of energy loss bins 
-  task->GetEnergyFitter().SetNEbins(300);
-  // Set whether to use increasing bin sizes 
-  task->GetEnergyFitter().SetUseIncreasingBins(true);
-  // Set whether to do fit the energy distributions 
-  task->GetEnergyFitter().SetDoFits(kFALSE);
-  // Set whether to make the correction object 
-  task->GetEnergyFitter().SetDoMakeObject(kFALSE);
-  // Set the low cut used for energy
-  task->GetEnergyFitter().SetLowCut(0.4);
-  // Set the number of bins to subtract from maximum of distributions
-  // to get the lower bound of the fit range
-  task->GetEnergyFitter().SetFitRangeBinWidth(4);
-  // Set the maximum number of landaus to try to fit (max 5)
-  task->GetEnergyFitter().SetNParticles(5);
-  // Set the minimum number of entries in the distribution before
-  // trying to fit to the data
-  task->GetEnergyFitter().SetMinEntries(1000);
-  // Set the low cut used for sharing - overrides settings in eloss fits
-  task->GetSharingFilter().SetLowCut(0.3);
-  // Set the number of xi's (width of landau peak) to stop at 
-  task->GetSharingFilter().SetNXi(1);
-  // Set the maximum number of particle to try to reconstruct 
-  task->GetDensityCalculator().SetMaxParticles(3);
-  // Set the lower multiplicity cut.  Overrides setting in energy loss fits.
-  task->GetDensityCalculator().SetMultCut(0.3); //was 0.3
-  // Whether to use the secondary map correction
-  task->GetCorrections().SetUseSecondaryMap(true);
-  // Whether to use the vertex bias correction
-  task->GetCorrections().SetUseVertexBias(false);
-  // Whether to use the vertex bias correction
-  task->GetCorrections().SetUseAcceptance(true);
-  // Whether to use the merging efficiency correction 
-  task->GetCorrections().SetUseMergingEfficiency(false);
-  // Set the number of extra bins (beyond the secondary map border) 
-  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.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 
-  // task->GetEventInspector().SetDebug(4);
-  // --- Set limits on fits the energy -------------------------------
-  // Maximum relative error on parameters 
-  AliFMDCorrELossFit::ELossFit::fgMaxRelError = .12;
-  // Least weight to use 
-  AliFMDCorrELossFit::ELossFit::fgLeastWeight = 1e-5;
-  // Maximum value of reduced chi^2 
-  AliFMDCorrELossFit::ELossFit::fgMaxChi2nu   = 5;
 
+  // --- Configure the task ------------------------------------------
+  gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWG2/FORWARD/analysis2", 
+                          gROOT->GetMacroPath()));
+  const char* config = gSystem->Which(gROOT->GetMacroPath(),
+                                     "ForwardAODConfig.C");
+  if (!config) 
+    Warning("AddTaskFMD", "ForwardAODConfig.C not found in %s",
+           gROOT->GetMacroPath());
+  else {
+    Info("AddTaskFMD", 
+        "Loading configuration of AliForwardMultiplicityTask from %s",
+        config);
+    gROOT->Macro(Form("%s((AliForwardMultiplicityBase*)%p)", config, task));
+    delete config;
+  }
   
   // --- Make the output container and connect it --------------------
   TString outputfile = AliAnalysisManager::GetCommonFileName();
index aa98bf8cb53a118e09f1fe14e2d2890c4d325b6b..bf6ead9258adce4456c6612f4cfa1a7bac20a51b 100644 (file)
@@ -90,6 +90,7 @@ void AliCentralMultiplicityTask::UserCreateOutputObjects()
   ah->AddBranch("AliAODCentralMult", &obj);
   
   fList = new TList();
+  fList->SetOwner();
   PostData(1,fList);
   
 }
index 02d10b11576082ee63887f80f3965b58fb23d174..0b7681cd36c6b98f2c602a2c2a305f3473722b6e 100644 (file)
@@ -24,9 +24,9 @@ AliFMDCorrector::AliFMDCorrector()
   : TNamed(), 
     fRingHistos(),
     fUseSecondaryMap(true),
-    fUseVertexBias(true),
+    fUseVertexBias(false),
     fUseAcceptance(true),
-    fUseMergingEfficiency(true),
+    fUseMergingEfficiency(false),
     fDebug(0)
 {
   // Constructor
@@ -37,9 +37,9 @@ AliFMDCorrector::AliFMDCorrector(const char* title)
   : TNamed("fmdCorrector", title), 
     fRingHistos(), 
     fUseSecondaryMap(true),
-    fUseVertexBias(true),
+    fUseVertexBias(false),
     fUseAcceptance(true),
-    fUseMergingEfficiency(true),
+    fUseMergingEfficiency(false),
     fDebug(0)
 {
   // Constructor 
index 038eea6f17b34db698727f2bf596171908deb961..7e0a8382088a881ece1ecb140fbe37ae9c57a5d3 100644 (file)
@@ -28,6 +28,7 @@ AliFMDDensityCalculator::AliFMDDensityCalculator()
     fWeightedSum(0),
     fCorrections(0),
     fMaxParticles(5),
+    fUsePoisson(false),
     fAccI(0),
     fAccO(0),
     fFMD1iMax(0),
@@ -52,6 +53,7 @@ AliFMDDensityCalculator::AliFMDDensityCalculator(const char* title)
     fWeightedSum(0),
     fCorrections(0),
     fMaxParticles(5),
+    fUsePoisson(false),
     fAccI(0),
     fAccO(0),
     fFMD1iMax(0),
@@ -100,6 +102,7 @@ AliFMDDensityCalculator::AliFMDDensityCalculator(const
     fWeightedSum(o.fWeightedSum),
     fCorrections(o.fCorrections),
     fMaxParticles(o.fMaxParticles),
+    fUsePoisson(o.fUsePoisson),
     fAccI(o.fAccI),
     fAccO(o.fAccO),
     fFMD1iMax(o.fFMD1iMax),
@@ -147,6 +150,7 @@ AliFMDDensityCalculator::operator=(const AliFMDDensityCalculator& o)
   fMultCut      = o.fMultCut;
   fDebug        = o.fDebug;
   fMaxParticles = o.fMaxParticles;
+  fUsePoisson   = o.fUsePoisson;
   fAccI         = o.fAccI;
   fAccO         = o.fAccO;
   fFMD1iMax     = o.fFMD1iMax;
@@ -165,13 +169,18 @@ AliFMDDensityCalculator::operator=(const AliFMDDensityCalculator& o)
 
 //____________________________________________________________________
 void
-AliFMDDensityCalculator::Init(const TAxis&)
+AliFMDDensityCalculator::Init(const TAxis& axis)
 {
   // Intialize this sub-algorithm 
   //
   // Parameters:
   //   etaAxis   Not used
   CacheMaxWeights();
+
+  TIter    next(&fRingHistos);
+  RingHistos* o = 0;
+  while ((o = static_cast<RingHistos*>(next())))
+    o->Init(axis);
 }
 
 //____________________________________________________________________
@@ -245,17 +254,24 @@ AliFMDDensityCalculator::Calculate(const AliESDFMD&        fmd,
       UShort_t    nt= (q == 0 ? 512 : 256);
       TH2D*       h = hists.Get(d,r);
       RingHistos* rh= GetRingHistos(d,r);
+      rh->fEmptyStrips->Reset();
+      rh->fTotalStrips->Reset();
 
       for (UShort_t s=0; s<ns; s++) { 
        for (UShort_t t=0; t<nt; t++) {
          Float_t mult = fmd.Multiplicity(d,r,s,t);
-         
-         if (mult == 0 || mult > 20) continue;
+         Float_t phi  = fmd.Phi(d,r,s,t) / 180 * TMath::Pi();
+         Float_t eta  = fmd.Eta(d,r,s,t);
 
-         Float_t phi = fmd.Phi(d,r,s,t) / 180 * TMath::Pi();
-         Float_t eta = fmd.Eta(d,r,s,t);
+         if (mult == AliESDFMD::kInvalidMult || mult > 20) continue;
+         rh->fTotalStrips->Fill(eta, phi);
          
-         Float_t n = NParticles(mult,d,r,s,t,vtxbin,eta,lowFlux);
+         if (mult == 0) { 
+           rh->fEmptyStrips->Fill(eta,phi);
+           continue;
+         }
+           
+         Float_t n   = NParticles(mult,d,r,s,t,vtxbin,eta,lowFlux);
          rh->fEvsN->Fill(mult,n);
          rh->fEtaVsN->Fill(eta, n);
 
@@ -266,10 +282,35 @@ AliFMDDensityCalculator::Calculate(const AliESDFMD&        fmd,
          rh->fEtaVsM->Fill(eta, n);
          rh->fCorr->Fill(eta, c);
 
+         if (n < .9) rh->fEmptyStrips->Fill(eta,phi);
+
          h->Fill(eta,phi,n);
          rh->fDensity->Fill(eta,phi,n);
        } // for t
       } // for s 
+
+      // Loop over poisson histograms 
+      for (Int_t ieta = 1; ieta <= h->GetNbinsX(); ieta++) { 
+       for (Int_t iphi = 1; iphi <= h->GetNbinsY(); iphi++) { 
+         Double_t eLossV   = h->GetBinContent(ieta, iphi);
+         // Double_t eLossE   = h->GetBinError(ieta, iphi);
+         
+         Double_t empty    = rh->fEmptyStrips->GetBinContent(ieta,iphi);
+         Double_t total    = rh->fTotalStrips->GetBinContent(ieta,iphi);
+
+         Double_t poissonV =  (total <= 0 || empty <= 0 ? 0 : 
+                               -TMath::Log(empty / total));
+         Double_t poissonE = TMath::Sqrt(poissonV);
+
+         rh->fELossVsPoisson->Fill(eLossV, poissonV);
+         rh->fEmptyVsTotal->Fill(total, empty);
+         if (fUsePoisson) {
+           h->SetBinContent(ieta,iphi,poissonV);
+           h->SetBinError(ieta,iphi,poissonE);
+         }
+       }
+      }
+       
     } // for q
   } // for d
   
@@ -404,25 +445,6 @@ AliFMDDensityCalculator::NParticles(Float_t  mult,
   fSumOfWeights->Fill(ret);
 
   return ret;
-#if 0
-  Float_t mpv  = pars->GetMPV(d,r,eta);
-  Float_t w    = pars->GetSigma(d,r,eta);
-  Float_t w2   = pars->Get2MIPWeight(d,r,eta);
-  Float_t w3   = pars->Get3MIPWeight(d,r,eta);
-  Float_t mpv2 = 2*mpv+2*w*TMath::Log(2);
-  Float_t mpv3 = 3*mpv+3*w*TMath::Log(3);
-  
-  Float_t sum  = (TMath::Landau(mult,mpv,w,kTRUE) +
-                 w2 * TMath::Landau(mult,mpv2,2*w,kTRUE) + 
-                 w3  * TMath::Landau(mult,mpv3,3*w,kTRUE));
-  Float_t wsum = (TMath::Landau(mult,mpv,w,kTRUE) +
-                 2*w2 * TMath::Landau(mult,mpv2,2*w,kTRUE) + 
-                 3*w3 * TMath::Landau(mult,mpv3,3*w,kTRUE));
-  
-  fWeightedSum->Fill(wsum);
-  fSumOfWeights->Fill(sum);
-  return (sum > 0) ? wsum / sum : 1;
-#endif
 }
 
 //_____________________________________________________________________
@@ -470,16 +492,6 @@ AliFMDDensityCalculator::Correction(UShort_t d,
       AliWarning(Form("Missing double hit correction for FMD%d%c",d,r));
     }
   }
-  
-#if 0
-  TH1F* deadCor = pars->GetFMDDeadCorrection(v);
-  if (deadCor) { 
-    Float_t deadC = deadCor->GetBinContent(deadCor->FindBin(eta));
-    if (deadC > 0) 
-      correction *= deadC; 
-  }
-#endif  
-
   return correction;
 }
 
@@ -577,59 +589,6 @@ AliFMDDensityCalculator::AcceptanceCorrection(Char_t r, UShort_t t) const
   //
   TH1D* acc = (r == 'I' || r == 'i' ? fAccI : fAccO);
   return acc->GetBinContent(t+1);
-
-#if 0
-  const Double_t ic1[] = { 4.9895, 15.3560 };
-  const Double_t ic2[] = { 1.8007, 17.2000 };
-  const Double_t oc1[] = { 4.2231, 26.6638 };
-  const Double_t oc2[] = { 1.8357, 27.9500 };
-  const Double_t* c1   = (r == 'I' ? ic1      : oc1);
-  const Double_t* c2   = (r == 'I' ? ic2      : oc2);
-  Double_t  minR       = (r == 'I' ?   4.5213 :  15.4);
-  Double_t  maxR       = (r == 'I' ?  17.2    :  28.0);
-  Int_t     nStrips    = (r == 'I' ? 512      : 256);
-  Int_t     nSec       = (r == 'I' ?  20      :  40);
-  Float_t   basearc    = 2 * TMath::Pi() / nSec;
-  Double_t  rad        = maxR - minR;
-  Float_t   segment    = rad / nStrips;
-  Float_t   radius     = minR + t * segment;
-
-  // Old method - calculate full strip area and take ratio to extended
-  // strip area
-  Float_t   baselen    = basearc * radius;
-  Float_t   slope      = (c1[1] - c2[1]) / (c1[0] - c2[0]);
-  Float_t   constant   = (c2[1] * c1[0] - c2[0] * c1[1]) / (c1[0]-c2[0]);
-  Float_t   d          = (TMath::Power(TMath::Abs(radius*slope),2) + 
-                         TMath::Power(radius,2) - TMath::Power(constant,2));
-
-  // If below corners return 1
-  if (d >= 0) return 1;
-  Float_t   x         = ((-TMath::Sqrt(d) - slope * constant) / 
-                        (1+TMath::Power(slope, 2)));
-  Float_t   y         = slope*x + constant;
-
-  // If x is larger than corner x or y less than corner y, we have full
-  // length strip
-  if(x >= c1[0] || y <= c1[1]) return 1;
-
-  //One sector since theta is by definition half-hybrid
-  Float_t   theta     = TMath::ATan2(x,y);
-  Float_t   arclen    = radius * theta;
-  
-  // Calculate the area of a strip with no cut
-  Float_t   basearea1 = 0.5 * baselen * TMath::Power(radius,2);
-  Float_t   basearea2 = 0.5 * baselen * TMath::Power((radius-segment),2);
-  Float_t   basearea  = basearea1 - basearea2;
-
-  // Calculate the area of a strip with cut
-  Float_t   area1     = 0.5 * arclen * TMath::Power(radius,2);
-  Float_t   area2     = 0.5 * arclen * TMath::Power((radius-segment),2);
-  Float_t   area      = area1 - area2;
-  
-  // Acceptance is ratio 
-  return area/basearea;
-#endif
 }
 
 //____________________________________________________________________
@@ -734,7 +693,11 @@ AliFMDDensityCalculator::RingHistos::RingHistos()
     fEtaVsN(0),
     fEtaVsM(0),
     fCorr(0),
-    fDensity(0)
+    fDensity(0),
+    fELossVsPoisson(0),
+    fTotalStrips(0),
+    fEmptyStrips(0),
+    fEmptyVsTotal(0)
 {
   // 
   // Default CTOR
@@ -749,7 +712,11 @@ AliFMDDensityCalculator::RingHistos::RingHistos(UShort_t d, Char_t r)
     fEtaVsN(0),
     fEtaVsM(0),
     fCorr(0),
-    fDensity(0)
+    fDensity(0),
+    fELossVsPoisson(0),
+    fTotalStrips(0),
+    fEmptyStrips(0),
+    fEmptyVsTotal(0)
 {
   // 
   // Constructor
@@ -810,6 +777,23 @@ AliFMDDensityCalculator::RingHistos::RingHistos(UShort_t d, Char_t r)
   fDensity->SetXTitle("#eta");
   fDensity->SetYTitle("#phi [radians]");
   fDensity->SetZTitle("Inclusive N_{ch} density");
+
+  fELossVsPoisson = new TH2D(Form("%s_eloss_vs_poisson", fName.Data()),
+                            Form("N_{ch} from energy loss vs from Poission %s",
+                                 fName.Data()), 100, 0, 20, 100, 0, 20);
+  fELossVsPoisson->SetDirectory(0);
+  fELossVsPoisson->SetXTitle("N_{ch} from #DeltaE");
+  fELossVsPoisson->SetYTitle("N_{ch} from Poisson");
+  fELossVsPoisson->SetZTitle("Correlation");
+
+  fEmptyVsTotal = new TH2D(Form("%s_empty_vs_total", fName.Data()), 
+                          Form("# of empty strips vs. total @ # strips in %s", 
+                               fName.Data()), 21, -.5, 20.5, 21, -0.5, 20.5);
+  fEmptyVsTotal->SetDirectory(0);
+  fEmptyVsTotal->SetXTitle("Total # strips");
+  fEmptyVsTotal->SetYTitle("Empty # strips");
+  fEmptyVsTotal->SetZTitle("Correlation");
+
 }
 //____________________________________________________________________
 AliFMDDensityCalculator::RingHistos::RingHistos(const RingHistos& o)
@@ -819,7 +803,11 @@ AliFMDDensityCalculator::RingHistos::RingHistos(const RingHistos& o)
     fEtaVsN(o.fEtaVsN),
     fEtaVsM(o.fEtaVsM),
     fCorr(o.fCorr),
-    fDensity(o.fDensity)
+    fDensity(o.fDensity),
+    fELossVsPoisson(o.fELossVsPoisson),
+    fTotalStrips(o.fTotalStrips),
+    fEmptyStrips(o.fEmptyStrips),
+    fEmptyVsTotal(o.fEmptyVsTotal)
 {
   // 
   // Copy constructor 
@@ -844,19 +832,25 @@ AliFMDDensityCalculator::RingHistos::operator=(const RingHistos& o)
   //
   AliForwardUtil::RingHistos::operator=(o);
   
-  if (fEvsN)    delete  fEvsN;
-  if (fEvsM)    delete  fEvsM;
-  if (fEtaVsN)  delete  fEtaVsN;
-  if (fEtaVsM)  delete  fEtaVsM;
-  if (fCorr)    delete  fCorr;
-  if (fDensity) delete fDensity;
+  if (fEvsN)           delete fEvsN;
+  if (fEvsM)           delete fEvsM;
+  if (fEtaVsN)         delete fEtaVsN;
+  if (fEtaVsM)         delete fEtaVsM;
+  if (fCorr)           delete fCorr;
+  if (fDensity)        delete fDensity;
+  if (fELossVsPoisson) delete fELossVsPoisson;
+  if (fTotalStrips)    delete fTotalStrips;
+  if (fEmptyStrips)    delete fEmptyStrips;
   
-  fEvsN    = static_cast<TH2D*>(o.fEvsN->Clone());
-  fEvsM    = static_cast<TH2D*>(o.fEvsM->Clone());
-  fEtaVsN  = static_cast<TProfile*>(o.fEtaVsN->Clone());
-  fEtaVsM  = static_cast<TProfile*>(o.fEtaVsM->Clone());
-  fCorr    = static_cast<TProfile*>(o.fCorr->Clone());
-  fDensity = static_cast<TH2D*>(o.fDensity->Clone());
+  fEvsN           = static_cast<TH2D*>(o.fEvsN->Clone());
+  fEvsM           = static_cast<TH2D*>(o.fEvsM->Clone());
+  fEtaVsN         = static_cast<TProfile*>(o.fEtaVsN->Clone());
+  fEtaVsM         = static_cast<TProfile*>(o.fEtaVsM->Clone());
+  fCorr           = static_cast<TProfile*>(o.fCorr->Clone());
+  fDensity        = static_cast<TH2D*>(o.fDensity->Clone());
+  fELossVsPoisson = static_cast<TH2D*>(o.fELossVsPoisson);
+  fTotalStrips    = static_cast<TH2D*>(o.fTotalStrips);
+  fEmptyStrips    = static_cast<TH2D*>(o.fEmptyStrips);
   
   return *this;
 }
@@ -866,12 +860,47 @@ AliFMDDensityCalculator::RingHistos::~RingHistos()
   // 
   // Destructor 
   //
-  if (fEvsN)    delete fEvsN;
-  if (fEvsM)    delete fEvsM;
-  if (fEtaVsN)  delete fEtaVsN;
-  if (fEtaVsM)  delete fEtaVsM;
-  if (fCorr)    delete fCorr;
-  if (fDensity) delete fDensity;
+  if (fEvsN)           delete fEvsN;
+  if (fEvsM)           delete fEvsM;
+  if (fEtaVsN)         delete fEtaVsN;
+  if (fEtaVsM)         delete fEtaVsM;
+  if (fCorr)           delete fCorr;
+  if (fDensity)        delete fDensity;
+  if (fELossVsPoisson) delete fELossVsPoisson;
+  if (fTotalStrips)    delete fTotalStrips;
+  if (fEmptyStrips)    delete fEmptyStrips;
+}
+
+  //____________________________________________________________________
+void
+AliFMDDensityCalculator::RingHistos::Init(const TAxis& eAxis)
+{
+  if (fTotalStrips) delete fTotalStrips;
+  if (fEmptyStrips) delete fEmptyStrips;
+
+  fTotalStrips = new TH2D(Form("total%s", fName.Data()), 
+                         Form("Total number of strips in %s", fName.Data()), 
+                         eAxis.GetNbins(), 
+                         eAxis.GetXmin(),
+                         eAxis.GetXmax(), 
+                         (fRing == 'I' || fRing == 'i' ? 20 : 40), 
+                         0, 2*TMath::Pi());
+  fEmptyStrips = new TH2D(Form("empty%s", fName.Data()), 
+                         Form("Empty number of strips in %s", fName.Data()), 
+                         eAxis.GetNbins(), 
+                         eAxis.GetXmin(), 
+                         eAxis.GetXmax(), 
+                         (fRing == 'I' || fRing == 'i' ? 20 : 40), 
+                         0, 2*TMath::Pi());
+  fTotalStrips->SetDirectory(0);
+  fEmptyStrips->SetDirectory(0);
+  fTotalStrips->SetXTitle("#eta");
+  fEmptyStrips->SetXTitle("#eta");
+  fTotalStrips->SetYTitle("#varphi [radians]");
+  fEmptyStrips->SetYTitle("#varphi [radians]");
+  fTotalStrips->Sumw2();
+  fEmptyStrips->Sumw2();
+
 }
 
 //____________________________________________________________________
@@ -891,6 +920,8 @@ AliFMDDensityCalculator::RingHistos::Output(TList* dir)
   d->Add(fEtaVsM);
   d->Add(fCorr);
   d->Add(fDensity);
+  d->Add(fELossVsPoisson);
+  d->Add(fEmptyVsTotal);
 }
 
 //____________________________________________________________________
index 3bed8ea48943af62bb2b5a32242d06eb933facbd..108ff44ac9f8b5fc5e42aab775210f567735d9e4 100644 (file)
@@ -106,6 +106,21 @@ public:
    * @param m 
    */
   void SetMaxParticles(UShort_t m) { fMaxParticles = m; }  
+  /** 
+   * Set whether to use poisson statistics to estimate the 
+   * number of particles that has hit within a region.  If this is true, 
+   * then the average charge particle density is given by 
+   * @f[
+   *  \lamda = -\log\left(\frac{N_e}{N_t}\right)
+   * @f]
+   * where $N_e$ is the number of strips within the region that has no
+   * hits over threshold, and $N_t$ is the total number of strips in the 
+   * region/ 
+   * 
+   * @param u Whether to use poisson statistics to estimate the 
+   * number of particles that has hit within a region.
+   */
+  void SetUsePoisson(Bool_t u) { fUsePoisson = u; }
   /** 
    * Set the lower multiplicity cut.  This overrides the setting in
    * the energy loss fits.
@@ -255,6 +270,12 @@ protected:
      * Destructor 
      */
     ~RingHistos();
+    /** 
+     * Initialize the object 
+     * 
+     * @param eAxis 
+     */
+    void Init(const TAxis& eAxis);
     /** 
      * Make output 
      * 
@@ -274,6 +295,10 @@ protected:
     TProfile* fEtaVsM;       // Average corrected Nch vs eta
     TProfile* fCorr;         // Average correction vs eta
     TH2D*     fDensity;      // Distribution inclusive Nch
+    TH2D*     fELossVsPoisson; // Correlation of energy loss vs Poisson N_ch
+    TH2D*     fTotalStrips;  //! Total number of strips in a region
+    TH2D*     fEmptyStrips;  //! Total number of strips in a region
+    TH2D*     fEmptyVsTotal; // # of empty strips vs total number of strips 
     ClassDef(RingHistos,1);
   };
   /** 
@@ -292,6 +317,7 @@ protected:
   TH1D*    fWeightedSum;   //  Histogram
   TH1D*    fCorrections;   //  Histogram
   UShort_t fMaxParticles;  //  Maximum particle weight to use 
+  Bool_t   fUsePoisson;    //  If true, then use poisson statistics 
   TH1D*    fAccI;          //  Acceptance correction for inner rings
   TH1D*    fAccO;          //  Acceptance correction for outer rings
   TArrayI  fFMD1iMax;      //  Array of max weights 
@@ -301,7 +327,7 @@ protected:
   TArrayI  fFMD3oMax;      //  Array of max weights 
   Int_t    fDebug;         //  Debug level 
 
-  ClassDef(AliFMDDensityCalculator,1); // Calculate Nch density 
+  ClassDef(AliFMDDensityCalculator,2); // Calculate Nch density 
 };
 
 #endif
index 732b0005d1c8c1945337e1880bf1d58a536eebfc..395aa21508245aa2b8e761ae573075a8cbd0f119 100644 (file)
@@ -30,9 +30,9 @@ namespace {
 AliFMDEnergyFitter::AliFMDEnergyFitter()
   : TNamed(), 
     fRingHistos(),
-    fLowCut(0.3),
-    fNParticles(3),
-    fMinEntries(100),
+    fLowCut(0.4),
+    fNParticles(5),
+    fMinEntries(1000),
     fFitRangeBinWidth(4),
     fDoFits(false),
     fDoMakeObject(false),
@@ -55,9 +55,9 @@ AliFMDEnergyFitter::AliFMDEnergyFitter()
 AliFMDEnergyFitter::AliFMDEnergyFitter(const char* title)
   : TNamed("fmdEnergyFitter", title), 
     fRingHistos(), 
-    fLowCut(0.3),
-    fNParticles(3),
-    fMinEntries(100),
+    fLowCut(0.4),
+    fNParticles(5),
+    fMinEntries(1000),
     fFitRangeBinWidth(4),
     fDoFits(false),
     fDoMakeObject(false),
index a83788a9a6e2d664c6c91f5049923d3cbde6b6a0..9f76c527e9ec512b7d8358d02c0edf0814a82fa8 100644 (file)
@@ -46,7 +46,7 @@ AliFMDEventInspector::AliFMDEventInspector()
     fHWords(0),
     fHCent(0),
     fLowFluxCut(1000),
-    fMaxVzErr(0.1),
+    fMaxVzErr(0.2),
     fList(0),
     fEnergy(0),
     fField(999), 
@@ -68,7 +68,7 @@ AliFMDEventInspector::AliFMDEventInspector(const char* name)
     fHWords(0),
     fHCent(0),
     fLowFluxCut(1000),
-    fMaxVzErr(0.1),
+    fMaxVzErr(0.2),
     fList(0),
     fEnergy(0),
     fField(999), 
@@ -344,11 +344,12 @@ AliFMDEventInspector::Process(const AliESDEvent* event,
   cent = -10;
   AliCentrality* centObj = const_cast<AliESDEvent*>(event)->GetCentrality();
   if (centObj) {
-    AliInfo(Form("Got centrality object %p with quality %d", centObj, centObj->GetQuality()));
-    centObj->Print();
+    // AliInfo(Form("Got centrality object %p with quality %d", 
+    //              centObj, centObj->GetQuality()));
+    // centObj->Print();
     cent = centObj->GetCentralityPercentileUnchecked("V0M");  
   }
-  AliInfo(Form("Centrality is %f", cent));
+  // AliInfo(Form("Centrality is %f", cent));
   fHCent->Fill(cent);
 
   // Check the FMD ESD data 
index 58690194b177c1a2304ec4df4bb328d68419a057..e6b030e27a81dcc5229abb1b936b248f2c38c1ce 100644 (file)
@@ -43,8 +43,11 @@ public:
    */
   AliFMDHistCollector(const char* title)
     : TNamed("fmdHistCollector", title), 
-      fNCutBins(1), fCorrectionCut(0.5), 
-      fFirstBins(1), fLastBins(1), fDebug(0)
+      fNCutBins(2), 
+      fCorrectionCut(0.5), 
+      fFirstBins(1), 
+      fLastBins(1), 
+      fDebug(0)
   {}
   /** 
    * Copy constructor 
index 008e1b492d28cffdce322771259c64fead122c1b..5401f2ff00f3c5b60c029897fd2b828b6f1013f9 100644 (file)
@@ -51,6 +51,7 @@ Bool_t AliForwardFlowBase::LoopAODFMD(AliAODEvent* AODevent)
   if (!AODFMult) return kFALSE;
   if (!AODCheck(AODFMult)) return kFALSE;
 
+  // Memory leak starts here!
   TH2D* hdNdetadphi = static_cast<TH2D*>(AODFMult->GetHistogram().Clone("d2ndetadphi"));
   TH1D* dNdphi = (TH1D*)fList->FindObject("hdNdphiSE");
   TH2D* dNdetadphi = (TH2D*)fList->FindObject("hdNdetadphiSE");
@@ -83,6 +84,7 @@ Bool_t AliForwardFlowBase::LoopAODSPD(AliAODEvent* AODevent)
   AliAODCentralMult* AODCMult = static_cast<AliAODCentralMult*>(AODevent->FindListObject("CentralClusters"));
   if (!AODCMult) return kFALSE;
 
+  // Memory leak starts here! 
   TH2D* hdNdetadphi = static_cast<TH2D*>(AODCMult->GetHistogram().Clone("central"));
   TH1D* dNdphi = (TH1D*)fList->FindObject("hdNdphiSE");
   TH2D* dNdetadphi = (TH2D*)fList->FindObject("hdNdetadphiSE");
index c86d01b453b649aaa75dcd94eab8fd7401023bb1..0738f15e48577439b50df7f2316b57573087ce7e 100644 (file)
  * Class used to handle the input from AODs and put it into histograms
  * the Forward Flow tasks can run on.
  * 
+ * @todo Rename class to something like AliForwardFlowUtil or the like
+ * - but no 'Base' in the name - misleading - one thinks its a base
+ * class
+ *
  * @ingroup pwg2_forward_tasks
  */
 class AliForwardFlowBase : public TNamed
@@ -20,21 +24,21 @@ public:
   /**
    * Constructor
    */
-   AliForwardFlowBase();
+  AliForwardFlowBase();
   /*
    * Constructor
    *
    * @param fList list of histograms for flow analysis
    */
-   AliForwardFlowBase(TList* fList);
+  AliForwardFlowBase(TList* fList);
   /*
    * Copy constructor
    *
    * @param o Object to copy from
    */
-    AliForwardFlowBase(const AliForwardFlowBase& o) : TNamed(),
-       fList(o.fList) {}
- /** 
+  AliForwardFlowBase(const AliForwardFlowBase& o) : TNamed(),
+                                                   fList(o.fList) {}
 /** 
    * Assignment operator 
    * 
    * @param o Object to assign from 
@@ -42,47 +46,52 @@ public:
    * @return Reference to this object 
    */
   AliForwardFlowBase& operator=(const AliForwardFlowBase&) { return *this; }
-  /*
+  /**
    * Check that AOD event meet trigger requirements
    */
-   Bool_t  AODCheck(AliAODForwardMult* aodfm);
-  /*
+  Bool_t  AODCheck(AliAODForwardMult* aodfm);
+  /**
    * Loop over AliAODForwardMult object and fill flow histograms
    */
-   Bool_t  LoopAODFMD(AliAODEvent* AODevent);
+  Bool_t  LoopAODFMD(AliAODEvent* AODevent);
   /*
    * Loop over AliAODForwardCentral object and fill flow histograms
    */
-   Bool_t  LoopAODSPD(AliAODEvent* AODevent);
-  /*
-   * Loop over AliAODForwardMult and Central object and fill flow histograms
+  Bool_t  LoopAODSPD(AliAODEvent* AODevent);
+  /**
+   * Loop over AliAODForwardMult and Central object and fill flow
+   * histograms
    */
-   Bool_t  LoopAODFMDandSPD(AliAODEvent* AODevent);
-  /*
+  Bool_t  LoopAODFMDandSPD(AliAODEvent* AODevent);
+  /**
    * Loop over AliAODMCParticle branch object and fill flow histograms
    */
-   Bool_t  LoopAODmc(AliAODEvent* AODevent);
-  /*
-   * Loop over AliAODForwardMult object and fill flow histograms from track refs
+  Bool_t  LoopAODmc(AliAODEvent* AODevent);
+  /**
+   * Loop over AliAODForwardMult object and fill flow histograms from
+   * track refs
    */
-   Bool_t  LoopAODtrrefHits(AliAODEvent* AODevent);
-  /*
-   * Loop over AliAODMCParticle branch object, add pt dep. flow and fill flow histograms
+  Bool_t  LoopAODtrrefHits(AliAODEvent* AODevent);
+  /**
+   * Loop over AliAODMCParticle branch object, add pt dep. flow and
+   * fill flow histograms
    */
-   Bool_t  LoopMCaddptFlow(AliAODEvent* AODevent);
-  /*
-   * Loop over AliAODMCParticle branch object, add pit dep. flow and fill flow histograms
+  Bool_t  LoopMCaddptFlow(AliAODEvent* AODevent);
+  /**
+   * Loop over AliAODMCParticle branch object, add pit dep. flow and
+   * fill flow histograms
    */
-   Bool_t  LoopMCaddpdgFlow(AliAODEvent* AODevent);
-  /*
-   * Loop over AliAODMCParticle branch object, add eta dep flow and fill flow histograms
+  Bool_t  LoopMCaddpdgFlow(AliAODEvent* AODevent);
+  /**
+   * Loop over AliAODMCParticle branch object, add eta dep flow and
+   * fill flow histograms
    */
-   Bool_t  LoopMCaddetaFlow(AliAODEvent* AODevent);
+  Bool_t  LoopMCaddetaFlow(AliAODEvent* AODevent);
 
-private:  
+protected:
   TList*       fList;  // List of flow histograms
   
-  ClassDef(AliForwardFlowBase, 0); 
+  ClassDef(AliForwardFlowBase, 1); 
 };
  
 #endif
index 099c6fc24aa74890124e25f831f3d48d0eef561a..770a529fb346b03b8f53fa363569169f45ca085e 100644 (file)
 #include "AliAODForwardMult.h"
 
 ClassImp(AliForwardFlowTaskQC)
+#if 0
+; // For emacs 
+#endif
 
 AliForwardFlowTaskQC::AliForwardFlowTaskQC()
-: fDebug(0),           // Debug flag
-  fOutputList(0),      // Output list
-  fAOD(0),             // AOD input event
-  fMC(kFALSE),         // MC flag
-  fEtaBins(20)         // # of eta bins in histograms
+  : fDebug(0),                 // Debug flag
+    fOutputList(0),    // Output list
+    fAOD(0),           // AOD input event
+    fMC(kFALSE),               // MC flag
+    fEtaBins(20)               // # of eta bins in histograms
 {
   // 
   // Default constructor
@@ -222,8 +225,8 @@ void AliForwardFlowTaskQC::UserExec(Option_t */*option*/)
   AliForwardFlowBase* common = new AliForwardFlowBase(fOutputList);
 
   if (!common->LoopAODFMD(fAOD)) return;
-//  else if (!common->LoopAODSPD(fAOD)) return;
-//  if (!common->LoopAODFMDandSPD(fAOD)) return;
+  //  else if (!common->LoopAODSPD(fAOD)) return;
+  //  if (!common->LoopAODFMDandSPD(fAOD)) return;
 
   // Run analysis
   for (Int_t n = 1; n <= 4; n++) {
@@ -242,7 +245,8 @@ void AliForwardFlowTaskQC::UserExec(Option_t */*option*/)
 
 }
 //_____________________________________________________________________
-void AliForwardFlowTaskQC::CumulantsMethod(TString type = "", Int_t harmonic = 2)
+void AliForwardFlowTaskQC::CumulantsMethod(TString type = "", 
+                                          Int_t harmonic = 2)
 {
   // 
   // Calculate the Q cumulant of order n
@@ -616,7 +620,7 @@ void AliForwardFlowTaskQC::Terminate(Option_t */*option*/)
         d4 += 4.*AvgSin2Psi1Phi2*AvgCos2Phi*AvgSin2Phi;
         d4 += 4.*Avg2p*(TMath::Power(AvgCos2Phi, 2) + TMath::Power(AvgSin2Phi, 2)); 
         d4 -= 6.*(TMath::Power(AvgCos2Phi, 2) - TMath::Power(AvgSin2Phi, 2))
-              *(AvgCos2Psi*AvgCos2Phi-AvgSin2Psi*AvgSin2Phi); 
+         *(AvgCos2Psi*AvgCos2Phi-AvgSin2Psi*AvgSin2Phi); 
         d4 -= 12.*AvgCos2Phi*AvgSin2Phi*(AvgSin2Psi*AvgCos2Phi+AvgCos2Psi*AvgSin2Phi); 
  
         vTwo4diff = - d4 / TMath::Power(-c4, 0.75);
@@ -657,30 +661,30 @@ void AliForwardFlowTaskQC::Terminate(Option_t */*option*/)
         Cov2p2np = CovXY(W2pW4pavg2pavg4p, W2pW4p, Avg2p*Avg4p, W2p, W4p);
  
         // Numbers on the side reference term number in paper (cite needed) loosely
-/*1*/ vTwo4diffErr =  TMath::Power(2.*Avg2*Avg2*Avg2p - 3.*Avg2*Avg4p + 2.*Avg4*Avg2p, 2) 
-                      * TMath::Power(sqrtW2sq, 2) * sAvg2sq / (W2*W2);
-/*2*/ vTwo4diffErr += 9. * TMath::Power(2.*Avg2*Avg2p - Avg4p, 2) * TMath::Power(sqrtW4sq, 2)
-                      * sAvg4sq / (16. * W4*W4);
-/*3*/ vTwo4diffErr += 4. * Avg2*Avg2 * TMath::Power(2.*Avg2*Avg2 - Avg4, 2) * TMath::Power(sqrtW2psq, 2)
-                      * sAvg2psq / (W2p*W2p);
-/*4*/ vTwo4diffErr += TMath::Power(2.*Avg2*Avg2 - Avg4, 2) * TMath::Power(sqrtW4psq, 2) * sAvg4psq
-                      / (W4p*W4p);
-/*5*/ vTwo4diffErr -= 1.5 * (2.*Avg2*Avg2p - Avg4p) * (2.*Avg2*Avg2*Avg2p - 3.*Avg2*Avg4p + 2.*Avg4*Avg2p)
-                      * W2W4 * Cov24 / (W2*W4);
-/*6*/ vTwo4diffErr -= 4. * Avg2 * (2.*Avg2*Avg2 - Avg4) 
-                      * (2.*Avg2*Avg2*Avg2p - 3.*Avg2*Avg4p + 2.*Avg4*Avg2p)
-                      * W2W2p * Cov22p / (W2 * W2p);
-/*7*/ vTwo4diffErr += 2. * (2.*Avg2*Avg2 - Avg4)
-                      * (2.*Avg2*Avg2*Avg2p - 3.*Avg2*Avg4p + 2.*Avg4*Avg2p)
-                      * W2W4p * Cov24p / (W2 * W4p);
-/*8*/ vTwo4diffErr += 3.*Avg2*(2.*Avg2*Avg2 - Avg4)*(2.*Avg2*Avg2p - Avg4p)
-                      * W4W2p * Cov42p / (W4*W2p);
-/*9*/ vTwo4diffErr -= 1.5 * (2.*Avg2*Avg2 - Avg4)*(2.*Avg2*Avg2p - Avg4p)
-                      * W4W4p * Cov44p / (W4 * W4p);
-/*10*/vTwo4diffErr -= 4.*Avg2*TMath::Power(2.*Avg2*Avg2 - Avg4, 2)
-                      * W2pW4p * Cov2p2np / (W2p * W4p);
-/*11*/vTwo4diffErr /= TMath::Power(2.*Avg2*Avg2 - Avg4, 3.5);
-      vTwo4diffErr = TMath::Sqrt(vTwo4diffErr);
+       /*1*/ vTwo4diffErr =  TMath::Power(2.*Avg2*Avg2*Avg2p - 3.*Avg2*Avg4p + 2.*Avg4*Avg2p, 2) 
+         * TMath::Power(sqrtW2sq, 2) * sAvg2sq / (W2*W2);
+       /*2*/ vTwo4diffErr += 9. * TMath::Power(2.*Avg2*Avg2p - Avg4p, 2) * TMath::Power(sqrtW4sq, 2)
+         * sAvg4sq / (16. * W4*W4);
+       /*3*/ vTwo4diffErr += 4. * Avg2*Avg2 * TMath::Power(2.*Avg2*Avg2 - Avg4, 2) * TMath::Power(sqrtW2psq, 2)
+         * sAvg2psq / (W2p*W2p);
+       /*4*/ vTwo4diffErr += TMath::Power(2.*Avg2*Avg2 - Avg4, 2) * TMath::Power(sqrtW4psq, 2) * sAvg4psq
+         / (W4p*W4p);
+       /*5*/ vTwo4diffErr -= 1.5 * (2.*Avg2*Avg2p - Avg4p) * (2.*Avg2*Avg2*Avg2p - 3.*Avg2*Avg4p + 2.*Avg4*Avg2p)
+         * W2W4 * Cov24 / (W2*W4);
+       /*6*/ vTwo4diffErr -= 4. * Avg2 * (2.*Avg2*Avg2 - Avg4) 
+         * (2.*Avg2*Avg2*Avg2p - 3.*Avg2*Avg4p + 2.*Avg4*Avg2p)
+         * W2W2p * Cov22p / (W2 * W2p);
+       /*7*/ vTwo4diffErr += 2. * (2.*Avg2*Avg2 - Avg4)
+         * (2.*Avg2*Avg2*Avg2p - 3.*Avg2*Avg4p + 2.*Avg4*Avg2p)
+         * W2W4p * Cov24p / (W2 * W4p);
+       /*8*/ vTwo4diffErr += 3.*Avg2*(2.*Avg2*Avg2 - Avg4)*(2.*Avg2*Avg2p - Avg4p)
+         * W4W2p * Cov42p / (W4*W2p);
+       /*9*/ vTwo4diffErr -= 1.5 * (2.*Avg2*Avg2 - Avg4)*(2.*Avg2*Avg2p - Avg4p)
+         * W4W4p * Cov44p / (W4 * W4p);
+       /*10*/vTwo4diffErr -= 4.*Avg2*TMath::Power(2.*Avg2*Avg2 - Avg4, 2)
+         * W2pW4p * Cov2p2np / (W2p * W4p);
+       /*11*/vTwo4diffErr /= TMath::Power(2.*Avg2*Avg2 - Avg4, 3.5);
+       vTwo4diffErr = TMath::Sqrt(vTwo4diffErr);
 
         cumulant4Hist->SetBinError(eta, vTwo4diffErr);
       } // End of eta loop
@@ -719,9 +723,9 @@ void AliForwardFlowTaskQC::ProcessPrimary()
   if (fAOD) {
     if (!common->LoopAODmc(fAOD)) return;
     if (!common->LoopAODtrrefHits(fAOD)) return;
-//    if (!common->LoopMCaddptFlow(fAOD)) return;
-//    if (!common->LoopMCaddpdgFlow(fAOD)) return;
-//    if (!common->LoopMCaddetaFlow(fAOD)) return;
+    //    if (!common->LoopMCaddptFlow(fAOD)) return;
+    //    if (!common->LoopMCaddpdgFlow(fAOD)) return;
+    //    if (!common->LoopMCaddetaFlow(fAOD)) return;
   }
 
   // Run analysis on MC truth
@@ -738,7 +742,8 @@ void AliForwardFlowTaskQC::ProcessPrimary()
 
 }
 //_____________________________________________________________________
-Double_t AliForwardFlowTaskQC::VarSQ(Double_t wxx2, Double_t x, Double_t wx, Double_t wxx, Double_t sqrtwx2)
+Double_t AliForwardFlowTaskQC::VarSQ(Double_t wxx2, Double_t x, Double_t wx, 
+                                    Double_t wxx, Double_t sqrtwx2)
 {
   //
   // Small function to compute the variance squared - used by Terminte()
@@ -752,7 +757,8 @@ Double_t AliForwardFlowTaskQC::VarSQ(Double_t wxx2, Double_t x, Double_t wx, Dou
   return sx;
 }
 //_____________________________________________________________________
-Double_t AliForwardFlowTaskQC::CovXY(Double_t wxwyxy, Double_t wxwy, Double_t xy, Double_t wx, Double_t wy)
+Double_t AliForwardFlowTaskQC::CovXY(Double_t wxwyxy, Double_t wxwy, 
+                                    Double_t xy, Double_t wx, Double_t wy)
 {
   //
   // Small function to compute the covariance between two numbers
index 882057f86ed47d6853c19daea5db0042fcc37b94..70532aca88a5c27c9bce7a6c00541701f72f9ba5 100644 (file)
@@ -17,6 +17,8 @@
  *
  * @ingroup pwg2_forward_tasks
  *
+ * @todo Add centrality stuff
+ *
  */
 class AliForwardFlowTaskQC : public AliAnalysisTaskSE
 {
@@ -24,23 +26,23 @@ public:
   /** 
    * Constructor 
    */
-   AliForwardFlowTaskQC();
+  AliForwardFlowTaskQC();
   /** 
    * Constructor
    * 
    * @param name Name of task 
    */
-   AliForwardFlowTaskQC(const char* name);
+  AliForwardFlowTaskQC(const char* name);
   /**
    * Destructor
    */
-   virtual ~AliForwardFlowTaskQC() {;}
+  virtual ~AliForwardFlowTaskQC() {;}
   /** 
    * Copy constructor 
    * 
    * @param o Object to copy from 
    */
-   AliForwardFlowTaskQC(const AliForwardFlowTaskQC& o);
+  AliForwardFlowTaskQC(const AliForwardFlowTaskQC& o);
   /** 
    * Assignment operator 
    * 
@@ -48,8 +50,8 @@ public:
    * 
    * @return Reference to this object 
    */
-   AliForwardFlowTaskQC& operator=(const AliForwardFlowTaskQC&) { return *this; }
-   /** 
+  AliForwardFlowTaskQC& operator=(const AliForwardFlowTaskQC&) { return *this; }
+  /** 
    * @{ 
    * @name Interface methods 
    */
@@ -57,81 +59,87 @@ public:
    * Create output objects 
    * 
    */
-   virtual void CreateOutputObjects();
+  virtual void CreateOutputObjects();
   /**
    * Initialize the task
    *
    */
-   virtual void Init() {}
+  virtual void Init() {}
   /** 
    * Process each event 
    *
    * @param option Not used
    */  
-   virtual void UserExec(Option_t *option);
-  /**
-   * Calculate Q cumulant
-   * Parameters:
-   *  type: Determines which histograms should be used
-   *        - "" = data histograms
-   *        - "TrRef" = track reference histograms
-   *        - "MC" = MC truth histograms
-   *  harmonic: Which harmonic to calculate
-   */
-   void CumulantsMethod(TString type, Int_t harmonic);
+  virtual void UserExec(Option_t *option);
   /** 
    * End of job
    * 
    * @param option Not used 
    */
-   virtual void Terminate(Option_t *option);
-  /*
-   * if MC information is available do analysis on Monte Carlo truth and track references
-   *
-   */
-   void ProcessPrimary();
+  virtual void Terminate(Option_t *option);
   /*
    * Returns the outputlist
    *
    */
-   TList* GetOutputList() { return fOutputList; }
+  TList* GetOutputList() { return fOutputList; }
   /* 
    * Set Number of eta bins to be used in flow analysis
    *
    */
-   void SetUseNEtaBins(Int_t nbins) { fEtaBins = nbins; }
+  void SetUseNEtaBins(Int_t nbins) { fEtaBins = nbins; }
   /* 
    * Set which harmonics to calculate
    * v_{1} to v_{4} is available and calculated as default
    *
    */
-   void SetDoHarmonics(Bool_t v1 = kTRUE, Bool_t v2 = kTRUE, Bool_t v3 = kTRUE, Bool_t v4 = kTRUE)
-    { fv[1] = v1; fv[2] = v2; fv[3] = v3; fv[4] = v4; }
+  void SetDoHarmonics(Bool_t v1 = kTRUE, Bool_t v2 = kTRUE, 
+                     Bool_t v3 = kTRUE, Bool_t v4 = kTRUE) { 
+    fv[1] = v1; fv[2] = v2; fv[3] = v3; fv[4] = v4; }
   /** 
    * @} 
    */
-   virtual void SetDebugLevel(Int_t level) { fDebug = level; }
+  virtual void SetDebugLevel(Int_t level) { fDebug = level; }
   
-private:
+protected:
   /*
-   * Caclulate the variance of x squared - used to finalize calculations in Terminate()
+   * if MC information is available do analysis on Monte Carlo truth
+   * and track references
    *
    */
-   Double_t VarSQ(Double_t wxx2, Double_t x, Double_t wx, Double_t wxx, Double_t sqrtwx2);
-  /*
-   * Caclulate the covariance between x and y - used to finalize calculations in Terminate()
+  void ProcessPrimary();
+  /**
+   * Calculate Q cumulant
+   * Parameters:
+   *  type: Determines which histograms should be used
+   *        - "" = data histograms
+   *        - "TrRef" = track reference histograms
+   *        - "MC" = MC truth histograms
+   *  harmonic: Which harmonic to calculate
+   */
+  void CumulantsMethod(TString type, Int_t harmonic);
+  /**
+   * Caclulate the variance of x squared - used to finalize
+   * calculations in Terminate()
+   *
+   */
+  Double_t VarSQ(Double_t wxx2, Double_t x, Double_t wx, 
+                Double_t wxx, Double_t sqrtwx2);
+  /**
+   * Caclulate the covariance between x and y - used to finalize
+   * calculations in Terminate()
    *
    */
-   Double_t CovXY(Double_t wxwyxy, Double_t wxwy, Double_t XY, Double_t wx, Double_t wy);
+  Double_t CovXY(Double_t wxwyxy, Double_t wxwy, Double_t XY, 
+                Double_t wx, Double_t wy);
 
-   Int_t          fDebug;        //  Debug flag
-   TList*         fOutputList;   //  Output list
-   AliAODEvent*   fAOD;          //  AOD event
-   Bool_t         fMC;           //  Is MC flags
-   Int_t          fEtaBins;      //  Number of eta bins in histograms
-   Bool_t         fv[5];         //  Calculate v_{n} flag
+  Int_t          fDebug;        //  Debug flag
+  TList*         fOutputList;   //  Output list
+  AliAODEvent*   fAOD;          //  AOD event
+  Bool_t         fMC;           //  Is MC flags
+  Int_t          fEtaBins;      //  Number of eta bins in histograms
+  Bool_t         fv[5];         //  Calculate v_{n} flag
    
-   ClassDef(AliForwardFlowTaskQC, 1); // Analysis task for FMD analysis
+  ClassDef(AliForwardFlowTaskQC, 1); // Analysis task for FMD analysis
 };
  
 #endif
index 12ec76b7aec6efaa4d32d383f7da99797bd4fc4c..29a9e7d6c3751277bfc8559472bf126f68476c32 100644 (file)
@@ -32,7 +32,7 @@
 //====================================================================
 AliForwardMultiplicityBase::AliForwardMultiplicityBase(const char* name) 
   : AliAnalysisTaskSE(name), 
-    fEnableLowFlux(true), 
+    fEnableLowFlux(false), 
     fFirstEvent(true),
     fCorrManager(0)
 {
index 1c7df91e3590ec2c5de205bfc66b45d326c9b108..b2cbd2b20b812d55c5f1a664f506b3277e2694ee 100644 (file)
@@ -181,6 +181,7 @@ AliForwardMultiplicityTask::UserCreateOutputObjects()
   // 
   //
   fList = new TList;
+  fList->SetOwner();
 
   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
   AliAODHandler*      ah = 
diff --git a/PWG2/FORWARD/analysis2/ForwardAODConfig.C b/PWG2/FORWARD/analysis2/ForwardAODConfig.C
new file mode 100644 (file)
index 0000000..0c6585e
--- /dev/null
@@ -0,0 +1,83 @@
+/**
+ * Configuration script for forward multiplicity task.  
+ *
+ * You can copy this to your working directory or to some other
+ * directory in up-front the ROOT macro path, and edit it to suit your
+ * needs.
+ * 
+ */
+void
+ForwardAODConfig(AliForwardMultiplicityBase* task)
+{
+  if (!task) return;
+
+  Info("ForwardAODConfig", "Setting up task %s (%p)", task->GetName(), task);
+  // Whether to enable low flux specific code 
+  task->SetEnableLowFlux(kFALSE);
+  // Set the number of SPD tracklets for which we consider the event a
+  // low flux event
+  task->GetEventInspector().SetLowFluxCut(1000); 
+  // Set the maximum error on v_z [cm]
+  task->GetEventInspector().SetMaxVzErr(0.2);
+  // Set the eta axis to use - note, this overrides whatever is used
+  // by the rest of the algorithms - but only for the energy fitter
+  // algorithm. 
+  task->GetEnergyFitter().SetEtaAxis(200, -4, 6);
+  // Set maximum energy loss to consider 
+  task->GetEnergyFitter().SetMaxE(10); 
+  // Set number of energy loss bins 
+  task->GetEnergyFitter().SetNEbins(300);
+  // Set whether to use increasing bin sizes 
+  task->GetEnergyFitter().SetUseIncreasingBins(true);
+  // Set whether to do fit the energy distributions 
+  task->GetEnergyFitter().SetDoFits(kFALSE);
+  // Set whether to make the correction object 
+  task->GetEnergyFitter().SetDoMakeObject(kFALSE);
+  // Set the low cut used for energy
+  task->GetEnergyFitter().SetLowCut(0.4);
+  // Set the number of bins to subtract from maximum of distributions
+  // to get the lower bound of the fit range
+  task->GetEnergyFitter().SetFitRangeBinWidth(4);
+  // Set the maximum number of landaus to try to fit (max 5)
+  task->GetEnergyFitter().SetNParticles(5);
+  // Set the minimum number of entries in the distribution before
+  // trying to fit to the data
+  task->GetEnergyFitter().SetMinEntries(1000);
+  // Set the low cut used for sharing - overrides settings in eloss fits
+  task->GetSharingFilter().SetLowCut(0.3);
+  // Set the number of xi's (width of landau peak) to stop at 
+  task->GetSharingFilter().SetNXi(1);
+  // Set the maximum number of particle to try to reconstruct 
+  task->GetDensityCalculator().SetMaxParticles(3);
+  // Wet whether to use poisson statistics to estimate N_ch
+  task->GetDensityCalculator().SetUsePoisson(false);
+  // Set the lower multiplicity cut.  Overrides setting in energy loss fits.
+  task->GetDensityCalculator().SetMultCut(0.3); //was 0.3
+  // Whether to use the secondary map correction
+  task->GetCorrections().SetUseSecondaryMap(true);
+  // Whether to use the vertex bias correction
+  task->GetCorrections().SetUseVertexBias(false);
+  // Whether to use the vertex bias correction
+  task->GetCorrections().SetUseAcceptance(true);
+  // Whether to use the merging efficiency correction 
+  task->GetCorrections().SetUseMergingEfficiency(false);
+  // Set the number of extra bins (beyond the secondary map border) 
+  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.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 
+  // task->GetEventInspector().SetDebug(4);
+  // --- Set limits on fits the energy -------------------------------
+  // Maximum relative error on parameters 
+  AliFMDCorrELossFit::ELossFit::fgMaxRelError = .12;
+  // Least weight to use 
+  AliFMDCorrELossFit::ELossFit::fgLeastWeight = 1e-5;
+  // Maximum value of reduced chi^2 
+  AliFMDCorrELossFit::ELossFit::fgMaxChi2nu   = 20;
+}
+//
+// EOF
+//
index 95a15fc513efb0074a564eeb78e85df5494c170d..fe9960b7f7ae06a4a5db41026d788a222957adde 100644 (file)
@@ -89,14 +89,14 @@ void MakeAOD(const char* esddir,
   // --- AOD output handler ------------------------------------------
   AliAODHandler* aodHandler   = new AliAODHandler();
   mgr->SetOutputEventHandler(aodHandler);
-  aodHandler->SetOutputFileName("AliAODs.root");
+  aodHandler->SetOutputFileName("AliAOD.root");
 
   // --- Add tasks ---------------------------------------------------
   // Physics selection 
   gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
   AddTaskPhysicsSelection(mc, kTRUE, kTRUE);
 
-#if 1
+#if 0
   // Centrality 
   if (!proof) {
     gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/Compile.C");
diff --git a/PWG2/FORWARD/analysis2/MakeFlow.C b/PWG2/FORWARD/analysis2/MakeFlow.C
new file mode 100644 (file)
index 0000000..1f401fb
--- /dev/null
@@ -0,0 +1,71 @@
+/**
+ * Script to analyse AOD input for flow 
+ * 
+ * @par Inputs: 
+ * - 
+ * 
+ * @par Outputs: 
+ * - 
+ * 
+ */
+void MakeFlow(TString path      = "", 
+             Bool_t  recursive = true,
+             Int_t   nevents   = 100, 
+             TString type      = "", 
+             Int_t   etabins)
+{
+  Bool_t proof = false;
+
+  // --- Load libs -------------------------------------------------
+  gROOT->Macro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadLibs.C");
+
+  // --- Check for proof mode, and possibly upload pars --------------
+  if (proof> 0) { 
+    gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadPars.C");
+    if (!LoadPars(proof)) { 
+      Error("MakeAOD", "Failed to load PARs");
+      return;
+    }
+  }
+
+  // --- Add to chain either ESD or AOD ----------------------------
+  gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/MakeChain.C");
+  TChain* chain = MakeChain("AOD", path, recursive);
+  // If 0 or less events is select, choose all 
+  if (nevents <= 0) nevents = chain->GetEntries();
+
+  // --- Initiate the event handlers -------------------------------
+  AliAnalysisManager *mgr  = new AliAnalysisManager("Forward Flow", 
+                                                   "Flow in forward region");
+  mgr->SetUseProgressBar(kTRUE);
+
+  // --- AOD input handler -------------------------------------------
+  AliAODInputHandler *aodInputHandler = new AliAODInputHandler();
+  mgr->SetInputEventHandler(aodInputHandler);      
+
+  // --- Create output handler ---------------------------------------
+  AliAODHandler* aodOut = new AliAODHandler();
+  mgr->SetOutputEventHandler(aodOut);
+  aodOut->SetOutputFileName("AliAODs.Flow.root");
+
+  // --- Add the tasks ---------------------------------------------
+  gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/AddTaskForwardFlow.C");
+  AddTaskForwardFlow(type.Data(), etabins);
+
+  // --- Run the analysis --------------------------------------------
+  TStopwatch t;
+  if (!mgr->InitAnalysis()) {
+    Error("MakeFlow", "Failed to initialize analysis train!");
+    return;
+  }
+  mgr->PrintStatus();
+
+  t.Start();
+  if (nevents == 0) mgr->StartAnalysis("local", chain);
+  if (nevents != 0) mgr->StartAnalysis("local", chain, nevents);
+  t.Stop();
+  t.Print();
+}
+//
+// EOF
+//
diff --git a/PWG2/FORWARD/analysis2/RunManagerFlow.C b/PWG2/FORWARD/analysis2/RunManagerFlow.C
deleted file mode 100644 (file)
index ec42540..0000000
+++ /dev/null
@@ -1,76 +0,0 @@
-void RunManagerFlow(TString path = "", TString file ="AliAODs.root", Int_t nevents = 100, TString type="", Int_t etabins)
-{
-
-  // --- Load libs -------------------------------------------------
-  
-  gSystem->Load("libVMC");
-  gSystem->Load("libGeom");
-  gSystem->Load("libXMLIO");
-  gSystem->Load("libTree");
-  
-  gSystem->Load("libSTEERBase");
-  
-  gSystem->Load("libESD") ;
-  gSystem->Load("libAOD") ;
-  gSystem->Load("libANALYSIS");
-  gSystem->Load("libANALYSISalice");
-  
-  gSystem->Load("libPhysics");
-  gSystem->Load("libPWG0base");
-  gSystem->Load("libPWG0dep");
-  gSystem->Load("libPWG2forward");
-  gSystem->Load("libPWG2forward2");
-
-  // --- Add to chain either ESD or AOD ----------------------------
-
-  TChain* chain = 0;
-
-  TChain* chainAOD = new TChain("aodTree");
-  if (path.Sizeof() == 1) { // Is it a single AOD file?
-    chainAOD->Add(Form("%s", file.Data()));
-    chain = chainAOD;
-  }
-  if (path.Sizeof() > 1) {  // More AOD files, assumes AliAODs1.root numbering
-    for (Int_t i=1;i<=9;i++) {
-      if (Form("%s/%s%d", path.Data(), file.Data(), i))
-        chainAOD->Add(Form("%s/AliAODs%d.root", path.Data(), i));
-    }
-    chain = chainAOD;
-  }
-
-  // --- Initiate the event handlers -------------------------------
-
-  AliAnalysisManager *mgr  = new AliAnalysisManager("Analysis Train", "A test setup for the analysis train");
-  mgr->SetUseProgressBar(kTRUE);
-
-  // AOD input handler
-  AliAODInputHandler* aodHandler = new AliAODInputHandler();
-  mgr->SetInputEventHandler(aodHandler);
-
-  // Create output handler
-  AliAODHandler* aodOut = new AliAODHandler();
-  mgr->SetOutputEventHandler(aodOut);
-  aodOut->SetOutputFileName("AliAODs.Flow.root");
-
-  // --- Add the tasks ---------------------------------------------
-  
-  gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/AddTaskForwardFlow.C");
-  AddTaskForwardFlow(type.Data(), etabins);
-
-  // --- Run the analysis ------------------------------------------
-
-  TStopwatch t;
-  if (mgr->InitAnalysis()) {
-    mgr->PrintStatus();
-    t.Start();
-    if (nevents == 0) mgr->StartAnalysis("local", chain);
-    if (nevents != 0) mgr->StartAnalysis("local", chain, nevents);
-    t.Stop();
-    t.Print();
-  }  
-
-  // ---------------------------------------------------------------
-}
-//
-// EOF
-//