[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 "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"
#include "AliTrackReference.h"
#include "AliFMDStripIndex.h"

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

//_____________________________________________________________________
AliFMDBackgroundCorrection::AliFMDInputBG::AliFMDInputBG(Bool_t hits_not_trackref) : 
  AliFMDInput("galice.root"),
  fPrimaryArray(),
  fHitArray(),
  fHitMap(),
  fPrim(0),
  fHits(0),
  fZvtxCut(0),
  fNvtxBins(0),
  fPrevTrack(-1),
  fPrevDetector(-1),
  fPrevRing('Q'),
  fPrevSec(-1),
  fNbinsEta(100)
{
  AddLoad(kTracks); 
  if(hits_not_trackref)
    AddLoad(kHits);
  else
    AddLoad(kTrackRefs);
  AddLoad(kKinematics); 
  AddLoad(kHeader);
  
  
}

//_____________________________________________________________________

void AliFMDBackgroundCorrection::GenerateBackgroundCorrection(Bool_t from_hits,
							      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/OCDB");
    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(from_hits);
  
  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/OCDB");
    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*/) 
{
  
  if(!h)
    return kTRUE;
  Bool_t retval = ProcessEvent(h->Detector(),
			       h->Ring(),
			       h->Sector(),
			       h->Strip(),
			       h->Track(),
			       h->Q());
  
  return retval;
}
//_____________________________________________________________________
Bool_t 
AliFMDBackgroundCorrection::AliFMDInputBG::ProcessTrackRef(AliTrackReference* tr, 
							   TParticle* p) 
{
  if(!tr)
    return kTRUE;
  UShort_t det,sec,strip;
  Char_t   ring;
  AliFMDStripIndex::Unpack(tr->UserId(),det,ring,sec,strip);
  Int_t         nTrack  = tr->GetTrack();
  TDatabasePDG* pdgDB   = TDatabasePDG::Instance();
  TParticlePDG* pdgPart = pdgDB->GetParticle(p->GetPdgCode());
  Float_t       charge  = (pdgPart ? pdgPart->Charge() : 0);
  Bool_t        retval  = ProcessEvent(det,ring,sec,strip,nTrack,charge);
  return retval;
  
}
//_____________________________________________________________________
Bool_t 
AliFMDBackgroundCorrection::AliFMDInputBG::ProcessEvent(UShort_t det,
							Char_t   ring, 
							UShort_t sec, 
							UShort_t strip,
							Int_t    nTrack,
							Float_t  charge)
{
  
  
  
  if(charge !=  0 && 
     ((nTrack != fPrevTrack) || 
      (det != fPrevDetector) || 
      (ring != fPrevRing)    ||
      (sec != fPrevSec))) {
    fHitMap.operator()(det,ring,sec,strip) = 1;
    fHits++;
    
  }
  
  fPrevDetector = det;
  fPrevRing     = ring;
  fPrevSec      = sec;
  fPrevTrack    = nTrack;
  
  return kTRUE;

}

//_____________________________________________________________________
Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::Init() 
{
  fPrimaryArray.SetOwner();
  fPrimaryArray.SetName("FMD_primary");
  
  fPrimaryMapInner.SetBins(fNbinsEta, -6,6, 20, 0, 2*TMath::Pi());
  fPrimaryMapOuter.SetBins(fNbinsEta, -6,6, 40, 0, 2*TMath::Pi());
  fPrimaryMapInner.SetName("fPrimaryMapInner");
  fPrimaryMapInner.SetName("fPrimaryMapOuter");
  
  fPrimaryMapInner.Sumw2();
  fPrimaryMapOuter.Sumw2();
  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();
  
  
  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   = p->Eta();   
    
    // std::cout<<-1*TMath::Log(TMath::Tan(0.5*p->Theta()))<<std::endl;
    
    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 && charge !=0) {
      
      
      
      fPrim++;
      
      fPrimaryMapInner.Fill(eta,phi);
      fPrimaryMapOuter.Fill(eta,phi);
    }
  }
  
  return retVal;
}
//_____________________________________________________________________
Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::End()  {
  
  Bool_t retval = AliFMDInput::End();
  
  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;
  //Primaries
  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->Add(&fPrimaryMapInner);
  hPrimaryOuter->Add(&fPrimaryMapOuter);
  
  //Hits
  for(UShort_t det=1;det<=3;det++) {
    Int_t nRings = (det==1 ? 1 : 2);
    for (UShort_t ir = 0; ir < nRings; ir++) {
      Char_t   ring = (ir == 0 ? 'I' : 'O');
      UShort_t nsec = (ir == 0 ? 20  : 40);
      UShort_t nstr = (ir == 0 ? 512 : 256);
      
      for(UShort_t sec =0; sec < nsec;  sec++) {
	
	for(UShort_t strip = 0; strip < nstr; strip++) {
	  
	  if(fHitMap.operator()(det,ring,sec,strip) > 0) {
	  
	    Double_t x,y,z;
	    AliFMDGeometry* fmdgeo = AliFMDGeometry::Instance();
	    fmdgeo->Detector2XYZ(det,ring,sec,strip,x,y,z);
	    
	    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(y,x);
	    if(phi<0) phi   = phi+2*TMath::Pi();
	    Float_t   r     = TMath::Sqrt(TMath::Power(x,2)+TMath::Power(y,2));
	    Float_t   theta = TMath::ATan2(r,z-vertex.At(2));
	    Float_t   eta   = -1*TMath::Log(TMath::Tan(0.5*theta));
	    hHits->Fill(eta,phi);
	  }
	  
	}
      }
    }
  }
  
  fPrimaryMapInner.Reset();
  fPrimaryMapOuter.Reset();
  fHitMap.Reset(0);
  
  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"
df0d5480	45	#include "AliGenEventHeader.h"
	46	#include "AliHeader.h"
	47	#include "TFile.h"
	48	#include "TAxis.h"
	49	#include "AliCDBId.h"
	50	#include "AliCDBMetaData.h"
	51	#include "TROOT.h"
	52	#include "AliFMDParameters.h"
	53	#include "AliLog.h"
26c8b5ff	54	#include "TList.h"
8dc823cc	55	#include "AliFMDAnaParameters.h"
b82e76e0	56	#include "AliFMDAnaCalibBackgroundCorrection.h"
08d168d9	57	#include "AliTrackReference.h"
08d168d9	58	#include "AliFMDStripIndex.h"
df0d5480	59
	60	ClassImp(AliFMDBackgroundCorrection)
	61	//_____________________________________________________________________
	62	AliFMDBackgroundCorrection::AliFMDBackgroundCorrection() :
	63	TNamed(),
ca3af5d9	64	fCorrectionArray(),
ca3af5d9	65	fPrimaryList()
df0d5480	66	{}
	67
	68	//_____________________________________________________________________
08d168d9	69	AliFMDBackgroundCorrection::AliFMDInputBG::AliFMDInputBG(Bool_t hits_not_trackref) :
5d5bc6be	70	AliFMDInput("galice.root"),
8dc823cc	71	fPrimaryArray(),
8dc823cc	72	fHitArray(),
627ba74b	73	fHitMap(),
df0d5480	74	fPrim(0),
	75	fHits(0),
	76	fZvtxCut(0),
	77	fNvtxBins(0),
	78	fPrevTrack(-1),
	79	fPrevDetector(-1),
	80	fPrevRing('Q'),
5d5bc6be	81	fPrevSec(-1),
df0d5480	82	fNbinsEta(100)
	83	{
	84	AddLoad(kTracks);
08d168d9	85	if(hits_not_trackref)
	86	AddLoad(kHits);
	87	else
	88	AddLoad(kTrackRefs);
df0d5480	89	AddLoad(kKinematics);
df0d5480	90	AddLoad(kHeader);
627ba74b	91
627ba74b	92
df0d5480	93	}
	94
	95	//_____________________________________________________________________
	96
08d168d9	97	void AliFMDBackgroundCorrection::GenerateBackgroundCorrection(Bool_t from_hits,
08d168d9	98	const Int_t nvtxbins,
df0d5480	99	Float_t zvtxcut,
ca3af5d9	100	const Int_t nBinsEta,
df0d5480	101	Bool_t storeInAlien,
	102	Int_t runNo,
	103	Int_t endRunNo,
	104	const Char_t* filename,
	105	Bool_t simulate,
109e97cd	106	Int_t nEvents,
	107	Bool_t inFile,
	108	const Char_t* infilename) {
df0d5480	109
ca3af5d9	110	//TGrid::Connect("alien:",0,0,"t");
df0d5480	111	if(simulate)
	112	Simulate(nEvents);
	113	else {
ca3af5d9	114	//AliCDBManager::Instance()->SetDefaultStorage("alien://Folder=/alice/data/2008/LHC08d/OCDB/");
162637e4	115	AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
df0d5480	116	AliCDBManager::Instance()->SetRun(runNo);
	117
	118	#if defined(__CINT__)
	119	gSystem->Load("liblhapdf");
	120	gSystem->Load("libEGPythia6");
	121	gSystem->Load("libpythia6");
	122	gSystem->Load("libAliPythia6");
	123	gSystem->Load("libgeant321");
	124	#endif
ca3af5d9	125	//
df0d5480	126	}
	127
	128	//Setting up the geometry
	129	//-----------------------------------------------
	130	if (AliGeomManager::GetGeometry() == NULL)
	131	AliGeomManager::LoadGeometry();
	132
	133	AliFMDGeometry* geo = AliFMDGeometry::Instance();
	134	geo->Init();
	135	geo->InitTransformations();
ca3af5d9	136
df0d5480	137	AliInfo("Processing hits and primaries ");
df0d5480	138
08d168d9	139	AliFMDInputBG input(from_hits);
df0d5480	140
109e97cd	141	if(!inFile) {
	142
	143	input.SetVtxCutZ(zvtxcut);
	144	input.SetNvtxBins(nvtxbins);
	145	input.SetNbinsEta(nBinsEta);
	146	input.Run();
	147	}
df0d5480	148
109e97cd	149	AliInfo(Form("Found %d primaries and %d hits.", input.GetNprim(),input.GetNhits()));
	150	TObjArray* hitArray ;
	151	TObjArray* primaryArray;
	152	if(inFile) {
	153	TFile* infile = TFile::Open(infilename);
	154	hitArray = new TObjArray();
	155	primaryArray = new TObjArray();
	156
	157	for(Int_t det =1; det<=3;det++)
	158	{
	159	TObjArray* detArrayHits = new TObjArray();
	160	detArrayHits->SetName(Form("FMD%d",det));
	161	hitArray->AddAtAndExpand(detArrayHits,det);
	162	Int_t nRings = (det==1 ? 1 : 2);
	163	for(Int_t ring = 0;ring<nRings;ring++)
	164	{
	165
	166	Char_t ringChar = (ring == 0 ? 'I' : 'O');
	167	TObjArray* vtxArrayHits = new TObjArray();
	168	vtxArrayHits->SetName(Form("FMD%d%c",det,ringChar));
	169	detArrayHits->AddAtAndExpand(vtxArrayHits,ring);
	170	for(Int_t v=0; v<nvtxbins;v++)
	171	{
	172
	173	TH2F* hHits = (TH2F*)infile->Get(Form("hHits_FMD%d%c_vtx%d",det,ringChar,v));
	174
	175
	176	vtxArrayHits->AddAtAndExpand(hHits,v);
	177
	178	}
	179	}
	180
	181	}
	182
	183	for(Int_t iring = 0; iring<2;iring++) {
	184	Char_t ringChar = (iring == 0 ? 'I' : 'O');
	185	TObjArray* ringArray = new TObjArray();
	186	ringArray->SetName(Form("FMD_%c",ringChar));
	187	primaryArray->AddAtAndExpand(ringArray,iring);
	188	for(Int_t v=0; v<nvtxbins;v++) {
	189
	190	TH2F* hPrimary = (TH2F*)infile->Get(Form("hPrimary_FMD_%c_vtx%d",ringChar,v));
	191	ringArray->AddAtAndExpand(hPrimary,v);
	192	}
	193	}
	194
	195
	196	}
	197	else {
	198	hitArray = input.GetHits();
	199	primaryArray = input.GetPrimaries();
	200	}
df0d5480	201	fCorrectionArray.SetName("FMD_bg_correction");
df0d5480	202	fCorrectionArray.SetOwner();
ca3af5d9	203
ca3af5d9	204	TList* primaryList = new TList();
26c8b5ff	205	primaryList->SetName("primaries");
ca3af5d9	206
26c8b5ff	207	TList* hitList = new TList();
	208	hitList->SetName("hits");
	209	TList* corrList = new TList();
	210	corrList->SetName("corrections");
109e97cd	211
b82e76e0	212	AliFMDAnaCalibBackgroundCorrection* background = new AliFMDAnaCalibBackgroundCorrection();
b82e76e0	213
df0d5480	214	for(Int_t det= 1; det <=3; det++) {
	215	Int_t nRings = (det==1 ? 1 : 2);
	216
	217	TObjArray* detArrayCorrection = new TObjArray();
	218	detArrayCorrection->SetName(Form("FMD%d",det));
	219	fCorrectionArray.AddAtAndExpand(detArrayCorrection,det);
	220
	221
	222	for(Int_t iring = 0; iring<nRings; iring++) {
	223	TObjArray* primRingArray = (TObjArray*)primaryArray->At(iring);
	224	Char_t ringChar = (iring == 0 ? 'I' : 'O');
	225	TObjArray* vtxArrayCorrection = new TObjArray();
	226	vtxArrayCorrection->SetName(Form("FMD%d%c",det,ringChar));
	227	detArrayCorrection->AddAtAndExpand(vtxArrayCorrection,iring);
	228
	229	for(Int_t vertexBin=0;vertexBin<nvtxbins;vertexBin++) {
	230	TObjArray* detArray = (TObjArray*)hitArray->At(det);
	231	TObjArray* vtxArray = (TObjArray*)detArray->At(iring);
	232	TH2F* hHits = (TH2F*)vtxArray->At(vertexBin);
26c8b5ff	233	hitList->Add(hHits);
df0d5480	234	TH2F* hPrimary = (TH2F*)primRingArray->At(vertexBin);
26c8b5ff	235	primaryList->Add(hPrimary);
df0d5480	236	TH2F* hCorrection = (TH2F*)hHits->Clone(Form("FMD%d%c_vtxbin_%d_correction",det,ringChar,vertexBin));
ca3af5d9	237	hCorrection->Divide(hPrimary);
26c8b5ff	238	corrList->Add(hCorrection);
8dc823cc	239	hCorrection->SetTitle(hCorrection->GetName());
df0d5480	240	vtxArrayCorrection->AddAtAndExpand(hCorrection,vertexBin);
b82e76e0	241	background->SetBgCorrection(det,ringChar,vertexBin,hCorrection);
df0d5480	242	}
	243
	244	}
	245	}
	246
	247	TAxis refAxis(nvtxbins,-1*zvtxcut,zvtxcut);
109e97cd	248	if(inFile) {
	249	TFile* infile = TFile::Open(infilename);
	250	TAxis* refaxis = (TAxis*)infile->Get("vertexbins");
	251	background->SetRefAxis(refaxis);
	252
	253	}
	254	else background->SetRefAxis(&refAxis);
df0d5480	255
8dc823cc	256	TFile* fout = new TFile(filename,"RECREATE");
df0d5480	257	refAxis.Write("vertexbins");
26c8b5ff	258
	259	hitList->Write();
	260	primaryList->Write();
	261	corrList->Write();
ca3af5d9	262
8dc823cc	263	TObjArray* container = new TObjArray();
	264	container->SetOwner();
	265	container->AddAtAndExpand(&refAxis,0);
	266	container->AddAtAndExpand(&fCorrectionArray,1);
	267	container->AddAtAndExpand(hitArray,2);
	268	container->AddAtAndExpand(primaryArray,3);
b82e76e0	269
b82e76e0	270
df0d5480	271	if(storeInAlien) {
df0d5480	272	AliCDBManager* cdb = AliCDBManager::Instance();
162637e4	273	cdb->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
8dc823cc	274	AliCDBId id(AliFMDAnaParameters::GetBackgroundPath(),runNo,endRunNo);
df0d5480	275
	276	AliCDBMetaData* meta = new AliCDBMetaData;
	277	meta->SetResponsible(gSystem->GetUserInfo()->fRealName.Data());
	278	meta->SetAliRootVersion(gROOT->GetVersion());
	279	meta->SetBeamPeriod(1);
	280	meta->SetComment("Background Correction for FMD");
b82e76e0	281	meta->SetProperty("key1", background );
b82e76e0	282	cdb->Put(background, id, meta);
df0d5480	283
	284	}
	285
8dc823cc	286	fout->Close();
8dc823cc	287
df0d5480	288
ca3af5d9	289	}
ca3af5d9	290
df0d5480	291	//_____________________________________________________________________
	292	void AliFMDBackgroundCorrection::Simulate(Int_t nEvents) {
	293
	294	AliSimulation sim ;
	295	sim.SetRunNumber(0);
	296	TGrid::Connect("alien:",0,0,"t");
	297	sim.SetDefaultStorage("alien://Folder=/alice/data/2008/LHC08d/OCDB/");
	298	sim.SetConfigFile("Config.C");
	299	sim.SetRunQA("FMD:");
	300	TStopwatch timer;
	301	timer.Start();
	302	sim.RunSimulation(nEvents);
	303	timer.Stop();
	304	timer.Print();
	305
	306	}
	307
	308	//_____________________________________________________________________
faf80567	309	Bool_t
	310	AliFMDBackgroundCorrection::AliFMDInputBG::ProcessHit(AliFMDHit* h,
	311	TParticle* /p/)
	312	{
df0d5480	313
08d168d9	314	if(!h)
08d168d9	315	return kTRUE;
faf80567	316	Bool_t retval = ProcessEvent(h->Detector(),
faf80567	317	h->Ring(),
98accb9e	318	h->Sector(),
98accb9e	319	h->Strip(),
faf80567	320	h->Track(),
faf80567	321	h->Q());
08d168d9	322
	323	return retval;
	324	}
	325	//_____________________________________________________________________
faf80567	326	Bool_t
	327	AliFMDBackgroundCorrection::AliFMDInputBG::ProcessTrackRef(AliTrackReference* tr,
	328	TParticle* p)
	329	{
08d168d9	330	if(!tr)
	331	return kTRUE;
	332	UShort_t det,sec,strip;
	333	Char_t ring;
	334	AliFMDStripIndex::Unpack(tr->UserId(),det,ring,sec,strip);
faf80567	335	Int_t nTrack = tr->GetTrack();
faf80567	336	TDatabasePDG* pdgDB = TDatabasePDG::Instance();
08d168d9	337	TParticlePDG* pdgPart = pdgDB->GetParticle(p->GetPdgCode());
faf80567	338	Float_t charge = (pdgPart ? pdgPart->Charge() : 0);
faf80567	339	Bool_t retval = ProcessEvent(det,ring,sec,strip,nTrack,charge);
08d168d9	340	return retval;
	341
	342	}
	343	//_____________________________________________________________________
faf80567	344	Bool_t
	345	AliFMDBackgroundCorrection::AliFMDInputBG::ProcessEvent(UShort_t det,
	346	Char_t ring,
	347	UShort_t sec,
	348	UShort_t strip,
	349	Int_t nTrack,
	350	Float_t charge)
08d168d9	351	{
08d168d9	352
df0d5480	353
df0d5480	354
faf80567	355	if(charge != 0 &&
627ba74b	356	((nTrack != fPrevTrack) \|\|
	357	(det != fPrevDetector) \|\|
	358	(ring != fPrevRing) \|\|
	359	(sec != fPrevSec))) {
	360	fHitMap.operator()(det,ring,sec,strip) = 1;
08d168d9	361	fHits++;
5d5bc6be	362
08d168d9	363	}
ca3af5d9	364
5d5bc6be	365	fPrevDetector = det;
	366	fPrevRing = ring;
	367	fPrevSec = sec;
	368	fPrevTrack = nTrack;
	369
df0d5480	370	return kTRUE;
08d168d9	371
df0d5480	372	}
	373
	374	//_____________________________________________________________________
faf80567	375	Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::Init()
faf80567	376	{
df0d5480	377	fPrimaryArray.SetOwner();
	378	fPrimaryArray.SetName("FMD_primary");
	379
627ba74b	380	fPrimaryMapInner.SetBins(fNbinsEta, -6,6, 20, 0, 2*TMath::Pi());
	381	fPrimaryMapOuter.SetBins(fNbinsEta, -6,6, 40, 0, 2*TMath::Pi());
	382	fPrimaryMapInner.SetName("fPrimaryMapInner");
	383	fPrimaryMapInner.SetName("fPrimaryMapOuter");
	384
	385	fPrimaryMapInner.Sumw2();
	386	fPrimaryMapOuter.Sumw2();
df0d5480	387	for(Int_t iring = 0; iring<2;iring++) {
	388	Char_t ringChar = (iring == 0 ? 'I' : 'O');
	389	TObjArray* ringArray = new TObjArray();
	390	ringArray->SetName(Form("FMD_%c",ringChar));
	391	fPrimaryArray.AddAtAndExpand(ringArray,iring);
	392	Int_t nSec = (iring == 1 ? 40 : 20);
	393	for(Int_t v=0; v<fNvtxBins;v++) {
	394
	395	TH2F* hPrimary = new TH2F(Form("hPrimary_FMD_%c_vtx%d",ringChar,v),
	396	Form("hPrimary_FMD_%c_vtx%d",ringChar,v),
	397	fNbinsEta, -6,6, nSec, 0,2*TMath::Pi());
b82e76e0	398	hPrimary->Sumw2();
df0d5480	399	ringArray->AddAtAndExpand(hPrimary,v);
	400	}
	401	}
	402
	403
	404	fHitArray.SetOwner();
	405	fHitArray.SetName("FMD_hits");
	406
faf80567	407	for(Int_t det =1; det<=3;det++) {
	408	TObjArray* detArrayHits = new TObjArray();
	409	detArrayHits->SetName(Form("FMD%d",det));
	410	fHitArray.AddAtAndExpand(detArrayHits,det);
	411	Int_t nRings = (det==1 ? 1 : 2);
	412	for(Int_t ring = 0;ring<nRings;ring++) {
	413	Int_t nSec = (ring == 1 ? 40 : 20);
	414	Char_t ringChar = (ring == 0 ? 'I' : 'O');
	415	TObjArray* vtxArrayHits = new TObjArray();
	416	vtxArrayHits->SetName(Form("FMD%d%c",det,ringChar));
	417	detArrayHits->AddAtAndExpand(vtxArrayHits,ring);
	418	for(Int_t v=0; v<fNvtxBins;v++) {
	419	TH2F* hHits = new TH2F(Form("hHits_FMD%d%c_vtx%d", det,ringChar,v),
	420	Form("hHits_FMD%d%c_vtx%d", det,ringChar,v),
	421	fNbinsEta, -6,6, nSec, 0, 2*TMath::Pi());
	422	hHits->Sumw2();
	423	vtxArrayHits->AddAtAndExpand(hHits,v);
	424
	425	}
df0d5480	426	}
faf80567	427	}
df0d5480	428
df0d5480	429	AliFMDInput::Init();
df0d5480	430
df0d5480	431	return kTRUE;
df0d5480	432	}
ca3af5d9	433	//_____________________________________________________________________
ca3af5d9	434
faf80567	435	Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::Begin(Int_t event )
faf80567	436	{
df0d5480	437
faf80567	438	Bool_t retVal = AliFMDInput::Begin(event);
	439	AliStack* partStack = fLoader->Stack();
	440	Int_t nTracks = partStack->GetNtrack();
ca3af5d9	441
faf80567	442
faf80567	443	for(Int_t j=0;j<nTracks;j++) {
5d5bc6be	444	TParticle* p = partStack->Particle(j);
	445	TDatabasePDG* pdgDB = TDatabasePDG::Instance();
	446	TParticlePDG* pdgPart = pdgDB->GetParticle(p->GetPdgCode());
	447	Float_t charge = (pdgPart ? pdgPart->Charge() : 0);
	448	Float_t phi = TMath::ATan2(p->Py(),p->Px());
faf80567	449
faf80567	450	if(phi<0) phi = phi+2*TMath::Pi();
627ba74b	451	// if(p->Theta() == 0) continue;
	452	Float_t eta = p->Eta();
	453
	454	// std::cout<<-1TMath::Log(TMath::Tan(0.5p->Theta()))<<std::endl;
faf80567	455
	456	Bool_t primary = partStack->IsPhysicalPrimary(j);
	457	//(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);
	458	if(primary && charge !=0) {
	459
627ba74b	460
627ba74b	461
faf80567	462	fPrim++;
faf80567	463
627ba74b	464	fPrimaryMapInner.Fill(eta,phi);
627ba74b	465	fPrimaryMapOuter.Fill(eta,phi);
faf80567	466	}
faf80567	467	}
109e97cd	468
ca3af5d9	469	return retVal;
	470	}
	471	//_____________________________________________________________________
627ba74b	472	Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::End() {
	473
	474	Bool_t retval = AliFMDInput::End();
	475
	476	AliGenEventHeader* genHeader = fLoader->GetHeader()->GenEventHeader();
	477	TArrayF vertex;
	478	genHeader->PrimaryVertex(vertex);
	479
	480	if(TMath::Abs(vertex.At(2)) > fZvtxCut)
	481	return kTRUE;
	482
	483	Double_t delta = 2*fZvtxCut/fNvtxBins;
	484	Double_t vertexBinDouble = (vertex.At(2) + fZvtxCut) / delta;
	485	Int_t vertexBin = (Int_t)vertexBinDouble;
	486	//Primaries
	487	TObjArray* innerArray = (TObjArray*)fPrimaryArray.At(0);
	488	TObjArray* outerArray = (TObjArray*)fPrimaryArray.At(1);
	489
	490	TH2F* hPrimaryInner = (TH2F*)innerArray->At(vertexBin);
	491	TH2F* hPrimaryOuter = (TH2F*)outerArray->At(vertexBin);
	492
	493	hPrimaryInner->Add(&fPrimaryMapInner);
	494	hPrimaryOuter->Add(&fPrimaryMapOuter);
	495
	496	//Hits
	497	for(UShort_t det=1;det<=3;det++) {
	498	Int_t nRings = (det==1 ? 1 : 2);
	499	for (UShort_t ir = 0; ir < nRings; ir++) {
	500	Char_t ring = (ir == 0 ? 'I' : 'O');
	501	UShort_t nsec = (ir == 0 ? 20 : 40);
	502	UShort_t nstr = (ir == 0 ? 512 : 256);
	503
	504	for(UShort_t sec =0; sec < nsec; sec++) {
	505
	506	for(UShort_t strip = 0; strip < nstr; strip++) {
	507
	508	if(fHitMap.operator()(det,ring,sec,strip) > 0) {
	509
	510	Double_t x,y,z;
	511	AliFMDGeometry* fmdgeo = AliFMDGeometry::Instance();
	512	fmdgeo->Detector2XYZ(det,ring,sec,strip,x,y,z);
	513
	514	Int_t iring = (ring == 'I' ? 0 : 1);
	515
	516	TObjArray* detArray = (TObjArray*)fHitArray.At(det);
	517	TObjArray* vtxArray = (TObjArray*)detArray->At(iring);
	518	TH2F* hHits = (TH2F*)vtxArray->At(vertexBin);
	519
	520	Float_t phi = TMath::ATan2(y,x);
	521	if(phi<0) phi = phi+2*TMath::Pi();
	522	Float_t r = TMath::Sqrt(TMath::Power(x,2)+TMath::Power(y,2));
	523	Float_t theta = TMath::ATan2(r,z-vertex.At(2));
	524	Float_t eta = -1TMath::Log(TMath::Tan(0.5theta));
	525	hHits->Fill(eta,phi);
	526	}
	527
	528	}
	529	}
	530	}
	531	}
	532
	533	fPrimaryMapInner.Reset();
	534	fPrimaryMapOuter.Reset();
	535	fHitMap.Reset(0);
536
537	return retval;
538	}
539
ca3af5d9	540	//
	541	// EOF
	542	//
df0d5480	543
df0d5480	544