[u/mrichter/AliRoot.git] / FMD / analysis / AliFMDAnalysisTaskBackgroundCorrection.cxx

 
#include <TROOT.h>
#include <TSystem.h>
#include <TInterpreter.h>
#include <TChain.h>
#include <TFile.h>
#include <TList.h>
#include <iostream>
#include "TH2F.h"
#include "AliFMDAnalysisTaskBackgroundCorrection.h"
#include "AliAnalysisManager.h"
#include "AliESDFMD.h"
#include "AliESDEvent.h"
#include "AliAODEvent.h"
#include "AliAODHandler.h"
#include "AliMCEventHandler.h"
#include "AliStack.h"
#include "AliLog.h"
#include "AliESDVertex.h"
#include "TMath.h"
#include "AliFMDAnaParameters.h"
//#include "AliFMDGeometry.h"

ClassImp(AliFMDAnalysisTaskBackgroundCorrection)


AliFMDAnalysisTaskBackgroundCorrection::AliFMDAnalysisTaskBackgroundCorrection()
: fDebug(0),
  fOutputList(0),
  fInputList(0),
  fHitList(0),
  fArray(0),
  fInputArray(0),
  fVertexString(0x0),
  fNevents(),
  fStandalone(kTRUE),
  fOutputVertexString(0)
{
  // Default constructor
  DefineInput (0, TList::Class());
  DefineOutput(0, TList::Class());
  DefineOutput(1, TObjString::Class());
}
//_____________________________________________________________________
AliFMDAnalysisTaskBackgroundCorrection::AliFMDAnalysisTaskBackgroundCorrection(const char* name, Bool_t SE):
    AliAnalysisTask(name, "Density"),
    fDebug(0),
    fOutputList(0),
    fInputList(0),
    fHitList(0),
    fArray(),
    fInputArray(0),
    fVertexString(0x0),
    fNevents(),
    fStandalone(kTRUE),
    fOutputVertexString(0)
{
  fStandalone = SE;
  if(fStandalone) {
    DefineInput (0, TList::Class());
    DefineOutput(0, TList::Class());
    DefineOutput(1, TObjString::Class());
  }
}
//_____________________________________________________________________
void AliFMDAnalysisTaskBackgroundCorrection::CreateOutputObjects()
{
  AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  
  fArray.SetName("FMD");
  fArray.SetOwner();
  
  if(!fOutputList)
    fOutputList = new TList();
  fOutputList->SetName("BackgroundCorrectedPerEvent");
  if(!fHitList)
    fHitList = new TList();
  fHitList->SetName("HitsList");
  
  
  TH2F* hMult = 0;
  TH2F* hHits = 0;
  Int_t nVtxbins = pars->GetNvtxBins();
  
  for(Int_t det =1; det<=3;det++)
    {
      TObjArray* detArray = new TObjArray();
      detArray->SetName(Form("FMD%d",det));
      fArray.AddAtAndExpand(detArray,det);
      Int_t nRings = (det==1 ? 1 : 2);
      for(Int_t ring = 0;ring<nRings;ring++)
	{
	  Char_t ringChar = (ring == 0 ? 'I' : 'O');
	  Int_t  nSec     = (ring == 0 ? 20 : 40);
	  
	  TObjArray* vtxArray = new TObjArray();
	  vtxArray->SetName(Form("FMD%d%c",det,ringChar));
	  detArray->AddAtAndExpand(vtxArray,ring);
	  for(Int_t i = 0; i< nVtxbins; i++) {
	    TH2F* hBg = pars->GetBackgroundCorrection(det, ringChar, i);
	    hMult  = new TH2F(Form("mult_FMD%d%c_vtxbin%d",det,ringChar,i),Form("mult_FMD%d%c_vtxbin%d",det,ringChar,i),
			      hBg->GetNbinsX(),
			      hBg->GetXaxis()->GetXmin(),
			      hBg->GetXaxis()->GetXmax(),
			      nSec, 0, 2*TMath::Pi());
	    hMult->Sumw2();
	    fOutputList->Add(hMult);
	    hHits  = new TH2F(Form("hits_FMD%d%c_vtxbin%d",det,ringChar,i),Form("hits_FMD%d%c_vtxbin%d",det,ringChar,i),
			      hBg->GetNbinsX(),
			      hBg->GetXaxis()->GetXmin(),
			      hBg->GetXaxis()->GetXmax(),
			      nSec, 0, 2*TMath::Pi());
	    
	    hHits->Sumw2();
	    fHitList->Add(hHits);
	    fOutputList->Add(hHits);
	    vtxArray->AddAtAndExpand(hMult,i);
	    
	  }
	} 
    }
  if(fStandalone) {
    fOutputVertexString = new TObjString();
  }
  fOutputList->Add(fOutputVertexString);
  
  
}
//_____________________________________________________________________
void AliFMDAnalysisTaskBackgroundCorrection::ConnectInputData(Option_t */*option*/)
{
  if(fStandalone) {
    fInputList   = (TList*)GetInputData(0);
    
  }
}
//_____________________________________________________________________
void AliFMDAnalysisTaskBackgroundCorrection::Exec(Option_t */*option*/)
{
  AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  
  fInputArray   = (TObjArray*)fInputList->At(0);
  fVertexString = (TObjString*)fInputList->At(1);
  
  
  Int_t vtxbin   = fVertexString->GetString().Atoi();
  fOutputVertexString->SetString(Form("%d",vtxbin));
  
  //fNevents.operator[](vtxbin)++;
  fNevents.Fill(vtxbin);
  
  //Reset everything
  for(UShort_t det=1;det<=3;det++) {
    TObjArray* detArray = (TObjArray*)fArray.At(det);
    Int_t nRings = (det==1 ? 1 : 2);
    for (UShort_t ir = 0; ir < nRings; ir++) {
      TObjArray* vtxArray = (TObjArray*)detArray->At(ir);
      
      TH2F* hMult   = (TH2F*)vtxArray->At(vtxbin); 
      hMult->Reset();
    }
    
  }
  
  
  for(UShort_t det=1;det<=3;det++) {
    TObjArray* detInputArray = (TObjArray*)fInputArray->At(det);
    TObjArray* detArray = (TObjArray*)fArray.At(det);
    Int_t nRings = (det==1 ? 1 : 2);
    for (UShort_t ir = 0; ir < nRings; ir++) {
      Char_t ringChar = (ir == 0 ? 'I' : 'O');
      TObjArray* vtxInputArray = (TObjArray*)detInputArray->At(ir);
      TObjArray* vtxArray = (TObjArray*)detArray->At(ir);
      TH2F* hMultTotal = (TH2F*)vtxArray->At(vtxbin);
      TH2F* hMultInput = (TH2F*)vtxInputArray->At(vtxbin);
      TH2F* hHits      = (TH2F*)fOutputList->FindObject(Form("hits_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
      
      hHits->Add(hMultInput);
      TH2F* hBg        = pars->GetBackgroundCorrection(det, ringChar, vtxbin);
      
      //TH2F* hTmp       = (TH2F*)hMultInput->Clone("hMult_from_event");
      hMultTotal->Add(hMultInput);
      
      hMultTotal->Divide(hBg);//,"B");
      /*for(Int_t i = 1; i<=hTmp->GetNbinsX();i++) {
	for(Int_t j = 1; j<=hTmp->GetNbinsY();j++) {
	  Float_t mult = hTmp->GetBinContent(i,j);
	  if(mult == 0) continue;
	  Float_t correction = hBg->GetBinContent(i,j);
	  
	  Float_t multcor = mult;
	  if(correction != 0)
	    multcor = multcor/correction;
	  else
	    std::cout<<"Warning! No correction for bin "<<i<<" , "<<j<<std::endl;
	  
	  hTmp->SetBinContent(i,j,multcor);
	}
	}*/
      
      
      // delete hTmp;
      
    }
  }
  if(fStandalone) {
    PostData(0, fOutputList); 
    PostData(1, fOutputVertexString);
  }
  
}
//_____________________________________________________________________
void AliFMDAnalysisTaskBackgroundCorrection::Terminate(Option_t */*option*/) {
  
  /*
  AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  
  Int_t nVtxbins = pars->GetNvtxBins();
  
  for(UShort_t det=1;det<=3;det++) {
    TObjArray* detArray = (TObjArray*)fArray.At(det);
    Int_t nRings = (det==1 ? 1 : 2);
    for (UShort_t ir = 0; ir < nRings; ir++) {
      TObjArray* vtxArray = (TObjArray*)detArray->At(ir);
      for(Int_t i =0; i<nVtxbins; i++) {
	TH2F* hMultTotal = (TH2F*)vtxArray->At(i);
	if(fNevents.At(i))
	  hMultTotal->Scale(1/(Float_t)fNevents.At(i));
      }
    }
  }
  */
}
//_____________________________________________________________________
//
//
// EOF
Commit	Line	Data
3bb122c7	1
	2	#include <TROOT.h>
	3	#include <TSystem.h>
	4	#include <TInterpreter.h>
	5	#include <TChain.h>
	6	#include <TFile.h>
	7	#include <TList.h>
	8	#include <iostream>
	9	#include "TH2F.h"
	10	#include "AliFMDAnalysisTaskBackgroundCorrection.h"
	11	#include "AliAnalysisManager.h"
	12	#include "AliESDFMD.h"
	13	#include "AliESDEvent.h"
	14	#include "AliAODEvent.h"
	15	#include "AliAODHandler.h"
	16	#include "AliMCEventHandler.h"
	17	#include "AliStack.h"
7c3e5162	18	#include "AliLog.h"
3bb122c7	19	#include "AliESDVertex.h"
	20	#include "TMath.h"
	21	#include "AliFMDAnaParameters.h"
78f6f750	22	//#include "AliFMDGeometry.h"
3bb122c7	23
	24	ClassImp(AliFMDAnalysisTaskBackgroundCorrection)
	25
	26
	27	AliFMDAnalysisTaskBackgroundCorrection::AliFMDAnalysisTaskBackgroundCorrection()
	28	: fDebug(0),
3bb122c7	29	fOutputList(0),
7c3e5162	30	fInputList(0),
cc066cb9	31	fHitList(0),
3bb122c7	32	fArray(0),
8dc823cc	33	fInputArray(0),
c78bc12b	34	fVertexString(0x0),
7c3e5162	35	fNevents(),
	36	fStandalone(kTRUE),
	37	fOutputVertexString(0)
3bb122c7	38	{
	39	// Default constructor
	40	DefineInput (0, TList::Class());
8dc7c4c2	41	DefineOutput(0, TList::Class());
7c3e5162	42	DefineOutput(1, TObjString::Class());
3bb122c7	43	}
8dc823cc	44	//_____________________________________________________________________
7c3e5162	45	AliFMDAnalysisTaskBackgroundCorrection::AliFMDAnalysisTaskBackgroundCorrection(const char* name, Bool_t SE):
3bb122c7	46	AliAnalysisTask(name, "Density"),
3bb122c7	47	fDebug(0),
7c3e5162	48	fOutputList(0),
7c3e5162	49	fInputList(0),
cc066cb9	50	fHitList(0),
8dc7c4c2	51	fArray(),
8dc823cc	52	fInputArray(0),
c78bc12b	53	fVertexString(0x0),
7c3e5162	54	fNevents(),
	55	fStandalone(kTRUE),
	56	fOutputVertexString(0)
3bb122c7	57	{
7c3e5162	58	fStandalone = SE;
	59	if(fStandalone) {
	60	DefineInput (0, TList::Class());
	61	DefineOutput(0, TList::Class());
	62	DefineOutput(1, TObjString::Class());
	63	}
3bb122c7	64	}
8dc823cc	65	//_____________________________________________________________________
3bb122c7	66	void AliFMDAnalysisTaskBackgroundCorrection::CreateOutputObjects()
	67	{
	68	AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
3bb122c7	69
8dc7c4c2	70	fArray.SetName("FMD");
8dc7c4c2	71	fArray.SetOwner();
3bb122c7	72
7c3e5162	73	if(!fOutputList)
	74	fOutputList = new TList();
	75	fOutputList->SetName("BackgroundCorrectedPerEvent");
cc066cb9	76	if(!fHitList)
	77	fHitList = new TList();
	78	fHitList->SetName("HitsList");
7c3e5162	79
7c3e5162	80
3bb122c7	81	TH2F* hMult = 0;
bb8a464f	82	TH2F* hHits = 0;
3bb122c7	83	Int_t nVtxbins = pars->GetNvtxBins();
	84
	85	for(Int_t det =1; det<=3;det++)
	86	{
	87	TObjArray* detArray = new TObjArray();
	88	detArray->SetName(Form("FMD%d",det));
8dc7c4c2	89	fArray.AddAtAndExpand(detArray,det);
3bb122c7	90	Int_t nRings = (det==1 ? 1 : 2);
	91	for(Int_t ring = 0;ring<nRings;ring++)
	92	{
	93	Char_t ringChar = (ring == 0 ? 'I' : 'O');
	94	Int_t nSec = (ring == 0 ? 20 : 40);
	95
	96	TObjArray* vtxArray = new TObjArray();
	97	vtxArray->SetName(Form("FMD%d%c",det,ringChar));
	98	detArray->AddAtAndExpand(vtxArray,ring);
	99	for(Int_t i = 0; i< nVtxbins; i++) {
	100	TH2F* hBg = pars->GetBackgroundCorrection(det, ringChar, i);
	101	hMult = new TH2F(Form("mult_FMD%d%c_vtxbin%d",det,ringChar,i),Form("mult_FMD%d%c_vtxbin%d",det,ringChar,i),
	102	hBg->GetNbinsX(),
	103	hBg->GetXaxis()->GetXmin(),
	104	hBg->GetXaxis()->GetXmax(),
	105	nSec, 0, 2*TMath::Pi());
8dc823cc	106	hMult->Sumw2();
7c3e5162	107	fOutputList->Add(hMult);
bb8a464f	108	hHits = new TH2F(Form("hits_FMD%d%c_vtxbin%d",det,ringChar,i),Form("hits_FMD%d%c_vtxbin%d",det,ringChar,i),
	109	hBg->GetNbinsX(),
	110	hBg->GetXaxis()->GetXmin(),
	111	hBg->GetXaxis()->GetXmax(),
	112	nSec, 0, 2*TMath::Pi());
	113
	114	hHits->Sumw2();
cc066cb9	115	fHitList->Add(hHits);
bb8a464f	116	fOutputList->Add(hHits);
3bb122c7	117	vtxArray->AddAtAndExpand(hMult,i);
	118
	119	}
	120	}
	121	}
7c3e5162	122	if(fStandalone) {
	123	fOutputVertexString = new TObjString();
	124	}
	125	fOutputList->Add(fOutputVertexString);
	126
3bb122c7	127
bb8a464f	128
3bb122c7	129	}
8dc823cc	130	//_____________________________________________________________________
3bb122c7	131	void AliFMDAnalysisTaskBackgroundCorrection::ConnectInputData(Option_t /option*/)
3bb122c7	132	{
7c3e5162	133	if(fStandalone) {
	134	fInputList = (TList*)GetInputData(0);
	135
	136	}
3bb122c7	137	}
8dc823cc	138	//_____________________________________________________________________
3bb122c7	139	void AliFMDAnalysisTaskBackgroundCorrection::Exec(Option_t /option*/)
	140	{
	141	AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
	142
7c3e5162	143	fInputArray = (TObjArray*)fInputList->At(0);
	144	fVertexString = (TObjString*)fInputList->At(1);
	145
	146
8dc823cc	147	Int_t vtxbin = fVertexString->GetString().Atoi();
7c3e5162	148	fOutputVertexString->SetString(Form("%d",vtxbin));
8dc7c4c2	149
875e75b6	150	//fNevents.operator[](vtxbin)++;
875e75b6	151	fNevents.Fill(vtxbin);
3bb122c7	152
7c3e5162	153	//Reset everything
	154	for(UShort_t det=1;det<=3;det++) {
	155	TObjArray* detArray = (TObjArray*)fArray.At(det);
	156	Int_t nRings = (det==1 ? 1 : 2);
	157	for (UShort_t ir = 0; ir < nRings; ir++) {
	158	TObjArray* vtxArray = (TObjArray*)detArray->At(ir);
	159
	160	TH2F* hMult = (TH2F*)vtxArray->At(vtxbin);
	161	hMult->Reset();
	162	}
	163
	164	}
	165
	166
	167
3bb122c7	168	for(UShort_t det=1;det<=3;det++) {
3bb122c7	169	TObjArray* detInputArray = (TObjArray*)fInputArray->At(det);
8dc7c4c2	170	TObjArray* detArray = (TObjArray*)fArray.At(det);
3bb122c7	171	Int_t nRings = (det==1 ? 1 : 2);
	172	for (UShort_t ir = 0; ir < nRings; ir++) {
	173	Char_t ringChar = (ir == 0 ? 'I' : 'O');
	174	TObjArray* vtxInputArray = (TObjArray*)detInputArray->At(ir);
	175	TObjArray* vtxArray = (TObjArray*)detArray->At(ir);
8dc823cc	176	TH2F* hMultTotal = (TH2F*)vtxArray->At(vtxbin);
7c3e5162	177	TH2F* hMultInput = (TH2F*)vtxInputArray->At(vtxbin);
bb8a464f	178	TH2F* hHits = (TH2F*)fOutputList->FindObject(Form("hits_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
	179
	180	hHits->Add(hMultInput);
8dc823cc	181	TH2F* hBg = pars->GetBackgroundCorrection(det, ringChar, vtxbin);
7c3e5162	182
78f6f750	183	//TH2F* hTmp = (TH2F*)hMultInput->Clone("hMult_from_event");
78f6f750	184	hMultTotal->Add(hMultInput);
bb8a464f	185
78f6f750	186	hMultTotal->Divide(hBg);//,"B");
6289b3e8	187	/*for(Int_t i = 1; i<=hTmp->GetNbinsX();i++) {
	188	for(Int_t j = 1; j<=hTmp->GetNbinsY();j++) {
	189	Float_t mult = hTmp->GetBinContent(i,j);
	190	if(mult == 0) continue;
	191	Float_t correction = hBg->GetBinContent(i,j);
	192
	193	Float_t multcor = mult;
	194	if(correction != 0)
	195	multcor = multcor/correction;
	196	else
	197	std::cout<<"Warning! No correction for bin "<<i<<" , "<<j<<std::endl;
	198
	199	hTmp->SetBinContent(i,j,multcor);
	200	}
	201	}*/
78f6f750	202
	203
	204
	205	// delete hTmp;
8dc823cc	206
3bb122c7	207	}
3bb122c7	208	}
7c3e5162	209	if(fStandalone) {
	210	PostData(0, fOutputList);
	211	PostData(1, fOutputVertexString);
	212	}
3bb122c7	213
3bb122c7	214	}
8dc823cc	215	//_____________________________________________________________________
	216	void AliFMDAnalysisTaskBackgroundCorrection::Terminate(Option_t /option*/) {
	217
875e75b6	218	/*
8dc823cc	219	AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
	220
	221	Int_t nVtxbins = pars->GetNvtxBins();
	222
	223	for(UShort_t det=1;det<=3;det++) {
8dc7c4c2	224	TObjArray* detArray = (TObjArray*)fArray.At(det);
8dc823cc	225	Int_t nRings = (det==1 ? 1 : 2);
	226	for (UShort_t ir = 0; ir < nRings; ir++) {
	227	TObjArray* vtxArray = (TObjArray*)detArray->At(ir);
	228	for(Int_t i =0; i<nVtxbins; i++) {
	229	TH2F* hMultTotal = (TH2F*)vtxArray->At(i);
8dc7c4c2	230	if(fNevents.At(i))
8dc7c4c2	231	hMultTotal->Scale(1/(Float_t)fNevents.At(i));
8dc823cc	232	}
	233	}
	234	}
875e75b6	235	*/
8dc823cc	236	}
	237	//_____________________________________________________________________
	238	//
	239	//
	240	// EOF