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

#include "AliFMDAnalysisTaskGenerateCorrection.h"
#include "AliESDEvent.h"
#include "iostream"
#include "AliESDFMD.h"
#include "TH2F.h"
#include "AliTrackReference.h"
#include "AliStack.h"
#include "AliFMDAnaParameters.h"
#include "AliFMDStripIndex.h"
#include "AliStack.h"
#include "AliMCParticle.h"
#include "AliMCEvent.h"
//#include "AliFMDGeometry.h"
#include "TArray.h"
#include "AliGenEventHeader.h"
#include "AliHeader.h"
#include "AliFMDAnaCalibBackgroundCorrection.h"
#include "AliFMDAnaCalibEventSelectionEfficiency.h"
//#include "AliCDBManager.h"
//#include "AliCDBId.h"
//#include "AliCDBMetaData.h"
#include "TSystem.h"
#include "TROOT.h"
#include "TAxis.h"
ClassImp(AliFMDAnalysisTaskGenerateCorrection)

//_____________________________________________________________________
AliFMDAnalysisTaskGenerateCorrection::AliFMDAnalysisTaskGenerateCorrection():
AliAnalysisTaskSE(),
  fListOfHits(), 
  fListOfPrimaries(),
  fListOfCorrection(),
  fVertexBins(),
  fLastTrackByStrip(0),
  fHitsByStrip(0),
  fZvtxCut(10),
  fNvtxBins(10),
  fNbinsEta(200),
  fBackground(0),
  fEventSelectionEff(0)
{
  // Default constructor
}
//_____________________________________________________________________
AliFMDAnalysisTaskGenerateCorrection::AliFMDAnalysisTaskGenerateCorrection(const char* name):
  AliAnalysisTaskSE(name),
  fListOfHits(), 
  fListOfPrimaries(),
  fListOfCorrection(),
  fVertexBins(),
  fLastTrackByStrip(0),
  fHitsByStrip(0),
  fZvtxCut(10),
  fNvtxBins(10),
  fNbinsEta(200),
  fBackground(0),
  fEventSelectionEff(0)
{
 
  DefineOutput(1, TList::Class());
  DefineOutput(2, TList::Class());
  DefineOutput(3, TH1F::Class());
  DefineOutput(4, TList::Class());
}
//_____________________________________________________________________
void AliFMDAnalysisTaskGenerateCorrection::UserCreateOutputObjects()
{
// Create the output containers
//
  
  std::cout<<"Creating output objects"<<std::endl;
  for(Int_t iring = 0; iring<2;iring++) {
    Char_t ringChar = (iring == 0 ? 'I' : 'O');
    Int_t nSec = (iring == 1 ? 40 : 20);
    for(Int_t v=0; v<fNvtxBins;v++) {

      TH2F* hPrimary       = new TH2F(Form("hPrimary_FMD_%c_vtx%d",ringChar,v),
				      Form("hPrimary_FMD_%c_vtx%d",ringChar,v),
				      fNbinsEta, -6,6, nSec, 0,2*TMath::Pi());
      hPrimary->Sumw2();
      fListOfPrimaries.Add(hPrimary);
    }
  }
  
  
  
  for(Int_t det =1; det<=3;det++) {
    Int_t nRings = (det==1 ? 1 : 2);
    for(Int_t ring = 0;ring<nRings;ring++) {
      Int_t nSec = (ring == 1 ? 40 : 20);
      Char_t ringChar = (ring == 0 ? 'I' : 'O');
      TH1F* doubleHits = new TH1F(Form("DoubleHits_FMD%d%c",det,ringChar),
				  Form("DoubleHits_FMD%d%c",det,ringChar),fNbinsEta, -6,6);
      TH1F* allHits = new TH1F(Form("allHits_FMD%d%c",det,ringChar),
			       Form("allHits_FMD%d%c",det,ringChar), fNbinsEta, -6,6);
      
      doubleHits->Sumw2();
      allHits->Sumw2();
      fListOfHits.Add(allHits);
      fListOfHits.Add(doubleHits);
	
      for(Int_t v=0; v<fNvtxBins;v++) {
	TH2F* hHits = new TH2F(Form("hHits_FMD%d%c_vtx%d", det,ringChar,v),
			       Form("hHits_FMD%d%c_vtx%d", det,ringChar,v),
			       fNbinsEta, -6,6, nSec, 0, 2*TMath::Pi());
	hHits->Sumw2();
	fListOfHits.Add(hHits);
	
      } 
    }
  }
  
  TH1F* hEventsSelected  = new TH1F("EventsSelected","EventsSelected",fNvtxBins,0,fNvtxBins);
  TH1F* hEventsAll    = new TH1F("EventsAll","EventsAll",fNvtxBins,0,fNvtxBins);
  //  TH1F* hTriggered    = new TH1F("Triggered","Triggered",fNvtxBins,0,fNvtxBins);
  // TH1F* hTriggeredAll = new TH1F("TriggeredAll","TriggeredAll",fNvtxBins,0,fNvtxBins);
  hEventsSelected->Sumw2();
  hEventsAll->Sumw2();
  fListOfHits.Add(hEventsSelected);
  fListOfPrimaries.Add(hEventsAll);
  
  // fListOfHits.Add(hTriggered);
  //  fListOfPrimaries.Add(hTriggeredAll);
  
  fVertexBins.SetName("VertexBins");
  fVertexBins.GetXaxis()->Set(fNvtxBins,-1*fZvtxCut,fZvtxCut);
  
}
//_____________________________________________________________________
void AliFMDAnalysisTaskGenerateCorrection::Init()
{
  fLastTrackByStrip.Reset(-1);
  
  
}
//_____________________________________________________________________
void AliFMDAnalysisTaskGenerateCorrection::UserExec(Option_t */*option*/)
{
  
  fLastTrackByStrip.Reset(-1);
  fHitsByStrip.Reset(0);
  AliMCEvent* mcevent = MCEvent();
  
  AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  
  
  AliESDEvent* esdevent = (AliESDEvent*)InputEvent();
  Double_t esdvertex[3];
  pars->GetVertex(esdevent,esdvertex);
  
  AliMCParticle* particle = 0;
  AliStack* stack = mcevent->Stack();
  
  UShort_t det,sec,strip;
  Char_t   ring;
  
  Int_t nTracks                = mcevent->GetNumberOfTracks();
  AliHeader* header            = mcevent->Header();
  AliGenEventHeader* genHeader = header->GenEventHeader();
  
  TArrayF vertex;
  genHeader->PrimaryVertex(vertex);
  
  if(TMath::Abs(vertex.At(2)) > fZvtxCut)
    return;
  
  Double_t delta           = 2*fZvtxCut/fNvtxBins;
  Double_t vertexBinDouble = (vertex.At(2) + fZvtxCut) / delta;
  Int_t    vertexBin       = (Int_t)vertexBinDouble;
  
  // Vertex determination correction
  TH1F* hEventsSelected    = (TH1F*)fListOfHits.FindObject("EventsSelected");
  TH1F* hEventsAll         = (TH1F*)fListOfPrimaries.FindObject("EventsAll");
  // TH1F* hTriggered      = (TH1F*)fListOfHits.FindObject("Triggered");
  //  TH1F* hTriggeredAll   = (TH1F*)fListOfPrimaries.FindObject("TriggeredAll");
  
  Bool_t vtxFound = kTRUE;
  if(esdvertex[0] == 0 && esdvertex[1] == 0 && esdvertex[2] == 0)
    vtxFound = kFALSE;
    
  Bool_t isTriggered = pars->IsEventTriggered(esdevent);
  
  if(vtxFound && isTriggered) {
    //hTriggered->Fill(vertexBin);
    hEventsSelected->Fill(vertexBin);
  }
  
  //if(isTriggered)
  
  hEventsAll->Fill(vertexBin);
  //if(vtxFound)
  //   hTriggeredAll->Fill(vertexBin);
  
  for(Int_t i = 0 ;i<nTracks;i++) {
    particle = mcevent->GetTrack(i);
    
    if(!particle)
      continue;
        
    if(stack->IsPhysicalPrimary(i) && particle->Charge() != 0) {
      
      
      TH2F* hPrimaryInner = (TH2F*)fListOfPrimaries.FindObject( Form("hPrimary_FMD_%c_vtx%d",'I',vertexBin));
      TH2F* hPrimaryOuter = (TH2F*)fListOfPrimaries.FindObject( Form("hPrimary_FMD_%c_vtx%d",'O',vertexBin));
      hPrimaryInner->Fill(particle->Eta(),particle->Phi());
      hPrimaryOuter->Fill(particle->Eta(),particle->Phi());      
    }
    
    for(Int_t j=0; j<particle->GetNumberOfTrackReferences();j++) {
      
      AliTrackReference* ref = particle->GetTrackReference(j);
      
      if(ref->DetectorId() != AliTrackReference::kFMD)
	continue;
      AliFMDStripIndex::Unpack(ref->UserId(),det,ring,sec,strip);
      Float_t thisStripTrack = fLastTrackByStrip(det,ring,sec,strip);
      if(particle->Charge() != 0 && i != thisStripTrack ) {
	
	Float_t phi = pars->GetPhiFromSector(det,ring,sec);
	Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex.At(2));
	
	TH2F* hHits = (TH2F*)fListOfHits.FindObject(Form("hHits_FMD%d%c_vtx%d", det,ring,vertexBin));
	hHits->Fill(eta,phi);
	Float_t nstrips = (ring =='O' ? 256 : 512);
	fHitsByStrip(det,ring,sec,strip) +=1;
	TH1F* allHits = (TH1F*)fListOfHits.FindObject(Form("allHits_FMD%d%c",det,ring));
	TH1F* doubleHits = (TH1F*)fListOfHits.FindObject(Form("DoubleHits_FMD%d%c",det,ring));
	
	if(fHitsByStrip(det,ring,sec,strip) == 1)
	  allHits->Fill(eta);
	
	doubleHits->Fill(eta);
		
	fLastTrackByStrip(det,ring,sec,strip) = (Float_t)i;
	if(strip >0)
	  fLastTrackByStrip(det,ring,sec,strip-1) = (Float_t)i;
	if(strip < (nstrips - 1))
	  fLastTrackByStrip(det,ring,sec,strip+1) = (Float_t)i;
      }
    }

  }
    
  	
  PostData(1, &fListOfHits);
  PostData(2, &fListOfPrimaries);
  PostData(3, &fVertexBins);
}
//_____________________________________________________________________
void AliFMDAnalysisTaskGenerateCorrection::Terminate(Option_t */*option*/)
{
  /*  TH1F* allHits = (TH1F*)fListOfHits.FindObject("allHits");
  TH1F* doubleHits = (TH1F*)fListOfHits.FindObject("DoubleHits");
  
  doubleHits->Divide(allHits);
  GenerateCorrection();
  PostData(1, &fListOfHits);
  PostData(4, &fListOfCorrection);*/
  
}
//_____________________________________________________________________
void AliFMDAnalysisTaskGenerateCorrection::GenerateCorrection() {
  
  fBackground         = new AliFMDAnaCalibBackgroundCorrection();
  fEventSelectionEff  = new AliFMDAnaCalibEventSelectionEfficiency();
  
  //TH1F* hTriggered      = (TH1F*)fListOfHits.FindObject("Triggered");
  //TH1F* hTriggeredAll   = (TH1F*)fListOfPrimaries.FindObject("TriggeredAll");
  
  TH1F* hEventsSelected    = (TH1F*)fListOfHits.FindObject("EventsSelected");
  TH1F* hEventsAll         = (TH1F*)fListOfPrimaries.FindObject("EventsAll");
  
  //hEventsSelected->Divide(hEventsSelected,hEventsAll,1,1,"B");
  
  for(Int_t i = 1; i<=hEventsSelected->GetNbinsX(); i++) {
    if(hEventsSelected->GetBinContent(i) == 0 )
      continue;
    Float_t a    = hEventsSelected->GetBinContent(i);
    Float_t da   = hEventsSelected->GetBinError(i);
    Float_t sum  = hEventsAll->GetBinContent(i);
    Float_t dsum = hEventsAll->GetBinError(i);
    Float_t b    = sum-a;
    Float_t db   = TMath::Sqrt(TMath::Power(da,2) + TMath::Power(dsum,2));
    
    Float_t cor  = a / sum;
    Float_t ecor = TMath::Sqrt(TMath::Power(b*da,2) + TMath::Power(a*db,2)) / TMath::Power(sum,2);
    
    hEventsSelected->SetBinContent(i,cor);
    hEventsSelected->SetBinError(i,ecor);
    
  }
  
  fEventSelectionEff->SetCorrection(hEventsSelected);
  
  for(Int_t det= 1; det <=3; det++) {
    Int_t nRings = (det==1 ? 1 : 2);
    
    for(Int_t iring = 0; iring<nRings; iring++) {
      Char_t ring = (iring == 0 ? 'I' : 'O');
      TH1F* allHits = (TH1F*)fListOfHits.FindObject(Form("allHits_FMD%d%c",det,ring));
      TH1F* doubleHits = (TH1F*)fListOfHits.FindObject(Form("DoubleHits_FMD%d%c",det,ring));
      fBackground->SetDoubleHitCorrection(det,ring,doubleHits);
      doubleHits->Divide(allHits);
      for(Int_t vertexBin=0;vertexBin<fNvtxBins  ;vertexBin++) {
	TH2F* hHits          = (TH2F*)fListOfHits.FindObject(Form("hHits_FMD%d%c_vtx%d", det,ring,vertexBin));
	TH2F* hPrimary  = (TH2F*)fListOfPrimaries.FindObject( Form("hPrimary_FMD_%c_vtx%d",ring,vertexBin));
	TH2F* hCorrection = (TH2F*)hHits->Clone(Form("FMD%d%c_vtxbin_%d_correction",det,ring,vertexBin));
	hCorrection->Divide(hPrimary);
	/*for(Int_t i = 1; i<=hCorrection->GetNbinsX();i++)  {
	  for(Int_t j=1; j<=hCorrection->GetNbinsY();j++) {
	    if(hCorrection()->GetBinContent(i,j) == 0)
	      continue;
	    Float_t a = h 
	    
	    
	  }
	}
	*/
	
	hCorrection->SetTitle(hCorrection->GetName());
	fListOfCorrection.Add(hCorrection);
	fBackground->SetBgCorrection(det,ring,vertexBin,hCorrection);
	
	
      }
      
    }
  }
  TAxis refAxis(fNvtxBins,-1*fZvtxCut,fZvtxCut);
  fBackground->SetRefAxis(&refAxis);

}
//_____________________________________________________________________
void AliFMDAnalysisTaskGenerateCorrection::ReadFromFile(const Char_t* filename, Bool_t storeInOCDB, Int_t /*runNo*/) {
  
  TFile infile(filename);
  TH1F* hVertex = (TH1F*)infile.Get("VertexBins");
  fZvtxCut = hVertex->GetXaxis()->GetXmax();
  fNvtxBins = hVertex->GetXaxis()->GetNbins();
  fVertexBins.SetName("VertexBins");
  fVertexBins.GetXaxis()->Set(fNvtxBins,-1*fZvtxCut,fZvtxCut);
  
  TList* listOfHits = (TList*)infile.Get("Hits");
  TList* listOfPrim = (TList*)infile.Get("Primaries");
  
  TH1F* hEventsSelected    = (TH1F*)listOfHits->FindObject("EventsSelected");
  TH1F* hEventsAll         = (TH1F*)listOfPrim->FindObject("EventsAll");
  fListOfHits.Add(hEventsSelected);
  fListOfPrimaries.Add(hEventsAll);
  
  for(Int_t det =1; det<=3;det++)
      {
	Int_t nRings = (det==1 ? 1 : 2);
	for(Int_t ring = 0;ring<nRings;ring++)
	  {
	    Char_t ringChar = (ring == 0 ? 'I' : 'O');
	    TH1F* allHits = (TH1F*)listOfHits->FindObject(Form("allHits_FMD%d%c",det,ringChar));
	    TH1F* doubleHits = (TH1F*)listOfHits->FindObject(Form("DoubleHits_FMD%d%c",det,ringChar));
	    fListOfHits.Add(allHits);
	    fListOfHits.Add(doubleHits);
	    for(Int_t v=0; v<fNvtxBins;v++)
	      {
		
		TH2F* hHits          = (TH2F*)listOfHits->FindObject(Form("hHits_FMD%d%c_vtx%d", det,ringChar,v));
		fListOfHits.Add(hHits);
	      }
	  }
      }
  for(Int_t iring = 0; iring<2;iring++) {
    Char_t ringChar = (iring == 0 ? 'I' : 'O');
    for(Int_t v=0; v<fNvtxBins;v++) {
      
      TH2F* hPrimary       = (TH2F*)listOfPrim->FindObject( Form("hPrimary_FMD_%c_vtx%d",ringChar,v));
      fListOfPrimaries.Add(hPrimary);
      
    }
  }
  GenerateCorrection();
  
  TFile fout("backgroundFromFile.root","recreate");
  fListOfHits.Write();
  fListOfPrimaries.Write();
  fListOfCorrection.Write();
  fVertexBins.Write();
  
  fout.Close();
  AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  if(storeInOCDB) {
    TFile fbg(pars->GetPath(pars->GetBackgroundID()),"RECREATE");
    fBackground->Write(AliFMDAnaParameters::GetBackgroundID());
    fbg.Close();
    TFile feselect(pars->GetPath(pars->GetEventSelectionEffID()),"RECREATE");
    fEventSelectionEff->Write(AliFMDAnaParameters::GetEventSelectionEffID());
    feselect.Close();
    
  }
  
  
  
  
}
//_____________________________________________________________________
//
// EOF
//
Commit	Line	Data
d31dce33	1	#include "AliFMDAnalysisTaskGenerateCorrection.h"
	2	#include "AliESDEvent.h"
	3	#include "iostream"
	4	#include "AliESDFMD.h"
	5	#include "TH2F.h"
	6	#include "AliTrackReference.h"
	7	#include "AliStack.h"
	8	#include "AliFMDAnaParameters.h"
	9	#include "AliFMDStripIndex.h"
	10	#include "AliStack.h"
	11	#include "AliMCParticle.h"
	12	#include "AliMCEvent.h"
	13	//#include "AliFMDGeometry.h"
	14	#include "TArray.h"
	15	#include "AliGenEventHeader.h"
	16	#include "AliHeader.h"
	17	#include "AliFMDAnaCalibBackgroundCorrection.h"
	18	#include "AliFMDAnaCalibEventSelectionEfficiency.h"
	19	//#include "AliCDBManager.h"
	20	//#include "AliCDBId.h"
	21	//#include "AliCDBMetaData.h"
	22	#include "TSystem.h"
	23	#include "TROOT.h"
	24	#include "TAxis.h"
	25	ClassImp(AliFMDAnalysisTaskGenerateCorrection)
	26
	27	//_____________________________________________________________________
	28	AliFMDAnalysisTaskGenerateCorrection::AliFMDAnalysisTaskGenerateCorrection():
	29	AliAnalysisTaskSE(),
	30	fListOfHits(),
	31	fListOfPrimaries(),
	32	fListOfCorrection(),
	33	fVertexBins(),
	34	fLastTrackByStrip(0),
	35	fHitsByStrip(0),
	36	fZvtxCut(10),
	37	fNvtxBins(10),
	38	fNbinsEta(200),
9f55be54	39	fBackground(0),
9f55be54	40	fEventSelectionEff(0)
d31dce33	41	{
	42	// Default constructor
	43	}
	44	//_____________________________________________________________________
	45	AliFMDAnalysisTaskGenerateCorrection::AliFMDAnalysisTaskGenerateCorrection(const char* name):
	46	AliAnalysisTaskSE(name),
	47	fListOfHits(),
	48	fListOfPrimaries(),
	49	fListOfCorrection(),
	50	fVertexBins(),
	51	fLastTrackByStrip(0),
	52	fHitsByStrip(0),
	53	fZvtxCut(10),
	54	fNvtxBins(10),
	55	fNbinsEta(200),
9f55be54	56	fBackground(0),
9f55be54	57	fEventSelectionEff(0)
d31dce33	58	{
	59
	60	DefineOutput(1, TList::Class());
	61	DefineOutput(2, TList::Class());
	62	DefineOutput(3, TH1F::Class());
	63	DefineOutput(4, TList::Class());
	64	}
	65	//_____________________________________________________________________
	66	void AliFMDAnalysisTaskGenerateCorrection::UserCreateOutputObjects()
	67	{
	68	// Create the output containers
	69	//
	70
	71	std::cout<<"Creating output objects"<<std::endl;
	72	for(Int_t iring = 0; iring<2;iring++) {
	73	Char_t ringChar = (iring == 0 ? 'I' : 'O');
	74	Int_t nSec = (iring == 1 ? 40 : 20);
	75	for(Int_t v=0; v<fNvtxBins;v++) {
	76
	77	TH2F* hPrimary = new TH2F(Form("hPrimary_FMD_%c_vtx%d",ringChar,v),
	78	Form("hPrimary_FMD_%c_vtx%d",ringChar,v),
	79	fNbinsEta, -6,6, nSec, 0,2*TMath::Pi());
	80	hPrimary->Sumw2();
	81	fListOfPrimaries.Add(hPrimary);
	82	}
	83	}
	84
	85
	86
	87	for(Int_t det =1; det<=3;det++) {
	88	Int_t nRings = (det==1 ? 1 : 2);
	89	for(Int_t ring = 0;ring<nRings;ring++) {
	90	Int_t nSec = (ring == 1 ? 40 : 20);
	91	Char_t ringChar = (ring == 0 ? 'I' : 'O');
	92	TH1F* doubleHits = new TH1F(Form("DoubleHits_FMD%d%c",det,ringChar),
	93	Form("DoubleHits_FMD%d%c",det,ringChar),fNbinsEta, -6,6);
	94	TH1F* allHits = new TH1F(Form("allHits_FMD%d%c",det,ringChar),
	95	Form("allHits_FMD%d%c",det,ringChar), fNbinsEta, -6,6);
	96
	97	doubleHits->Sumw2();
	98	allHits->Sumw2();
	99	fListOfHits.Add(allHits);
	100	fListOfHits.Add(doubleHits);
	101
	102	for(Int_t v=0; v<fNvtxBins;v++) {
	103	TH2F* hHits = new TH2F(Form("hHits_FMD%d%c_vtx%d", det,ringChar,v),
	104	Form("hHits_FMD%d%c_vtx%d", det,ringChar,v),
	105	fNbinsEta, -6,6, nSec, 0, 2*TMath::Pi());
	106	hHits->Sumw2();
	107	fListOfHits.Add(hHits);
	108
	109	}
	110	}
	111	}
	112
	113	TH1F* hEventsSelected = new TH1F("EventsSelected","EventsSelected",fNvtxBins,0,fNvtxBins);
	114	TH1F* hEventsAll = new TH1F("EventsAll","EventsAll",fNvtxBins,0,fNvtxBins);
	115	// TH1F* hTriggered = new TH1F("Triggered","Triggered",fNvtxBins,0,fNvtxBins);
	116	// TH1F* hTriggeredAll = new TH1F("TriggeredAll","TriggeredAll",fNvtxBins,0,fNvtxBins);
0b0a4ae5	117	hEventsSelected->Sumw2();
0b0a4ae5	118	hEventsAll->Sumw2();
d31dce33	119	fListOfHits.Add(hEventsSelected);
	120	fListOfPrimaries.Add(hEventsAll);
	121
	122	// fListOfHits.Add(hTriggered);
	123	// fListOfPrimaries.Add(hTriggeredAll);
	124
	125	fVertexBins.SetName("VertexBins");
	126	fVertexBins.GetXaxis()->Set(fNvtxBins,-1*fZvtxCut,fZvtxCut);
	127
	128	}
	129	//_____________________________________________________________________
	130	void AliFMDAnalysisTaskGenerateCorrection::Init()
	131	{
	132	fLastTrackByStrip.Reset(-1);
	133
	134
	135	}
	136	//_____________________________________________________________________
	137	void AliFMDAnalysisTaskGenerateCorrection::UserExec(Option_t /option*/)
	138	{
	139
	140	fLastTrackByStrip.Reset(-1);
	141	fHitsByStrip.Reset(0);
	142	AliMCEvent* mcevent = MCEvent();
	143
	144	AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
	145
	146
	147	AliESDEvent* esdevent = (AliESDEvent*)InputEvent();
	148	Double_t esdvertex[3];
	149	pars->GetVertex(esdevent,esdvertex);
	150
	151	AliMCParticle* particle = 0;
	152	AliStack* stack = mcevent->Stack();
	153
	154	UShort_t det,sec,strip;
	155	Char_t ring;
	156
	157	Int_t nTracks = mcevent->GetNumberOfTracks();
	158	AliHeader* header = mcevent->Header();
	159	AliGenEventHeader* genHeader = header->GenEventHeader();
	160
	161	TArrayF vertex;
	162	genHeader->PrimaryVertex(vertex);
	163
	164	if(TMath::Abs(vertex.At(2)) > fZvtxCut)
	165	return;
	166
	167	Double_t delta = 2*fZvtxCut/fNvtxBins;
	168	Double_t vertexBinDouble = (vertex.At(2) + fZvtxCut) / delta;
	169	Int_t vertexBin = (Int_t)vertexBinDouble;
	170
	171	// Vertex determination correction
	172	TH1F* hEventsSelected = (TH1F*)fListOfHits.FindObject("EventsSelected");
	173	TH1F* hEventsAll = (TH1F*)fListOfPrimaries.FindObject("EventsAll");
	174	// TH1F* hTriggered = (TH1F*)fListOfHits.FindObject("Triggered");
	175	// TH1F* hTriggeredAll = (TH1F*)fListOfPrimaries.FindObject("TriggeredAll");
	176
	177	Bool_t vtxFound = kTRUE;
	178	if(esdvertex[0] == 0 && esdvertex[1] == 0 && esdvertex[2] == 0)
	179	vtxFound = kFALSE;
	180
	181	Bool_t isTriggered = pars->IsEventTriggered(esdevent);
	182
183	if(vtxFound && isTriggered) {
184	//hTriggered->Fill(vertexBin);
185	hEventsSelected->Fill(vertexBin);
186	}
187
188	//if(isTriggered)
189
190	hEventsAll->Fill(vertexBin);
191	//if(vtxFound)
192	// hTriggeredAll->Fill(vertexBin);
193
194	for(Int_t i = 0 ;i<nTracks;i++) {
195	particle = mcevent->GetTrack(i);
196
197	if(!particle)
198	continue;
199
200	if(stack->IsPhysicalPrimary(i) && particle->Charge() != 0) {
201
202
203	TH2F* hPrimaryInner = (TH2F*)fListOfPrimaries.FindObject( Form("hPrimary_FMD_%c_vtx%d",'I',vertexBin));
204	TH2F* hPrimaryOuter = (TH2F*)fListOfPrimaries.FindObject( Form("hPrimary_FMD_%c_vtx%d",'O',vertexBin));
205	hPrimaryInner->Fill(particle->Eta(),particle->Phi());
206	hPrimaryOuter->Fill(particle->Eta(),particle->Phi());
207	}
208
209	for(Int_t j=0; j<particle->GetNumberOfTrackReferences();j++) {
210
211	AliTrackReference* ref = particle->GetTrackReference(j);
212
213	if(ref->DetectorId() != AliTrackReference::kFMD)
214	continue;
215	AliFMDStripIndex::Unpack(ref->UserId(),det,ring,sec,strip);
216	Float_t thisStripTrack = fLastTrackByStrip(det,ring,sec,strip);
217	if(particle->Charge() != 0 && i != thisStripTrack ) {
218
219	Float_t phi = pars->GetPhiFromSector(det,ring,sec);
220	Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex.At(2));
221
222	TH2F* hHits = (TH2F*)fListOfHits.FindObject(Form("hHits_FMD%d%c_vtx%d", det,ring,vertexBin));
223	hHits->Fill(eta,phi);
224	Float_t nstrips = (ring =='O' ? 256 : 512);
225	fHitsByStrip(det,ring,sec,strip) +=1;
226	TH1F* allHits = (TH1F*)fListOfHits.FindObject(Form("allHits_FMD%d%c",det,ring));
227	TH1F* doubleHits = (TH1F*)fListOfHits.FindObject(Form("DoubleHits_FMD%d%c",det,ring));
228
229	if(fHitsByStrip(det,ring,sec,strip) == 1)
230	allHits->Fill(eta);
231
232	doubleHits->Fill(eta);
233
234	fLastTrackByStrip(det,ring,sec,strip) = (Float_t)i;
235	if(strip >0)
236	fLastTrackByStrip(det,ring,sec,strip-1) = (Float_t)i;
237	if(strip < (nstrips - 1))
238	fLastTrackByStrip(det,ring,sec,strip+1) = (Float_t)i;
239	}
240	}
241
242	}
243
244
245	PostData(1, &fListOfHits);
246	PostData(2, &fListOfPrimaries);
247	PostData(3, &fVertexBins);
248	}
249	//_____________________________________________________________________
250	void AliFMDAnalysisTaskGenerateCorrection::Terminate(Option_t /option*/)
251	{
252	/* TH1F* allHits = (TH1F*)fListOfHits.FindObject("allHits");
253	TH1F* doubleHits = (TH1F*)fListOfHits.FindObject("DoubleHits");
254
255	doubleHits->Divide(allHits);
256	GenerateCorrection();
257	PostData(1, &fListOfHits);
258	PostData(4, &fListOfCorrection);*/
259
260	}
261	//_____________________________________________________________________
262	void AliFMDAnalysisTaskGenerateCorrection::GenerateCorrection() {
263
264	fBackground = new AliFMDAnaCalibBackgroundCorrection();
265	fEventSelectionEff = new AliFMDAnaCalibEventSelectionEfficiency();
266
267	//TH1F* hTriggered = (TH1F*)fListOfHits.FindObject("Triggered");
268	//TH1F* hTriggeredAll = (TH1F*)fListOfPrimaries.FindObject("TriggeredAll");
269
270	TH1F* hEventsSelected = (TH1F*)fListOfHits.FindObject("EventsSelected");
271	TH1F* hEventsAll = (TH1F*)fListOfPrimaries.FindObject("EventsAll");
272
0b0a4ae5	273	//hEventsSelected->Divide(hEventsSelected,hEventsAll,1,1,"B");
	274
	275	for(Int_t i = 1; i<=hEventsSelected->GetNbinsX(); i++) {
	276	if(hEventsSelected->GetBinContent(i) == 0 )
	277	continue;
	278	Float_t a = hEventsSelected->GetBinContent(i);
	279	Float_t da = hEventsSelected->GetBinError(i);
	280	Float_t sum = hEventsAll->GetBinContent(i);
	281	Float_t dsum = hEventsAll->GetBinError(i);
	282	Float_t b = sum-a;
	283	Float_t db = TMath::Sqrt(TMath::Power(da,2) + TMath::Power(dsum,2));
	284
	285	Float_t cor = a / sum;
	286	Float_t ecor = TMath::Sqrt(TMath::Power(bda,2) + TMath::Power(adb,2)) / TMath::Power(sum,2);
	287
	288	hEventsSelected->SetBinContent(i,cor);
	289	hEventsSelected->SetBinError(i,ecor);
	290
	291	}
d31dce33	292
	293	fEventSelectionEff->SetCorrection(hEventsSelected);
	294
	295	for(Int_t det= 1; det <=3; det++) {
	296	Int_t nRings = (det==1 ? 1 : 2);
	297
	298	for(Int_t iring = 0; iring<nRings; iring++) {
	299	Char_t ring = (iring == 0 ? 'I' : 'O');
	300	TH1F* allHits = (TH1F*)fListOfHits.FindObject(Form("allHits_FMD%d%c",det,ring));
	301	TH1F* doubleHits = (TH1F*)fListOfHits.FindObject(Form("DoubleHits_FMD%d%c",det,ring));
	302	fBackground->SetDoubleHitCorrection(det,ring,doubleHits);
	303	doubleHits->Divide(allHits);
	304	for(Int_t vertexBin=0;vertexBin<fNvtxBins ;vertexBin++) {
	305	TH2F* hHits = (TH2F*)fListOfHits.FindObject(Form("hHits_FMD%d%c_vtx%d", det,ring,vertexBin));
	306	TH2F* hPrimary = (TH2F*)fListOfPrimaries.FindObject( Form("hPrimary_FMD_%c_vtx%d",ring,vertexBin));
	307	TH2F* hCorrection = (TH2F*)hHits->Clone(Form("FMD%d%c_vtxbin_%d_correction",det,ring,vertexBin));
	308	hCorrection->Divide(hPrimary);
0b0a4ae5	309	/*for(Int_t i = 1; i<=hCorrection->GetNbinsX();i++) {
	310	for(Int_t j=1; j<=hCorrection->GetNbinsY();j++) {
	311	if(hCorrection()->GetBinContent(i,j) == 0)
	312	continue;
	313	Float_t a = h
	314
	315
	316	}
	317	}
	318	*/
	319
d31dce33	320	hCorrection->SetTitle(hCorrection->GetName());
	321	fListOfCorrection.Add(hCorrection);
	322	fBackground->SetBgCorrection(det,ring,vertexBin,hCorrection);
	323
	324
	325	}
	326
	327	}
	328	}
	329	TAxis refAxis(fNvtxBins,-1*fZvtxCut,fZvtxCut);
	330	fBackground->SetRefAxis(&refAxis);
	331
	332	}
	333	//_____________________________________________________________________
	334	void AliFMDAnalysisTaskGenerateCorrection::ReadFromFile(const Char_t* filename, Bool_t storeInOCDB, Int_t /runNo/) {
0b0a4ae5	335
d31dce33	336	TFile infile(filename);
	337	TH1F* hVertex = (TH1F*)infile.Get("VertexBins");
	338	fZvtxCut = hVertex->GetXaxis()->GetXmax();
	339	fNvtxBins = hVertex->GetXaxis()->GetNbins();
	340	fVertexBins.SetName("VertexBins");
	341	fVertexBins.GetXaxis()->Set(fNvtxBins,-1*fZvtxCut,fZvtxCut);
	342
	343	TList* listOfHits = (TList*)infile.Get("Hits");
	344	TList* listOfPrim = (TList*)infile.Get("Primaries");
	345
	346	TH1F* hEventsSelected = (TH1F*)listOfHits->FindObject("EventsSelected");
	347	TH1F* hEventsAll = (TH1F*)listOfPrim->FindObject("EventsAll");
	348	fListOfHits.Add(hEventsSelected);
	349	fListOfPrimaries.Add(hEventsAll);
	350
	351	for(Int_t det =1; det<=3;det++)
	352	{
	353	Int_t nRings = (det==1 ? 1 : 2);
	354	for(Int_t ring = 0;ring<nRings;ring++)
	355	{
	356	Char_t ringChar = (ring == 0 ? 'I' : 'O');
	357	TH1F* allHits = (TH1F*)listOfHits->FindObject(Form("allHits_FMD%d%c",det,ringChar));
	358	TH1F* doubleHits = (TH1F*)listOfHits->FindObject(Form("DoubleHits_FMD%d%c",det,ringChar));
	359	fListOfHits.Add(allHits);
	360	fListOfHits.Add(doubleHits);
	361	for(Int_t v=0; v<fNvtxBins;v++)
	362	{
	363
	364	TH2F* hHits = (TH2F*)listOfHits->FindObject(Form("hHits_FMD%d%c_vtx%d", det,ringChar,v));
	365	fListOfHits.Add(hHits);
	366	}
	367	}
	368	}
	369	for(Int_t iring = 0; iring<2;iring++) {
	370	Char_t ringChar = (iring == 0 ? 'I' : 'O');
	371	for(Int_t v=0; v<fNvtxBins;v++) {
	372
	373	TH2F* hPrimary = (TH2F*)listOfPrim->FindObject( Form("hPrimary_FMD_%c_vtx%d",ringChar,v));
	374	fListOfPrimaries.Add(hPrimary);
	375
	376	}
	377	}
	378	GenerateCorrection();
	379
	380	TFile fout("backgroundFromFile.root","recreate");
	381	fListOfHits.Write();
	382	fListOfPrimaries.Write();
	383	fListOfCorrection.Write();
	384	fVertexBins.Write();
	385
	386	fout.Close();
0b0a4ae5	387	AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
d31dce33	388	if(storeInOCDB) {
0b0a4ae5	389	TFile fbg(pars->GetPath(pars->GetBackgroundID()),"RECREATE");
d31dce33	390	fBackground->Write(AliFMDAnaParameters::GetBackgroundID());
d31dce33	391	fbg.Close();
0b0a4ae5	392	TFile feselect(pars->GetPath(pars->GetEventSelectionEffID()),"RECREATE");
d31dce33	393	fEventSelectionEff->Write(AliFMDAnaParameters::GetEventSelectionEffID());
	394	feselect.Close();
	395
	396	}
	397
	398
	399
	400
	401	}
	402	//_____________________________________________________________________
	403	//
	404	// EOF
	405	//