]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - FMD/AliFMDBackgroundCorrection.cxx
remove dependency to aliroot libraries, access of ESDEvent object through abstract...
[u/mrichter/AliRoot.git] / FMD / AliFMDBackgroundCorrection.cxx
index d2dee003b58595662b7b9140203c7eb8ee71a752..89123f8794a79fa9acc8ab9c3f2d54a63b333874 100644 (file)
  * 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
 //
@@ -42,7 +43,6 @@
 #include "AliFMDHit.h"
 #include "AliLoader.h" 
 #include "AliFMD.h"
-#include "TH2F.h"
 #include "AliGenEventHeader.h"
 #include "AliHeader.h"
 #include "TFile.h"
@@ -55,6 +55,8 @@
 #include "TList.h"
 #include "AliFMDAnaParameters.h"
 #include "AliFMDAnaCalibBackgroundCorrection.h"
+#include "AliTrackReference.h"
+#include "AliFMDStripIndex.h"
 
 ClassImp(AliFMDBackgroundCorrection)
 //_____________________________________________________________________
@@ -65,43 +67,56 @@ AliFMDBackgroundCorrection::AliFMDBackgroundCorrection() :
 {} 
 
 //_____________________________________________________________________
-AliFMDBackgroundCorrection::AliFMDInputBG::AliFMDInputBG() : 
-  AliFMDInput(),
-  fPrimaryArray(),
-  fHitArray(),
-  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(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) {
-  
+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");
+    AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
     AliCDBManager::Instance()->SetRun(runNo);
     
 #if defined(__CINT__)
@@ -125,16 +140,68 @@ void AliFMDBackgroundCorrection::GenerateBackgroundCorrection(const Int_t nvtxbi
     
   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) {
   
-  TObjArray* hitArray = input.GetHits();
-  TObjArray* primaryArray = input.GetPrimaries();
+    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();
    
@@ -145,8 +212,9 @@ void AliFMDBackgroundCorrection::GenerateBackgroundCorrection(const Int_t nvtxbi
   hitList->SetName("hits");
   TList* corrList    = new TList();
   corrList->SetName("corrections");
-
-  AliFMDAnaCalibBackgroundCorrection* background = new AliFMDAnaCalibBackgroundCorrection();
+  
+  AliFMDAnaCalibBackgroundCorrection* background = 
+    new AliFMDAnaCalibBackgroundCorrection();
   
   for(Int_t det= 1; det <=3; det++) {
     Int_t nRings = (det==1 ? 1 : 2);
@@ -157,8 +225,9 @@ void AliFMDBackgroundCorrection::GenerateBackgroundCorrection(const Int_t nvtxbi
     
     
     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);
@@ -168,22 +237,28 @@ void AliFMDBackgroundCorrection::GenerateBackgroundCorrection(const Int_t nvtxbi
        TObjArray* vtxArray  = (TObjArray*)detArray->At(iring);
        TH2F* hHits          = (TH2F*)vtxArray->At(vertexBin);
        hitList->Add(hHits);
-       TH2F* hPrimary  = (TH2F*)primRingArray->At(vertexBin);
+       TH2F* hPrimary  = static_cast<TH2F*>(primRingArray->At(vertexBin));
        primaryList->Add(hPrimary);
-       TH2F* hCorrection = (TH2F*)hHits->Clone(Form("FMD%d%c_vtxbin_%d_correction",det,ringChar,vertexBin));
+       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);
-  
-  background->SetRefAxis(&refAxis);
+  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");
@@ -202,7 +277,7 @@ void AliFMDBackgroundCorrection::GenerateBackgroundCorrection(const Int_t nvtxbi
   
   if(storeInAlien) {
     AliCDBManager* cdb = AliCDBManager::Instance();
-    cdb->SetDefaultStorage("local://$ALICE_ROOT");
+    cdb->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
     AliCDBId      id(AliFMDAnaParameters::GetBackgroundPath(),runNo,endRunNo);
     
     AliCDBMetaData* meta = new AliCDBMetaData;                         
@@ -226,6 +301,7 @@ 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:");
@@ -238,53 +314,110 @@ void AliFMDBackgroundCorrection::Simulate(Int_t nEvents) {
 }
 
 //_____________________________________________________________________
-Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::ProcessHit( AliFMDHit* h, TParticle* p) {
-  
-  AliGenEventHeader* genHeader = fLoader->GetHeader()->GenEventHeader();
-  TArrayF vertex;
-  genHeader->PrimaryVertex(vertex);
+Bool_t 
+AliFMDBackgroundCorrection::AliFMDInputBG::ProcessHit(AliFMDHit* h, 
+                                                     TParticle* /*p*/) 
+{
   
-  if(TMath::Abs(vertex.At(2)) > fZvtxCut) 
+  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;
   
-  Double_t delta = 2*fZvtxCut/fNvtxBins;
-  Double_t vertexBinDouble = (vertex.At(2) + fZvtxCut) / delta;
-  Int_t vertexBin = (Int_t)vertexBinDouble;
+}
+//_____________________________________________________________________
+Bool_t 
+AliFMDBackgroundCorrection::AliFMDInputBG::ProcessEvent(UShort_t det,
+                                                       Char_t   ring, 
+                                                       UShort_t sec, 
+                                                       UShort_t strip,
+                                                       Int_t    nTrack,
+                                                       Float_t  charge)
+{
   
-  Int_t i = h->Track();
+  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;
   
-  if(h)
-    if(h->Q() !=  0 && (i != fPrevTrack || h->Detector() != fPrevDetector || h->Ring() != fPrevRing)) {
-      
-      Int_t det = h->Detector();
-      Char_t ring = h->Ring();
-      Int_t iring = (ring == 'I' ? 0 : 1);
-      
-      TObjArray* detArray  = (TObjArray*)fHitArray.At(det);
-      TObjArray* vtxArray  = (TObjArray*)detArray->At(iring);
-      TH2F* hHits          = (TH2F*)vtxArray->At(vertexBin);
-      
-      Float_t phi   = TMath::ATan2(h->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++;
-      fPrevDetector = det;
-      fPrevRing     = ring;
-    }
-   
-  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();
@@ -305,49 +438,72 @@ Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::Init() {
   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);
-                     
-           } 
-       }
-      
+  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 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;
 
-  Bool_t retVal = AliFMDInput::Begin(event); 
+    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()  {
   
-  AliStack* partStack=fLoader->Stack();
+  Bool_t retval = AliFMDInput::End();
   
-  Int_t nTracks=partStack->GetNtrack();
   AliGenEventHeader* genHeader = fLoader->GetHeader()->GenEventHeader();
   TArrayF vertex;
   genHeader->PrimaryVertex(vertex);
@@ -355,39 +511,64 @@ Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::Begin(Int_t event ) {
   if(TMath::Abs(vertex.At(2)) > fZvtxCut) 
     return kTRUE;
   
-  Double_t delta = 2*fZvtxCut/fNvtxBins;
+  Double_t delta           = 2*fZvtxCut/fNvtxBins;
   Double_t vertexBinDouble = (vertex.At(2) + fZvtxCut) / delta;
-  Int_t vertexBin = (Int_t)vertexBinDouble;
-  
-  for(Int_t j=0;j<nTracks;j++)
-       {
-         TParticle* p  = partStack->Particle(j);
-         TDatabasePDG* pdgDB = TDatabasePDG::Instance();
-         TParticlePDG* pdgPart = pdgDB->GetParticle(p->GetPdgCode());
-         Float_t charge = (pdgPart ? pdgPart->Charge() : 0);
-         Float_t phi = TMath::ATan2(p->Py(),p->Px());
-         
-         if(phi<0)
-           phi = phi+2*TMath::Pi();
-         if(p->Theta() == 0) continue;
-             Float_t eta   = -1*TMath::Log(TMath::Tan(0.5*p->Theta()));
+  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++) {
          
-         Bool_t primary = (charge!=0)&&(TMath::Abs(p->Vx() - vertex.At(0))<0.01)&&(TMath::Abs(p->Vy() - vertex.At(1))<0.01)&&(TMath::Abs(p->Vz() - vertex.At(2))<0.01);
-         if(primary) {
-           fPrim++;
-           TObjArray* innerArray = (TObjArray*)fPrimaryArray.At(0);
-           TObjArray* outerArray = (TObjArray*)fPrimaryArray.At(1);
+         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);
            
-           TH2F* hPrimaryInner  = (TH2F*)innerArray->At(vertexBin);
-           TH2F* hPrimaryOuter  = (TH2F*)outerArray->At(vertexBin);
-           hPrimaryInner->Fill(eta,phi);
-           hPrimaryOuter->Fill(eta,phi);
+           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));
          }
+         
        }
-
-  return retVal;
+      }
+    }
+  }
+  
+  fPrimaryMapInner.Reset();
+  fPrimaryMapOuter.Reset();
+  fHitMap.Reset(0);
+  fLastTrackByStrip.Reset(-1);
+  
+  return retval;
 }
-//_____________________________________________________________________
+
 //
 // EOF
 //