[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(),
    fPrimaryMapInner(),
    fPrimaryMapOuter(),
    fHitMap(0),           // nDetector=0 means 51200 slots
    fLastTrackByStrip(0), // nDetector=0 means 51200 slots
    fPrim(0),
    fHits(0),
    fZvtxCut(0),
    fNvtxBins(0),
    fPrevTrack(-1),
    fPrevDetector(-1),
    fPrevRing('Q'),
    fPrevSec(-1),
    fNbinsEta(100)
{
  if(hits_not_trackref) {
    AddLoad(kHits);
    AddLoad(kTracks); 
  }
  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 = 
	    static_cast<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 = 
	  static_cast<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      = 
	static_cast<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  = static_cast<TH2F*>(primRingArray->At(vertexBin));
	primaryList->Add(hPrimary);
	TH2F* hCorrection = 
	  static_cast<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");
  // FIXME: Do not use a hard-coded path!
  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());
  
  // std::cout<<"Length: "<<h->Length()<<std::endl;
  // std::cout<<"Is stopped "<<h->IsStop()<<"    "<<kTRUE<<std::endl;
  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) return kTRUE;
  
  /*  if(charge !=  0 && 
     ((nTrack != fPrevTrack) || 
      (det != fPrevDetector))){ || 
      (ring != fPrevRing))){    ||
				  (sec != fPrevSec))) {*/
  /*Float_t nstrips = (ring =='O' ? 256 : 512);
  
  
  Float_t prevStripTrack = -1;
  if(strip >0)
    prevStripTrack = fLastTrackByStrip.operator()(det,ring,sec,strip-1);
  
  Float_t nextStripTrack = -1;
  if(strip < (nstrips - 1))
    nextStripTrack = fLastTrackByStrip.operator()(det,ring,sec,strip+1);
  */
  
  Float_t thisStripTrack = fLastTrackByStrip(det,ring,sec,strip);
  // if(nTrack == fLastTrackByStrip.operator()(det,ring,sec,strip))
  //  std::cout<<"Track # "<<nTrack<<"  failed the cut in "<<det<<"   "<<ring<<"   "<<sec<<"  "<<strip<<std::endl;
  // else
  //  std::cout<<"Track # "<<nTrack<<"  passed the cut in "<<det<<"   "<<ring<<"   "<<sec<<"  "<<strip<<std::endl;
  if(nTrack == thisStripTrack) return kTRUE;
    
  fHitMap(det,ring,sec,strip) += 1;
  fHits++;
  Float_t nstrips = (ring =='O' ? 256 : 512);
  fLastTrackByStrip(det,ring,sec,strip)     = Float_t(nTrack);
  if(strip >0)
    fLastTrackByStrip(det,ring,sec,strip-1) = Float_t(nTrack);
  if(strip < (nstrips - 1))
    fLastTrackByStrip(det,ring,sec,strip+1) = Float_t(nTrack);
  
  fPrevDetector = det;
  fPrevRing     = ring;
  fPrevSec      = sec;
  fPrevTrack    = nTrack;
  
  return kTRUE;

}

//_____________________________________________________________________
Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::Init() 
{
  fLastTrackByStrip.Reset(-1);
  
  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);
    if (charge == 0) continue;

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