[u/mrichter/AliRoot.git] / FMD / AliFMDBackgroundCorrection.cxx

/**************************************************************************
 * Copyright(c) 2004, 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.                  *
 **************************************************************************/

// Thil class computes background corrections for the FMD. The correction is computed 
// in eta,phi cells and the objects stored can be put into alien to use with analysis.
// It is based on the AliFMDInput class that is used to loop over hits and primaries.
//
// Author: Hans Hjersing Dalsgaard, NBI, hans.dalsgaard@cern.ch
//
//

#include "AliSimulation.h"
#include "TStopwatch.h"
#include "iostream"
#include "TGrid.h"
#include "AliRunLoader.h"
#include "AliGeomManager.h"
#include "AliFMDGeometry.h"
#include "AliStack.h"
#include "TParticle.h"
#include "TDatabasePDG.h"
#include "TParticlePDG.h"
#include "AliRun.h"
#include "AliFMDBackgroundCorrection.h"
#include "TSystem.h"
#include "AliCDBManager.h"
#include "TTree.h"
#include "TClonesArray.h"
#include "TBranch.h"
#include "AliFMDHit.h"
#include "AliLoader.h" 
#include "AliFMD.h"
#include "TH2F.h"
#include "AliGenEventHeader.h"
#include "AliHeader.h"
#include "TFile.h"
#include "TAxis.h"
#include "AliCDBId.h"
#include "AliCDBMetaData.h"
#include "TROOT.h"
#include "AliFMDParameters.h"
#include "AliLog.h"
#include "TList.h"
#include "AliFMDAnaParameters.h"
#include "AliFMDAnaCalibBackgroundCorrection.h"

ClassImp(AliFMDBackgroundCorrection)
//_____________________________________________________________________
AliFMDBackgroundCorrection::AliFMDBackgroundCorrection() : 
  TNamed(),
  fCorrectionArray(),
  fPrimaryList()
{} 

//_____________________________________________________________________
AliFMDBackgroundCorrection::AliFMDInputBG::AliFMDInputBG() : 
  AliFMDInput(),
  fPrimaryArray(),
  fHitArray(),
  fPrim(0),
  fHits(0),
  fZvtxCut(0),
  fNvtxBins(0),
  fPrevTrack(-1),
  fPrevDetector(-1),
  fPrevRing('Q'),
  fNbinsEta(100)
{
  AddLoad(kTracks); 
  AddLoad(kHits); 
  AddLoad(kKinematics); 
  AddLoad(kHeader);
}

//_____________________________________________________________________

void AliFMDBackgroundCorrection::GenerateBackgroundCorrection(const Int_t nvtxbins,
							      Float_t zvtxcut, 
							      const Int_t nBinsEta, 
							      Bool_t storeInAlien, 
							      Int_t runNo,
							      Int_t endRunNo, 
							      const Char_t* filename, 
							      Bool_t simulate,
							      Int_t nEvents,
							      Bool_t inFile,
							      const Char_t* infilename) {
  
  //TGrid::Connect("alien:",0,0,"t");
  if(simulate)
    Simulate(nEvents);
  else {
    //AliCDBManager::Instance()->SetDefaultStorage("alien://Folder=/alice/data/2008/LHC08d/OCDB/");
    AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT");
    AliCDBManager::Instance()->SetRun(runNo);
    
#if defined(__CINT__)
    gSystem->Load("liblhapdf");
    gSystem->Load("libEGPythia6");
    gSystem->Load("libpythia6");
    gSystem->Load("libAliPythia6");
    gSystem->Load("libgeant321");
#endif
      // 
  }  
  
  //Setting up the geometry
  //-----------------------------------------------
  if (AliGeomManager::GetGeometry() == NULL)
    AliGeomManager::LoadGeometry();
  
  AliFMDGeometry* geo = AliFMDGeometry::Instance();
  geo->Init();
  geo->InitTransformations();
    
  AliInfo("Processing hits and primaries ");
  
  AliFMDInputBG input;
  
  if(!inFile) {
  
    input.SetVtxCutZ(zvtxcut);
    input.SetNvtxBins(nvtxbins);
    input.SetNbinsEta(nBinsEta);
    input.Run();
  }
  
  AliInfo(Form("Found %d primaries and %d hits.", input.GetNprim(),input.GetNhits()));
  TObjArray* hitArray ;
  TObjArray* primaryArray;
  if(inFile) {
    TFile* infile = TFile::Open(infilename);
    hitArray     = new TObjArray();
    primaryArray = new TObjArray();
    
    for(Int_t det =1; det<=3;det++)
      {
	TObjArray* detArrayHits = new TObjArray();
	detArrayHits->SetName(Form("FMD%d",det));
	hitArray->AddAtAndExpand(detArrayHits,det);
	Int_t nRings = (det==1 ? 1 : 2);
	for(Int_t ring = 0;ring<nRings;ring++)
	  {
	    
	    Char_t ringChar = (ring == 0 ? 'I' : 'O');
	    TObjArray* vtxArrayHits = new TObjArray();
	    vtxArrayHits->SetName(Form("FMD%d%c",det,ringChar));
	    detArrayHits->AddAtAndExpand(vtxArrayHits,ring);
	    for(Int_t v=0; v<nvtxbins;v++)
	      {
	      
		TH2F* hHits          = (TH2F*)infile->Get(Form("hHits_FMD%d%c_vtx%d",det,ringChar,v));
		
		
	      vtxArrayHits->AddAtAndExpand(hHits,v);
	      	      
	    } 
	}
      
    }
    
    for(Int_t iring = 0; iring<2;iring++) {
      Char_t ringChar = (iring == 0 ? 'I' : 'O');
      TObjArray* ringArray = new TObjArray();
      ringArray->SetName(Form("FMD_%c",ringChar));
      primaryArray->AddAtAndExpand(ringArray,iring);
      for(Int_t v=0; v<nvtxbins;v++) {
	
	TH2F* hPrimary       = (TH2F*)infile->Get(Form("hPrimary_FMD_%c_vtx%d",ringChar,v));
	ringArray->AddAtAndExpand(hPrimary,v);
      }
    }
    
    
  }
  else {
    hitArray     = input.GetHits();
    primaryArray = input.GetPrimaries();
  }
  fCorrectionArray.SetName("FMD_bg_correction");
  fCorrectionArray.SetOwner();
   
  TList* primaryList     = new TList();
  primaryList->SetName("primaries");
  
  TList* hitList     = new TList();
  hitList->SetName("hits");
  TList* corrList    = new TList();
  corrList->SetName("corrections");
  
  AliFMDAnaCalibBackgroundCorrection* background = new AliFMDAnaCalibBackgroundCorrection();
  
  for(Int_t det= 1; det <=3; det++) {
    Int_t nRings = (det==1 ? 1 : 2);
    
    TObjArray* detArrayCorrection = new TObjArray();
    detArrayCorrection->SetName(Form("FMD%d",det));
    fCorrectionArray.AddAtAndExpand(detArrayCorrection,det);
    
    
    for(Int_t iring = 0; iring<nRings; iring++) {
      TObjArray* primRingArray = (TObjArray*)primaryArray->At(iring);
      Char_t ringChar = (iring == 0 ? 'I' : 'O');
      TObjArray* vtxArrayCorrection = new TObjArray();
      vtxArrayCorrection->SetName(Form("FMD%d%c",det,ringChar));
      detArrayCorrection->AddAtAndExpand(vtxArrayCorrection,iring);
      
      for(Int_t vertexBin=0;vertexBin<nvtxbins;vertexBin++) {
	TObjArray* detArray  = (TObjArray*)hitArray->At(det);
	TObjArray* vtxArray  = (TObjArray*)detArray->At(iring);
	TH2F* hHits          = (TH2F*)vtxArray->At(vertexBin);
	hitList->Add(hHits);
	TH2F* hPrimary  = (TH2F*)primRingArray->At(vertexBin);
	primaryList->Add(hPrimary);
	TH2F* hCorrection = (TH2F*)hHits->Clone(Form("FMD%d%c_vtxbin_%d_correction",det,ringChar,vertexBin));
	hCorrection->Divide(hPrimary);
	corrList->Add(hCorrection);
	hCorrection->SetTitle(hCorrection->GetName());
	vtxArrayCorrection->AddAtAndExpand(hCorrection,vertexBin);
	background->SetBgCorrection(det,ringChar,vertexBin,hCorrection);
      }
      
    }
  }
  
  TAxis refAxis(nvtxbins,-1*zvtxcut,zvtxcut);
  if(inFile) {
    TFile* infile = TFile::Open(infilename);
    TAxis* refaxis = (TAxis*)infile->Get("vertexbins");
    background->SetRefAxis(refaxis);
      
  }
  else background->SetRefAxis(&refAxis);
  
  TFile*  fout = new TFile(filename,"RECREATE");
  refAxis.Write("vertexbins");
  
  hitList->Write();
  primaryList->Write();
  corrList->Write();
   
  TObjArray* container = new TObjArray();
  container->SetOwner();
  container->AddAtAndExpand(&refAxis,0);
  container->AddAtAndExpand(&fCorrectionArray,1);
  container->AddAtAndExpand(hitArray,2);
  container->AddAtAndExpand(primaryArray,3);
  
  
  if(storeInAlien) {
    AliCDBManager* cdb = AliCDBManager::Instance();
    cdb->SetDefaultStorage("local://$ALICE_ROOT");
    AliCDBId      id(AliFMDAnaParameters::GetBackgroundPath(),runNo,endRunNo);
    
    AliCDBMetaData* meta = new AliCDBMetaData;				
    meta->SetResponsible(gSystem->GetUserInfo()->fRealName.Data());	
    meta->SetAliRootVersion(gROOT->GetVersion());			
    meta->SetBeamPeriod(1);						
    meta->SetComment("Background Correction for FMD");
    meta->SetProperty("key1", background );
    cdb->Put(background, id, meta);
    
  }
  
  fout->Close();
 
  
  }

//_____________________________________________________________________
void AliFMDBackgroundCorrection::Simulate(Int_t nEvents) {
  
  AliSimulation sim ; 
  sim.SetRunNumber(0);
  TGrid::Connect("alien:",0,0,"t");
  sim.SetDefaultStorage("alien://Folder=/alice/data/2008/LHC08d/OCDB/");
  sim.SetConfigFile("Config.C");
  sim.SetRunQA("FMD:");
  TStopwatch timer;
  timer.Start();
  sim.RunSimulation(nEvents);    
  timer.Stop();
  timer.Print();
  
}

//_____________________________________________________________________
Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::ProcessHit( AliFMDHit* h, TParticle* p) {
  
  AliGenEventHeader* genHeader = fLoader->GetHeader()->GenEventHeader();
  TArrayF vertex;
  genHeader->PrimaryVertex(vertex);
  
  if(TMath::Abs(vertex.At(2)) > fZvtxCut) 
    return kTRUE;
  
  Double_t delta = 2*fZvtxCut/fNvtxBins;
  Double_t vertexBinDouble = (vertex.At(2) + fZvtxCut) / delta;
  Int_t vertexBin = (Int_t)vertexBinDouble;
  
  Int_t i = h->Track();
  
  if(h)
    if(h->Q() !=  0 && (i != fPrevTrack || h->Detector() != fPrevDetector || h->Ring() != fPrevRing)) {
      
      Int_t det = h->Detector();
      Char_t ring = h->Ring();
      Int_t iring = (ring == 'I' ? 0 : 1);
      
      TObjArray* detArray  = (TObjArray*)fHitArray.At(det);
      TObjArray* vtxArray  = (TObjArray*)detArray->At(iring);
      TH2F* hHits          = (TH2F*)vtxArray->At(vertexBin);
      
      Float_t phi   = TMath::ATan2(h->Y(),h->X());
      if(phi<0)
	phi = phi+2*TMath::Pi();
      Float_t theta = TMath::ATan2(TMath::Sqrt(TMath::Power(h->X(),2)+TMath::Power(h->Y(),2)),h->Z()+vertex.At(2));
      Float_t eta   = -1*TMath::Log(TMath::Tan(0.5*theta));
      hHits->Fill(eta,phi);
      fHits++;
      fPrevDetector = det;
      fPrevRing     = ring;
    }
   
  fPrevTrack = i;
  return kTRUE;
}

//_____________________________________________________________________
Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::Init() {
  
  fPrimaryArray.SetOwner();
  fPrimaryArray.SetName("FMD_primary");
  
  for(Int_t iring = 0; iring<2;iring++) {
    Char_t ringChar = (iring == 0 ? 'I' : 'O');
    TObjArray* ringArray = new TObjArray();
    ringArray->SetName(Form("FMD_%c",ringChar));
    fPrimaryArray.AddAtAndExpand(ringArray,iring);
    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();
      ringArray->AddAtAndExpand(hPrimary,v);
    }
  }
  
  
  fHitArray.SetOwner();
  fHitArray.SetName("FMD_hits");
   
  for(Int_t det =1; det<=3;det++)
    {
      TObjArray* detArrayHits = new TObjArray();
      detArrayHits->SetName(Form("FMD%d",det));
      fHitArray.AddAtAndExpand(detArrayHits,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');
	  TObjArray* vtxArrayHits = new TObjArray();
	  vtxArrayHits->SetName(Form("FMD%d%c",det,ringChar));
	  detArrayHits->AddAtAndExpand(vtxArrayHits,ring);
	  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();
	      vtxArrayHits->AddAtAndExpand(hHits,v);
	      	      
	    } 
	}
      
    }


  AliFMDInput::Init();
  
  
  return kTRUE;
}
//_____________________________________________________________________

Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::Begin(Int_t event ) {

  Bool_t retVal = AliFMDInput::Begin(event); 
  
  AliStack* partStack=fLoader->Stack();
  
  Int_t nTracks=partStack->GetNtrack();
  AliGenEventHeader* genHeader = fLoader->GetHeader()->GenEventHeader();
  TArrayF vertex;
  genHeader->PrimaryVertex(vertex);
  
  if(TMath::Abs(vertex.At(2)) > fZvtxCut) 
    return kTRUE;
  
  Double_t delta = 2*fZvtxCut/fNvtxBins;
  Double_t vertexBinDouble = (vertex.At(2) + fZvtxCut) / delta;
  Int_t vertexBin = (Int_t)vertexBinDouble;
  
  for(Int_t j=0;j<nTracks;j++)
	{
	  TParticle* p  = partStack->Particle(j);
	  //TDatabasePDG* pdgDB = TDatabasePDG::Instance();
	  //TParticlePDG* pdgPart = pdgDB->GetParticle(p->GetPdgCode());
	  //Float_t charge = (pdgPart ? pdgPart->Charge() : 0);
	  Float_t phi = TMath::ATan2(p->Py(),p->Px());
	  
	  if(phi<0)
	    phi = phi+2*TMath::Pi();
	  if(p->Theta() == 0) continue;
	      Float_t eta   = -1*TMath::Log(TMath::Tan(0.5*p->Theta()));
	  
	      Bool_t primary = partStack->IsPhysicalPrimary(j);
	      //(charge!=0)&&(TMath::Abs(p->Vx() - vertex.At(0))<0.01)&&(TMath::Abs(p->Vy() - vertex.At(1))<0.01)&&(TMath::Abs(p->Vz() - vertex.At(2))<0.01);
	      if(primary) {
		fPrim++;
		TObjArray* innerArray = (TObjArray*)fPrimaryArray.At(0);
		TObjArray* outerArray = (TObjArray*)fPrimaryArray.At(1);
		
		TH2F* hPrimaryInner  = (TH2F*)innerArray->At(vertexBin);
		TH2F* hPrimaryOuter  = (TH2F*)outerArray->At(vertexBin);
		hPrimaryInner->Fill(eta,phi);
		hPrimaryOuter->Fill(eta,phi);
	      }
	}
  
  return retVal;
}
//_____________________________________________________________________
//
// EOF
//
Commit	Line	Data
df0d5480	1	/**************************************************************************
	2	* Copyright(c) 2004, 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	// Thil class computes background corrections for the FMD. The correction is computed
	17	// in eta,phi cells and the objects stored can be put into alien to use with analysis.
	18	// It is based on the AliFMDInput class that is used to loop over hits and primaries.
	19	//
	20	// Author: Hans Hjersing Dalsgaard, NBI, hans.dalsgaard@cern.ch
	21	//
	22	//
	23
	24	#include "AliSimulation.h"
	25	#include "TStopwatch.h"
	26	#include "iostream"
	27	#include "TGrid.h"
	28	#include "AliRunLoader.h"
	29	#include "AliGeomManager.h"
	30	#include "AliFMDGeometry.h"
	31	#include "AliStack.h"
	32	#include "TParticle.h"
	33	#include "TDatabasePDG.h"
	34	#include "TParticlePDG.h"
	35	#include "AliRun.h"
	36	#include "AliFMDBackgroundCorrection.h"
	37	#include "TSystem.h"
	38	#include "AliCDBManager.h"
	39	#include "TTree.h"
	40	#include "TClonesArray.h"
	41	#include "TBranch.h"
	42	#include "AliFMDHit.h"
	43	#include "AliLoader.h"
	44	#include "AliFMD.h"
	45	#include "TH2F.h"
	46	#include "AliGenEventHeader.h"
	47	#include "AliHeader.h"
	48	#include "TFile.h"
	49	#include "TAxis.h"
	50	#include "AliCDBId.h"
	51	#include "AliCDBMetaData.h"
	52	#include "TROOT.h"
	53	#include "AliFMDParameters.h"
	54	#include "AliLog.h"
26c8b5ff	55	#include "TList.h"
8dc823cc	56	#include "AliFMDAnaParameters.h"
b82e76e0	57	#include "AliFMDAnaCalibBackgroundCorrection.h"
df0d5480	58
	59	ClassImp(AliFMDBackgroundCorrection)
	60	//_____________________________________________________________________
	61	AliFMDBackgroundCorrection::AliFMDBackgroundCorrection() :
	62	TNamed(),
ca3af5d9	63	fCorrectionArray(),
ca3af5d9	64	fPrimaryList()
df0d5480	65	{}
	66
	67	//_____________________________________________________________________
	68	AliFMDBackgroundCorrection::AliFMDInputBG::AliFMDInputBG() :
	69	AliFMDInput(),
8dc823cc	70	fPrimaryArray(),
8dc823cc	71	fHitArray(),
df0d5480	72	fPrim(0),
	73	fHits(0),
	74	fZvtxCut(0),
	75	fNvtxBins(0),
	76	fPrevTrack(-1),
	77	fPrevDetector(-1),
	78	fPrevRing('Q'),
	79	fNbinsEta(100)
	80	{
	81	AddLoad(kTracks);
	82	AddLoad(kHits);
	83	AddLoad(kKinematics);
	84	AddLoad(kHeader);
	85	}
	86
	87	//_____________________________________________________________________
	88
ca3af5d9	89	void AliFMDBackgroundCorrection::GenerateBackgroundCorrection(const Int_t nvtxbins,
df0d5480	90	Float_t zvtxcut,
ca3af5d9	91	const Int_t nBinsEta,
df0d5480	92	Bool_t storeInAlien,
	93	Int_t runNo,
	94	Int_t endRunNo,
	95	const Char_t* filename,
	96	Bool_t simulate,
109e97cd	97	Int_t nEvents,
	98	Bool_t inFile,
	99	const Char_t* infilename) {
df0d5480	100
ca3af5d9	101	//TGrid::Connect("alien:",0,0,"t");
df0d5480	102	if(simulate)
	103	Simulate(nEvents);
	104	else {
ca3af5d9	105	//AliCDBManager::Instance()->SetDefaultStorage("alien://Folder=/alice/data/2008/LHC08d/OCDB/");
ca3af5d9	106	AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT");
df0d5480	107	AliCDBManager::Instance()->SetRun(runNo);
	108
	109	#if defined(__CINT__)
	110	gSystem->Load("liblhapdf");
	111	gSystem->Load("libEGPythia6");
	112	gSystem->Load("libpythia6");
	113	gSystem->Load("libAliPythia6");
	114	gSystem->Load("libgeant321");
	115	#endif
ca3af5d9	116	//
df0d5480	117	}
	118
	119	//Setting up the geometry
	120	//-----------------------------------------------
	121	if (AliGeomManager::GetGeometry() == NULL)
	122	AliGeomManager::LoadGeometry();
	123
	124	AliFMDGeometry* geo = AliFMDGeometry::Instance();
	125	geo->Init();
	126	geo->InitTransformations();
ca3af5d9	127
df0d5480	128	AliInfo("Processing hits and primaries ");
	129
	130	AliFMDInputBG input;
df0d5480	131
109e97cd	132	if(!inFile) {
	133
	134	input.SetVtxCutZ(zvtxcut);
	135	input.SetNvtxBins(nvtxbins);
	136	input.SetNbinsEta(nBinsEta);
	137	input.Run();
	138	}
df0d5480	139
109e97cd	140	AliInfo(Form("Found %d primaries and %d hits.", input.GetNprim(),input.GetNhits()));
	141	TObjArray* hitArray ;
	142	TObjArray* primaryArray;
	143	if(inFile) {
	144	TFile* infile = TFile::Open(infilename);
	145	hitArray = new TObjArray();
	146	primaryArray = new TObjArray();
	147
	148	for(Int_t det =1; det<=3;det++)
	149	{
	150	TObjArray* detArrayHits = new TObjArray();
	151	detArrayHits->SetName(Form("FMD%d",det));
	152	hitArray->AddAtAndExpand(detArrayHits,det);
	153	Int_t nRings = (det==1 ? 1 : 2);
	154	for(Int_t ring = 0;ring<nRings;ring++)
	155	{
	156
	157	Char_t ringChar = (ring == 0 ? 'I' : 'O');
	158	TObjArray* vtxArrayHits = new TObjArray();
	159	vtxArrayHits->SetName(Form("FMD%d%c",det,ringChar));
	160	detArrayHits->AddAtAndExpand(vtxArrayHits,ring);
	161	for(Int_t v=0; v<nvtxbins;v++)
	162	{
	163
	164	TH2F* hHits = (TH2F*)infile->Get(Form("hHits_FMD%d%c_vtx%d",det,ringChar,v));
	165
	166
	167	vtxArrayHits->AddAtAndExpand(hHits,v);
	168
	169	}
	170	}
	171
	172	}
	173
	174	for(Int_t iring = 0; iring<2;iring++) {
	175	Char_t ringChar = (iring == 0 ? 'I' : 'O');
	176	TObjArray* ringArray = new TObjArray();
	177	ringArray->SetName(Form("FMD_%c",ringChar));
	178	primaryArray->AddAtAndExpand(ringArray,iring);
	179	for(Int_t v=0; v<nvtxbins;v++) {
	180
	181	TH2F* hPrimary = (TH2F*)infile->Get(Form("hPrimary_FMD_%c_vtx%d",ringChar,v));
	182	ringArray->AddAtAndExpand(hPrimary,v);
	183	}
	184	}
	185
	186
	187	}
	188	else {
	189	hitArray = input.GetHits();
	190	primaryArray = input.GetPrimaries();
	191	}
df0d5480	192	fCorrectionArray.SetName("FMD_bg_correction");
df0d5480	193	fCorrectionArray.SetOwner();
ca3af5d9	194
ca3af5d9	195	TList* primaryList = new TList();
26c8b5ff	196	primaryList->SetName("primaries");
ca3af5d9	197
26c8b5ff	198	TList* hitList = new TList();
	199	hitList->SetName("hits");
	200	TList* corrList = new TList();
	201	corrList->SetName("corrections");
109e97cd	202
b82e76e0	203	AliFMDAnaCalibBackgroundCorrection* background = new AliFMDAnaCalibBackgroundCorrection();
b82e76e0	204
df0d5480	205	for(Int_t det= 1; det <=3; det++) {
	206	Int_t nRings = (det==1 ? 1 : 2);
	207
	208	TObjArray* detArrayCorrection = new TObjArray();
	209	detArrayCorrection->SetName(Form("FMD%d",det));
	210	fCorrectionArray.AddAtAndExpand(detArrayCorrection,det);
	211
	212
	213	for(Int_t iring = 0; iring<nRings; iring++) {
	214	TObjArray* primRingArray = (TObjArray*)primaryArray->At(iring);
	215	Char_t ringChar = (iring == 0 ? 'I' : 'O');
	216	TObjArray* vtxArrayCorrection = new TObjArray();
	217	vtxArrayCorrection->SetName(Form("FMD%d%c",det,ringChar));
	218	detArrayCorrection->AddAtAndExpand(vtxArrayCorrection,iring);
	219
	220	for(Int_t vertexBin=0;vertexBin<nvtxbins;vertexBin++) {
	221	TObjArray* detArray = (TObjArray*)hitArray->At(det);
	222	TObjArray* vtxArray = (TObjArray*)detArray->At(iring);
	223	TH2F* hHits = (TH2F*)vtxArray->At(vertexBin);
26c8b5ff	224	hitList->Add(hHits);
df0d5480	225	TH2F* hPrimary = (TH2F*)primRingArray->At(vertexBin);
26c8b5ff	226	primaryList->Add(hPrimary);
df0d5480	227	TH2F* hCorrection = (TH2F*)hHits->Clone(Form("FMD%d%c_vtxbin_%d_correction",det,ringChar,vertexBin));
ca3af5d9	228	hCorrection->Divide(hPrimary);
26c8b5ff	229	corrList->Add(hCorrection);
8dc823cc	230	hCorrection->SetTitle(hCorrection->GetName());
df0d5480	231	vtxArrayCorrection->AddAtAndExpand(hCorrection,vertexBin);
b82e76e0	232	background->SetBgCorrection(det,ringChar,vertexBin,hCorrection);
df0d5480	233	}
	234
	235	}
	236	}
	237
	238	TAxis refAxis(nvtxbins,-1*zvtxcut,zvtxcut);
109e97cd	239	if(inFile) {
	240	TFile* infile = TFile::Open(infilename);
	241	TAxis* refaxis = (TAxis*)infile->Get("vertexbins");
	242	background->SetRefAxis(refaxis);
	243
	244	}
	245	else background->SetRefAxis(&refAxis);
df0d5480	246
8dc823cc	247	TFile* fout = new TFile(filename,"RECREATE");
df0d5480	248	refAxis.Write("vertexbins");
26c8b5ff	249
	250	hitList->Write();
	251	primaryList->Write();
	252	corrList->Write();
ca3af5d9	253
8dc823cc	254	TObjArray* container = new TObjArray();
	255	container->SetOwner();
	256	container->AddAtAndExpand(&refAxis,0);
	257	container->AddAtAndExpand(&fCorrectionArray,1);
	258	container->AddAtAndExpand(hitArray,2);
	259	container->AddAtAndExpand(primaryArray,3);
b82e76e0	260
b82e76e0	261
df0d5480	262	if(storeInAlien) {
	263	AliCDBManager* cdb = AliCDBManager::Instance();
	264	cdb->SetDefaultStorage("local://$ALICE_ROOT");
8dc823cc	265	AliCDBId id(AliFMDAnaParameters::GetBackgroundPath(),runNo,endRunNo);
df0d5480	266
	267	AliCDBMetaData* meta = new AliCDBMetaData;
	268	meta->SetResponsible(gSystem->GetUserInfo()->fRealName.Data());
	269	meta->SetAliRootVersion(gROOT->GetVersion());
	270	meta->SetBeamPeriod(1);
	271	meta->SetComment("Background Correction for FMD");
b82e76e0	272	meta->SetProperty("key1", background );
b82e76e0	273	cdb->Put(background, id, meta);
df0d5480	274
	275	}
	276
8dc823cc	277	fout->Close();
8dc823cc	278
df0d5480	279
ca3af5d9	280	}
ca3af5d9	281
df0d5480	282	//_____________________________________________________________________
	283	void AliFMDBackgroundCorrection::Simulate(Int_t nEvents) {
	284
	285	AliSimulation sim ;
	286	sim.SetRunNumber(0);
	287	TGrid::Connect("alien:",0,0,"t");
	288	sim.SetDefaultStorage("alien://Folder=/alice/data/2008/LHC08d/OCDB/");
	289	sim.SetConfigFile("Config.C");
	290	sim.SetRunQA("FMD:");
	291	TStopwatch timer;
	292	timer.Start();
	293	sim.RunSimulation(nEvents);
	294	timer.Stop();
	295	timer.Print();
	296
	297	}
	298
	299	//_____________________________________________________________________
ca3af5d9	300	Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::ProcessHit( AliFMDHit* h, TParticle* p) {
df0d5480	301
df0d5480	302	AliGenEventHeader* genHeader = fLoader->GetHeader()->GenEventHeader();
	303	TArrayF vertex;
	304	genHeader->PrimaryVertex(vertex);
	305
	306	if(TMath::Abs(vertex.At(2)) > fZvtxCut)
	307	return kTRUE;
	308
	309	Double_t delta = 2*fZvtxCut/fNvtxBins;
	310	Double_t vertexBinDouble = (vertex.At(2) + fZvtxCut) / delta;
	311	Int_t vertexBin = (Int_t)vertexBinDouble;
	312
ca3af5d9	313	Int_t i = h->Track();
ca3af5d9	314
df0d5480	315	if(h)
	316	if(h->Q() != 0 && (i != fPrevTrack \|\| h->Detector() != fPrevDetector \|\| h->Ring() != fPrevRing)) {
	317
	318	Int_t det = h->Detector();
	319	Char_t ring = h->Ring();
df0d5480	320	Int_t iring = (ring == 'I' ? 0 : 1);
	321
	322	TObjArray* detArray = (TObjArray*)fHitArray.At(det);
	323	TObjArray* vtxArray = (TObjArray*)detArray->At(iring);
	324	TH2F* hHits = (TH2F*)vtxArray->At(vertexBin);
	325
109e97cd	326	Float_t phi = TMath::ATan2(h->Y(),h->X());
df0d5480	327	if(phi<0)
	328	phi = phi+2*TMath::Pi();
	329	Float_t theta = TMath::ATan2(TMath::Sqrt(TMath::Power(h->X(),2)+TMath::Power(h->Y(),2)),h->Z()+vertex.At(2));
	330	Float_t eta = -1TMath::Log(TMath::Tan(0.5theta));
	331	hHits->Fill(eta,phi);
	332	fHits++;
df0d5480	333	fPrevDetector = det;
	334	fPrevRing = ring;
	335	}
ca3af5d9	336
df0d5480	337	fPrevTrack = i;
	338	return kTRUE;
	339	}
	340
	341	//_____________________________________________________________________
	342	Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::Init() {
	343
	344	fPrimaryArray.SetOwner();
	345	fPrimaryArray.SetName("FMD_primary");
	346
	347	for(Int_t iring = 0; iring<2;iring++) {
	348	Char_t ringChar = (iring == 0 ? 'I' : 'O');
	349	TObjArray* ringArray = new TObjArray();
	350	ringArray->SetName(Form("FMD_%c",ringChar));
	351	fPrimaryArray.AddAtAndExpand(ringArray,iring);
	352	Int_t nSec = (iring == 1 ? 40 : 20);
	353	for(Int_t v=0; v<fNvtxBins;v++) {
	354
	355	TH2F* hPrimary = new TH2F(Form("hPrimary_FMD_%c_vtx%d",ringChar,v),
	356	Form("hPrimary_FMD_%c_vtx%d",ringChar,v),
	357	fNbinsEta, -6,6, nSec, 0,2*TMath::Pi());
b82e76e0	358	hPrimary->Sumw2();
df0d5480	359	ringArray->AddAtAndExpand(hPrimary,v);
	360	}
	361	}
	362
	363
	364	fHitArray.SetOwner();
	365	fHitArray.SetName("FMD_hits");
	366
	367	for(Int_t det =1; det<=3;det++)
	368	{
	369	TObjArray* detArrayHits = new TObjArray();
	370	detArrayHits->SetName(Form("FMD%d",det));
	371	fHitArray.AddAtAndExpand(detArrayHits,det);
	372	Int_t nRings = (det==1 ? 1 : 2);
	373	for(Int_t ring = 0;ring<nRings;ring++)
	374	{
	375	Int_t nSec = (ring == 1 ? 40 : 20);
	376	Char_t ringChar = (ring == 0 ? 'I' : 'O');
	377	TObjArray* vtxArrayHits = new TObjArray();
	378	vtxArrayHits->SetName(Form("FMD%d%c",det,ringChar));
	379	detArrayHits->AddAtAndExpand(vtxArrayHits,ring);
	380	for(Int_t v=0; v<fNvtxBins;v++)
	381	{
	382
	383	TH2F* hHits = new TH2F(Form("hHits_FMD%d%c_vtx%d",det,ringChar,v),
	384	Form("hHits_FMD%d%c_vtx%d",det,ringChar,v),
	385	fNbinsEta, -6,6, nSec, 0, 2*TMath::Pi());
b82e76e0	386	hHits->Sumw2();
df0d5480	387	vtxArrayHits->AddAtAndExpand(hHits,v);
	388
	389	}
	390	}
	391
	392	}
	393
	394
	395
	396	AliFMDInput::Init();
	397
	398
	399	return kTRUE;
	400	}
ca3af5d9	401	//_____________________________________________________________________
	402
	403	Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::Begin(Int_t event ) {
df0d5480	404
ca3af5d9	405	Bool_t retVal = AliFMDInput::Begin(event);
	406
	407	AliStack* partStack=fLoader->Stack();
	408
	409	Int_t nTracks=partStack->GetNtrack();
	410	AliGenEventHeader* genHeader = fLoader->GetHeader()->GenEventHeader();
	411	TArrayF vertex;
	412	genHeader->PrimaryVertex(vertex);
	413
	414	if(TMath::Abs(vertex.At(2)) > fZvtxCut)
	415	return kTRUE;
	416
	417	Double_t delta = 2*fZvtxCut/fNvtxBins;
	418	Double_t vertexBinDouble = (vertex.At(2) + fZvtxCut) / delta;
	419	Int_t vertexBin = (Int_t)vertexBinDouble;
	420
	421	for(Int_t j=0;j<nTracks;j++)
	422	{
	423	TParticle* p = partStack->Particle(j);
109e97cd	424	//TDatabasePDG* pdgDB = TDatabasePDG::Instance();
	425	//TParticlePDG* pdgPart = pdgDB->GetParticle(p->GetPdgCode());
	426	//Float_t charge = (pdgPart ? pdgPart->Charge() : 0);
ca3af5d9	427	Float_t phi = TMath::ATan2(p->Py(),p->Px());
	428
	429	if(phi<0)
	430	phi = phi+2*TMath::Pi();
	431	if(p->Theta() == 0) continue;
	432	Float_t eta = -1TMath::Log(TMath::Tan(0.5p->Theta()));
	433
109e97cd	434	Bool_t primary = partStack->IsPhysicalPrimary(j);
	435	//(charge!=0)&&(TMath::Abs(p->Vx() - vertex.At(0))<0.01)&&(TMath::Abs(p->Vy() - vertex.At(1))<0.01)&&(TMath::Abs(p->Vz() - vertex.At(2))<0.01);
	436	if(primary) {
	437	fPrim++;
	438	TObjArray* innerArray = (TObjArray*)fPrimaryArray.At(0);
	439	TObjArray* outerArray = (TObjArray*)fPrimaryArray.At(1);
	440
	441	TH2F* hPrimaryInner = (TH2F*)innerArray->At(vertexBin);
	442	TH2F* hPrimaryOuter = (TH2F*)outerArray->At(vertexBin);
	443	hPrimaryInner->Fill(eta,phi);
	444	hPrimaryOuter->Fill(eta,phi);
	445	}
ca3af5d9	446	}
109e97cd	447
ca3af5d9	448	return retVal;
	449	}
	450	//_____________________________________________________________________
	451	//
	452	// EOF
	453	//
df0d5480	454
df0d5480	455