#include "AliFMDHit.h"
#include "AliLoader.h"
#include "AliFMD.h"
-#include "TH2F.h"
#include "AliGenEventHeader.h"
#include "AliHeader.h"
#include "TFile.h"
//_____________________________________________________________________
AliFMDBackgroundCorrection::AliFMDInputBG::AliFMDInputBG(Bool_t hits_not_trackref) :
- AliFMDInput(),
+ AliFMDInput("galice.root"),
fPrimaryArray(),
fHitArray(),
+ fPrimaryMapInner(),
+ fPrimaryMapOuter(),
+ fHitMap(),
+ fLastTrackByStrip(),
fPrim(0),
fHits(0),
fZvtxCut(0),
fPrevTrack(-1),
fPrevDetector(-1),
fPrevRing('Q'),
+ 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);
+
+
}
//_____________________________________________________________________
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__)
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;
return kTRUE;
Bool_t retval = ProcessEvent(h->Detector(),
h->Ring(),
- h->Sector()
- ,h->Strip(),
+ 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)
+AliFMDBackgroundCorrection::AliFMDInputBG::ProcessTrackRef(AliTrackReference* tr, TParticle* p)
{
if(!tr)
return kTRUE;
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)) {
- 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;
}
- 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();
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);
- 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());
+ 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()));
+ // 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
//