* 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.
+// 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 "AliFMDHit.h"
#include "AliLoader.h"
#include "AliFMD.h"
-#include "TH2F.h"
#include "AliGenEventHeader.h"
#include "AliHeader.h"
#include "TFile.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()
+ fCorrectionArray(),
+ fPrimaryList()
{}
//_____________________________________________________________________
-AliFMDBackgroundCorrection::AliFMDInputBG::AliFMDInputBG() :
- AliFMDInput(),
- fPrim(0),
- fHits(0),
- fZvtxCut(0),
- fNvtxBins(0),
- fPrevTrack(-1),
- fPrevDetector(-1),
- fPrevRing('Q'),
- fNbinsEta(100)
+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)
{
- AddLoad(kTracks);
- AddLoad(kHits);
+ if(hits_not_trackref) {
+ AddLoad(kHits);
+ AddLoad(kTracks);
+ }
+ else
+ AddLoad(kTrackRefs);
AddLoad(kKinematics);
AddLoad(kHeader);
}
//_____________________________________________________________________
-void AliFMDBackgroundCorrection::GenerateBackgroundCorrection(Int_t nvtxbins,
- Float_t zvtxcut,
- Int_t nBinsEta,
- Bool_t storeInAlien,
- Int_t runNo,
- Int_t endRunNo,
- const Char_t* filename,
- Bool_t simulate,
- Int_t nEvents) {
-
- TGrid::Connect("alien:",0,0,"t");
+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("alien://Folder=/alice/data/2008/LHC08d/OCDB/");
+ AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
AliCDBManager::Instance()->SetRun(runNo);
#if defined(__CINT__)
gSystem->Load("libAliPythia6");
gSystem->Load("libgeant321");
#endif
+ //
}
//Setting up the geometry
AliFMDGeometry* geo = AliFMDGeometry::Instance();
geo->Init();
geo->InitTransformations();
-
-
-
+
AliInfo("Processing hits and primaries ");
- AliFMDInputBG input;
- input.SetVtxCutZ(zvtxcut);
- input.SetNvtxBins(nvtxbins);
- input.SetNbinsEta(nBinsEta);
- input.Run();
+ AliFMDInputBG input(from_hits);
- AliInfo(Form("Found %d primaries and %d hits.", input.GetNprim(),input.GetNhits()));
+ if(!inFile) {
+
+ input.SetVtxCutZ(zvtxcut);
+ input.SetNvtxBins(nvtxbins);
+ input.SetNbinsEta(nBinsEta);
+ input.Run();
+ }
- TObjArray* hitArray = input.GetHits();
- TObjArray* primaryArray = input.GetPrimaries();
+ 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);
for(Int_t iring = 0; iring<nRings; iring++) {
- TObjArray* primRingArray = (TObjArray*)primaryArray->At(iring);
- Char_t ringChar = (iring == 0 ? 'I' : 'O');
+ 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);
TObjArray* detArray = (TObjArray*)hitArray->At(det);
TObjArray* vtxArray = (TObjArray*)detArray->At(iring);
TH2F* hHits = (TH2F*)vtxArray->At(vertexBin);
- TH2F* hPrimary = (TH2F*)primRingArray->At(vertexBin);
- TH2F* hCorrection = (TH2F*)hHits->Clone(Form("FMD%d%c_vtxbin_%d_correction",det,ringChar,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);
- TFile fout(filename,"RECREATE");
- refAxis.Write("vertexbins");
- fout.mkdir("Hits");
- fout.cd("Hits");
- hitArray->Write();
- fout.mkdir("Primaries");
- fout.cd("Primaries");
- primaryArray->Write();
- fout.mkdir("Correction");
- fout.cd("Correction");
- fCorrectionArray.Write();
if(storeInAlien) {
AliCDBManager* cdb = AliCDBManager::Instance();
- cdb->SetDefaultStorage("local://$ALICE_ROOT");
- AliCDBId id("FMD/AnalysisCalib/Background",runNo,endRunNo);
+ 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", &fout );
- cdb->Put(&fout, id, meta);
+ meta->SetProperty("key1", background );
+ cdb->Put(background, id, meta);
}
-
- fout.Close();
+ 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:");
}
//_____________________________________________________________________
-Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::ProcessTrack(Int_t i , TParticle* p, AliFMDHit* h) {
+Bool_t
+AliFMDBackgroundCorrection::AliFMDInputBG::ProcessHit(AliFMDHit* h,
+ TParticle* /*p*/)
+{
-
- // if(!h)
- // return kTRUE;
- //if(!p)
- // return kTRUE;
-
- // if(h) {
- // if(h->Detector() == 3)
- // std::cout<<"Track "<<i<<" hit "<<"FMD3"<<h->Ring()<<std::endl;
-
- //}
+ 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;
- AliGenEventHeader* genHeader = fLoader->GetHeader()->GenEventHeader();
- TArrayF vertex;
- genHeader->PrimaryVertex(vertex);
+}
+//_____________________________________________________________________
+Bool_t
+AliFMDBackgroundCorrection::AliFMDInputBG::ProcessEvent(UShort_t det,
+ Char_t ring,
+ UShort_t sec,
+ UShort_t strip,
+ Int_t nTrack,
+ Float_t charge)
+{
- if(TMath::Abs(vertex.At(2)) > fZvtxCut)
- return kTRUE;
+ if (charge == 0) return kTRUE;
- Double_t delta = 2*fZvtxCut/fNvtxBins;
- Double_t vertexBinDouble = (vertex.At(2) + fZvtxCut) / delta;
- Int_t vertexBin = (Int_t)vertexBinDouble;
+ /* if(charge != 0 &&
+ ((nTrack != fPrevTrack) ||
+ (det != fPrevDetector))){ ||
+ (ring != fPrevRing))){ ||
+ (sec != fPrevSec))) {*/
+ /*Float_t nstrips = (ring =='O' ? 256 : 512);
- if(h)
- if(h->Q() != 0 && (i != fPrevTrack || h->Detector() != fPrevDetector || h->Ring() != fPrevRing)) {
-
- Int_t det = h->Detector();
- Char_t ring = h->Ring();
- // if(h->Detector() == 3)
- // std::cout<<"Detected a hit in " <<"FMD"<<det<<ring<<" ! "<<std::endl;
-
- Int_t iring = (ring == 'I' ? 0 : 1);
-
- TObjArray* detArray = (TObjArray*)fHitArray.At(det);
- TObjArray* vtxArray = (TObjArray*)detArray->At(iring);
- TH2F* hHits = (TH2F*)vtxArray->At(vertexBin);
-
- Float_t phi = TMath::ATan2(h->Py(),h->Px());
- if(phi<0)
- phi = phi+2*TMath::Pi();
- Float_t theta = TMath::ATan2(TMath::Sqrt(TMath::Power(h->X(),2)+TMath::Power(h->Y(),2)),h->Z()+vertex.At(2));
- Float_t eta = -1*TMath::Log(TMath::Tan(0.5*theta));
- hHits->Fill(eta,phi);
- fHits++;
- //fPrevTrack = i;
- fPrevDetector = det;
- fPrevRing = ring;
- }
- if(p && i != fPrevTrack) {
-
- 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();
-
-
- Bool_t primary = (charge!=0)&&(TMath::Abs(p->Vx() - vertex.At(0))<0.1)&&(TMath::Abs(p->Vy() - vertex.At(1))<0.1)&&(TMath::Abs(p->Vz() - vertex.At(2))<0.1);
- if(primary) {
- fPrim++;
- TObjArray* innerArray = (TObjArray*)fPrimaryArray.At(0);
- TObjArray* outerArray = (TObjArray*)fPrimaryArray.At(1);
-
- TH2F* hPrimaryInner = (TH2F*)innerArray->At(vertexBin);
- TH2F* hPrimaryOuter = (TH2F*)outerArray->At(vertexBin);
- hPrimaryInner->Fill(p->Eta(),phi);
- hPrimaryOuter->Fill(p->Eta(),phi);
- }
- }
+ 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;
- fPrevTrack = i;
return kTRUE;
+
}
//_____________________________________________________________________
-Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::Init() {
+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();
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());
- vtxArrayHits->AddAtAndExpand(hHits,v);
-
- }
- }
-
+ 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();
- AliFMDInput::Init();
+ 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() {
- return kTRUE;
+ 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
+//
+// EOF
+//