[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(),
  fLastTrackByStrip(),
  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          = (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());
  
  // 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 && 
     ((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.operator()(det,ring,0,0);
  // 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(charge != 0 && nTrack != thisStripTrack) {
    
    fHitMap.operator()(det,ring,sec,strip) += 1;
    fHits++;
    Float_t nstrips = (ring =='O' ? 256 : 512);
    fLastTrackByStrip.operator()(det,ring,sec,strip) = (Float_t)nTrack;
    if(strip >0)
      fLastTrackByStrip.operator()(det,ring,sec,strip-1) = (Float_t)nTrack;
    if(strip < (nstrips - 1))
      fLastTrackByStrip.operator()(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);
    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) {
	    //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  = (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,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
	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(),
5ea1e0a9	73	fPrimaryMapInner(),
5ea1e0a9	74	fPrimaryMapOuter(),
627ba74b	75	fHitMap(),
6289b3e8	76	fLastTrackByStrip(),
df0d5480	77	fPrim(0),
	78	fHits(0),
	79	fZvtxCut(0),
	80	fNvtxBins(0),
	81	fPrevTrack(-1),
	82	fPrevDetector(-1),
	83	fPrevRing('Q'),
5d5bc6be	84	fPrevSec(-1),
df0d5480	85	fNbinsEta(100)
df0d5480	86	{
6289b3e8	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);
df0d5480	95	AddLoad(kHeader);
627ba74b	96
627ba74b	97
df0d5480	98	}
	99
	100	//_____________________________________________________________________
	101
08d168d9	102	void AliFMDBackgroundCorrection::GenerateBackgroundCorrection(Bool_t from_hits,
08d168d9	103	const Int_t nvtxbins,
df0d5480	104	Float_t zvtxcut,
ca3af5d9	105	const Int_t nBinsEta,
df0d5480	106	Bool_t storeInAlien,
	107	Int_t runNo,
	108	Int_t endRunNo,
	109	const Char_t* filename,
	110	Bool_t simulate,
109e97cd	111	Int_t nEvents,
	112	Bool_t inFile,
	113	const Char_t* infilename) {
df0d5480	114
ca3af5d9	115	//TGrid::Connect("alien:",0,0,"t");
df0d5480	116	if(simulate)
	117	Simulate(nEvents);
	118	else {
ca3af5d9	119	//AliCDBManager::Instance()->SetDefaultStorage("alien://Folder=/alice/data/2008/LHC08d/OCDB/");
162637e4	120	AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
df0d5480	121	AliCDBManager::Instance()->SetRun(runNo);
	122
	123	#if defined(__CINT__)
	124	gSystem->Load("liblhapdf");
	125	gSystem->Load("libEGPythia6");
	126	gSystem->Load("libpythia6");
	127	gSystem->Load("libAliPythia6");
	128	gSystem->Load("libgeant321");
	129	#endif
ca3af5d9	130	//
df0d5480	131	}
	132
	133	//Setting up the geometry
	134	//-----------------------------------------------
	135	if (AliGeomManager::GetGeometry() == NULL)
	136	AliGeomManager::LoadGeometry();
	137
	138	AliFMDGeometry* geo = AliFMDGeometry::Instance();
	139	geo->Init();
	140	geo->InitTransformations();
ca3af5d9	141
df0d5480	142	AliInfo("Processing hits and primaries ");
df0d5480	143
08d168d9	144	AliFMDInputBG input(from_hits);
df0d5480	145
109e97cd	146	if(!inFile) {
	147
	148	input.SetVtxCutZ(zvtxcut);
	149	input.SetNvtxBins(nvtxbins);
	150	input.SetNbinsEta(nBinsEta);
	151	input.Run();
	152	}
df0d5480	153
109e97cd	154	AliInfo(Form("Found %d primaries and %d hits.", input.GetNprim(),input.GetNhits()));
	155	TObjArray* hitArray ;
	156	TObjArray* primaryArray;
	157	if(inFile) {
	158	TFile* infile = TFile::Open(infilename);
	159	hitArray = new TObjArray();
	160	primaryArray = new TObjArray();
	161
	162	for(Int_t det =1; det<=3;det++)
	163	{
	164	TObjArray* detArrayHits = new TObjArray();
	165	detArrayHits->SetName(Form("FMD%d",det));
	166	hitArray->AddAtAndExpand(detArrayHits,det);
	167	Int_t nRings = (det==1 ? 1 : 2);
	168	for(Int_t ring = 0;ring<nRings;ring++)
	169	{
	170
	171	Char_t ringChar = (ring == 0 ? 'I' : 'O');
	172	TObjArray* vtxArrayHits = new TObjArray();
	173	vtxArrayHits->SetName(Form("FMD%d%c",det,ringChar));
	174	detArrayHits->AddAtAndExpand(vtxArrayHits,ring);
	175	for(Int_t v=0; v<nvtxbins;v++)
	176	{
	177
	178	TH2F* hHits = (TH2F*)infile->Get(Form("hHits_FMD%d%c_vtx%d",det,ringChar,v));
	179
	180
	181	vtxArrayHits->AddAtAndExpand(hHits,v);
	182
	183	}
	184	}
	185
	186	}
	187
	188	for(Int_t iring = 0; iring<2;iring++) {
	189	Char_t ringChar = (iring == 0 ? 'I' : 'O');
	190	TObjArray* ringArray = new TObjArray();
	191	ringArray->SetName(Form("FMD_%c",ringChar));
	192	primaryArray->AddAtAndExpand(ringArray,iring);
	193	for(Int_t v=0; v<nvtxbins;v++) {
	194
	195	TH2F* hPrimary = (TH2F*)infile->Get(Form("hPrimary_FMD_%c_vtx%d",ringChar,v));
	196	ringArray->AddAtAndExpand(hPrimary,v);
	197	}
	198	}
	199
	200
	201	}
	202	else {
	203	hitArray = input.GetHits();
	204	primaryArray = input.GetPrimaries();
	205	}
df0d5480	206	fCorrectionArray.SetName("FMD_bg_correction");
df0d5480	207	fCorrectionArray.SetOwner();
ca3af5d9	208
ca3af5d9	209	TList* primaryList = new TList();
26c8b5ff	210	primaryList->SetName("primaries");
ca3af5d9	211
26c8b5ff	212	TList* hitList = new TList();
	213	hitList->SetName("hits");
	214	TList* corrList = new TList();
	215	corrList->SetName("corrections");
109e97cd	216
b82e76e0	217	AliFMDAnaCalibBackgroundCorrection* background = 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++) {
	228	TObjArray* primRingArray = (TObjArray*)primaryArray->At(iring);
	229	Char_t ringChar = (iring == 0 ? 'I' : 'O');
	230	TObjArray* vtxArrayCorrection = new TObjArray();
	231	vtxArrayCorrection->SetName(Form("FMD%d%c",det,ringChar));
	232	detArrayCorrection->AddAtAndExpand(vtxArrayCorrection,iring);
	233
	234	for(Int_t vertexBin=0;vertexBin<nvtxbins;vertexBin++) {
	235	TObjArray* detArray = (TObjArray*)hitArray->At(det);
	236	TObjArray* vtxArray = (TObjArray*)detArray->At(iring);
	237	TH2F* hHits = (TH2F*)vtxArray->At(vertexBin);
26c8b5ff	238	hitList->Add(hHits);
df0d5480	239	TH2F* hPrimary = (TH2F*)primRingArray->At(vertexBin);
26c8b5ff	240	primaryList->Add(hPrimary);
df0d5480	241	TH2F* hCorrection = (TH2F*)hHits->Clone(Form("FMD%d%c_vtxbin_%d_correction",det,ringChar,vertexBin));
ca3af5d9	242	hCorrection->Divide(hPrimary);
26c8b5ff	243	corrList->Add(hCorrection);
8dc823cc	244	hCorrection->SetTitle(hCorrection->GetName());
df0d5480	245	vtxArrayCorrection->AddAtAndExpand(hCorrection,vertexBin);
b82e76e0	246	background->SetBgCorrection(det,ringChar,vertexBin,hCorrection);
df0d5480	247	}
	248
	249	}
	250	}
	251
	252	TAxis refAxis(nvtxbins,-1*zvtxcut,zvtxcut);
109e97cd	253	if(inFile) {
	254	TFile* infile = TFile::Open(infilename);
	255	TAxis* refaxis = (TAxis*)infile->Get("vertexbins");
	256	background->SetRefAxis(refaxis);
	257
	258	}
	259	else background->SetRefAxis(&refAxis);
df0d5480	260
8dc823cc	261	TFile* fout = new TFile(filename,"RECREATE");
df0d5480	262	refAxis.Write("vertexbins");
26c8b5ff	263
	264	hitList->Write();
	265	primaryList->Write();
	266	corrList->Write();
ca3af5d9	267
8dc823cc	268	TObjArray* container = new TObjArray();
	269	container->SetOwner();
	270	container->AddAtAndExpand(&refAxis,0);
	271	container->AddAtAndExpand(&fCorrectionArray,1);
	272	container->AddAtAndExpand(hitArray,2);
	273	container->AddAtAndExpand(primaryArray,3);
b82e76e0	274
b82e76e0	275
df0d5480	276	if(storeInAlien) {
df0d5480	277	AliCDBManager* cdb = AliCDBManager::Instance();
162637e4	278	cdb->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
8dc823cc	279	AliCDBId id(AliFMDAnaParameters::GetBackgroundPath(),runNo,endRunNo);
df0d5480	280
	281	AliCDBMetaData* meta = new AliCDBMetaData;
	282	meta->SetResponsible(gSystem->GetUserInfo()->fRealName.Data());
	283	meta->SetAliRootVersion(gROOT->GetVersion());
	284	meta->SetBeamPeriod(1);
	285	meta->SetComment("Background Correction for FMD");
b82e76e0	286	meta->SetProperty("key1", background );
b82e76e0	287	cdb->Put(background, id, meta);
df0d5480	288
	289	}
	290
8dc823cc	291	fout->Close();
8dc823cc	292
df0d5480	293
ca3af5d9	294	}
ca3af5d9	295
df0d5480	296	//_____________________________________________________________________
	297	void AliFMDBackgroundCorrection::Simulate(Int_t nEvents) {
	298
	299	AliSimulation sim ;
	300	sim.SetRunNumber(0);
	301	TGrid::Connect("alien:",0,0,"t");
	302	sim.SetDefaultStorage("alien://Folder=/alice/data/2008/LHC08d/OCDB/");
	303	sim.SetConfigFile("Config.C");
	304	sim.SetRunQA("FMD:");
	305	TStopwatch timer;
	306	timer.Start();
	307	sim.RunSimulation(nEvents);
	308	timer.Stop();
	309	timer.Print();
	310
	311	}
	312
	313	//_____________________________________________________________________
faf80567	314	Bool_t
	315	AliFMDBackgroundCorrection::AliFMDInputBG::ProcessHit(AliFMDHit* h,
	316	TParticle* /p/)
	317	{
df0d5480	318
08d168d9	319	if(!h)
08d168d9	320	return kTRUE;
faf80567	321	Bool_t retval = ProcessEvent(h->Detector(),
faf80567	322	h->Ring(),
98accb9e	323	h->Sector(),
98accb9e	324	h->Strip(),
faf80567	325	h->Track(),
faf80567	326	h->Q());
08d168d9	327
6289b3e8	328	// std::cout<<"Length: "<<h->Length()<<std::endl;
6289b3e8	329	// std::cout<<"Is stopped "<<h->IsStop()<<" "<<kTRUE<<std::endl;
08d168d9	330	return retval;
	331	}
	332	//_____________________________________________________________________
faf80567	333	Bool_t
6289b3e8	334	AliFMDBackgroundCorrection::AliFMDInputBG::ProcessTrackRef(AliTrackReference* tr, TParticle* p)
faf80567	335	{
08d168d9	336	if(!tr)
	337	return kTRUE;
	338	UShort_t det,sec,strip;
	339	Char_t ring;
	340	AliFMDStripIndex::Unpack(tr->UserId(),det,ring,sec,strip);
faf80567	341	Int_t nTrack = tr->GetTrack();
faf80567	342	TDatabasePDG* pdgDB = TDatabasePDG::Instance();
08d168d9	343	TParticlePDG* pdgPart = pdgDB->GetParticle(p->GetPdgCode());
faf80567	344	Float_t charge = (pdgPart ? pdgPart->Charge() : 0);
faf80567	345	Bool_t retval = ProcessEvent(det,ring,sec,strip,nTrack,charge);
08d168d9	346	return retval;
	347
	348	}
	349	//_____________________________________________________________________
faf80567	350	Bool_t
	351	AliFMDBackgroundCorrection::AliFMDInputBG::ProcessEvent(UShort_t det,
	352	Char_t ring,
	353	UShort_t sec,
	354	UShort_t strip,
	355	Int_t nTrack,
	356	Float_t charge)
08d168d9	357	{
08d168d9	358
df0d5480	359
df0d5480	360
6289b3e8	361	/* if(charge != 0 &&
627ba74b	362	((nTrack != fPrevTrack) \|\|
6289b3e8	363	(det != fPrevDetector))){ \|\|
	364	(ring != fPrevRing))){ \|\|
	365	(sec != fPrevSec))) {*/
	366	/*Float_t nstrips = (ring =='O' ? 256 : 512);
	367
	368
	369	Float_t prevStripTrack = -1;
	370	if(strip >0)
	371	prevStripTrack = fLastTrackByStrip.operator()(det,ring,sec,strip-1);
	372
	373	Float_t nextStripTrack = -1;
	374	if(strip < (nstrips - 1))
	375	nextStripTrack = fLastTrackByStrip.operator()(det,ring,sec,strip+1);
	376	*/
	377
	378	Float_t thisStripTrack = fLastTrackByStrip.operator()(det,ring,0,0);
	379	// if(nTrack == fLastTrackByStrip.operator()(det,ring,sec,strip))
	380	// std::cout<<"Track # "<<nTrack<<" failed the cut in "<<det<<" "<<ring<<" "<<sec<<" "<<strip<<std::endl;
	381	// else
	382	// std::cout<<"Track # "<<nTrack<<" passed the cut in "<<det<<" "<<ring<<" "<<sec<<" "<<strip<<std::endl;
	383	if(charge != 0 && nTrack != thisStripTrack) {
	384
	385	fHitMap.operator()(det,ring,sec,strip) += 1;
08d168d9	386	fHits++;
6289b3e8	387	Float_t nstrips = (ring =='O' ? 256 : 512);
	388	fLastTrackByStrip.operator()(det,ring,sec,strip) = (Float_t)nTrack;
	389	if(strip >0)
	390	fLastTrackByStrip.operator()(det,ring,sec,strip-1) = (Float_t)nTrack;
	391	if(strip < (nstrips - 1))
	392	fLastTrackByStrip.operator()(det,ring,sec,strip+1) = (Float_t)nTrack;
5d5bc6be	393
6289b3e8	394	fPrevDetector = det;
	395	fPrevRing = ring;
	396	fPrevSec = sec;
	397	fPrevTrack = nTrack;
08d168d9	398	}
ca3af5d9	399
6289b3e8	400
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);
	482	Float_t phi = TMath::ATan2(p->Py(),p->Px());
faf80567	483
faf80567	484	if(phi<0) phi = phi+2*TMath::Pi();
627ba74b	485	// if(p->Theta() == 0) continue;
	486	Float_t eta = p->Eta();
	487
	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);
	492	if(primary && charge !=0) {
	493
	494	fPrim++;
faf80567	495
627ba74b	496	fPrimaryMapInner.Fill(eta,phi);
627ba74b	497	fPrimaryMapOuter.Fill(eta,phi);
faf80567	498	}
faf80567	499	}
109e97cd	500
ca3af5d9	501	return retVal;
	502	}
	503	//_____________________________________________________________________
627ba74b	504	Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::End() {
	505
	506	Bool_t retval = AliFMDInput::End();
	507
	508	AliGenEventHeader* genHeader = fLoader->GetHeader()->GenEventHeader();
	509	TArrayF vertex;
	510	genHeader->PrimaryVertex(vertex);
	511
	512	if(TMath::Abs(vertex.At(2)) > fZvtxCut)
	513	return kTRUE;
	514
	515	Double_t delta = 2*fZvtxCut/fNvtxBins;
	516	Double_t vertexBinDouble = (vertex.At(2) + fZvtxCut) / delta;
	517	Int_t vertexBin = (Int_t)vertexBinDouble;
	518	//Primaries
	519	TObjArray* innerArray = (TObjArray*)fPrimaryArray.At(0);
	520	TObjArray* outerArray = (TObjArray*)fPrimaryArray.At(1);
	521
	522	TH2F* hPrimaryInner = (TH2F*)innerArray->At(vertexBin);
	523	TH2F* hPrimaryOuter = (TH2F*)outerArray->At(vertexBin);
	524
	525	hPrimaryInner->Add(&fPrimaryMapInner);
	526	hPrimaryOuter->Add(&fPrimaryMapOuter);
	527
	528	//Hits
	529	for(UShort_t det=1;det<=3;det++) {
	530	Int_t nRings = (det==1 ? 1 : 2);
	531	for (UShort_t ir = 0; ir < nRings; ir++) {
	532	Char_t ring = (ir == 0 ? 'I' : 'O');
	533	UShort_t nsec = (ir == 0 ? 20 : 40);
	534	UShort_t nstr = (ir == 0 ? 512 : 256);
	535
	536	for(UShort_t sec =0; sec < nsec; sec++) {
	537
	538	for(UShort_t strip = 0; strip < nstr; strip++) {
	539
	540	if(fHitMap.operator()(det,ring,sec,strip) > 0) {
6289b3e8	541	//std::cout<<fHitMap.operator()(det,ring,sec,strip)<<std::endl;
627ba74b	542	Double_t x,y,z;
	543	AliFMDGeometry* fmdgeo = AliFMDGeometry::Instance();
	544	fmdgeo->Detector2XYZ(det,ring,sec,strip,x,y,z);
	545
	546	Int_t iring = (ring == 'I' ? 0 : 1);
	547
	548	TObjArray* detArray = (TObjArray*)fHitArray.At(det);
	549	TObjArray* vtxArray = (TObjArray*)detArray->At(iring);
	550	TH2F* hHits = (TH2F*)vtxArray->At(vertexBin);
	551
	552	Float_t phi = TMath::ATan2(y,x);
	553	if(phi<0) phi = phi+2*TMath::Pi();
	554	Float_t r = TMath::Sqrt(TMath::Power(x,2)+TMath::Power(y,2));
	555	Float_t theta = TMath::ATan2(r,z-vertex.At(2));
	556	Float_t eta = -1TMath::Log(TMath::Tan(0.5theta));
6289b3e8	557	hHits->Fill(eta,phi,fHitMap.operator()(det,ring,sec,strip));
627ba74b	558	}
	559
	560	}
	561	}
	562	}
	563	}
	564
	565	fPrimaryMapInner.Reset();
	566	fPrimaryMapOuter.Reset();
	567	fHitMap.Reset(0);
6289b3e8	568	fLastTrackByStrip.Reset(-1);
627ba74b	569
	570	return retval;
	571	}
	572
ca3af5d9	573	//
	574	// EOF
	575	//
df0d5480	576
df0d5480	577