]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - FMD/AliFMDBackgroundCorrection.cxx
Always export to FES the regional configuration file
[u/mrichter/AliRoot.git] / FMD / AliFMDBackgroundCorrection.cxx
index d313b539ce9615403c546cfd83c3e098b8814091..dfc413afa5105bca8b7475381550b3a8565e1e55 100644 (file)
@@ -42,7 +42,6 @@
 #include "AliFMDHit.h"
 #include "AliLoader.h" 
 #include "AliFMD.h"
-#include "TH2F.h"
 #include "AliGenEventHeader.h"
 #include "AliHeader.h"
 #include "TFile.h"
@@ -71,6 +70,10 @@ AliFMDBackgroundCorrection::AliFMDInputBG::AliFMDInputBG(Bool_t hits_not_trackre
   AliFMDInput("galice.root"),
   fPrimaryArray(),
   fHitArray(),
+  fPrimaryMapInner(),
+  fPrimaryMapOuter(),
+  fHitMap(),
+  fLastTrackByStrip(),
   fPrim(0),
   fHits(0),
   fZvtxCut(0),
@@ -81,13 +84,17 @@ AliFMDBackgroundCorrection::AliFMDInputBG::AliFMDInputBG(Bool_t hits_not_trackre
   fPrevSec(-1),
   fNbinsEta(100)
 {
-  AddLoad(kTracks); 
-  if(hits_not_trackref)
+
+  if(hits_not_trackref) {
     AddLoad(kHits);
+    AddLoad(kTracks); 
+  }
   else
     AddLoad(kTrackRefs);
   AddLoad(kKinematics); 
   AddLoad(kHeader);
+  
+  
 }
 
 //_____________________________________________________________________
@@ -110,7 +117,7 @@ void AliFMDBackgroundCorrection::GenerateBackgroundCorrection(Bool_t from_hits,
     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__)
@@ -268,7 +275,7 @@ void AliFMDBackgroundCorrection::GenerateBackgroundCorrection(Bool_t from_hits,
   
   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;                         
@@ -318,12 +325,13 @@ AliFMDBackgroundCorrection::AliFMDInputBG::ProcessHit(AliFMDHit* h,
                               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) 
+AliFMDBackgroundCorrection::AliFMDInputBG::ProcessTrackRef(AliTrackReference* tr, TParticle* p) 
 {
   if(!tr)
     return kTRUE;
@@ -348,46 +356,48 @@ AliFMDBackgroundCorrection::AliFMDInputBG::ProcessEvent(UShort_t det,
                                                        Float_t  charge)
 {
   
-  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);
-  
-  if(charge !=  0 && 
-     (nTrack != fPrevTrack || 
-      det != fPrevDetector || 
-      ring != fPrevRing) ||
-     sec != fPrevSec) {
-    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);
+  /*  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) {
     
-    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) += 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;
   }
   
-  fPrevDetector = det;
-  fPrevRing     = ring;
-  fPrevSec      = sec;
-  fPrevTrack    = nTrack;
+  
   
   return kTRUE;
 
@@ -396,9 +406,18 @@ AliFMDBackgroundCorrection::AliFMDInputBG::ProcessEvent(UShort_t det,
 //_____________________________________________________________________
 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();
@@ -453,16 +472,7 @@ Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::Begin(Int_t event )
   Bool_t             retVal    = AliFMDInput::Begin(event); 
   AliStack*          partStack = fLoader->Stack();
   Int_t              nTracks   = partStack->GetNtrack();
-  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;
   
   for(Int_t j=0;j<nTracks;j++) {
     TParticle* p           = partStack->Particle(j);
@@ -472,27 +482,94 @@ Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::Begin(Int_t event )
     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()));
+    // 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++;
-      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(eta,phi);
-      hPrimaryOuter->Fill(eta,phi);
+      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
 //