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

/**************************************************************************
 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

//*************************************************************************
// Class AliNormalizationCounter
// Class to store the informations relevant for the normalization in the 
// barrel for each run
// Authors: G. Ortona, ortona@to.infn.it
// D. Caffarri, davide.caffarri@pd.to.infn.it
// with many thanks to P. Pillot
/////////////////////////////////////////////////////////////

#include "AliLog.h"
#include "AliNormalizationCounter.h"
#include <AliESDEvent.h>
#include <AliESDtrack.h>
#include <AliAODEvent.h>
#include <AliVParticle.h>
#include <AliTriggerAnalysis.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TList.h>
#include <TString.h>
#include <TCanvas.h>
#include <AliPhysicsSelection.h>
#include <AliMultiplicity.h>

ClassImp(AliNormalizationCounter)

//____________________________________________
AliNormalizationCounter::AliNormalizationCounter(): 
TNamed(),
fCounters(),
fESD(kFALSE),
fMultiplicity(kFALSE),
fMultiplicityEtaRange(1.0),
fHistTrackFilterEvMult(0),
fHistTrackAnaEvMult(0),
fHistTrackFilterSpdMult(0),
fHistTrackAnaSpdMult(0)
{
  // empty constructor
}

//__________________________________________________				
AliNormalizationCounter::AliNormalizationCounter(const char *name): 
TNamed(name,name),
fCounters(name),
fESD(kFALSE),
fMultiplicity(kFALSE),
fMultiplicityEtaRange(1.0),
fHistTrackFilterEvMult(0),
fHistTrackAnaEvMult(0),
fHistTrackFilterSpdMult(0),
fHistTrackAnaSpdMult(0)
{
  ;
}

//______________________________________________
AliNormalizationCounter::~AliNormalizationCounter()
{
  //destructor
  if(fHistTrackFilterEvMult){
    delete fHistTrackFilterEvMult;
    fHistTrackFilterEvMult =0;
  }
  if(fHistTrackAnaEvMult){
    delete fHistTrackAnaEvMult;
    fHistTrackAnaEvMult=0;
  }
  if(fHistTrackFilterSpdMult){
    delete fHistTrackFilterSpdMult;
    fHistTrackFilterSpdMult=0;
  }
  if(fHistTrackAnaSpdMult){
    delete fHistTrackAnaSpdMult;
    fHistTrackAnaSpdMult=0;
  }
}

//______________________________________________
void AliNormalizationCounter::Init()
{
  //variables initialization
  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)");
  if(fMultiplicity)  fCounters.AddRubric("Multiplicity", 5000);
  fCounters.AddRubric("Run", 1000000);
  fCounters.Init();
  fHistTrackFilterEvMult=new TH2F("FiltCandidvsTracksinEv","FiltCandidvsTracksinEv",10000,-0.5,9999.5,200,-0.5,199.5);
  fHistTrackFilterEvMult->GetYaxis()->SetTitle("NCandidates");
  fHistTrackFilterEvMult->GetXaxis()->SetTitle("NTracksinEvent");
  fHistTrackAnaEvMult=new TH2F("AnaCandidvsTracksinEv","AnaCandidvsTracksinEv",10000,-0.5,9999.5,100,-0.5,99.5);
  fHistTrackAnaEvMult->GetYaxis()->SetTitle("NCandidates");
  fHistTrackAnaEvMult->GetXaxis()->SetTitle("NTracksinEvent");
  fHistTrackFilterSpdMult=new TH2F("FilterCandidvsSpdMult","FilterCandidvsSpdMult",5000,-0.5,4999.5,200,-0.5,199.5);
  fHistTrackFilterSpdMult->GetYaxis()->SetTitle("NCandidates");
  fHistTrackFilterSpdMult->GetXaxis()->SetTitle("NSPDTracklets");
  fHistTrackAnaSpdMult=new TH2F("AnaCandidvsSpdMult","AnaCandidvsSpdMult",5000,-0.5,4999.5,100,-0.5,99.5);
  fHistTrackAnaSpdMult->GetYaxis()->SetTitle("NCandidates");
  fHistTrackAnaSpdMult->GetXaxis()->SetTitle("NSPDTracklets");
}

//______________________________________________
Long64_t AliNormalizationCounter::Merge(TCollection* list){
  if (!list) return 0;
  if (list->IsEmpty()) return 0;//(Long64_t)fCounters.Merge(list);

  TIter next(list);
  const TObject* obj = 0x0;
  while ((obj = next())) {
    
    // check that "obj" is an object of the class AliNormalizationCounter
    const AliNormalizationCounter* counter = dynamic_cast<const AliNormalizationCounter*>(obj);
    if (!counter) {
      AliError(Form("object named %s is not AliNormalizationCounter! Skipping it.", counter->GetName()));
      continue;
    }

    Add(counter);

  }
  
  return (Long64_t)1;//(Long64_t)fCounters->GetEntries();
}
//_______________________________________
void AliNormalizationCounter::Add(const AliNormalizationCounter *norm){
  fCounters.Add(&(norm->fCounters));
  fHistTrackFilterEvMult->Add(norm->fHistTrackFilterEvMult);
  fHistTrackAnaEvMult->Add(norm->fHistTrackAnaEvMult);
  fHistTrackFilterSpdMult->Add(norm->fHistTrackFilterSpdMult);
  fHistTrackAnaSpdMult->Add(norm->fHistTrackAnaSpdMult);
}
//_______________________________________
/*
Stores the variables used for normalization as function of run number
returns kTRUE if the event is to be counted for normalization
(pass event selection cuts OR has no primary vertex)
 */
void AliNormalizationCounter::StoreEvent(AliVEvent *event,AliRDHFCuts *rdCut,Bool_t mc, Int_t multiplicity){
  //

  Bool_t isEventSelected = rdCut->IsEventSelected(event);

  // events not passing physics selection. do nothing
  if(rdCut->IsEventRejectedDuePhysicsSelection()) return;

  Bool_t v0A=kFALSE; 
  Bool_t v0B=kFALSE;
  Bool_t flag03=kFALSE;
  Bool_t flagPV=kFALSE;

  //Run Number
  Int_t runNumber = event->GetRunNumber();
 
  // Evaluate the multiplicity
  if(multiplicity==-9999) Multiplicity(event);

  //Find CINT1B
  AliESDEvent *eventESD = (AliESDEvent*)event;
  if(!eventESD){AliError("ESD event not available");return;}
  if(mc&&event->GetEventType() != 0)return;
  //event must be either physics or MC
  if(!(event->GetEventType() == 7||event->GetEventType() == 0))return;
  
  if(fMultiplicity) 
    fCounters.Count(Form("Event:triggered/Run:%d/Multiplicity:%d",runNumber,multiplicity));
  else 
    fCounters.Count(Form("Event:triggered/Run:%d",runNumber));

  //Find V0AND
  AliTriggerAnalysis trAn; /// Trigger Analysis
  v0B = trAn.IsOfflineTriggerFired(eventESD , AliTriggerAnalysis::kV0C);
  v0A = trAn.IsOfflineTriggerFired(eventESD , AliTriggerAnalysis::kV0A);
  if(v0A&&v0B){
    if(fMultiplicity) 
      fCounters.Count(Form("Event:V0AND/Run:%d/Multiplicity:%d",runNumber,multiplicity));
    else
      fCounters.Count(Form("Event:V0AND/Run:%d",runNumber));
  }
  
  //FindPrimary vertex  
  // AliVVertex *vtrc =  (AliVVertex*)event->GetPrimaryVertex();
  // if(vtrc && vtrc->GetNContributors()>0){
  //   fCounters.Count(Form("Event:PrimaryV/Run:%d",runNumber));
  //   flagPV=kTRUE;
  // }

  //trigger
  AliAODEvent *eventAOD = (AliAODEvent*)event;
  TString trigclass=eventAOD->GetFiredTriggerClasses();
  if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD")||trigclass.Contains("C0SMH-B-NOPF-ALL")){
    if(fMultiplicity) 
      fCounters.Count(Form("Event:PbPbC0SMH-B-NOPF-ALLNOTRD/Run:%d/Multiplicity:%d",runNumber,multiplicity));
    else 
      fCounters.Count(Form("Event:PbPbC0SMH-B-NOPF-ALLNOTRD/Run:%d",runNumber));
  }

  //FindPrimary vertex  
  if(isEventSelected){
    if(fMultiplicity) 
      fCounters.Count(Form("Event:PrimaryV/Run:%d/Multiplicity:%d",runNumber,multiplicity));
    else
      fCounters.Count(Form("Event:PrimaryV/Run:%d",runNumber));
    flagPV=kTRUE;
  }else{
    if(rdCut->GetWhyRejection()==0){
      if(fMultiplicity) 
	fCounters.Count(Form("Event:noPrimaryV/Run:%d/Multiplicity:%d",runNumber,multiplicity));
      else
	fCounters.Count(Form("Event:noPrimaryV/Run:%d",runNumber));
    }
    //find good vtx outside range
    if(rdCut->GetWhyRejection()==6){
      if(fMultiplicity) {
	fCounters.Count(Form("Event:zvtxGT10/Run:%d/Multiplicity:%d",runNumber,multiplicity));
	fCounters.Count(Form("Event:PrimaryV/Run:%d/Multiplicity:%d",runNumber,multiplicity));
      } else {
	fCounters.Count(Form("Event:zvtxGT10/Run:%d",runNumber));
	fCounters.Count(Form("Event:PrimaryV/Run:%d",runNumber));
      }
      flagPV=kTRUE;
    }
    if(rdCut->GetWhyRejection()==1){
      if(fMultiplicity) 
	fCounters.Count(Form("Event:PileUp/Run:%d/Multiplicity:%d",runNumber,multiplicity));
      else
	fCounters.Count(Form("Event:PileUp/Run:%d",runNumber));
    }
  }
  //to be counted for normalization
  if(rdCut->CountEventForNormalization()){
    if(fMultiplicity) 
      fCounters.Count(Form("Event:countForNorm/Run:%d/Multiplicity:%d",runNumber,multiplicity));
    else 
      fCounters.Count(Form("Event:countForNorm/Run:%d",runNumber));
  }


  //Find Candle
  Int_t trkEntries = (Int_t)event->GetNumberOfTracks();
  for(Int_t i=0;i<trkEntries&&!flag03;i++){
    AliAODTrack *track=(AliAODTrack*)event->GetTrack(i);
    if((track->Pt()>0.3)&&(!flag03)){
      if(fMultiplicity) 
	fCounters.Count(Form("Event:Candles0.3/Run:%d/Multiplicity:%d",runNumber,multiplicity));
      else
	fCounters.Count(Form("Event:Candles0.3/Run:%d",runNumber));
      flag03=kTRUE;
      break;
    }
  }
  
  if(!(v0A&&v0B)&&(flag03)){ 
    if(fMultiplicity) 
      fCounters.Count(Form("Event:!V0A&Candle03/Run:%d/Multiplicity:%d",runNumber,multiplicity));
    else 
      fCounters.Count(Form("Event:!V0A&Candle03/Run:%d",runNumber));
  }
  if(!(v0A&&v0B)&&flagPV){
    if(fMultiplicity) 
      fCounters.Count(Form("Event:!V0A&PrimaryV/Run:%d/Multiplicity:%d",runNumber,multiplicity));
    else
      fCounters.Count(Form("Event:!V0A&PrimaryV/Run:%d",runNumber));
  }
  
  return;
}
//_____________________________________________________________________
void AliNormalizationCounter::StoreCandidates(AliVEvent *event,Int_t nCand,Bool_t flagFilter){
  
  Int_t ntracks=event->GetNumberOfTracks();
  if(flagFilter)fHistTrackFilterEvMult->Fill(ntracks,nCand);
  else fHistTrackAnaEvMult->Fill(ntracks,nCand);
  Int_t nSPD=0;
  if(fESD){
    AliESDEvent *ESDevent=(AliESDEvent*)event;
    const AliMultiplicity *alimult = ESDevent->GetMultiplicity();
    nSPD = alimult->GetNumberOfTracklets();

  }else{
    AliAODEvent *aodEvent =(AliAODEvent*)event;
    AliAODTracklets *trklets=aodEvent->GetTracklets();
    nSPD = trklets->GetNumberOfTracklets();
  }
  if(flagFilter)fHistTrackFilterSpdMult->Fill(nSPD,nCand);
  else fHistTrackAnaSpdMult->Fill(nSPD,nCand);
  
  Int_t runNumber = event->GetRunNumber();
  Int_t multiplicity = Multiplicity(event);
  if(nCand==0)return;
  if(flagFilter){
    if(fMultiplicity) 
      fCounters.Count(Form("Event:Candid(Filter)/Run:%d/Multiplicity:%d",runNumber,multiplicity));
    else 
      fCounters.Count(Form("Event:Candid(Filter)/Run:%d",runNumber));
    for(Int_t i=0;i<nCand;i++){ 
      if(fMultiplicity) 
	fCounters.Count(Form("Event:NCandid(Filter)/Run:%d/Multiplicity:%d",runNumber,multiplicity));
      else 
	fCounters.Count(Form("Event:NCandid(Filter)/Run:%d",runNumber));
    }
  }else{
    if(fMultiplicity) 
      fCounters.Count(Form("Event:Candid(Analysis)/Run:%d/Multiplicity:%d",runNumber,multiplicity));
    else
      fCounters.Count(Form("Event:Candid(Analysis)/Run:%d",runNumber));
    for(Int_t i=0;i<nCand;i++){ 
      if(fMultiplicity) 
	fCounters.Count(Form("Event:NCandid(Analysis)/Run:%d/Multiplicity:%d",runNumber,multiplicity));
      else
	fCounters.Count(Form("Event:NCandid(Analysis)/Run:%d",runNumber));
    }
  }
  return;
}
//_______________________________________________________________________
TH1D* AliNormalizationCounter::DrawAgainstRuns(TString candle,Bool_t drawHist){
  //
  fCounters.SortRubric("Run");
  TString selection;
  selection.Form("event:%s",candle.Data());
  TH1D* histoneD = fCounters.Get("run",selection.Data());

  histoneD->Sumw2();
  if(drawHist)histoneD->DrawClone();
  return histoneD;
}
//___________________________________________________________________________
TH1D* AliNormalizationCounter::DrawRatio(TString candle1,TString candle2){
  //
  fCounters.SortRubric("Run");
  TString name;

  name.Form("%s/%s",candle1.Data(),candle2.Data());
  TH1D* num=DrawAgainstRuns(candle1.Data(),kFALSE);
  TH1D* den=DrawAgainstRuns(candle2.Data(),kFALSE);

  den->SetTitle(candle2.Data());
  den->SetName(candle2.Data());
  num->Divide(num,den,1,1,"B");
  num->SetTitle(name.Data());
  num->SetName(name.Data());
  num->DrawClone();
  return num;
}
//___________________________________________________________________________
void AliNormalizationCounter::PrintRubrics(){
  fCounters.PrintKeyWords();
}
//___________________________________________________________________________
Double_t AliNormalizationCounter::GetSum(TString candle){
  TString selection="event:";
  selection.Append(candle);
  return fCounters.GetSum(selection.Data());
}
//___________________________________________________________________________
TH2F* AliNormalizationCounter::GetHist(Bool_t filtercuts,Bool_t spdtracklets,Bool_t drawHist){
  if(filtercuts){
    if(spdtracklets){
      if(drawHist)fHistTrackFilterSpdMult->DrawCopy("LEGO2Z 0");
      return fHistTrackFilterSpdMult;
    }else{
      if(drawHist)fHistTrackFilterEvMult->DrawCopy("LEGO2Z 0");
      return fHistTrackFilterEvMult;
    }
  }else{
    if(spdtracklets){
      if(drawHist)fHistTrackAnaSpdMult->DrawCopy("LEGO2Z 0");
      return fHistTrackAnaSpdMult;
    }else{
      if(drawHist)fHistTrackAnaEvMult->DrawCopy("LEGO2Z 0");
      return fHistTrackAnaEvMult;
    }
  }
}
//___________________________________________________________________________
Double_t AliNormalizationCounter::GetNEventsForNorm(){
  Double_t noVtxzGT10=GetSum("noPrimaryV")*GetSum("zvtxGT10")/GetSum("PrimaryV");
  return GetSum("countForNorm")-noVtxzGT10;
}
//___________________________________________________________________________
Double_t AliNormalizationCounter::GetNEventsForNorm(Int_t runnumber){
  TString listofruns = fCounters.GetKeyWords("RUN");
  if(!listofruns.Contains(Form("%d",runnumber))){
    printf("WARNING: %d is not a valid run number\n",runnumber);
    fCounters.Print("Run","",kTRUE);
    return 0.;
  }
  TString suffix;suffix.Form("/RUN:%d",runnumber);
  TString zvtx;zvtx.Form("zvtxGT10%s",suffix.Data());
  TString noPV;noPV.Form("noPrimaryV%s",suffix.Data());
  TString pV;pV.Form("PrimaryV%s",suffix.Data());
  TString tbc;tbc.Form("countForNorm%s",suffix.Data());
  Double_t noVtxzGT10=GetSum(noPV.Data())*GetSum(zvtx.Data())/GetSum(pV.Data());
  return GetSum(tbc.Data())-noVtxzGT10;
}

//___________________________________________________________________________
Double_t AliNormalizationCounter::GetNEventsForNorm(Int_t minmultiplicity, Int_t maxmultiplicity){

  if(!fMultiplicity) {
    AliInfo("Sorry, you didn't activate the multiplicity in the counter!");
    return 0.;
  }

  TString listofruns = fCounters.GetKeyWords("Multiplicity");

  Int_t nmultbins = maxmultiplicity - minmultiplicity;
  Double_t sumnoPV=0., sumZvtx=0., sumPv=0., sumEvtNorm=0.;
  for (Int_t ibin=0; ibin<=nmultbins; ibin++) {
    //    cout << " Looking at bin "<< ibin+minmultiplicity<<endl;
    if(!listofruns.Contains(Form("%d",ibin+minmultiplicity))){
      AliInfo(Form("WARNING: %d is not a valid multiplicity number. \n",ibin+minmultiplicity));
      continue;
    }
    TString suffix;suffix.Form("/Multiplicity:%d",ibin+minmultiplicity);
    TString zvtx;zvtx.Form("zvtxGT10%s",suffix.Data());
    TString noPV;noPV.Form("noPrimaryV%s",suffix.Data());
    TString pV;pV.Form("PrimaryV%s",suffix.Data());
    TString tbc;tbc.Form("countForNorm%s",suffix.Data());
    sumnoPV += GetSum(noPV.Data());
    sumZvtx += GetSum(zvtx.Data());
    sumPv += GetSum(pV.Data());
    sumEvtNorm += GetSum(tbc.Data());
  }
  Double_t noVtxzGT10 = sumPv>0. ? sumnoPV * sumZvtx / sumPv : 0.;
  return sumEvtNorm - noVtxzGT10;
}

//___________________________________________________________________________
TH1D* AliNormalizationCounter::DrawNEventsForNorm(Bool_t drawRatio){
  //usare algebra histos
  fCounters.SortRubric("Run");
  TString selection;

  selection.Form("event:noPrimaryV");
  TH1D* hnoPrimV = fCounters.Get("run",selection.Data());
  hnoPrimV->Sumw2();

  selection.Form("event:zvtxGT10");
  TH1D*  hzvtx= fCounters.Get("run",selection.Data());
  hzvtx->Sumw2();

  selection.Form("event:PrimaryV");
  TH1D* hPrimV = fCounters.Get("run",selection.Data());
  hPrimV->Sumw2();

  hzvtx->Multiply(hnoPrimV);
  hzvtx->Divide(hPrimV);

  selection.Form("event:countForNorm");
  TH1D* hCountForNorm = fCounters.Get("run",selection.Data());
  hCountForNorm->Sumw2();

  hCountForNorm->Add(hzvtx,-1.);

  if(drawRatio){
    selection.Form("event:triggered");
    TH1D* htriggered = fCounters.Get("run",selection.Data());
    htriggered->Sumw2();
    hCountForNorm->Divide(htriggered);
  }

  hCountForNorm->DrawClone();
  return hCountForNorm;
}

//___________________________________________________________________________
Int_t AliNormalizationCounter::Multiplicity(AliVEvent* event){

  Int_t multiplicity = 0;
  AliAODEvent *eventAOD = (AliAODEvent*)event;
  AliAODTracklets * aodTracklets = (AliAODTracklets*)eventAOD->GetTracklets();
  Int_t ntracklets = (Int_t)aodTracklets->GetNumberOfTracklets();
  for(Int_t i=0;i<ntracklets; i++){
    Double_t theta = aodTracklets->GetTheta(i);
    Double_t eta = -TMath::Log( TMath::Tan(theta/2.) ); // check the formula
    if(TMath::Abs(eta)<fMultiplicityEtaRange){ // set the proper cut on eta
      multiplicity++;
    }
  }

  return multiplicity;
}
Commit	Line	Data
8eea8388	1	/**************************************************************************
	2	* Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	/* $Id$ */
	17
	18	//*************************************************************************
	19	// Class AliNormalizationCounter
	20	// Class to store the informations relevant for the normalization in the
	21	// barrel for each run
	22	// Authors: G. Ortona, ortona@to.infn.it
	23	// D. Caffarri, davide.caffarri@pd.to.infn.it
	24	// with many thanks to P. Pillot
	25	/////////////////////////////////////////////////////////////
	26
	27	#include "AliLog.h"
	28	#include "AliNormalizationCounter.h"
	29	#include <AliESDEvent.h>
	30	#include <AliESDtrack.h>
	31	#include <AliAODEvent.h>
	32	#include <AliVParticle.h>
	33	#include <AliTriggerAnalysis.h>
	34	#include <TH1F.h>
	35	#include <TH2F.h>
	36	#include <TList.h>
	37	#include <TString.h>
	38	#include <TCanvas.h>
	39	#include <AliPhysicsSelection.h>
	40	#include <AliMultiplicity.h>
	41
	42	ClassImp(AliNormalizationCounter)
	43
	44	//____________________________________________
	45	AliNormalizationCounter::AliNormalizationCounter():
	46	TNamed(),
	47	fCounters(),
	48	fESD(kFALSE),
	49	fMultiplicity(kFALSE),
	50	fMultiplicityEtaRange(1.0),
	51	fHistTrackFilterEvMult(0),
	52	fHistTrackAnaEvMult(0),
	53	fHistTrackFilterSpdMult(0),
	54	fHistTrackAnaSpdMult(0)
	55	{
	56	// empty constructor
	57	}
	58
	59	//__________________________________________________
	60	AliNormalizationCounter::AliNormalizationCounter(const char *name):
	61	TNamed(name,name),
	62	fCounters(name),
	63	fESD(kFALSE),
	64	fMultiplicity(kFALSE),
65	fMultiplicityEtaRange(1.0),
66	fHistTrackFilterEvMult(0),
67	fHistTrackAnaEvMult(0),
68	fHistTrackFilterSpdMult(0),
69	fHistTrackAnaSpdMult(0)
70	{
71	;
72	}
73
74	//______________________________________________
75	AliNormalizationCounter::~AliNormalizationCounter()
76	{
77	//destructor
78	if(fHistTrackFilterEvMult){
79	delete fHistTrackFilterEvMult;
80	fHistTrackFilterEvMult =0;
81	}
82	if(fHistTrackAnaEvMult){
83	delete fHistTrackAnaEvMult;
84	fHistTrackAnaEvMult=0;
85	}
86	if(fHistTrackFilterSpdMult){
87	delete fHistTrackFilterSpdMult;
88	fHistTrackFilterSpdMult=0;
89	}
90	if(fHistTrackAnaSpdMult){
91	delete fHistTrackAnaSpdMult;
92	fHistTrackAnaSpdMult=0;
93	}
94	}
95
96	//______________________________________________
97	void AliNormalizationCounter::Init()
98	{
99	//variables initialization
100	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)");
101	if(fMultiplicity) fCounters.AddRubric("Multiplicity", 5000);
102	fCounters.AddRubric("Run", 1000000);
103	fCounters.Init();
104	fHistTrackFilterEvMult=new TH2F("FiltCandidvsTracksinEv","FiltCandidvsTracksinEv",10000,-0.5,9999.5,200,-0.5,199.5);
105	fHistTrackFilterEvMult->GetYaxis()->SetTitle("NCandidates");
106	fHistTrackFilterEvMult->GetXaxis()->SetTitle("NTracksinEvent");
107	fHistTrackAnaEvMult=new TH2F("AnaCandidvsTracksinEv","AnaCandidvsTracksinEv",10000,-0.5,9999.5,100,-0.5,99.5);
108	fHistTrackAnaEvMult->GetYaxis()->SetTitle("NCandidates");
109	fHistTrackAnaEvMult->GetXaxis()->SetTitle("NTracksinEvent");
110	fHistTrackFilterSpdMult=new TH2F("FilterCandidvsSpdMult","FilterCandidvsSpdMult",5000,-0.5,4999.5,200,-0.5,199.5);
111	fHistTrackFilterSpdMult->GetYaxis()->SetTitle("NCandidates");
112	fHistTrackFilterSpdMult->GetXaxis()->SetTitle("NSPDTracklets");
113	fHistTrackAnaSpdMult=new TH2F("AnaCandidvsSpdMult","AnaCandidvsSpdMult",5000,-0.5,4999.5,100,-0.5,99.5);
114	fHistTrackAnaSpdMult->GetYaxis()->SetTitle("NCandidates");
115	fHistTrackAnaSpdMult->GetXaxis()->SetTitle("NSPDTracklets");
116	}
117
118	//______________________________________________
119	Long64_t AliNormalizationCounter::Merge(TCollection* list){
120	if (!list) return 0;
121	if (list->IsEmpty()) return 0;//(Long64_t)fCounters.Merge(list);
122
123	TIter next(list);
124	const TObject* obj = 0x0;
125	while ((obj = next())) {
126
127	// check that "obj" is an object of the class AliNormalizationCounter
128	const AliNormalizationCounter* counter = dynamic_cast<const AliNormalizationCounter*>(obj);
129	if (!counter) {
130	AliError(Form("object named %s is not AliNormalizationCounter! Skipping it.", counter->GetName()));
131	continue;
132	}
133
134	Add(counter);
135
136	}
137
138	return (Long64_t)1;//(Long64_t)fCounters->GetEntries();
139	}
140	//_______________________________________
141	void AliNormalizationCounter::Add(const AliNormalizationCounter *norm){
142	fCounters.Add(&(norm->fCounters));
143	fHistTrackFilterEvMult->Add(norm->fHistTrackFilterEvMult);
144	fHistTrackAnaEvMult->Add(norm->fHistTrackAnaEvMult);
145	fHistTrackFilterSpdMult->Add(norm->fHistTrackFilterSpdMult);
146	fHistTrackAnaSpdMult->Add(norm->fHistTrackAnaSpdMult);
147	}
148	//_______________________________________
149	/*
150	Stores the variables used for normalization as function of run number
151	returns kTRUE if the event is to be counted for normalization
152	(pass event selection cuts OR has no primary vertex)
153	*/
71d2489f	154	void AliNormalizationCounter::StoreEvent(AliVEvent event,AliRDHFCuts rdCut,Bool_t mc, Int_t multiplicity){
8eea8388	155	//
	156
	157	Bool_t isEventSelected = rdCut->IsEventSelected(event);
	158
	159	// events not passing physics selection. do nothing
	160	if(rdCut->IsEventRejectedDuePhysicsSelection()) return;
	161
	162	Bool_t v0A=kFALSE;
	163	Bool_t v0B=kFALSE;
	164	Bool_t flag03=kFALSE;
	165	Bool_t flagPV=kFALSE;
	166
	167	//Run Number
	168	Int_t runNumber = event->GetRunNumber();
	169
	170	// Evaluate the multiplicity
71d2489f	171	if(multiplicity==-9999) Multiplicity(event);
8eea8388	172
	173	//Find CINT1B
	174	AliESDEvent eventESD = (AliESDEvent)event;
	175	if(!eventESD){AliError("ESD event not available");return;}
	176	if(mc&&event->GetEventType() != 0)return;
	177	//event must be either physics or MC
	178	if(!(event->GetEventType() == 7\|\|event->GetEventType() == 0))return;
	179
	180	if(fMultiplicity)
	181	fCounters.Count(Form("Event:triggered/Run:%d/Multiplicity:%d",runNumber,multiplicity));
	182	else
	183	fCounters.Count(Form("Event:triggered/Run:%d",runNumber));
	184
	185	//Find V0AND
	186	AliTriggerAnalysis trAn; /// Trigger Analysis
	187	v0B = trAn.IsOfflineTriggerFired(eventESD , AliTriggerAnalysis::kV0C);
	188	v0A = trAn.IsOfflineTriggerFired(eventESD , AliTriggerAnalysis::kV0A);
	189	if(v0A&&v0B){
	190	if(fMultiplicity)
	191	fCounters.Count(Form("Event:V0AND/Run:%d/Multiplicity:%d",runNumber,multiplicity));
	192	else
	193	fCounters.Count(Form("Event:V0AND/Run:%d",runNumber));
	194	}
	195
	196	//FindPrimary vertex
	197	// AliVVertex vtrc = (AliVVertex)event->GetPrimaryVertex();
	198	// if(vtrc && vtrc->GetNContributors()>0){
	199	// fCounters.Count(Form("Event:PrimaryV/Run:%d",runNumber));
	200	// flagPV=kTRUE;
	201	// }
	202
	203	//trigger
	204	AliAODEvent eventAOD = (AliAODEvent)event;
	205	TString trigclass=eventAOD->GetFiredTriggerClasses();
	206	if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD")\|\|trigclass.Contains("C0SMH-B-NOPF-ALL")){
	207	if(fMultiplicity)
	208	fCounters.Count(Form("Event:PbPbC0SMH-B-NOPF-ALLNOTRD/Run:%d/Multiplicity:%d",runNumber,multiplicity));
	209	else
	210	fCounters.Count(Form("Event:PbPbC0SMH-B-NOPF-ALLNOTRD/Run:%d",runNumber));
	211	}
	212
	213	//FindPrimary vertex
	214	if(isEventSelected){
	215	if(fMultiplicity)
	216	fCounters.Count(Form("Event:PrimaryV/Run:%d/Multiplicity:%d",runNumber,multiplicity));
	217	else
	218	fCounters.Count(Form("Event:PrimaryV/Run:%d",runNumber));
	219	flagPV=kTRUE;
	220	}else{
	221	if(rdCut->GetWhyRejection()==0){
	222	if(fMultiplicity)
	223	fCounters.Count(Form("Event:noPrimaryV/Run:%d/Multiplicity:%d",runNumber,multiplicity));
	224	else
	225	fCounters.Count(Form("Event:noPrimaryV/Run:%d",runNumber));
	226	}
	227	//find good vtx outside range
	228	if(rdCut->GetWhyRejection()==6){
	229	if(fMultiplicity) {
	230	fCounters.Count(Form("Event:zvtxGT10/Run:%d/Multiplicity:%d",runNumber,multiplicity));
	231	fCounters.Count(Form("Event:PrimaryV/Run:%d/Multiplicity:%d",runNumber,multiplicity));
	232	} else {
	233	fCounters.Count(Form("Event:zvtxGT10/Run:%d",runNumber));
	234	fCounters.Count(Form("Event:PrimaryV/Run:%d",runNumber));
	235	}
236	flagPV=kTRUE;
237	}
238	if(rdCut->GetWhyRejection()==1){
239	if(fMultiplicity)
240	fCounters.Count(Form("Event:PileUp/Run:%d/Multiplicity:%d",runNumber,multiplicity));
241	else
242	fCounters.Count(Form("Event:PileUp/Run:%d",runNumber));
243	}
244	}
245	//to be counted for normalization
246	if(rdCut->CountEventForNormalization()){
247	if(fMultiplicity)
248	fCounters.Count(Form("Event:countForNorm/Run:%d/Multiplicity:%d",runNumber,multiplicity));
249	else
250	fCounters.Count(Form("Event:countForNorm/Run:%d",runNumber));
251	}
252
253
254	//Find Candle
255	Int_t trkEntries = (Int_t)event->GetNumberOfTracks();
256	for(Int_t i=0;i<trkEntries&&!flag03;i++){
257	AliAODTrack track=(AliAODTrack)event->GetTrack(i);
258	if((track->Pt()>0.3)&&(!flag03)){
259	if(fMultiplicity)
260	fCounters.Count(Form("Event:Candles0.3/Run:%d/Multiplicity:%d",runNumber,multiplicity));
261	else
262	fCounters.Count(Form("Event:Candles0.3/Run:%d",runNumber));
263	flag03=kTRUE;
264	break;
265	}
266	}
267
268	if(!(v0A&&v0B)&&(flag03)){
269	if(fMultiplicity)
270	fCounters.Count(Form("Event:!V0A&Candle03/Run:%d/Multiplicity:%d",runNumber,multiplicity));
271	else
272	fCounters.Count(Form("Event:!V0A&Candle03/Run:%d",runNumber));
273	}
274	if(!(v0A&&v0B)&&flagPV){
275	if(fMultiplicity)
276	fCounters.Count(Form("Event:!V0A&PrimaryV/Run:%d/Multiplicity:%d",runNumber,multiplicity));
277	else
278	fCounters.Count(Form("Event:!V0A&PrimaryV/Run:%d",runNumber));
279	}
280
281	return;
282	}
283	//_____________________________________________________________________
284	void AliNormalizationCounter::StoreCandidates(AliVEvent *event,Int_t nCand,Bool_t flagFilter){
285
286	Int_t ntracks=event->GetNumberOfTracks();
287	if(flagFilter)fHistTrackFilterEvMult->Fill(ntracks,nCand);
288	else fHistTrackAnaEvMult->Fill(ntracks,nCand);
289	Int_t nSPD=0;
290	if(fESD){
291	AliESDEvent ESDevent=(AliESDEvent)event;
292	const AliMultiplicity *alimult = ESDevent->GetMultiplicity();
293	nSPD = alimult->GetNumberOfTracklets();
294
295	}else{
296	AliAODEvent aodEvent =(AliAODEvent)event;
297	AliAODTracklets *trklets=aodEvent->GetTracklets();
298	nSPD = trklets->GetNumberOfTracklets();
299	}
300	if(flagFilter)fHistTrackFilterSpdMult->Fill(nSPD,nCand);
301	else fHistTrackAnaSpdMult->Fill(nSPD,nCand);
302
303	Int_t runNumber = event->GetRunNumber();
304	Int_t multiplicity = Multiplicity(event);
305	if(nCand==0)return;
306	if(flagFilter){
307	if(fMultiplicity)
308	fCounters.Count(Form("Event:Candid(Filter)/Run:%d/Multiplicity:%d",runNumber,multiplicity));
309	else
310	fCounters.Count(Form("Event:Candid(Filter)/Run:%d",runNumber));
311	for(Int_t i=0;i<nCand;i++){
312	if(fMultiplicity)
313	fCounters.Count(Form("Event:NCandid(Filter)/Run:%d/Multiplicity:%d",runNumber,multiplicity));
314	else
315	fCounters.Count(Form("Event:NCandid(Filter)/Run:%d",runNumber));
316	}
317	}else{
318	if(fMultiplicity)
319	fCounters.Count(Form("Event:Candid(Analysis)/Run:%d/Multiplicity:%d",runNumber,multiplicity));
320	else
321	fCounters.Count(Form("Event:Candid(Analysis)/Run:%d",runNumber));
322	for(Int_t i=0;i<nCand;i++){
323	if(fMultiplicity)
324	fCounters.Count(Form("Event:NCandid(Analysis)/Run:%d/Multiplicity:%d",runNumber,multiplicity));
325	else
326	fCounters.Count(Form("Event:NCandid(Analysis)/Run:%d",runNumber));
327	}
328	}
329	return;
330	}
331	//_______________________________________________________________________
332	TH1D* AliNormalizationCounter::DrawAgainstRuns(TString candle,Bool_t drawHist){
333	//
334	fCounters.SortRubric("Run");
335	TString selection;
336	selection.Form("event:%s",candle.Data());
337	TH1D* histoneD = fCounters.Get("run",selection.Data());
338
339	histoneD->Sumw2();
340	if(drawHist)histoneD->DrawClone();
341	return histoneD;
342	}
343	//___________________________________________________________________________
344	TH1D* AliNormalizationCounter::DrawRatio(TString candle1,TString candle2){
345	//
346	fCounters.SortRubric("Run");
347	TString name;
348
349	name.Form("%s/%s",candle1.Data(),candle2.Data());
350	TH1D* num=DrawAgainstRuns(candle1.Data(),kFALSE);
351	TH1D* den=DrawAgainstRuns(candle2.Data(),kFALSE);
352
353	den->SetTitle(candle2.Data());
354	den->SetName(candle2.Data());
355	num->Divide(num,den,1,1,"B");
356	num->SetTitle(name.Data());
357	num->SetName(name.Data());
358	num->DrawClone();
359	return num;
360	}
361	//___________________________________________________________________________
362	void AliNormalizationCounter::PrintRubrics(){
363	fCounters.PrintKeyWords();
364	}
365	//___________________________________________________________________________
366	Double_t AliNormalizationCounter::GetSum(TString candle){
367	TString selection="event:";
368	selection.Append(candle);
369	return fCounters.GetSum(selection.Data());
370	}
371	//___________________________________________________________________________
372	TH2F* AliNormalizationCounter::GetHist(Bool_t filtercuts,Bool_t spdtracklets,Bool_t drawHist){
373	if(filtercuts){
374	if(spdtracklets){
375	if(drawHist)fHistTrackFilterSpdMult->DrawCopy("LEGO2Z 0");
376	return fHistTrackFilterSpdMult;
377	}else{
378	if(drawHist)fHistTrackFilterEvMult->DrawCopy("LEGO2Z 0");
379	return fHistTrackFilterEvMult;
380	}
381	}else{
382	if(spdtracklets){
383	if(drawHist)fHistTrackAnaSpdMult->DrawCopy("LEGO2Z 0");
384	return fHistTrackAnaSpdMult;
385	}else{
386	if(drawHist)fHistTrackAnaEvMult->DrawCopy("LEGO2Z 0");
387	return fHistTrackAnaEvMult;
388	}
389	}
390	}
391	//___________________________________________________________________________
392	Double_t AliNormalizationCounter::GetNEventsForNorm(){
393	Double_t noVtxzGT10=GetSum("noPrimaryV")*GetSum("zvtxGT10")/GetSum("PrimaryV");
394	return GetSum("countForNorm")-noVtxzGT10;
395	}
396	//___________________________________________________________________________
397	Double_t AliNormalizationCounter::GetNEventsForNorm(Int_t runnumber){
398	TString listofruns = fCounters.GetKeyWords("RUN");
399	if(!listofruns.Contains(Form("%d",runnumber))){
400	printf("WARNING: %d is not a valid run number\n",runnumber);
401	fCounters.Print("Run","",kTRUE);
402	return 0.;
403	}
404	TString suffix;suffix.Form("/RUN:%d",runnumber);
405	TString zvtx;zvtx.Form("zvtxGT10%s",suffix.Data());
406	TString noPV;noPV.Form("noPrimaryV%s",suffix.Data());
407	TString pV;pV.Form("PrimaryV%s",suffix.Data());
408	TString tbc;tbc.Form("countForNorm%s",suffix.Data());
409	Double_t noVtxzGT10=GetSum(noPV.Data())*GetSum(zvtx.Data())/GetSum(pV.Data());
410	return GetSum(tbc.Data())-noVtxzGT10;
411	}
412
413	//___________________________________________________________________________
414	Double_t AliNormalizationCounter::GetNEventsForNorm(Int_t minmultiplicity, Int_t maxmultiplicity){
415
416	if(!fMultiplicity) {
417	AliInfo("Sorry, you didn't activate the multiplicity in the counter!");
418	return 0.;
419	}
420
421	TString listofruns = fCounters.GetKeyWords("Multiplicity");
422
423	Int_t nmultbins = maxmultiplicity - minmultiplicity;
424	Double_t sumnoPV=0., sumZvtx=0., sumPv=0., sumEvtNorm=0.;
425	for (Int_t ibin=0; ibin<=nmultbins; ibin++) {
426	// cout << " Looking at bin "<< ibin+minmultiplicity<<endl;
427	if(!listofruns.Contains(Form("%d",ibin+minmultiplicity))){
428	AliInfo(Form("WARNING: %d is not a valid multiplicity number. \n",ibin+minmultiplicity));
429	continue;
430	}
431	TString suffix;suffix.Form("/Multiplicity:%d",ibin+minmultiplicity);
432	TString zvtx;zvtx.Form("zvtxGT10%s",suffix.Data());
433	TString noPV;noPV.Form("noPrimaryV%s",suffix.Data());
434	TString pV;pV.Form("PrimaryV%s",suffix.Data());
435	TString tbc;tbc.Form("countForNorm%s",suffix.Data());
436	sumnoPV += GetSum(noPV.Data());
437	sumZvtx += GetSum(zvtx.Data());
438	sumPv += GetSum(pV.Data());
439	sumEvtNorm += GetSum(tbc.Data());
440	}
441	Double_t noVtxzGT10 = sumPv>0. ? sumnoPV * sumZvtx / sumPv : 0.;
442	return sumEvtNorm - noVtxzGT10;
443	}
444
445	//___________________________________________________________________________
446	TH1D* AliNormalizationCounter::DrawNEventsForNorm(Bool_t drawRatio){
447	//usare algebra histos
448	fCounters.SortRubric("Run");
449	TString selection;
450
451	selection.Form("event:noPrimaryV");
452	TH1D* hnoPrimV = fCounters.Get("run",selection.Data());
453	hnoPrimV->Sumw2();
454
455	selection.Form("event:zvtxGT10");
456	TH1D* hzvtx= fCounters.Get("run",selection.Data());
457	hzvtx->Sumw2();
458
459	selection.Form("event:PrimaryV");
460	TH1D* hPrimV = fCounters.Get("run",selection.Data());
461	hPrimV->Sumw2();
462
463	hzvtx->Multiply(hnoPrimV);
464	hzvtx->Divide(hPrimV);
465
466	selection.Form("event:countForNorm");
467	TH1D* hCountForNorm = fCounters.Get("run",selection.Data());
468	hCountForNorm->Sumw2();
469
470	hCountForNorm->Add(hzvtx,-1.);
471
472	if(drawRatio){
473	selection.Form("event:triggered");
474	TH1D* htriggered = fCounters.Get("run",selection.Data());
475	htriggered->Sumw2();
476	hCountForNorm->Divide(htriggered);
477	}
478
479	hCountForNorm->DrawClone();
480	return hCountForNorm;
481	}
482
483	//___________________________________________________________________________
484	Int_t AliNormalizationCounter::Multiplicity(AliVEvent* event){
485
486	Int_t multiplicity = 0;
487	AliAODEvent eventAOD = (AliAODEvent)event;
488	AliAODTracklets * aodTracklets = (AliAODTracklets*)eventAOD->GetTracklets();
489	Int_t ntracklets = (Int_t)aodTracklets->GetNumberOfTracklets();
490	for(Int_t i=0;i<ntracklets; i++){
491	Double_t theta = aodTracklets->GetTheta(i);
492	Double_t eta = -TMath::Log( TMath::Tan(theta/2.) ); // check the formula
493	if(TMath::Abs(eta)<fMultiplicityEtaRange){ // set the proper cut on eta
494	multiplicity++;
495	}
496	}
497
498	return multiplicity;
499	}