Add class AliHFEdebugTreeTask

[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliNormalizationCounter.cxx

/**************************************************************************\r
 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *\r
 *                                                                        *\r
 * Author: The ALICE Off-line Project.                                    *\r
 * Contributors are mentioned in the code where appropriate.              *\r
 *                                                                        *\r
 * Permission to use, copy, modify and distribute this software and its   *\r
 * documentation strictly for non-commercial purposes is hereby granted   *\r
 * without fee, provided that the above copyright notice appears in all   *\r
 * copies and that both the copyright notice and this permission notice   *\r
 * appear in the supporting documentation. The authors make no claims     *\r
 * about the suitability of this software for any purpose. It is          *\r
 * provided "as is" without express or implied warranty.                  *\r
 **************************************************************************/\r
\r
/* $Id$ */\r
\r
//*************************************************************************\r
// Class AliNormalizationCounter\r
// Class to store the informations relevant for the normalization in the \r
// barrel for each run\r
// Authors: G. Ortona, ortona@to.infn.it\r
// D. Caffarri, davide.caffarri@pd.to.infn.it\r
// with many thanks to P. Pillot\r
/////////////////////////////////////////////////////////////\r
\r
#include "AliLog.h"\r
#include "AliNormalizationCounter.h"\r
#include <AliESDEvent.h>\r
#include <AliESDtrack.h>\r
#include <AliAODEvent.h>\r
#include <AliVParticle.h>\r
#include <AliTriggerAnalysis.h>\r
#include <TH1F.h>\r
#include <TH2F.h>\r
#include <TList.h>\r
#include <TString.h>\r
#include <TCanvas.h>\r
#include <AliPhysicsSelection.h>\r
#include <AliMultiplicity.h>\r
\r
ClassImp(AliNormalizationCounter)\r
\r
//____________________________________________\r
AliNormalizationCounter::AliNormalizationCounter(): \r
TNamed(),\r
fCounters(),\r
fESD(kFALSE),\r
fMultiplicity(kFALSE),\r
fHistTrackFilterEvMult(0),\r
fHistTrackAnaEvMult(0),\r
fHistTrackFilterSpdMult(0),\r
fHistTrackAnaSpdMult(0)\r
{\r
  // empty constructor\r
}\r
\r
//__________________________________________________				\r
AliNormalizationCounter::AliNormalizationCounter(const char *name): \r
TNamed(name,name),\r
fCounters(name),\r
fESD(kFALSE),\r
fMultiplicity(kFALSE),\r
fHistTrackFilterEvMult(0),\r
fHistTrackAnaEvMult(0),\r
fHistTrackFilterSpdMult(0),\r
fHistTrackAnaSpdMult(0)\r
{\r
  ;\r
}\r
\r
//______________________________________________\r
AliNormalizationCounter::~AliNormalizationCounter()\r
{\r
  //destructor\r
  if(fHistTrackFilterEvMult){\r
    delete fHistTrackFilterEvMult;\r
    fHistTrackFilterEvMult =0;\r
  }\r
  if(fHistTrackAnaEvMult){\r
    delete fHistTrackAnaEvMult;\r
    fHistTrackAnaEvMult=0;\r
  }\r
  if(fHistTrackFilterSpdMult){\r
    delete fHistTrackFilterSpdMult;\r
    fHistTrackFilterSpdMult=0;\r
  }\r
  if(fHistTrackAnaSpdMult){\r
    delete fHistTrackAnaSpdMult;\r
    fHistTrackAnaSpdMult=0;\r
  }\r
}\r
\r
//______________________________________________\r
void AliNormalizationCounter::Init()\r
{\r
  //variables initialization\r
  fCounters.AddRubric("Event","triggered/V0AND/PileUp/PbPbC0SMH-B-NOPF-ALLNOTRD/Candles0.3/PrimaryV/countForNorm/noPrimaryV/zvtxGT10/!V0A&Candle03/!V0A&PrimaryV/Candid(Filter)/Candid(Analysis)/NCandid(Filter)/NCandid(Analysis)");\r
  if(fMultiplicity)  fCounters.AddRubric("Multiplicity", 5000);\r
  fCounters.AddRubric("Run", 1000000);\r
  fCounters.Init();\r
  fHistTrackFilterEvMult=new TH2F("FiltCandidvsTracksinEv","FiltCandidvsTracksinEv",10000,-0.5,9999.5,200,-0.5,199.5);\r
  fHistTrackFilterEvMult->GetYaxis()->SetTitle("NCandidates");\r
  fHistTrackFilterEvMult->GetXaxis()->SetTitle("NTracksinEvent");\r
  fHistTrackAnaEvMult=new TH2F("AnaCandidvsTracksinEv","AnaCandidvsTracksinEv",10000,-0.5,9999.5,100,-0.5,99.5);\r
  fHistTrackAnaEvMult->GetYaxis()->SetTitle("NCandidates");\r
  fHistTrackAnaEvMult->GetXaxis()->SetTitle("NTracksinEvent");\r
  fHistTrackFilterSpdMult=new TH2F("FilterCandidvsSpdMult","FilterCandidvsSpdMult",5000,-0.5,4999.5,200,-0.5,199.5);\r
  fHistTrackFilterSpdMult->GetYaxis()->SetTitle("NCandidates");\r
  fHistTrackFilterSpdMult->GetXaxis()->SetTitle("NSPDTracklets");\r
  fHistTrackAnaSpdMult=new TH2F("AnaCandidvsSpdMult","AnaCandidvsSpdMult",5000,-0.5,4999.5,100,-0.5,99.5);\r
  fHistTrackAnaSpdMult->GetYaxis()->SetTitle("NCandidates");\r
  fHistTrackAnaSpdMult->GetXaxis()->SetTitle("NSPDTracklets");\r
}\r
\r
//______________________________________________\r
Long64_t AliNormalizationCounter::Merge(TCollection* list){\r
  if (!list) return 0;\r
  if (list->IsEmpty()) return 0;//(Long64_t)fCounters.Merge(list);\r
\r
  TIter next(list);\r
  const TObject* obj = 0x0;\r
  while ((obj = next())) {\r
    \r
    // check that "obj" is an object of the class AliNormalizationCounter\r
    const AliNormalizationCounter* counter = dynamic_cast<const AliNormalizationCounter*>(obj);\r
    if (!counter) {\r
      AliError(Form("object named %s is not AliNormalizationCounter! Skipping it.", counter->GetName()));\r
      continue;\r
    }\r
\r
    Add(counter);\r
\r
  }\r
  \r
  return (Long64_t)1;//(Long64_t)fCounters->GetEntries();\r
}\r
//_______________________________________\r
void AliNormalizationCounter::Add(const AliNormalizationCounter *norm){\r
  fCounters.Add(&(norm->fCounters));\r
  fHistTrackFilterEvMult->Add(norm->fHistTrackFilterEvMult);\r
  fHistTrackAnaEvMult->Add(norm->fHistTrackAnaEvMult);\r
  fHistTrackFilterSpdMult->Add(norm->fHistTrackFilterSpdMult);\r
  fHistTrackAnaSpdMult->Add(norm->fHistTrackAnaSpdMult);\r
}\r
//_______________________________________\r
/*\r
Stores the variables used for normalization as function of run number\r
returns kTRUE if the event is to be counted for normalization\r
(pass event selection cuts OR has no primary vertex)\r
 */\r
void AliNormalizationCounter::StoreEvent(AliVEvent *event,AliRDHFCuts *rdCut,Bool_t mc){\r
  //\r
\r
  Bool_t isEventSelected = rdCut->IsEventSelected(event);\r
\r
  // events not passing physics selection. do nothing\r
  if(rdCut->IsEventRejectedDuePhysicsSelection()) return;\r
\r
  Bool_t v0A=kFALSE; \r
  Bool_t v0B=kFALSE;\r
  Bool_t flag03=kFALSE;\r
  Bool_t flagPV=kFALSE;\r
\r
  //Run Number\r
  Int_t runNumber = event->GetRunNumber();\r
 \r
  // Evaluate the multiplicity\r
  Int_t multiplicity = Multiplicity(event);\r
\r
  //Find CINT1B\r
  AliESDEvent *eventESD = (AliESDEvent*)event;\r
  if(!eventESD){AliError("ESD event not available");return;}\r
  if(mc&&event->GetEventType() != 0)return;\r
  //event must be either physics or MC\r
  if(!(event->GetEventType() == 7||event->GetEventType() == 0))return;\r
  \r
  if(fMultiplicity) \r
    fCounters.Count(Form("Event:triggered/Run:%d/Multiplicity:%d",runNumber,multiplicity));\r
  else \r
    fCounters.Count(Form("Event:triggered/Run:%d",runNumber));\r
\r
  //Find V0AND\r
  AliTriggerAnalysis trAn; /// Trigger Analysis\r
  v0B = trAn.IsOfflineTriggerFired(eventESD , AliTriggerAnalysis::kV0C);\r
  v0A = trAn.IsOfflineTriggerFired(eventESD , AliTriggerAnalysis::kV0A);\r
  if(v0A&&v0B){\r
    if(fMultiplicity) \r
      fCounters.Count(Form("Event:V0AND/Run:%d/Multiplicity:%d",runNumber,multiplicity));\r
    else\r
      fCounters.Count(Form("Event:V0AND/Run:%d",runNumber));\r
  }\r
  \r
  //FindPrimary vertex  \r
  // AliVVertex *vtrc =  (AliVVertex*)event->GetPrimaryVertex();\r
  // if(vtrc && vtrc->GetNContributors()>0){\r
  //   fCounters.Count(Form("Event:PrimaryV/Run:%d",runNumber));\r
  //   flagPV=kTRUE;\r
  // }\r
\r
  //trigger\r
  AliAODEvent *eventAOD = (AliAODEvent*)event;\r
  TString trigclass=eventAOD->GetFiredTriggerClasses();\r
  if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD")||trigclass.Contains("C0SMH-B-NOPF-ALL")){\r
    if(fMultiplicity) \r
      fCounters.Count(Form("Event:PbPbC0SMH-B-NOPF-ALLNOTRD/Run:%d/Multiplicity:%d",runNumber,multiplicity));\r
    else \r
      fCounters.Count(Form("Event:PbPbC0SMH-B-NOPF-ALLNOTRD/Run:%d",runNumber));\r
  }\r
\r
  //FindPrimary vertex  \r
  if(isEventSelected){\r
    if(fMultiplicity) \r
      fCounters.Count(Form("Event:PrimaryV/Run:%d/Multiplicity:%d",runNumber,multiplicity));\r
    else\r
      fCounters.Count(Form("Event:PrimaryV/Run:%d",runNumber));\r
    flagPV=kTRUE;\r
  }else{\r
    if(rdCut->GetWhyRejection()==0){\r
      if(fMultiplicity) \r
	fCounters.Count(Form("Event:noPrimaryV/Run:%d/Multiplicity:%d",runNumber,multiplicity));\r
      else\r
	fCounters.Count(Form("Event:noPrimaryV/Run:%d",runNumber));\r
    }\r
    //find good vtx outside range\r
    if(rdCut->GetWhyRejection()==6){\r
      if(fMultiplicity) {\r
	fCounters.Count(Form("Event:zvtxGT10/Run:%d/Multiplicity:%d",runNumber,multiplicity));\r
	fCounters.Count(Form("Event:PrimaryV/Run:%d/Multiplicity:%d",runNumber,multiplicity));\r
      } else {\r
	fCounters.Count(Form("Event:zvtxGT10/Run:%d",runNumber));\r
	fCounters.Count(Form("Event:PrimaryV/Run:%d",runNumber));\r
      }\r
      flagPV=kTRUE;\r
    }\r
    if(rdCut->GetWhyRejection()==1){\r
      if(fMultiplicity) \r
	fCounters.Count(Form("Event:PileUp/Run:%d/Multiplicity:%d",runNumber,multiplicity));\r
      else\r
	fCounters.Count(Form("Event:PileUp/Run:%d",runNumber));\r
    }\r
  }\r
  //to be counted for normalization\r
  if(rdCut->CountEventForNormalization()){\r
    if(fMultiplicity) \r
      fCounters.Count(Form("Event:countForNorm/Run:%d/Multiplicity:%d",runNumber,multiplicity));\r
    else \r
      fCounters.Count(Form("Event:countForNorm/Run:%d",runNumber));\r
  }\r
\r
\r
  //Find Candle\r
  Int_t trkEntries = (Int_t)event->GetNumberOfTracks();\r
  for(Int_t i=0;i<trkEntries&&!flag03;i++){\r
    AliAODTrack *track=(AliAODTrack*)event->GetTrack(i);\r
    if((track->Pt()>0.3)&&(!flag03)){\r
      if(fMultiplicity) \r
	fCounters.Count(Form("Event:Candles0.3/Run:%d/Multiplicity:%d",runNumber,multiplicity));\r
      else\r
	fCounters.Count(Form("Event:Candles0.3/Run:%d",runNumber));\r
      flag03=kTRUE;\r
      break;\r
    }\r
  }\r
  \r
  if(!(v0A&&v0B)&&(flag03)){ \r
    if(fMultiplicity) \r
      fCounters.Count(Form("Event:!V0A&Candle03/Run:%d/Multiplicity:%d",runNumber,multiplicity));\r
    else \r
      fCounters.Count(Form("Event:!V0A&Candle03/Run:%d",runNumber));\r
  }\r
  if(!(v0A&&v0B)&&flagPV){\r
    if(fMultiplicity) \r
      fCounters.Count(Form("Event:!V0A&PrimaryV/Run:%d/Multiplicity:%d",runNumber,multiplicity));\r
    else\r
      fCounters.Count(Form("Event:!V0A&PrimaryV/Run:%d",runNumber));\r
  }\r
  \r
  return;\r
}\r
//_____________________________________________________________________\r
void AliNormalizationCounter::StoreCandidates(AliVEvent *event,Int_t nCand,Bool_t flagFilter){\r
  \r
  Int_t ntracks=event->GetNumberOfTracks();\r
  if(flagFilter)fHistTrackFilterEvMult->Fill(ntracks,nCand);\r
  else fHistTrackAnaEvMult->Fill(ntracks,nCand);\r
  Int_t nSPD=0;\r
  if(fESD){\r
    AliESDEvent *ESDevent=(AliESDEvent*)event;\r
    const AliMultiplicity *alimult = ESDevent->GetMultiplicity();\r
    nSPD = alimult->GetNumberOfTracklets();\r
\r
  }else{\r
    AliAODEvent *aodEvent =(AliAODEvent*)event;\r
    AliAODTracklets *trklets=aodEvent->GetTracklets();\r
    nSPD = trklets->GetNumberOfTracklets();\r
  }\r
  if(flagFilter)fHistTrackFilterSpdMult->Fill(nSPD,nCand);\r
  else fHistTrackAnaSpdMult->Fill(nSPD,nCand);\r
  \r
  Int_t runNumber = event->GetRunNumber();\r
  Int_t multiplicity = Multiplicity(event);\r
  if(nCand==0)return;\r
  if(flagFilter){\r
    if(fMultiplicity) \r
      fCounters.Count(Form("Event:Candid(Filter)/Run:%d/Multiplicity:%d",runNumber,multiplicity));\r
    else \r
      fCounters.Count(Form("Event:Candid(Filter)/Run:%d",runNumber));\r
    for(Int_t i=0;i<nCand;i++){ \r
      if(fMultiplicity) \r
	fCounters.Count(Form("Event:NCandid(Filter)/Run:%d/Multiplicity:%d",runNumber,multiplicity));\r
      else \r
	fCounters.Count(Form("Event:NCandid(Filter)/Run:%d",runNumber));\r
    }\r
  }else{\r
    if(fMultiplicity) \r
      fCounters.Count(Form("Event:Candid(Analysis)/Run:%d/Multiplicity:%d",runNumber,multiplicity));\r
    else\r
      fCounters.Count(Form("Event:Candid(Analysis)/Run:%d",runNumber));\r
    for(Int_t i=0;i<nCand;i++){ \r
      if(fMultiplicity) \r
	fCounters.Count(Form("Event:NCandid(Analysis)/Run:%d/Multiplicity:%d",runNumber,multiplicity));\r
      else\r
	fCounters.Count(Form("Event:NCandid(Analysis)/Run:%d",runNumber));\r
    }\r
  }\r
  return;\r
}\r
//_______________________________________________________________________\r
TH1D* AliNormalizationCounter::DrawAgainstRuns(TString candle,Bool_t drawHist){\r
  //\r
  fCounters.SortRubric("Run");\r
  TString selection;\r
  selection.Form("event:%s",candle.Data());\r
  TH1D* histoneD = fCounters.Get("run",selection.Data());\r
\r
  histoneD->Sumw2();\r
  if(drawHist)histoneD->DrawClone();\r
  return histoneD;\r
}\r
//___________________________________________________________________________\r
TH1D* AliNormalizationCounter::DrawRatio(TString candle1,TString candle2){\r
  //\r
  fCounters.SortRubric("Run");\r
  TString name;\r
\r
  name.Form("%s/%s",candle1.Data(),candle2.Data());\r
  TH1D* num=DrawAgainstRuns(candle1.Data(),kFALSE);\r
  TH1D* den=DrawAgainstRuns(candle2.Data(),kFALSE);\r
\r
  den->SetTitle(candle2.Data());\r
  den->SetName(candle2.Data());\r
  num->Divide(num,den,1,1,"B");\r
  num->SetTitle(name.Data());\r
  num->SetName(name.Data());\r
  num->DrawClone();\r
  return num;\r
}\r
//___________________________________________________________________________\r
void AliNormalizationCounter::PrintRubrics(){\r
  fCounters.PrintKeyWords();\r
}\r
//___________________________________________________________________________\r
Double_t AliNormalizationCounter::GetSum(TString candle){\r
  TString selection="event:";\r
  selection.Append(candle);\r
  return fCounters.GetSum(selection.Data());\r
}\r
//___________________________________________________________________________\r
TH2F* AliNormalizationCounter::GetHist(Bool_t filtercuts,Bool_t spdtracklets,Bool_t drawHist){\r
  if(filtercuts){\r
    if(spdtracklets){\r
      if(drawHist)fHistTrackFilterSpdMult->DrawCopy("LEGO2Z 0");\r
      return fHistTrackFilterSpdMult;\r
    }else{\r
      if(drawHist)fHistTrackFilterEvMult->DrawCopy("LEGO2Z 0");\r
      return fHistTrackFilterEvMult;\r
    }\r
  }else{\r
    if(spdtracklets){\r
      if(drawHist)fHistTrackAnaSpdMult->DrawCopy("LEGO2Z 0");\r
      return fHistTrackAnaSpdMult;\r
    }else{\r
      if(drawHist)fHistTrackAnaEvMult->DrawCopy("LEGO2Z 0");\r
      return fHistTrackAnaEvMult;\r
    }\r
  }\r
}\r
//___________________________________________________________________________\r
Double_t AliNormalizationCounter::GetNEventsForNorm(){\r
  Double_t noVtxzGT10=GetSum("noPrimaryV")*GetSum("zvtxGT10")/GetSum("PrimaryV");\r
  return GetSum("countForNorm")-noVtxzGT10;\r
}\r
//___________________________________________________________________________\r
Double_t AliNormalizationCounter::GetNEventsForNorm(Int_t runnumber){\r
  TString listofruns = fCounters.GetKeyWords("RUN");\r
  if(!listofruns.Contains(Form("%d",runnumber))){\r
    printf("WARNING: %d is not a valid run number\n",runnumber);\r
    fCounters.Print("Run","",kTRUE);\r
    return 0.;\r
  }\r
  TString suffix;suffix.Form("/RUN:%d",runnumber);\r
  TString zvtx;zvtx.Form("zvtxGT10%s",suffix.Data());\r
  TString noPV;noPV.Form("noPrimaryV%s",suffix.Data());\r
  TString pV;pV.Form("PrimaryV%s",suffix.Data());\r
  TString tbc;tbc.Form("countForNorm%s",suffix.Data());\r
  Double_t noVtxzGT10=GetSum(noPV.Data())*GetSum(zvtx.Data())/GetSum(pV.Data());\r
  return GetSum(tbc.Data())-noVtxzGT10;\r
}\r
\r
//___________________________________________________________________________\r
Double_t AliNormalizationCounter::GetNEventsForNorm(Int_t minmultiplicity, Int_t maxmultiplicity){\r
\r
  if(!fMultiplicity) {\r
    AliInfo("Sorry, you didn't activate the multiplicity in the counter!");\r
    return 0.;\r
  }\r
\r
  TString listofruns = fCounters.GetKeyWords("Multiplicity");\r
\r
  Int_t nmultbins = maxmultiplicity - minmultiplicity;\r
  Double_t sumnoPV=0., sumZvtx=0., sumPv=0., sumEvtNorm=0.;\r
  for (Int_t ibin=0; ibin<=nmultbins; ibin++) {\r
    //    cout << " Looking at bin "<< ibin+minmultiplicity<<endl;\r
    if(!listofruns.Contains(Form("%d",ibin+minmultiplicity))){\r
      AliInfo(Form("WARNING: %d is not a valid multiplicity number. \n",ibin+minmultiplicity));\r
      continue;\r
    }\r
    TString suffix;suffix.Form("/Multiplicity:%d",ibin+minmultiplicity);\r
    TString zvtx;zvtx.Form("zvtxGT10%s",suffix.Data());\r
    TString noPV;noPV.Form("noPrimaryV%s",suffix.Data());\r
    TString pV;pV.Form("PrimaryV%s",suffix.Data());\r
    TString tbc;tbc.Form("countForNorm%s",suffix.Data());\r
    sumnoPV += GetSum(noPV.Data());\r
    sumZvtx += GetSum(zvtx.Data());\r
    sumPv += GetSum(pV.Data());\r
    sumEvtNorm += GetSum(tbc.Data());\r
  }\r
  Double_t noVtxzGT10 = sumPv>0. ? sumnoPV * sumZvtx / sumPv : 0.;\r
  return sumEvtNorm - noVtxzGT10;\r
}\r
\r
//___________________________________________________________________________\r
TH1D* AliNormalizationCounter::DrawNEventsForNorm(Bool_t drawRatio){\r
  //usare algebra histos\r
  fCounters.SortRubric("Run");\r
  TString selection;\r
\r
  selection.Form("event:noPrimaryV");\r
  TH1D* hnoPrimV = fCounters.Get("run",selection.Data());\r
  hnoPrimV->Sumw2();\r
\r
  selection.Form("event:zvtxGT10");\r
  TH1D*  hzvtx= fCounters.Get("run",selection.Data());\r
  hzvtx->Sumw2();\r
\r
  selection.Form("event:PrimaryV");\r
  TH1D* hPrimV = fCounters.Get("run",selection.Data());\r
  hPrimV->Sumw2();\r
\r
  hzvtx->Multiply(hnoPrimV);\r
  hzvtx->Divide(hPrimV);\r
\r
  selection.Form("event:countForNorm");\r
  TH1D* hCountForNorm = fCounters.Get("run",selection.Data());\r
  hCountForNorm->Sumw2();\r
\r
  hCountForNorm->Add(hzvtx,-1.);\r
\r
  if(drawRatio){\r
    selection.Form("event:triggered");\r
    TH1D* htriggered = fCounters.Get("run",selection.Data());\r
    htriggered->Sumw2();\r
    hCountForNorm->Divide(htriggered);\r
  }\r
\r
  hCountForNorm->DrawClone();\r
  return hCountForNorm;\r
}\r
\r
//___________________________________________________________________________\r
Int_t AliNormalizationCounter::Multiplicity(AliVEvent* event){\r
\r
  Int_t multiplicity = 0;\r
  AliAODEvent *eventAOD = (AliAODEvent*)event;\r
  AliAODTracklets * aodTracklets = (AliAODTracklets*)eventAOD->GetTracklets();\r
  Int_t ntracklets = (Int_t)aodTracklets->GetNumberOfTracklets();\r
  for(Int_t i=0;i<ntracklets; i++){\r
    Double_t theta = aodTracklets->GetTheta(i);\r
    Double_t eta = -TMath::Log( TMath::Tan(theta/2.) ); // check the formula\r
    if(TMath::Abs(eta)<1.6){ // set the proper cut on eta\r
      multiplicity++;\r
    }\r
  }\r
\r
  return multiplicity;\r
}\r