+++ /dev/null
-// $Id$
-
-#include "AutoCorr.h"
-
-ClassImp(AutoCorr)
-
-Int_t AutoCorr::InitEventPools(Int_t depth,
- Int_t nMultBins, Double_t multbin[],
- Int_t nZvtxBins, Double_t zvtxbin[],
- Double_t /*ptMin*/, Double_t /*ptMax*/)
-{
- // First assign AutoCorr members
- fNMultBins = nMultBins;
- fNZvtxBins = nZvtxBins;
-
- for (Int_t iM=0; iM<nMultBins; iM++) {
- std::vector<EventPool*> evp;
- for (Int_t iZ=0; iZ<nZvtxBins; iZ++) {
- evp.push_back(new EventPool(depth,
- multbin[iM], multbin[iM+1],
- zvtxbin[iZ], zvtxbin[iZ+1]));
- }
- fEvPool.push_back(evp);
- }
-
- for (Int_t iM=0; iM<nMultBins; iM++) {
- for (Int_t iZ=0; iZ<nZvtxBins; iZ++) {
- fEvPool.at(iM).at(iZ)->SetMultBinIndex(iM);
- fEvPool.at(iM).at(iZ)->SetZvtxBinIndex(iZ);
- }
- }
-
- bool print_this = false;
- if (print_this) {
- cout << "fEvPool outer size: " << fEvPool.size() << endl;
- for (Int_t iM=0; iM<nMultBins; iM++) {
- for (Int_t iZ=0; iZ<nZvtxBins; iZ++) {
- if(fEvPool.at(iM).at(iZ)) {
- cout << "multiplicity bin: " << iM;
- cout << ", z-vertex bin: " << iZ;
- fEvPool.at(iM).at(iZ)->PrintInfo();
- }
- }
- }
- }
-
- return 0;
-}
-
-EventPool* AutoCorr::GetEventPool(Int_t iMult, Int_t iZvtx) const
-{
- if (iMult < 0 || iMult >= fNMultBins) return 0x0;
- if (iZvtx < 0 || iZvtx >= fNZvtxBins) return 0x0;
- return fEvPool.at(iMult).at(iZvtx);
-}
-
-Int_t AutoCorr::UpdatePools(Int_t iEvent, const MyHeader* ev, TClonesArray* trk)
-{
- for (Int_t iM=0; iM<fNMultBins; iM++) {
- for (Int_t iZ=0; iZ<fNZvtxBins; iZ++) {
- fEvPool.at(iM).at(iZ)->UpdatePool(iEvent, ev, trk);
- }
- }
- return 0;
-}
-
-Double_t AutoCorr::DeltaPhi(const MyPart &t1, const MyPart &t2,
- Double_t rangeMin, Double_t rangeMax) const
-{
- Double_t dphi = -999;
- Double_t pi = TMath::Pi();
- Double_t phia = t1.fPhi;
- Double_t phib = t2.fPhi;
-
- if (phia < 0) phia += 2*pi;
- else if (phia > 2*pi) phia -= 2*pi;
- if (phib < 0) phib += 2*pi;
- else if (phib > 2*pi) phib -= 2*pi;
- dphi = phib - phia;
- if (dphi < rangeMin) dphi += 2*pi;
- else if (dphi > rangeMax) dphi -= 2*pi;
-
- return dphi;
-}
-
-Double_t AutoCorr::DeltaEta(const MyPart &t1, const MyPart &t2) const
-{
- return t1.fEta - t2.fEta;
-}
-
-Bool_t AutoCorr::InBounds(Double_t val, Double_t min, Double_t max) const
-{
- if (val<min)
- return 0;
- if (val>max)
- return 0;
- return 1;
-}
-
-Bool_t AutoCorr::InBounds(Int_t val, Int_t min, Int_t max) const
-{
- if (val<min)
- return 0;
- if (val>max)
- return 0;
- return 1;
-}
-
-Bool_t AutoCorr::IsEventOk(const MyHeader &ev, Int_t minVc,
- Int_t maxNTracklets, Double_t zMin, Double_t zMax,
- Int_t trbit) const
-{
- Bool_t VcOk = ev.fVc >= minVc;
- Bool_t NTrackletsOK = ev.fNTracklets <= maxNTracklets;
- Bool_t zOk = InBounds(ev.fVz, zMin, zMax);
- Bool_t trOk = 1;
- if (trbit)
- trOk = (ev.fTrClassMask>>trbit&1);
- return (!ev.fIsPileupSPD && VcOk && NTrackletsOK && zOk && trOk);
-}
-
-Bool_t AutoCorr::IsTrackOk(const MyPart &t, Double_t etaMin, Double_t etaMax) const
-{
- return InBounds(t.fEta, etaMin, etaMax);
-}
-
-Bool_t AutoCorr::IsTrackOk(const MyPart &t, Double_t etaMin, Double_t etaMax,
- Double_t ptMin, Double_t ptMax) const
-{
- Bool_t etaOk = InBounds(t.fEta, etaMin, etaMax);
- Bool_t ptOk = InBounds(t.fPt, ptMin, ptMax);
- return etaOk && ptOk;
-}
-
-Bool_t AutoCorr::IsPairOk(const MyPart &t1, const MyPart &t2) const
-{
- Double_t deta = DeltaEta(t1, t2);
- Double_t dphi = DeltaPhi(t1, t2);
- Double_t dpmax = 0.03;
- Double_t demax = 0.03;
- Double_t dr = dphi*dphi/(dpmax*dpmax) + deta*deta/(demax*demax);
- return (dr > 1);
-}
-
-Bool_t AutoCorr::IsMixedPairOk(const MyPart &t1, const MyPart &t2) const
-{
- Double_t deta = DeltaEta(t1, t2);
- Double_t dphi = DeltaPhi(t1, t2);
- Double_t dpmax = 0.03;
- Double_t demax = 0.03;
- Double_t dr = dphi*dphi/(dpmax*dpmax) + deta*deta/(demax*demax);
- return (dr > 1);
-}
-
-Double_t AutoCorr::InvMass(const MyPart &p1, const MyPart &p2) const
-{
- Double_t px1 = p1.Px();
- Double_t py1 = p1.Py();
- Double_t pz1 = p1.Pz();
- Double_t px2 = p2.Px();
- Double_t py2 = p2.Py();
- Double_t pz2 = p2.Pz();
- Double_t pm1 = TMath::Sqrt(px1*px1+py1*py1+pz1*pz1);
- Double_t pm2 = TMath::Sqrt(px2*px2+py2*py2+pz1*pz2);
- Double_t p12 = px1*px2+py1*py2+pz1*pz2;
- Double_t m = TMath::Sqrt(pm1*pm2-p12);
- return m;
-}
+++ /dev/null
-// $Id$
-
-#ifndef AutoCorr_h
-#define AutoCorr_h
-
-#include <Rtypes.h>
-#include <Riostream.h>
-#include <TMath.h>
-#include <vector>
-#include <TRandom.h>
-#include "TreeClasses.h"
-#include "EventPool.h"
-
-class AutoCorr : public TObject
-{
-public:
- AutoCorr(){;}
- ~AutoCorr(){;}
-
- Double_t DeltaPhi(const MyPart &t1, const MyPart &t2,
- Double_t rangeMin = -TMath::Pi()/2,
- Double_t rangeMax = 3*TMath::Pi()/2) const;
- Double_t DeltaEta(const MyPart &t1, const MyPart &t2) const;
- Bool_t IsMixedPairOk(const MyPart &t1, const MyPart &t2) const;
- Bool_t IsPairOk(const MyPart &t1, const MyPart &t2) const;
- Bool_t IsTrackOk(const MyPart &t,
- Double_t etaMin, Double_t etaMax) const;
- Bool_t IsTrackOk(const MyPart &t,
- Double_t etaMin, Double_t etaMax,
- Double_t ptMin, Double_t ptMax) const;
- Double_t InvMass(const MyPart &p1, const MyPart &p2) const;
- Bool_t IsEventOk(const MyHeader &ev, Int_t minVc,
- Int_t maxNTracklets, Double_t zMin,
- Double_t zMax, Int_t trmask = -1) const;
- Bool_t InBounds(Double_t val, Double_t min, Double_t max) const;
- Bool_t InBounds(Int_t val, Int_t min, Int_t max) const;
-
- EventPool* GetEventPool(Int_t iMult, Int_t iZvtx) const;
- Int_t InitEventPools(Int_t depth,
- Int_t nmultbins, Double_t multbins[],
- Int_t nzvtxbins, Double_t zvtxbins[],
- Double_t ptMin, Double_t ptMax);
- Int_t UpdatePools(Int_t iEvent, const MyHeader* ev, TClonesArray* trk);
-
-protected:
- Int_t fNMultBins; // mult bins
- Int_t fNZvtxBins; // vertex bins
- std::vector<std::vector<EventPool*> > fEvPool; // pool in bins of [fMultBin][fZvtxBin]
-
- ClassDef(AutoCorr,1)
-};
-#endif
+++ /dev/null
-#if !defined(__CINT__) || defined(__MAKECINT__)
-#include <Riostream.h>
-#include <TFile.h>
-#include <TKey.h>
-#include <TTree.h>
-#include <TMath.h>
-#include <TChain.h>
-#include <TH2F.h>
-#include <TSystem.h>
-#include <TTimeStamp.h>
-#include <TStyle.h>
-#include <TLegend.h>
-#include <vector>
-#include "TreeClasses.h"
-#include "Utils.h"
-#endif
-
-Bool_t printPDF = 1;
-const Int_t nCuts = 3;
-enum eCuts {CUT0=0, CUT1, CUT2}; // may split this into ev and trk
-enum {EVT=0, TRK=1};
-Int_t fillCol[2][nCuts] = {{kBlue, kYellow, kSpring-5},
- {kAzure+1, kGreen, kRed}};
-Int_t lineCol=kBlack, mkrCol=kBlack, mkrSty=kFullCircle;
-Float_t mkrSize=1.0;
-TLegend* leg[2];
-TString cutDef[2][nCuts];
-TFile* outFile = new TFile("out.root", "recreate");
-const Double_t PI = TMath::Pi();
-Int_t nEvents = -1; // Run over whole chain if < 0
-MyHeader *event = 0;
-TClonesArray *tracks = 0;
-TBranch *hBranch = 0;
-TBranch *tBranch = 0;
-
-// Event
-vector<TH1F*> vhBx; // crossing times (ns)
-vector<TH1F*> vhNChips1;
-vector<TH1F*> vhNChips2;
-vector<TH1F*> vhNTracks;
-vector<TH1F*> vhNSelTracks; // no kink daughters, NClTPC > 20.
-vector<TH1F*> vhNTracklets;
-vector<TH1F*> vhVx;
-vector<TH1F*> vhVy;
-vector<TH1F*> vhVz; // from V0 (?)
-vector<TH1F*> vhVc; // from TPC (?)
-vector<TH1F*> vhIsPileupSPD;
-vector<TH1F*> vhNPileupSPD;
-vector<TH1F*> vhNPileup;
-vector<TH1F*> vhTrClassMask; // get names from constantin
-vector<TH1F*> vhTrCluster;
-vector<TH1F*> vhVxSPD;
-vector<TH1F*> vhVySPD;
-vector<TH1F*> vhVzSPD;
-vector<TH1F*> vhVcSPD;
-vector<TH1F*> vhVxTPC;
-vector<TH1F*> vhVyTPC;
-vector<TH1F*> vhVzTPC;
-vector<TH1F*> vhVcTPC;
-
-// Track
-vector<TH1F*> vhSt; // status
-vector<TH1F*> vhC; // charge
-vector<TH1F*> vhPt;
-vector<TH1F*> vhEta;
-vector<TH1F*> vhPhi;
-vector<TH1F*> vhNClTPC;
-vector<TH1F*> vhNClTPC1; // ?
-vector<TH1F*> vhNClTPCShared;
-vector<TH1F*> vhNClITS;
-vector<TH1F*> vhChi2TPC;
-vector<TH1F*> vhChi2TPC1;
-vector<TH1F*> vhChi2ITS;
-vector<TH1F*> vhD; // Transverse DCA or "impact parameter"
-vector<TH1F*> vhZ; // Z DCA
-vector<TH1F*> vhDTPC;
-vector<TH1F*> vhZTPC;
-
-vector<TH1F*> vhStBits;
-vector<TH1F*> vhL0Bits;
-vector<TH1F*> vhTrClBits;
-vector<TH1F*> vhChi2ndfTPC;
-vector<TH1F*> vhChi2ndfITS; // <---- Still need to make this
-
-vector<TH2F*> vhDZ;
-vector<TH2F*> vhDZTPC;
-
-void BookHistos();
-TObjArray* GetHistList(TFile& inFile, TString clname);
-void SetHistProps();
-void DrawHistos(TFile& inFile);
-void FillEventInfo(const MyHeader& ev, const Int_t k);
-void FillTrackInfo(const MyPart& trk, const Int_t k);
-
-void DataChecker(const char* inFileNames = "../rootfiles/res_LHC10c_09212010/merged_*.root")
-{
- // --- TTree/TChain input ---
- Noti *nme = new Noti;
- TChain *tree = new TChain("MyTree");
- tree->Add(inFileNames);
- tree->SetNotify(nme);
-
- Int_t nents = tree->GetEntries();
- cout << "Found " << nents << " entries!" << endl;
- if (nents<=0)
- return;
- if (nEvents>0 && nEvents<nents) {
- nents=nEvents;
- cout << "Using " << nents << " entries!" << endl;
- }
-
- BookHistos();
-
- // --- Event loop ---
- for (Int_t i=0;i<nents;++i) {
- Int_t li = tree->LoadTree(i);
- if (nme->Notified()) {
- hBranch = tree->GetBranch("header");
- hBranch->SetAddress(&event);
- tBranch = tree->GetBranch("parts");
- tBranch->SetAddress(&tracks);
- nme->Reset();
- // get bits from maps here and fill array of names
-#if 0
-TFile f("rootfiles/merged_run119037.root")
-f.ls()
-TrigClass_run119037_map->ls()
-TrigClass_run119037_map->FindObject("CINT1B-ABCE-NOPF-ALL")
-TPair *p = TrigClass_run119037_map->FindObject("CINT1B-ABCE-NOPF-ALL")
-p->Value()->GetName()
-#endif
- }
-
- hBranch->GetEntry(li);
- tBranch->GetEntry(li);
-
- Bool_t evCut0=1, evCut1=1, evCut2=1;
- Bool_t trCut0=1, trCut1=1, trCut2=1;
- evCut0 = 1;
- evCut1 = !event->fIsPileupSPD && fabs(event->fVz) < 10;
- evCut2 = event->fVc > 15;
-
- if (evCut0) {
- FillEventInfo(*event, CUT0);
- if (evCut1) {
- FillEventInfo(*event, CUT1);
- if (evCut2) {
- FillEventInfo(*event, CUT2);
- }
- }
- }
-
- // --- Track loop ---
- Int_t ntracks = tracks->GetEntries();
- for (Int_t t=0;t<ntracks;++t) {
- MyPart *trk = (MyPart*)tracks->At(t);
-
- trCut0 = evCut1;
- trCut1 = trk->fNClTPC > 0;
- trCut2 = trk->fNClTPC > 70;
-
- if (trCut0) {
- FillTrackInfo(*trk, CUT0);
- if (trCut1) {
- FillTrackInfo(*trk, CUT1);
- if (trCut2) {
- FillTrackInfo(*trk, CUT2);
- }
- }
- }
- }
- }
-
- outFile->Write();
- SetHistProps();
- DrawHistos(*outFile);
-
- return;
-}
-
-void BookHistos()
-{
- for (Int_t k=0; k<nCuts; k++) {
- // EVENT variables
- vhBx .push_back(new TH1F(Form("EVThBx_%d",k),"", 3000, 0, 3000));
- vhNChips1 .push_back(new TH1F(Form("EVThNChips1_%d",k),"", 150, 0, 150));
- vhNChips2 .push_back(new TH1F(Form("EVThNChips2_%d",k),"", 150, 0, 150));
- vhNTracks .push_back(new TH1F(Form("EVThNTracks_%d",k),"", 500, 0, 500));
- vhNSelTracks .push_back(new TH1F(Form("EVThNSelTracks_%d",k),"", 500, 0, 500));
- vhNTracklets .push_back(new TH1F(Form("EVThNTracklets_%d",k),"", 100, 0, 100));
- vhVx .push_back(new TH1F(Form("EVThVx_%d",k),"", 100, -0.1, 0.0));
- vhVy .push_back(new TH1F(Form("EVThVy_%d",k),"", 100, 0.15, 0.25));
- vhVz .push_back(new TH1F(Form("EVThVz_%d",k),"", 40, -20, 20));
- vhVc .push_back(new TH1F(Form("EVThVc_%d",k),"", 100, 0, 100));
- vhIsPileupSPD .push_back(new TH1F(Form("EVThIsPileupSPD_%d",k),"", 3, 0, 3));
- vhNPileupSPD .push_back(new TH1F(Form("EVThNPileupSPD_%d",k),"", 10, 0, 10));
- vhNPileup .push_back(new TH1F(Form("EVThNPileup_%d",k),"", 10, 0, 10));
- vhTrCluster .push_back(new TH1F(Form("EVThTrCluster_%d",k),"", 10, 0, 10));
- vhVxSPD .push_back(new TH1F(Form("EVThVxSPD_%d",k),"", 100, -0.5, 0.5));
- vhVySPD .push_back(new TH1F(Form("EVThVySPD_%d",k),"", 100, -0.5, 0.5));
- vhVzSPD .push_back(new TH1F(Form("EVThVzSPD_%d",k),"", 100, -20, 20));
- vhVcSPD .push_back(new TH1F(Form("EVThVcSPD_%d",k),"", 100, 0, 100));
- vhVxTPC .push_back(new TH1F(Form("EVThVxTPC_%d",k),"", 100, -0.01, 0.01));
- vhVyTPC .push_back(new TH1F(Form("EVThVyTPC_%d",k),"", 100, -0.01, 0.01));
- vhVzTPC .push_back(new TH1F(Form("EVThVzTPC_%d",k),"", 100, -5, 5));
- vhVcTPC .push_back(new TH1F(Form("EVThVcTPC_%d",k),"", 100, -1, 1));
- vhL0Bits .push_back(new TH1F(Form("EVThL0Bits_%d",k),"", 32, 0, 32));
- vhTrClBits .push_back(new TH1F(Form("EVThTrClBits_%d",k),"", 64, 0, 64));
-
- // TRACK variables
- vhC .push_back(new TH1F(Form("TRKhC_%d",k),"", 5, -2.5, 2.5));
- vhPt .push_back(new TH1F(Form("TRKhPt_%d",k),"", 200, 0, 20));
- vhEta .push_back(new TH1F(Form("TRKhEta_%d",k),"", 600, -3, 3));
- vhPhi .push_back(new TH1F(Form("TRKhPhi_%d",k),"", 100, -PI/6, 2*PI+PI/6));
- vhNClTPC .push_back(new TH1F(Form("TRKhNClTPC_%d",k),"", 200, 0, 200));
- vhNClTPC1 .push_back(new TH1F(Form("TRKhNClTPC1_%d",k),"", 200, 0, 200));
- vhNClTPCShared.push_back(new TH1F(Form("TRKhNClTPCShared_%d",k),"", 200, 0, 200));
- vhNClITS .push_back(new TH1F(Form("TRKhNClITS_%d",k),"", 20, 0, 20));
- vhChi2TPC .push_back(new TH1F(Form("TRKhChi2TPC_%d",k),"", 100, 0, 100));
- vhChi2ndfTPC .push_back(new TH1F(Form("TRKhChi2ndfTPC_%d",k),"", 100, 0, 100));
- vhChi2TPC1 .push_back(new TH1F(Form("TRKhChi2TPC1_%d",k),"", 100, 0, 100));
- vhChi2ITS .push_back(new TH1F(Form("TRKhChi2ITS_%d",k),"", 100, 0, 100));
- vhD .push_back(new TH1F(Form("TRKhD_%d",k),"", 200, -0.1, 0.1));
- vhZ .push_back(new TH1F(Form("TRKhZ_%d",k),"", 200, -0.1, 0.1));
- vhDTPC .push_back(new TH1F(Form("TRKhDTPC_%d",k),"", 200, -0.1, 0.1));
- vhZTPC .push_back(new TH1F(Form("TRKhZTPC_%d",k),"", 200, -0.1, 0.1));
- vhStBits .push_back(new TH1F(Form("TRKhStBits_%d",k),"", 32, 0, 32));
- vhDZ .push_back(new TH2F(Form("TRKhDZ_%d",k),"", 100,-0.1,0.1, 100,-0.1,0.1));
- vhDZTPC .push_back(new TH2F(Form("TRKhDZTPC_%d",k),"", 100,-0.1,0.1, 100,-0.1,0.1));
- }
-
- return;
-}
-
-void SetHistProps()
-{
- for (Int_t k=0; k<nCuts; k++) {
- vhBx .at(k)->SetTitle("Bx;bunch crossing time [ns];events");
- vhNChips1 .at(k)->SetTitle("NChips1;chips - SPD layer 1;events");
- vhNChips2 .at(k)->SetTitle("NChips2;chips - SPD layer 2;events");
- vhNTracks .at(k)->SetTitle("NTracks;reconstructed tracks;events");
- vhNSelTracks .at(k)->SetTitle("NSelTracks;selected tracks;events");
- vhNTracklets .at(k)->SetTitle("NTracklets;tracklets;events");
- vhVx .at(k)->SetTitle("x-vertex;v_{x} [cm];events");
- vhVy .at(k)->SetTitle("y-vertex;v_{y} [cm];events");
- vhVz .at(k)->SetTitle("z-vertex;v_{z} [cm];events");
- vhVc .at(k)->SetTitle("Vertex contributors;tracks;events");
- vhIsPileupSPD .at(k)->SetTitle("IsPileupSPD;IsPileupSPD;events");
- vhNPileupSPD .at(k)->SetTitle("NPileupSPD;pileup events;events");
- vhNPileup .at(k)->SetTitle("NPileup;pileup events;events");
- vhTrCluster .at(k)->SetTitle("Trigger cluster;trigger cluster;events");
- vhVxSPD .at(k)->SetTitle("SPD x-vertex;v_{x} [cm];events");
- vhVySPD .at(k)->SetTitle("SPD y-vertex;v_{y} [cm];events");
- vhVzSPD .at(k)->SetTitle("SPD z-vertex;v_{z} [cm];events");
- vhVcSPD .at(k)->SetTitle("SPD vertex contributors;tracks;events");
- vhVxTPC .at(k)->SetTitle("TPC x-vertex;v_{x} [cm];events");
- vhVyTPC .at(k)->SetTitle("TPC y-vertex;v_{y} [cm];events");
- vhVzTPC .at(k)->SetTitle("TPC z-vertex;v_{z} [cm];events");
- vhVcTPC .at(k)->SetTitle("TPC vertex contributors;tracks;events");
- vhL0Bits .at(k)->SetTitle("L0 trigger bits;;events");
- vhTrClBits .at(k)->SetTitle("Trigger class bits;;events");
- vhC .at(k)->SetTitle("Charge;charge;tracks");
- vhPt .at(k)->SetTitle("Pt;p_{T} [GeV/c];tracks");
- vhEta .at(k)->SetTitle("Eta;#eta;tracks");
- vhPhi .at(k)->SetTitle("Phi;#phi;tracks");
- vhNClTPC .at(k)->SetTitle("TPC clusters in track;TPC clusters;tracks");
- vhNClTPC1 .at(k)->SetTitle("TPC clusters in track (iter 1);TPC clusters;tracks");
- vhNClTPCShared.at(k)->SetTitle("NClTPCShared;ITS/TPC clusters;tracks");
- vhNClITS .at(k)->SetTitle("ITS clusters;clusters;tracks");
- vhChi2TPC .at(k)->SetTitle("TPC #chi^{2};total #chi^{2}_{TPC};tracks");
- vhChi2TPC1 .at(k)->SetTitle("TPC #chi^{2} (iter 1);total #chi^{2}_{TPC};tracks");
- vhChi2ITS .at(k)->SetTitle("ITS #chi^{2};#chi^{2}_{ITS};tracks");
- vhD .at(k)->SetTitle("Transverse DCA;x-y impact parameter [cm];tracks");
- vhZ .at(k)->SetTitle("Longitudinal DCA;z impact parameter [cm];tracks");
- vhDTPC .at(k)->SetTitle("Transverse TPC DCA;x-y impact parameter [cm];tracks");
- vhZTPC .at(k)->SetTitle("Longitudinal TPC DCA;z impact parameter [cm];tracks");
- vhStBits .at(k)->SetTitle("Track status flags;;tracks");
- vhDZ .at(k)->SetTitle("Impact parameter;x-y DCA [cm];z DCA [cm]");
- vhDZTPC .at(k)->SetTitle("TPC Impact parameter;x-y DCA [cm];z DCA [cm]");
- vhChi2ndfTPC .at(k)->SetTitle("#chi^{2} / N_{TPC clusters};#chi^{2}_{TPC} / N;tracks");
-
- for (Int_t n=0; n<32; ++n) {
- vhStBits.at(k)->GetXaxis()->SetBinLabel(n+1, flagLabel[n].Data());
- }
- }
-
- return;
-}
-
-void FillEventInfo(const MyHeader& ev, const Int_t k)
-{
- if (k>=nCuts) {
- cout << "ERROR in FillEventInfo(): bad cut index "
- << k << endl;
- return;
- }
- vhBx .at(k)->Fill(ev.fBx );
- vhNChips1 .at(k)->Fill(ev.fNChips1 );
- vhNChips2 .at(k)->Fill(ev.fNChips2 );
- vhNTracks .at(k)->Fill(ev.fNTracks );
- vhNSelTracks .at(k)->Fill(ev.fNSelTracks );
- vhNTracklets .at(k)->Fill(ev.fNTracklets );
- vhVx .at(k)->Fill(ev.fVx );
- vhVy .at(k)->Fill(ev.fVy );
- vhVz .at(k)->Fill(ev.fVz );
- vhVc .at(k)->Fill(ev.fVc );
- vhIsPileupSPD.at(k)->Fill(ev.fIsPileupSPD);
- vhNPileupSPD .at(k)->Fill(ev.fNPileupSPD );
- vhNPileup .at(k)->Fill(ev.fNPileup );
- vhTrCluster .at(k)->Fill(ev.fTrCluster );
- vhVxSPD .at(k)->Fill(ev.fVxSPD );
- vhVySPD .at(k)->Fill(ev.fVySPD );
- vhVzSPD .at(k)->Fill(ev.fVzSPD );
- vhVcSPD .at(k)->Fill(ev.fVcSPD );
- vhVxTPC .at(k)->Fill(ev.fVxTPC );
- vhVyTPC .at(k)->Fill(ev.fVyTPC );
- vhVzTPC .at(k)->Fill(ev.fVzTPC );
- vhVcTPC .at(k)->Fill(ev.fVcTPC );
-
- /*
- http://root.cern.ch/root/html/ListOfTypes.html
- UChar_t Unsigned Character 1 byte (unsigned char)
- UInt_t Unsigned integer 4 bytes (unsigned int)
- ULong64_t Portable unsigned long integer 8 bytes
- ULong_t Unsigned long integer 8 bytes (unsigned long)
- */
- for (UInt_t n=0; n<32; n++) {
- ULong64_t uno = 1;
- Bool_t bit = ev.fL0 & uno<<n;
- vhL0Bits.at(k)->Fill(n+0.5, bit);
- }
- for (ULong64_t n=0; n<64; n++) {
- ULong64_t uno = 1;
- Bool_t bit = ev.fTrClassMask & uno<<n;
- vhTrClBits.at(k)->Fill(n+0.5, bit);
- }
- return;
-}
-
-void FillTrackInfo(const MyPart& trk, const Int_t k)
-{
- // N.B. a pt > 0.4 cut is applied outside this fn.
- if (k>=nCuts) {
- cout << "ERROR in FillTrackInfo(): bad cut index "
- << k << endl;
- return;
- }
-
- Double_t chi2tpc = trk.fNClTPC > 0? trk.fChi2TPC/trk.fNClTPC : 999;
- vhC .at(k)->Fill(trk.fC);
- vhPt .at(k)->Fill(trk.fPt);
- vhEta .at(k)->Fill(trk.fEta);
- vhPhi .at(k)->Fill(trk.fPhi);
- vhNClTPC .at(k)->Fill(trk.fNClTPC);
- vhNClTPC1 .at(k)->Fill(trk.fNClTPC1);
- vhNClTPCShared.at(k)->Fill(trk.fNClTPCShared);
- vhNClITS .at(k)->Fill(trk.fNClITS);
- vhChi2TPC .at(k)->Fill(trk.fChi2TPC);
- vhChi2TPC1 .at(k)->Fill(trk.fChi2TPC1);
- vhChi2ITS .at(k)->Fill(trk.fChi2ITS);
- vhD .at(k)->Fill(trk.fD);
- vhZ .at(k)->Fill(trk.fZ);
- vhDTPC .at(k)->Fill(trk.fDTPC);
- vhZTPC .at(k)->Fill(trk.fZTPC);
- vhChi2ndfTPC .at(k)->Fill(chi2tpc);
- vhDZ .at(k)->Fill(trk.fD, trk.fZ);
- vhDZTPC .at(k)->Fill(trk.fDTPC, trk.fZTPC);
-
- // fSt is a ULong_t (8 bytes), but only the first 32 seem to make sense
- for (ULong_t n=0; n<32; ++n) {
- ULong64_t uno = 1;
- Bool_t bit = trk.fSt & uno<<n; // 0 or 1 * 2^n, downcasted
- vhStBits.at(k)->Fill(n+0.5, bit);
- }
-
- return;
-}
-
-TObjArray* GetHistList(TFile& inFile, TString clname)
-{
- TObjArray* hList = new TObjArray();
- TIter next(inFile.GetListOfKeys());
- TKey *key;
-
- while ((key=(TKey*)next())) {
- TString className(key->GetClassName());
- TString keyName(key->GetName());
- if (0) printf("%10s %20s\n", className.Data(), keyName.Data());
-
- if (className.Contains(clname) && clname=="TH1F") {
- // Transpose from how histos are listed in file
- if (keyName.EndsWith("_0")) {
- hList->Add((TH1F*)inFile.Get(keyName.Data()));
- for (Int_t k=1; k<nCuts; k++) {
- TString name(keyName.Replace(keyName.Length()-1, 1, Form("%d",k)));
- // printf("%20s %20s\n", keyName.Data(), name.Data());
- hList->Add((TH1F*)inFile.Get(name.Data()));
- }
- }
- }
- if (className.Contains(clname) && clname=="TH2F") {
- // Transpose from how histos are listed in file
- if (keyName.EndsWith("_0")) {
- hList->Add((TH2F*)inFile.Get(keyName.Data()));
- for (Int_t k=1; k<nCuts; k++) {
- TString name(keyName.Replace(keyName.Length()-1, 1, Form("%d",k)));
- // printf("%20s %20s\n", keyName.Data(), name.Data());
- hList->Add((TH2F*)inFile.Get(name.Data()));
- }
- }
- }
- }
- return hList;
-}
-
-void DrawHistos(TFile& inFile)
-{
- gStyle->SetTitleFontSize(0.06);
- gStyle->SetOptStat(0);
-
- TObjArray* cList = new TObjArray();
- TObjArray* h1List = GetHistList(inFile, "TH1F");
- TObjArray* h2List = GetHistList(inFile, "TH2F");
-
- PlotUtils::set_hist_props(h1List,lineCol,kNone,mkrCol,mkrSty,mkrSize);
-
- Int_t nCanv = 0, nPrinted = 0;
- TCanvas* c;
-
- for (int i=0; i<2; i++) {
- leg[i] = new TLegend(0.6, 0.85, 0.99, 0.99);
- leg[i]->SetFillColor(kNone);
- leg[i]->SetBorderSize(0);
- }
- cutDef[EVT][0] = "no cuts";
- cutDef[EVT][1] = "!fIsPileupSPD && fabs(fVz) < 10";
- cutDef[EVT][2] = "fVc > 15";
- cutDef[TRK][0] = "!fIsPileupSPD && fabs(fVz) < 10";
- cutDef[TRK][1] = "fNClTPC > 0";
- cutDef[TRK][2] = "fNClTPC > 70";
-
- // TH1Fs
- for (int i=0; i<h1List->GetEntries(); i++) {
-
- TH1F* h1 = (TH1F*)h1List->At(i);
- h1->SetStats(1);
- TString s(h1->GetName());
- TString options(h1->GetDrawOption());
-
- Int_t T = s.Contains("EVT")? 0 : 1;
-
- int nLegEntries = leg[T]->GetListOfPrimitives()->GetEntries();
-
- if (s.Contains("_0")) {
- h1->SetFillColor(fillCol[T][CUT0]);
- if (nLegEntries < nCuts)
- leg[T]->AddEntry(h1, cutDef[T][CUT0].Data(), "f");
- }
- if (s.Contains("_1")) {
- h1->SetFillColor(fillCol[T][CUT1]);
- options.Append("same");
- if (nLegEntries < nCuts)
- leg[T]->AddEntry(h1, cutDef[T][CUT1].Data(), "f");
- }
- if (s.Contains("_2")) {
- h1->SetFillColor(fillCol[T][CUT2]);
- options.Append("same");
- if (nLegEntries < nCuts)
- leg[T]->AddEntry(h1, cutDef[T][CUT2].Data(), "f");
- }
-
- if (!options.Contains("same")) {
- TString ylabel(h1->GetYaxis()->GetTitle());
- TString prefix("c");
- if (ylabel=="events") prefix = "ev";
- if (ylabel=="tracks") prefix = "tr";
- char* title = Form("%s_%s", prefix.Data(), h1->GetName());
- cList->Add(new TCanvas(Form("c%i", nCanv++), title, 1));
- }
- c = (TCanvas*)cList->Last();
- c->cd();
-
- if (s.Contains("hPt") ||
- s.Contains("NClTPC") ||
- s.Contains("Chi2") )
- c->SetLogy();
- if (s.Contains("hMult"))
- h1->SetStats(1);
- if (s.Contains("hStBits") ||
- s.Contains("hL0Bits") ||
- s.Contains("hTrClBits")) {
- // h1->SetFillColor(kRed);
- h1->SetBarWidth(0.9 * h1->GetBinWidth(1));
- options.Append("hbar");
- c->SetLeftMargin(0.2);
- }
- if (!options.Contains("hbar"))
- PlotUtils::make_nice_axes(c, h1, 1.5);
- if (s.Contains("hVx") || s.Contains("hVy") ||
- s.Contains("hBx") || s.Contains("hC") ||
- s.Contains("hOrbit") || s.Contains("hD") || s.Contains("hZ") ) {
- h1->GetXaxis()->SetNdivisions(5); // unclutter
-
- // if (s.Contains("_0"))
- // h1->GetYaxis()->SetRangeUser(0.0, h1->GetYaxis()->GetXmax());
- }
-
- printf("Drawing %20s opt: %s", s.Data(), options.Data());
- cout << endl;
- h1->Draw(options.Data());
- leg[T]->Draw();
-
- }
-
- // TH2Fs
- for (int i=0; i<h2List->GetEntries(); i++) {
-
- TH2F* h2 = (TH2F*)h2List->At(i);
- TString s(h2->GetName());
- TString options("colz");
-
- if (s.Contains("hDZ"))
- h2->GetXaxis()->SetNdivisions(5); // unclutter
-
- cout << "Drawing " << s.Data() << endl;
-
-
- if (!options.Contains("same")) {
- TString ylabel(h2->GetYaxis()->GetTitle());
- TString prefix("c");
- if (ylabel=="events") prefix = "ev";
- if (ylabel=="tracks") prefix = "tr";
- char* title = Form("%s_%s", prefix.Data(), h2->GetName());
- cList->Add(new TCanvas(Form("c%i", nCanv++), title, 1));
- }
- c = (TCanvas*)cList->Last();
- c->cd();
-
- PlotUtils::make_nice_axes(c, h2, 1.5);
-
- h2->Draw(options.Data());
- }
-
- if (printPDF) {
- TTimeStamp ts;
- UInt_t date = ts.GetDate();
- for (int i=0; i<cList->GetEntries(); i++) {
- TCanvas* c = (TCanvas*)cList->At(i);
- c->Print(Form("plots/pngs/%s_%i.png",
- c->GetTitle(), date));
- if (nPrinted==0) {
- c->Print(Form("plots/datacheck_%i.ps%s", date, "["));
- }
- c->Print(Form("plots/datacheck_%i.ps%s",
- date, i==cList->GetEntries()-1 ? "]" : ""));
- nPrinted++;
- }
- TString base = Form("plots/datacheck_%i", date);
- TString cmd = "ps2pdf " + base + ".ps " + base + ".pdf";
- gSystem->Exec(cmd.Data());
- }
-
- return;
-}
+++ /dev/null
-// $Id$
-
-#include "EventPool.h"
-
-ClassImp(EventPool)
-
-void EventPool::PrintInfo() const
-{
- cout << " --- --- --- " << endl;
- cout << Form("%20s: %d", "Pool capacity", fMixDepth) << endl;
- cout << Form("%20s: %d", "Current depth", Depth()) << endl;
- cout << Form("%20s: %d to %d", "Mult. range", fMultMin, fMultMax) << endl;
- cout << Form("%20s: %.1f to %.1f", "Z-vtx range", fZvtxMin, fZvtxMax) << endl;
-
- return;
-}
-
-Bool_t EventPool::EventMatchesBin(Int_t mult, Short_t zvtx) const
-{
- // N.B. Lower bin limit included; upper limit excluded.
-
- Bool_t multOK = (mult >= fMultMin && mult < fMultMax);
- Bool_t zvtxOK = (zvtx >= fZvtxMin && zvtx < fZvtxMax);
- return (multOK && zvtxOK);
-}
-
-Int_t EventPool::TracksInPool() const
-{
- Int_t ntrk=0;
- for (Int_t i=0; i<(Int_t)fEvents.size(); ++i) {
- ntrk += fNTracksInEvent.at(i);
- }
- return ntrk;
-}
-
-Int_t EventPool::SetEventMultRange(Int_t multMin, Int_t multMax)
-{
- fMultMin = multMin;
- fMultMax = multMax;
- return 0;
-}
-
-Int_t EventPool::SetEventZvtxRange(Int_t zvtxMin, Int_t zvtxMax)
-{
- fZvtxMin = zvtxMin;
- fZvtxMax = zvtxMax;
- return 0;
-}
-
-Int_t EventPool::UpdatePool(Int_t iEvent, const MyHeader *ev, TClonesArray *trk)
-{
- // Initialize at any chosen starting event
- if (!fTracks)
- fTracks = new TClonesArray("MyPart", 1000);
-
- fMult = trk->GetEntries();
- fZvtx = ev->fVz;
-
- if (!EventMatchesBin(fMult, fZvtx)) {
- fWasUpdated = false;
- return -1;
- }
-
- // Should see evsize = trsize (= fMixDepth once full).
- Int_t evsize = fEvents.size();
- Int_t ntsize = fNTracksInEvent.size();
-
- if (evsize != ntsize)
- cout << "WARNING: Event array and track counter array sizes do not match."
- << " evsize = " << evsize
- << " ntsize = " << ntsize
- << endl;
-
- Bool_t firstReachedCapacity = false;
- if (evsize == fMixDepth - 1) firstReachedCapacity = true;
-
- // Remove 0th element before appending this event
- if (evsize >= fMixDepth) {
- // TODO: find out - does popping delete the elements from memory?
- fNTracksInEvent.pop_front(); // remove first int
- fEvents.pop_front(); // remove first track array
- fEventIndex.pop_front();
- }
-
- // TODO: Clone() is expensive. Maybe a better way available?
- fTracks = (TClonesArray*) trk->Clone();
-
- fNTracksInEvent.push_back(ev->fNSelTracks);
- fEvents.push_back(fTracks);
- fEventIndex.push_back(iEvent);
-
- if (firstReachedCapacity) {
- cout << "Pool " << MultBinIndex() << ", " << ZvtxBinIndex()
- << " ready at event "<< iEvent;
- PrintInfo();
- }
-
- fWasUpdated = true;
-
- bool print = 1;
- if (fDebug && print) {
- cout << " Event " << fEventIndex.back();
- cout << " PoolDepth = " << Depth();
- cout << " NTracks = " << NTracksInCurrentEvent();
- cout << " TracksInPool = " << TracksInPool();
- }
-
- return 0;
-}
-
-MyPart* EventPool::GetRandomTrack() const
-{
- MyPart* trk = 0;
- TClonesArray* tca = 0;
- UInt_t ranEvt = gRandom->Integer(fEvents.size());
- tca = fEvents.at(ranEvt);
- UInt_t ranTrk = gRandom->Integer(tca->GetEntries());
- trk = (MyPart*)tca->At(ranTrk);
- return trk;
-}
-
-// Important!! This fn. will break if selective filling is implemented
-// (i.e. pool is not updated for every single event)
-// TODO: Implement internal counters to fix this.
-Int_t EventPool::NTracksInEvent(Int_t iEvent) const
-{
- Int_t n = -1;
- Int_t curEvent = fEventIndex.back();
- Int_t offset = curEvent - iEvent;
- // pos = position of iEvent in rolling buffer
- Int_t pos = fEventIndex.size() - offset - 1;
-
- if (offset==0) // (iEvent == curEvent)
- n = fNTracksInEvent.back();
- else if (offset < 0 || iEvent < 0) {
- n = 0;// cout << " No info for event " << iEvent << " ";
- }
- else if (offset > 0 && offset <= (int)fEventIndex.size()) {
- n = fNTracksInEvent.at(pos);
- }
- else
- cout << "Event info no longer in memory" << endl;
- return n;
-}
+++ /dev/null
-// $Id$
-
-#ifndef EventPool_h
-#define EventPool_h
-
-#include <vector>
-#include <deque>
-#include <Rtypes.h>
-#include <Riostream.h>
-#include <TClonesArray.h>
-#include <TFile.h>
-#include <TMath.h>
-#include <TRandom.h>
-#include <TSystem.h>
-#include "TreeClasses.h"
-
-class EventPool : public TObject
-{
-public:
- EventPool(Int_t d) : fTracks(0), fMixDepth(d), fMultMin(-999), fMultMax(+999),
- fZvtxMin(-999), fZvtxMax(+999), fMult(0), fZvtx(0),
- fWasUpdated(0), fMultBinIndex(0), fZvtxBinIndex(0), fDebug(0) {;}
- EventPool(Int_t d, Int_t multMin, Int_t multMax,
- Double_t zvtxMin, Double_t zvtxMax) : fTracks(0), fMixDepth(d), fMultMin(multMin), fMultMax(multMax),
- fZvtxMin(zvtxMin), fZvtxMax(zvtxMax), fMult(0), fZvtx(0),
- fWasUpdated(0), fMultBinIndex(0), fZvtxBinIndex(0), fDebug(0) {;}
- ~EventPool() {;}
-
- Bool_t EventMatchesBin(Int_t mult, Short_t zvtx) const;
- Bool_t IsPoolReady() const { return (Int_t)fEvents.size()==fMixDepth; }
- Int_t Depth() const { return fEvents.size(); }
- MyPart *GetRandomTrack() const;
- Int_t MultBinIndex() const { return fMultBinIndex; }
- Int_t NTracksInEvent(Int_t iEvent) const;
- Int_t NTracksInCurrentEvent() const { return fNTracksInEvent.back(); }
- void PrintInfo() const;
- Int_t TracksInPool() const;
- Bool_t WasUpdated() const { return fWasUpdated; }
- Int_t ZvtxBinIndex() const { return fZvtxBinIndex; }
- void SetDebug(Bool_t b) { fDebug = b; }
- Int_t SetEventMultRange(Int_t multMin, Int_t multMax);
- Int_t SetEventZvtxRange(Int_t zvtxMin, Int_t zvtxMax);
- void SetMultBinIndex(Int_t iM) { fMultBinIndex = iM; }
- void SetZvtxBinIndex(Int_t iZ) { fZvtxBinIndex = iZ; }
- Int_t UpdatePool(int iEvent, const MyHeader *ev, TClonesArray* trk);
-
-protected:
- TClonesArray *fTracks; // Copy of trk array. Refreshes each event.
- deque<TClonesArray*> fEvents; // The guts of the class. Holds TObjArrays of MyParts.
- deque<int> fNTracksInEvent; // Tracks in event
- deque<int> fEventIndex; // Original event index
- Int_t fMixDepth; // Number of evts. to mix with
- Int_t fMultMin, fMultMax; // Track multiplicity bin range
- Double_t fZvtxMin, fZvtxMax; // Event z-vertex bin range
- Int_t fMult; // Tracks in current event
- Short_t fZvtx; // Current z-vertex
- Bool_t fWasUpdated; // Evt. succesfully passed selection?
- Int_t fMultBinIndex; // Multiplicity bin
- Int_t fZvtxBinIndex; // Zvertex bin
- Int_t fDebug; // if 1 then debug on
-
- ClassDef(EventPool,2) // Event Pool class
-};
-#endif
+++ /dev/null
-// $Id$
-
-#include "EventPoolManager.h"
-
-ClassImp(GenericEventPool)
-
-void GenericEventPool::PrintInfo() const
-{
- cout << " --- --- --- " << endl;
- cout << Form("%20s: %d", "Pool capacity", fMixDepth) << endl;
- cout << Form("%20s: %d", "Current depth", Depth()) << endl;
- cout << Form("%20s: %d to %d", "Mult. range", fMultMin, fMultMax) << endl;
- cout << Form("%20s: %.1f to %.1f", "Z-vtx range", fZvtxMin, fZvtxMax) << endl;
-
- return;
-}
-
-Bool_t GenericEventPool::EventMatchesBin(Int_t mult, Double_t zvtx) const
-{
- // N.B. Lower bin limit included; upper limit excluded.
-
- Bool_t multOK = (mult >= fMultMin && mult < fMultMax);
- Bool_t zvtxOK = (zvtx >= fZvtxMin && zvtx < fZvtxMax);
- return (multOK && zvtxOK);
-}
-
-Int_t GenericEventPool::TracksInPool() const
-{
- Int_t ntrk=0;
- for (Int_t i=0; i<(Int_t)fEvents.size(); ++i) {
- ntrk += fNTracksInEvent.at(i);
- }
- return ntrk;
-}
-
-Int_t GenericEventPool::SetEventMultRange(Int_t multMin, Int_t multMax)
-{
- fMultMin = multMin;
- fMultMax = multMax;
- return 0;
-}
-
-Int_t GenericEventPool::SetEventZvtxRange(Double_t zvtxMin, Double_t zvtxMax)
-{
- fZvtxMin = zvtxMin;
- fZvtxMax = zvtxMax;
- return 0;
-}
-
-Int_t GenericEventPool::UpdatePool(Int_t iEvent, const MyHeader *ev, TObjArray *trk)
-{
- // Initialize at any chosen starting event
-
- Int_t mult = trk->GetEntries();
- Double_t zvtx = ev->fVz;
-
- if (!EventMatchesBin(mult, zvtx)) {
- fWasUpdated = false;
- return -1;
- }
-
- // Should see evsize = trsize (= fMixDepth once full).
- Int_t evsize = fEvents.size();
- Int_t ntsize = fNTracksInEvent.size();
-
- if (evsize != ntsize)
- cout << "WARNING: Event array and track counter array sizes do not match."
- << " evsize = " << evsize
- << " ntsize = " << ntsize
- << endl;
-
- Bool_t firstReachedCapacity = false;
- if (evsize == fMixDepth - 1)
- firstReachedCapacity = true;
-
- // remove 0th element before appending this event
- if (evsize >= fMixDepth) {
- TObjArray *fa = fEvents.front();
- fa->Delete();
- delete fa;
- fEvents.pop_front(); // remove first track array
- fNTracksInEvent.pop_front(); // remove first int
- fEventIndex.pop_front();
- }
-
- fNTracksInEvent.push_back(mult);
- fEvents.push_back(trk);
- fEventIndex.push_back(iEvent);
-
- if (firstReachedCapacity) {
- cout << "Pool " << MultBinIndex() << ", " << ZvtxBinIndex()
- << " ready at event "<< iEvent;
- PrintInfo();
- }
-
- fWasUpdated = true;
-
- bool print = 1;
- if (fDebug && print) {
- cout << " Event " << fEventIndex.back();
- cout << " PoolDepth = " << Depth();
- cout << " NTracks = " << NTracksInCurrentEvent();
- cout << " TracksInPool = " << TracksInPool();
- }
-
- return fEvents.size();
-}
-
-TObject* GenericEventPool::GetRandomTrack() const
-{
- UInt_t ranEvt = gRandom->Integer(fEvents.size());
- TObjArray *tca = fEvents.at(ranEvt);
- UInt_t ranTrk = gRandom->Integer(tca->GetEntries());
- TObject *trk = (TObject*)tca->At(ranTrk);
- return trk;
-}
-
-TObjArray* GenericEventPool::GetRandomEvent() const
-{
- UInt_t ranEvt = gRandom->Integer(fEvents.size());
- TObjArray *tca = fEvents.at(ranEvt);
- return tca;
-}
-
-Int_t GenericEventPool::NTracksInEvent(Int_t iEvent) const
-{
- Int_t n = -1;
- Int_t curEvent = fEventIndex.back();
- Int_t offset = curEvent - iEvent;
- Int_t pos = fEventIndex.size() - offset - 1;
-
- if (offset==0)
- n = fNTracksInEvent.back();
- else if (offset < 0 || iEvent < 0) {
- n = 0;
- }
- else if (offset > 0 && offset <= (int)fEventIndex.size()) {
- n = fNTracksInEvent.at(pos);
- }
- else
- cout << "Event info no longer in memory" << endl;
- return n;
-}
-
-ClassImp(EventPoolManager)
-
-Int_t EventPoolManager::InitEventPools(Int_t depth,
- Int_t nMultBins, Int_t *multbin,
- Int_t nZvtxBins, Double_t *zvtxbin)
-{
- // First assign EventPoolManager members
- fNMultBins = nMultBins;
- fNZvtxBins = nZvtxBins;
-
- for (Int_t iM=0; iM<nMultBins; iM++) {
- std::vector<GenericEventPool*> evp;
- for (Int_t iZ=0; iZ<nZvtxBins; iZ++) {
- evp.push_back(new GenericEventPool(depth,
- multbin[iM], multbin[iM+1],
- zvtxbin[iZ], zvtxbin[iZ+1] ));
- }
- fEvPool.push_back(evp);
- }
-
- for (Int_t iM=0; iM<nMultBins; iM++) {
- for (Int_t iZ=0; iZ<nZvtxBins; iZ++) {
- fEvPool.at(iM).at(iZ)->SetMultBinIndex(iM);
- fEvPool.at(iM).at(iZ)->SetZvtxBinIndex(iZ);
- }
- }
-
- if (0) {
- cout << "fEvPool outer size: " << fEvPool.size() << endl;
- for (Int_t iM=0; iM<nMultBins; iM++) {
- for (Int_t iZ=0; iZ<nZvtxBins; iZ++) {
- if(fEvPool.at(iM).at(iZ)) {
- cout << "multiplicity bin: " << iM;
- cout << ", z-vertex bin: " << iZ;
- fEvPool.at(iM).at(iZ)->PrintInfo();
- }
- }
- }
- }
-
- return fEvPool.size();
-}
-
-GenericEventPool *EventPoolManager::GetEventPool(Int_t iMult, Int_t iZvtx) const
-{
- if (iMult < 0 || iMult >= fNMultBins)
- return 0x0;
- if (iZvtx < 0 || iZvtx >= fNZvtxBins)
- return 0x0;
- return fEvPool.at(iMult).at(iZvtx);
-}
-
-Int_t EventPoolManager::UpdatePools(Int_t iEvent, const MyHeader* ev, TObjArray *trk)
-{
- for (Int_t iM=0; iM<fNMultBins; iM++) {
- for (Int_t iZ=0; iZ<fNZvtxBins; iZ++) {
- if (fEvPool.at(iM).at(iZ)->UpdatePool(iEvent, ev, trk)>-1)
- break;
- }
- }
- return 0;
-}
+++ /dev/null
-// $Id$
-
-#ifndef EventPoolManager_h
-#define EventPoolManager_h
-
-#include <vector>
-#include <deque>
-#include <Rtypes.h>
-#include <Riostream.h>
-#include <TClonesArray.h>
-#include <TFile.h>
-#include <TMath.h>
-#include <TRandom.h>
-#include <TSystem.h>
-#include "TreeClasses.h"
-
-class GenericEventPool : public TObject
-{
-public:
- GenericEventPool(Int_t d) :
- fMixDepth(d), fMultMin(-999), fMultMax(+999),
- fZvtxMin(-999), fZvtxMax(+999), fWasUpdated(0), fMultBinIndex(0),
- fZvtxBinIndex(0), fDebug(0) {;}
- GenericEventPool(Int_t d, Int_t multMin, Int_t multMax,
- Double_t zvtxMin, Double_t zvtxMax) :
- fMixDepth(d), fMultMin(multMin), fMultMax(multMax),
- fZvtxMin(zvtxMin), fZvtxMax(zvtxMax), fWasUpdated(0), fMultBinIndex(0),
- fZvtxBinIndex(0), fDebug(0) {;}
- ~GenericEventPool() {;}
-
- Bool_t EventMatchesBin(Int_t mult, Double_t zvtx) const;
- Bool_t IsReady() const { return (Int_t)fEvents.size()==fMixDepth; }
- Int_t Depth() const { return fEvents.size(); }
- TObject *GetRandomTrack() const;
- TObjArray *GetRandomEvent() const;
- Int_t MultBinIndex() const { return fMultBinIndex; }
- Int_t NTracksInEvent(Int_t iEvent) const;
- Int_t NTracksInCurrentEvent() const { return fNTracksInEvent.back(); }
- void PrintInfo() const;
- Int_t TracksInPool() const;
- Bool_t WasUpdated() const { return fWasUpdated; }
- Int_t ZvtxBinIndex() const { return fZvtxBinIndex; }
- void SetDebug(Bool_t b) { fDebug = b; }
- Int_t SetEventMultRange(Int_t multMin, Int_t multMax);
- Int_t SetEventZvtxRange(Double_t zvtxMin, Double_t zvtxMax);
- void SetMultBinIndex(Int_t iM) { fMultBinIndex = iM; }
- void SetZvtxBinIndex(Int_t iZ) { fZvtxBinIndex = iZ; }
- Int_t UpdatePool(int iEvent, const MyHeader *ev, TObjArray *trk);
-
-protected:
- deque<TObjArray*> fEvents; //Holds TObjArrays of MyTracklets
- deque<int> fNTracksInEvent; //Tracks in event
- deque<int> fEventIndex; //Original event index
- Int_t fMixDepth; //Number of evts. to mix with
- Int_t fMultMin, fMultMax; //Track multiplicity bin range
- Double_t fZvtxMin, fZvtxMax; //Event z-vertex bin range
- Bool_t fWasUpdated; //Evt. succesfully passed selection?
- Int_t fMultBinIndex; //Multiplicity bin
- Int_t fZvtxBinIndex; //Zvertex bin
- Int_t fDebug; //If 1 then debug on
-
- ClassDef(GenericEventPool,1) // Event pool class
-};
-
-class EventPoolManager : public TObject
-{
-public:
- EventPoolManager() : fNMultBins(0), fNZvtxBins(0) {;}
- ~EventPoolManager() {;}
-
- GenericEventPool *GetEventPool(Int_t iMult, Int_t iZvtx) const;
- Int_t InitEventPools(Int_t depth,
- Int_t nmultbins, Int_t *multbins,
- Int_t nzvtxbins, Double_t *zvtxbins);
- Int_t UpdatePools(Int_t iEvent, const MyHeader *ev, TObjArray *trk);
-
-protected:
- Int_t fNMultBins; // number mult bins
- Int_t fNZvtxBins; // number vertex bins
- std::vector<std::vector<GenericEventPool*> > fEvPool; // pool in bins of [fMultBin][fZvtxBin]
-
- ClassDef(EventPoolManager,1)
-};
-#endif
+++ /dev/null
-
-// Simple macro to plot AutoCorr output histos
-
-void Plot2pc(const Int_t trial = 6, Bool_t printPDF = false)
-{
- const char* histFile = Form("output/try%i/try%i.root", trial, trial);
-
- TObjArray* c1List = new TObjArray();
- TObjArray* c2List = new TObjArray();
- TObjArray* TH1List = new TObjArray();
- TObjArray* TH2List = new TObjArray();
-
- TFile* inFile = new TFile(histFile, "read");
-
- // TODO: Move this to another macro or a histo manager
- int nMultBinsRead = hMultBins->GetNbinsX();
- const int nMultBins = nMultBinsRead;
- int nPtBinsRead = hPtBins->GetNbinsX();
- const int nPtBins = nPtBinsRead;
-
- TString sMultBin[nMultBins];
- for (int i=0; i<nMultBins; i++) {
- double m1 = hMultBins->GetBinLowEdge(i+1);
- double m2 = m1 + hMultBins->GetBinWidth(i+1);
- sMultBin[i] = Form("%.0f - %.0f tracks", m1, m2);
- }
- TString sPtBin[nPtBins];
- for (int i=0; i<nPtBins; i++) {
- double pt1 = hPtBins->GetBinLowEdge(i+1);
- double pt2 = pt1 + hPtBins->GetBinWidth(i+1);
- sPtBin[i] = Form("%.2f - %.2f GeV/c", pt1, pt2);
- }
-
- gStyle->SetOptStat(0); // Turn on indiv. stats w/ TH1::SetStats(1)
-
-
- TIter next(inFile->GetListOfKeys());
- TKey *key;
- while ((key=(TKey*)next())) {
- TString className(key->GetClassName());
- TString keyName(key->GetName());
- printf("%10s %20s\n", className.Data(), keyName.Data());
-
- // Get TH1Fs
- if (className.Contains("TH1F")) {
- TH1F* h1 = (TH1F*)inFile->Get(keyName.Data());
- TH1List->Add(h1);
- Int_t i = TH1List->IndexOf(h1);
- c1List->Add(new TCanvas(Form("c1_%i",i),
- Form("c1_%s",keyName.Data()), 1));
- }
-
- // Get TH2Fs
- if (className.Contains("TH2F")) {
- TH2F* h2 = (TH2F*)inFile->Get(keyName.Data());
- TH2List->Add(h2);
- Int_t i = TH2List->IndexOf(h2);
- c2List->Add(new TCanvas(Form("c2_%i",i),
- Form("c2_%s",keyName.Data()), 1));
- }
-
-
- }
-
- // Give reasonable titles
- for (int i=0; i<TH2List->GetEntries(); i++) {
- for (int iM=0; iM<nMultBins; iM++) {
- for (int iPt=0; iPt<nPtBins; iPt++) {
- if (keyName.EndsWith(Form("_%i_%i", iM, iPt))) {
- TH2F* h = (TH2F*)TH2List->At(i);
- h->SetTitle(Form("%s, %s;#Delta#phi;#Delta#eta;",
- sMultBin[iM].Data(), sPtBin[iPt].Data()));
- }
- }
- }
- }
-
- cuts->Print();
-
- Int_t nPrinted = 0; // Number of canvases printed
-
- // Draw TH1s
- for (int i=0; i<c1List->GetEntries(); i++) {
- TCanvas* c = (TCanvas*)c1List->At(i);
- TH1F* h1 = (TH1F*)TH1List->At(i);
- TString s(h1->GetName());
- TString options("");
- Bool_t skipPrint = false;
-
- if (s.Contains("Bins")) {
- skipPrint = true;
- continue;
- }
- if (s.Contains("hPt") || s.Contains("hMass"))
- c->SetLogy();
- if (s.Contains("Sel") || s.Contains("hMassBg")) {
- h1->SetLineColor(kBlue);
- c = (TCanvas*)c1List->At(i-1);
- options = "same";
- skipPrint = true;
- }
- if (s.Contains("hMult"))
- h1->SetStats(1);
-
-
- c->cd();
- h1->Draw(options.Data());
-
-
- if (printPDF) {
- c->Print(Form("plots/pngs/%s_try%i.png", c->GetName(), trial));
- if (nPrinted==0)
- // "[" opens ps file but doesn't print c
- c->Print(Form("plots/try%i.ps%s", trial, "["));
-
- c->Print(Form("plots/try%i.ps%s", trial, ""));
- nPrinted++;
- }
- }
-
- // Draw TH2s
- for (int i=0; i<c2List->GetEntries(); i++) {
- TCanvas* c = (TCanvas*)c2List->At(i);
- TH2F* h2 = (TH2F*)TH2List->At(i);
- TString s(h2->GetName());
- TString options("");
- Bool_t skipPrint = false;
-
- if (s.Contains("Sig") ||
- s.Contains("Bkg") ||
- s.Contains("SB") ||
- s.Contains("hDR"))
- options = "surf1";
- if (s=="hMixEff")
- options = "colz";
-
-
- c->cd();
- h2->Draw(options.Data());
-
- if (printPDF) {
- c->Print(Form("plots/pngs/%s_try%i.png",
- c->GetName(), trial));
-
- c->Print(Form("plots/try%i.ps%s",
- trial, i==c2List->LastIndex()? "]" : ""));
-
- cout<<Form("plots/try%i.ps%s",
- trial, i==c2List->GetEntries()-1 ? "]" : "")<<endl;
- nPrinted++;
- }
- }
-
- TString base = Form("plots/try%i", trial);
- TString cmd = "ps2pdf " + base + ".ps " + base + ".pdf";
- gSystem->Exec(cmd.Data());
-
-
- return;
-}
-
+++ /dev/null
-// $Id$
-
-#include "TreeClasses.h"
-
-ClassImp(Noti)
-ClassImp(MyHeader)
-ClassImp(MyPart)
-ClassImp(MyTracklet)
-
+++ /dev/null
-// $Id$
-
-#ifndef TreeClasses_h
-#define TreeClasses_h
-
-#include <Riostream.h>
-#include <TTree.h>
-#include <TChain.h>
-#include <TCanvas.h>
-#include <TClonesArray.h>
-#include <TH1F.h>
-#include <TH2F.h>
-#include <TMath.h>
-
-class Noti: public TObject
-{
-public:
- Noti() : fc(0) {;}
- virtual ~Noti() {;}
- Bool_t Notify() { fc=1; return 1; }
- Bool_t Notified() const { return fc; }
- void Reset() { fc=0; }
-protected:
- Bool_t fc; //=1 when file changed
- ClassDef (Noti,0) // Use to be notified when file in chain changes
-};
-
-class MyHeader
-{
-public:
- MyHeader() : fRun(0), fOrbit(0), fTime(0), fPeriod(0), fBx(0), fL0(0), fL1(0), fL2(0),
- fNChips1(0), fNChips2(0), fNTracks(0), fNSelTracks(0), fNTracklets(0),
- fVx(0), fVy(0), fVz(0), fVc(-1), fIsPileupSPD(0), fNPileupSPD(0), fNPileup(0),
- fTrClassMask(0), fTrCluster(0), fEvNumberInFile(-1), fFileId(-1),
- fVxSPD(0), fVySPD(0), fVzSPD(0), fVcSPD(-1) {;}
- virtual ~MyHeader() {;}
- ULong64_t GetEventId() const {
- return (((ULong64_t)fPeriod << 36) |
- ((ULong64_t)fOrbit << 12) |
- (ULong64_t)fBx); }
- Bool_t IsL0Fired(Int_t bit) const { return fL0 & (1<<bit); }
- Bool_t IsTrClassFired(Int_t bit) const { return fTrClassMask & (1<<bit); }
-public:
- Int_t fRun; // run number
- UInt_t fOrbit; // orbit number
- UInt_t fTime; // timestamp from DAQ
- UInt_t fPeriod; // period number
- UShort_t fBx; // bunch crossing id
- UInt_t fL0; // l0 trigger bits
- UInt_t fL1; // l1 trigger bits
- UShort_t fL2; // l2 trigger bits
- Short_t fNChips1; // number of chips layer 1
- Short_t fNChips2; // number of chips layer 2
- Short_t fNTracks; // number of tracks before cuts
- Short_t fNSelTracks; // number of stored tracks
- Short_t fNTracklets; // number of tracklets
- Double_t fVx; //[0,0,16] global vertex x
- Double_t fVy; //[0,0,16] global vertex y
- Double_t fVz; //[0,0,16] global vertex z
- Double_t fVc; //[0,0,16] number of contributors to global vertex
- Bool_t fIsPileupSPD; // is pileup according to standard definition
- Char_t fNPileupSPD; // number of additional vertices in SPD
- Char_t fNPileup; // number of additional vertices in TPC
- ULong64_t fTrClassMask; // trigger class mask
- UChar_t fTrCluster; // trigger cluster mask
- Int_t fEvNumberInFile; // event number in ESD file
- Short_t fFileId; // file id to access path name
- Double_t fVxSPD; //[0,0,16] SPD vertex x
- Double_t fVySPD; //[0,0,16] SPD vertex y
- Double_t fVzSPD; //[0,0,16] SPD vertex z
- Double_t fVcSPD; //[0,0,16] number of contributors to SPD vertex
- Double_t fVxTPC; //[0,0,16] TPC vertex x
- Double_t fVyTPC; //[0,0,16] TPC vertex y
- Double_t fVzTPC; //[0,0,16] TPC vertex z
- Double_t fVcTPC; //[0,0,16] number of contributors to TPC vertex
-
- ClassDef(MyHeader,4) // My header class
-};
-
-class MyPart : public TObject
-{
-public:
- enum { /*from AliESDtrack.h */
- kITSin=0x0001,kITSout=0x0002,kITSrefit=0x0004,kITSpid=0x0008,
- kTPCin=0x0010,kTPCout=0x0020,kTPCrefit=0x0040,kTPCpid=0x0080,
- kTRDin=0x0100,kTRDout=0x0200,kTRDrefit=0x0400,kTRDpid=0x0800,
- kTOFin=0x1000,kTOFout=0x2000,kTOFrefit=0x4000,kTOFpid=0x8000,
- kHMPIDout=0x10000,kHMPIDpid=0x20000,
- kEMCALmatch=0x40000,
- kPHOSmatch=0x200000,
- kTRDbackup =0x80000,
- kTRDStop=0x20000000,
- kESDpid=0x40000000,
- kTIME=0x80000000,
- kGlobalMerge=0x08000000,
- kITSpureSA=0x10000000,
- kMultInV0=0x2000000,
- kMultSec=0x4000000
- };
-
- MyPart(ULong_t st=0, Char_t c=0, Double_t pt=0, Double_t eta=0, Double_t phi=0,
- Short_t ncltpc=0, Short_t ncltpc1=0, Short_t ncltpcs=0,Char_t nclits=0, Double_t chi2tpc=0,
- Double_t chi2tpc1=0, Double_t chi2its=0, Double_t d=0, Double_t z=0, Double_t dtpc=0,
- Double_t ztpc=0, Double_t sigma=0.) :
- TObject(), fSt(st), fC(c), fPt(pt), fEta(eta), fPhi(phi), fNClTPC(ncltpc), fNClTPC1(ncltpc1),
- fNClTPCShared(ncltpcs),fNClITS(nclits), fChi2TPC(chi2tpc), fChi2TPC1(chi2tpc1), fChi2ITS(chi2its),
- fD(d), fZ(z), fDTPC(dtpc), fZTPC(ztpc), fSigma(sigma) {;}
-
-
- Double_t Px() const { return fPt*TMath::Cos(fPhi); }
- Double_t Py() const { return fPt*TMath::Sin(fPhi); }
- Double_t Pz() const { return fPt*TMath::SinH(fEta); }
- Bool_t IsITSRefit() const { return (fSt&(ULong64_t)kITSrefit); }
- Bool_t IsTPCIn() const { return (fSt&(ULong64_t)kTPCin); }
- Bool_t IsTPCRefit() const { return (fSt&(ULong64_t)kITSrefit); }
-
-public:
- ULong64_t fSt; // status flag
- Double32_t fC; //[-1,1,2] charge
- Double32_t fPt; //[0,0,16] pt from constrained fit
- Double32_t fEta; //[0,0,10] eta from constrained fit
- Double32_t fPhi; //[0,0,10] phi from constrained fit
- Short_t fNClTPC; // number of clusters in TPC
- Short_t fNClTPC1; // number of clusters in TPC if standalone
- Short_t fNClTPCShared; // number of shared clusters in TPC
- Char_t fNClITS; // number of clusters in ITS
- Double32_t fChi2TPC; //[0,0,10] chi2 of TPC
- Double32_t fChi2TPC1; //[0,0,10] chi2 of TPC if standalone
- Double32_t fChi2ITS; //[0,0,10] chi2 of ITC
- Double32_t fD; //[0,0,16] transverse DCA
- Double32_t fZ; //[0,0,16] longitudinal DCA
- Double32_t fDTPC; //[0,0,16] transverse DCA TPC
- Double32_t fZTPC; //[0,0,16] longitudinal DCA TPC
- Double32_t fSigma; //[0,0,16] sigma to vertex
-
- ClassDef(MyPart,5) // My particle class in cylindrical coordinates
-};
-
-class MyTracklet : public TObject
-{
-public:
- MyTracklet(Double_t dphi=0, Double_t dth=0, Double_t eta=0, Double_t phi=0) :
- TObject(), fDPhi(dphi), fDTh(dth), fEta(eta), fPhi(phi) {;}
-
-public:
- Double32_t fDPhi; //[0,0,10] delta phi
- Double32_t fDTh; //[0,0,10] delta theta
- Double32_t fEta; //[0,0,10] eta
- Double32_t fPhi; //[0,0,10] phi
-
- ClassDef(MyTracklet,1) // My tracklet class
-};
-
-ULong_t flagValue[32] =
- {
- MyPart::kITSin,
- MyPart::kITSout,
- MyPart::kITSrefit,
- MyPart::kITSpid,
- MyPart::kTPCin,
- MyPart::kTPCout,
- MyPart::kTPCrefit,
- MyPart::kTPCpid,
- MyPart::kTRDin,
- MyPart::kTRDout,
- MyPart::kTRDrefit,
- MyPart::kTRDpid,
- MyPart::kTOFin,
- MyPart::kTOFout,
- MyPart::kTOFrefit,
- MyPart::kTOFpid,
- MyPart::kHMPIDout,
- MyPart::kHMPIDpid,
- MyPart::kEMCALmatch,
- MyPart::kTRDbackup,
- 0x0, // 20 missing
- MyPart::kPHOSmatch,
- 0x0, // 22-24 missing
- 0x0,
- 0x0,
- MyPart::kMultInV0,
- MyPart::kMultSec,
- MyPart::kGlobalMerge,
- MyPart::kITSpureSA,
- MyPart::kTRDStop,
- MyPart::kESDpid,
- MyPart::kTIME
- };
-
-TString flagLabel[32] =
- {
- "kITSin",
- "kITSout",
- "kITSrefit",
- "kITSpid",
- "kTPCin",
- "kTPCout",
- "kTPCrefit",
- "kTPCpid",
- "kTRDin",
- "kTRDout",
- "kTRDrefit",
- "kTRDpid",
- "kTOFin",
- "kTOFout",
- "kTOFrefit",
- "kTOFpid",
- "kHMPIDout",
- "kHMPIDpid",
- "kEMCALmatch",
- "kTRDbackup",
- "--",
- "kPHOSmatch",
- "--",
- "--",
- "--",
- "kMultInV0",
- "kMultSec",
- "kGlobalMerge",
- "kITSpureSA",
- "kTRDStop",
- "kESDpid",
- "kTIME"
- };
-#endif
+++ /dev/null
-// $Id$
-
-#include "Utils.h"
-#include "TMath.h"
-#include "TLine.h"
-#include "TList.h"
-
-using std::cout;
-using std::endl;
-
-ClassImp(PlotUtils)
-
-void PlotUtils::set_hist_props(TObjArray* arr,
- Int_t linecolor,
- Int_t fillcolor,
- Int_t markercolor,
- Int_t markerstyle,
- Double_t markersize)
-{
- TH1F* h = 0;
-
- for (Int_t i=0; i<arr->GetEntries(); i++) {
-
- if (arr->At(i)->InheritsFrom("TH1F")) {
- h = (TH1F*)arr->At(i);
- }
-
- h->SetLineColor(linecolor);
- h->SetFillColor(fillcolor);
- h->SetMarkerColor(markercolor);
- h->SetMarkerStyle(markerstyle);
- h->SetMarkerSize(markersize);
- }
-}
-
-void PlotUtils::set_hist_props(TH1F* h,
- Int_t linecolor,
- Int_t fillcolor,
- Int_t markercolor,
- Int_t markerstyle,
- Double_t markersize)
-{
- h->SetLineColor(linecolor);
- h->SetFillColor(fillcolor);
- h->SetMarkerColor(markercolor);
- h->SetMarkerStyle(markerstyle);
- h->SetMarkerSize(markersize);
-}
-
-void PlotUtils::set_tgraph_props(TGraphErrors* h,
- Int_t linecolor,
- Int_t markercolor,
- Int_t markerstyle,
- Double_t markersize)
-{
- h->SetLineColor(linecolor);
- h->SetMarkerColor(markercolor);
- h->SetMarkerStyle(markerstyle);
- h->SetMarkerSize(markersize);
- h->SetLineWidth(2);
-}
-
-void PlotUtils::set_tgraph_props(TGraphAsymmErrors* h,
- Int_t linecolor,
- Int_t markercolor,
- Int_t markerstyle,
- Double_t markersize)
-{
- h->SetLineColor(linecolor);
- h->SetMarkerColor(markercolor);
- h->SetMarkerStyle(markerstyle);
- h->SetMarkerSize(markersize);
- h->SetLineWidth(2);
-}
-
-void PlotUtils::shade(TF1 * fup, TF1 * flow)
-{
- TGraph * gr = new TGraph();
- gr->SetFillColor(fup->GetFillColor());
- gr->SetFillStyle(fup->GetFillStyle());
- Double_t xmin = fup->GetXmin();
- Double_t xmax = fup->GetXmax();
-
- //process top side function
- Int_t npx = flow->GetNpx();
- Int_t npoints=0;
- Double_t dx = (xmax-xmin)/npx;
- Double_t x = xmin+0.5*dx;
- Double_t y =0;
- while (x <= xmax) {
- Double_t yup = fup->Eval(x);
- Double_t ylow = flow->Eval(x);
-
- if ( yup > ylow )
- y=yup;
- else
- y=ylow;
- gr->SetPoint(npoints,x,y);
- npoints++;
- x += dx;
- }
- //process bottom side
- x = xmax-0.5*dx;
- while (x >= xmin) {
- Double_t ylow = flow->Eval(x);
- Double_t yup = fup->Eval(x);
- if ( yup < ylow )
- y=yup;
- else
- y=ylow;
- gr->SetPoint(npoints,x,y);
- npoints++;
- x -= dx;
- }
- gr->Draw("f");
-}
-
-void PlotUtils::draw_errorbox(TH1F* mid, TH1F* hi, TH1F* lo, Int_t color, Double_t dx)
-{
- Double_t x, y, dy_hi, dy_lo;
-
- for (Int_t j=1; j <= mid->GetNbinsX(); j++) {
-
- if ( ! mid->GetBinContent(j) ) continue;
-
- x = mid->GetBinCenter(j);
- dx = 0.08;
- y = mid->GetBinContent(j);
- dy_hi = fabs(hi->GetBinContent(j) - mid->GetBinContent(j) );
- dy_lo = fabs(lo->GetBinContent(j) - mid->GetBinContent(j));
-
- draw_box_at_point(x, dx, y, dy_hi, dy_lo, color);
- }
-}
-
-void PlotUtils::draw_errorbox(TGraph* mid,
- TGraph* hi,
- TGraph* lo, Int_t color, Double_t dx)
-{
- Double_t* x = mid->GetX();
- Double_t* y = mid->GetY();
- Double_t* y_hi= hi->GetY();
- Double_t* y_lo= lo->GetY();
-
- for (Int_t j=0; j < mid->GetN(); j++) {
-
- Double_t dy_hi = fabs(y_hi[j] - y[j]);
- Double_t dy_lo = fabs(y_lo[j] - y[j]);
-
- draw_box_at_point(x[j], dx, y[j], dy_hi, dy_lo, color);
- }
- return;
-}
-
-void PlotUtils::draw_errorbox(TGraphErrors* mid,
- TGraphErrors* hi,
- TGraphErrors* lo, Int_t color, Double_t dx)
-{
- Double_t* x = mid->GetX();
- Double_t* y = mid->GetY();
- Double_t* y_hi= hi->GetY();
- Double_t* y_lo= lo->GetY();
-
- for (Int_t j=0; j < mid->GetN(); j++) {
-
- Double_t dy_hi = fabs(y_hi[j] - y[j]);
- Double_t dy_lo = fabs(y_lo[j] - y[j]);
-
- draw_box_at_point(x[j], dx, y[j], dy_hi, dy_lo, color);
- }
- return;
-}
-
-void PlotUtils::draw_errorbox(TGraphAsymmErrors* mid,
- TGraphAsymmErrors* hi,
- TGraphAsymmErrors* lo, Int_t color, Double_t dx)
-{
- Double_t* x = mid->GetX();
- Double_t* y = mid->GetY();
- Double_t* y_hi= hi->GetY();
- Double_t* y_lo= lo->GetY();
-
- for (Int_t j=0; j < mid->GetN(); j++) {
-
- Double_t dy_hi = fabs(y_hi[j] - y[j]);
- Double_t dy_lo = fabs(y_lo[j] - y[j]);
-
- draw_box_at_point(x[j], dx, y[j], dy_hi, dy_lo, color);
- }
-}
-
-
-void PlotUtils::draw_box_at_point(Double_t x, Double_t dx, Double_t y,
- Double_t dy_plus, Double_t dy_minus,
- Int_t color)
-{
- TLine tl;
- tl.SetLineWidth(2);
- tl.SetLineColor(color);
-
- tl.DrawLine(x-dx, y-dy_minus, x+dx, y-dy_minus);
- tl.DrawLine(x-dx, y+dy_plus, x+dx, y+dy_plus);
- tl.DrawLine(x-dx, y+dy_plus, x-dx, y-dy_minus);
- tl.DrawLine(x+dx, y+dy_plus, x+dx, y-dy_minus);
-}
-
-// This function modifies the graphs so that 0, 1, 2 truly have the
-// mid, high, and low points at each x value.
-void PlotUtils::set_012(TH1F* h0,
- TH1F* h1,
- TH1F* h2)
-{
- for (Int_t j=0; j<h0->GetNbinsX(); j++) {
- Double_t y[3], dy[3];
- Int_t kmid=0, kmax=1, kmin=2;
-
- TH1F* h[3];
- h[0] = h0; // k=0: set to mid
- h[1] = h1; // k=1: set to high
- h[2] = h2; // k=2: set to low
-
- for (Int_t k=0; k<3; k++) {
- y[k] = h[k]->GetBinContent(j);
- dy[k] = h[k]->GetBinError(j);
- }
-
- kmid = midpos(y[0], y[1], y[2]);
- kmax = maxpos(y[0], y[1], y[2]);
- kmin = minpos(y[0], y[1], y[2]);
-
- h0->SetBinContent(j, y[kmid]);
- h1->SetBinContent(j, y[kmax]);
- h2->SetBinContent(j, y[kmin]);
-
- h0->SetBinError(j, dy[kmid]);
- h1->SetBinError(j, dy[kmax]);
- h2->SetBinError(j, dy[kmin]);
- } // j loop
-}
-
-// This function modifies the graphs so that 0, 1, 2 truly have the
-// mid, high, and low points at each x value.
-void PlotUtils::set_012(TGraphErrors* gr0,
- TGraphErrors* gr1,
- TGraphErrors* gr2)
-{
- // Assuming that gr0,1,2 have matching
- // x-values at each poInt_t j.
- for (Int_t j=0; j<gr0->GetN(); j++) {
-
- Double_t x, y[3], dy[3];
- Int_t kmid=0, kmax=1, kmin=2;
-
- TGraphErrors* g[3];
- g[0] = gr0; // k=0: set to mid
- g[1] = gr1; // k=1: set to high
- g[2] = gr2; // k=2: set to low
-
- for (Int_t k=0; k<3; k++) {
- Double_t ytmp=0;
- g[k]->GetPoint(j, x, ytmp);
- y[k] = ytmp;
- dy[k] = g[k]->GetErrorY(j);
- }
-
- kmid = midpos(y[0], y[1], y[2]);
- kmax = maxpos(y[0], y[1], y[2]);
- kmin = minpos(y[0], y[1], y[2]);
-
- gr0->SetPoint(j, x, y[kmid]);
- gr1->SetPoint(j, x, y[kmax]);
- gr2->SetPoint(j, x, y[kmin]);
-
- gr0->SetPointError(j, 0, dy[kmid]);
- gr1->SetPointError(j, 0, dy[kmax]);
- gr2->SetPointError(j, 0, dy[kmin]);
- } // j loop
-}
-
-// This function modifies the graphs so that 0, 1, 2 truly have the
-// mid, high, and low points at each x value.
-void PlotUtils::set_012(TGraphAsymmErrors* gr0,
- TGraphAsymmErrors* gr1,
- TGraphAsymmErrors* gr2)
-{
- // Assuming that gr0,1,2 have matching
- // x-values at each poInt_t j.
- for (Int_t j=0; j<gr0->GetN(); j++) {
-
- Double_t x, y[3], yhi[3], ylo[3];
- Int_t kmid=0, kmax=1, kmin=2;
-
- TGraphAsymmErrors* g[3];
- g[0] = gr0; // k=0: set to mid
- g[1] = gr1; // k=1: set to high
- g[2] = gr2; // k=2: set to low
-
- for (Int_t k=0; k<3; k++) {
- Double_t ytmp=0;
- g[k]->GetPoint(j, x, ytmp);
- y[k] = ytmp;
- yhi[k] = g[k]->GetErrorYhigh(j);
- ylo[k] = g[k]->GetErrorYlow(j);
- }
-
- kmid = midpos(y[0], y[1], y[2]);
- kmax = maxpos(y[0], y[1], y[2]);
- kmin = minpos(y[0], y[1], y[2]);
-
- gr0->SetPoint(j, x, y[kmid]);
- gr1->SetPoint(j, x, y[kmax]);
- gr2->SetPoint(j, x, y[kmin]);
-
- gr0->SetPointError(j, 0, 0, ylo[kmid], yhi[kmid]);
- gr1->SetPointError(j, 0, 0, ylo[kmax], yhi[kmax]);
- gr2->SetPointError(j, 0, 0, ylo[kmin], yhi[kmin]);
- } // j loop
-}
-
-Int_t PlotUtils::maxpos(Double_t x, Double_t y, Double_t z)
-{
- Double_t max = maximum(x,y,z);
- if(x==max) return 0;
- if(y==max) return 1;
- if(z==max) return 2;
-
- return -1;
-}
-
-Int_t PlotUtils::minpos(Double_t x, Double_t y, Double_t z)
-{
- Double_t min = minimum(x,y,z);
- if(x==min) return 0;
- if(y==min) return 1;
- if(z==min) return 2;
-
- return -1;
-}
-
-Int_t PlotUtils::midpos(Double_t x, Double_t y, Double_t z)
-{
- Int_t maxp = maxpos(x,y,z);
- Int_t minp = minpos(x,y,z);
-
- if (maxp==-1 || minp==-1){
- cout << "PlotUtils::midpos(): Bad parameters!"
- << endl;
- return -1;
- }
-
- if (maxp==minp) return 0;
- if ((maxp==1 && minp==2) || (maxp==2 && minp==1)) return 0;
- if ((maxp==0 && minp==2) || (maxp==2 && minp==0)) return 1;
- if ((maxp==0 && minp==1) || (maxp==1 && minp==0)) return 2;
-
- cout << "PlotUtils::midpos(): Did not return a middle position!"
- << endl;
- return -1;
-}
-
-Double_t PlotUtils::maximum(Double_t x, Double_t y, Double_t z)
-{
- Double_t max = (x > y ? x : y);
- max = (max > z ? max : z);
- return max;
-}
-
-Double_t PlotUtils::minimum(Double_t x, Double_t y, Double_t z)
-{
- Double_t min = (x < y ? x : y);
- min = (min < z ? min : z);
- return min;
-}
-
-void PlotUtils::offset_x(TGraphErrors* g, Double_t xoff)
-{
- Int_t npoints = g->GetN();
- Double_t* x = g->GetX();
- Double_t* y = g->GetY();
- for (Int_t j=0; j<npoints; j++){
- g->SetPoint(j, x[j] + xoff, y[j]);
- }
-}
-
-void PlotUtils::offset_x(TGraphAsymmErrors* g, Double_t xoff)
-{
- Int_t npoints = g->GetN();
- Double_t* x = g->GetX();
- Double_t* y = g->GetY();
- for (Int_t j=0; j<npoints; j++){
- g->SetPoint(j, x[j] + xoff, y[j]);
- }
-}
-
-void PlotUtils::padsetup(TCanvas* c, Int_t nx, Int_t ny, std::string opt,
- Double_t lmarg, Double_t rmarg,
- Double_t himarg, Double_t lomarg)
-{
- c->Divide(nx,ny,0,0);
-
- Int_t pad[100][100];
-
- for (Int_t x=1; x<=nx; x++) {
- for (Int_t y=1; y<=ny; y++) {
- pad[x][y] = nx * (y-1) + x;
- }
- }
-
- if (opt=="mergex" || opt=="MERGEX") {
- for (Int_t x=1; x<=nx; x++) {
- for (Int_t y=1; y<=ny; y++) {
- // left edge...
- if (x == 1)
- c->GetPad(pad[x][y])->SetLeftMargin(lmarg);
- else
- c->GetPad(pad[x][y])->SetLeftMargin(0);
-
- // right edge...
- if (x==nx)
- c->GetPad(pad[x][y])->SetRightMargin(rmarg);
- else
- c->GetPad(pad[x][y])->SetRightMargin(0);
- // top edge...
- if (y==1) c->GetPad(pad[x][y])->SetTopMargin(himarg);
- // bottom...
- if (y==ny) c->GetPad(pad[x][y])->SetBottomMargin(lomarg);
- }
- }
- }
- else {
- for (Int_t x=1; x<=nx; x++) {
- for (Int_t y=1; y<=ny; y++) {
- // left edge...
- c->GetPad(pad[x][y])->SetLeftMargin(lmarg);
- // right edge...
- c->GetPad(pad[x][y])->SetRightMargin(rmarg);
- // top edge...
- if (y==1) c->GetPad(pad[x][y])->SetTopMargin(himarg);
- // bottom...
- if (y==ny) c->GetPad(pad[x][y])->SetBottomMargin(lomarg);
- }
- }
- }
- for (Int_t i=1; i<=nx*ny; i++){
- c->GetPad(i)->SetTickx();
- c->GetPad(i)->SetTicky();
- // c->GetPad(i)->SetGridy();
- c->GetPad(i)->SetFrameLineWidth(2);
- }
-}
-
-void PlotUtils::set_ylimits(TObjArray* arr1, TObjArray* arr2,
- Double_t topspace, Double_t lowspace)
-{
- Int_t n1 = arr1->GetEntries();
- Int_t n2 = arr2->GetEntries();
- if (n1 != n2) {
- cout << __FILE__ << " " << __LINE__ <<": " << endl;
- cout << "Error: Arrays have unequal length: "
- << n1 << " vs " << n2 << "." << endl;
- return;
- }
-
- for (Int_t i=0; i<n1; i++) {
- TH1F* h1 = (TH1F*)arr1->At(i);
- TH1F* h2 = (TH1F*)arr2->At(i);
- set_ylimits(h1, h2, topspace, lowspace);
- }
-}
-
-
-void PlotUtils::set_ylimits(TH1F* h1, TH1F* h2,
- Double_t topspace, Double_t lowspace)
-{
- Int_t maxbin1 = h1->GetMaximumBin();
- Int_t maxbin2 = h2->GetMaximumBin();
- Int_t minbin1 = h1->GetMinimumBin();
- Int_t minbin2 = h2->GetMinimumBin();
-
- Double_t max1 = h1->GetBinContent(maxbin1) + h1->GetBinError(maxbin1);
- Double_t max2 = h2->GetBinContent(maxbin2) + h2->GetBinError(maxbin2);
-
- Double_t min1 = h1->GetBinContent(minbin1) - h1->GetBinError(minbin1);
- Double_t min2 = h2->GetBinContent(minbin2) - h2->GetBinError(minbin2);
-
- Double_t max = max1 > max2 ? max1 : max2;
- Double_t min = min1 < min2 ? min1 : min2;
-
- if (topspace == 0) topspace = 0.1*max;
- if (lowspace == 0) lowspace = 0.1*max;
-
- Double_t upperSpace = (max - min)*topspace;
- Double_t lowerSpace = (max - min)*lowspace;
-
- h1->GetYaxis()->SetRangeUser(min - lowerSpace, max + upperSpace);
- h2->GetYaxis()->SetRangeUser(min - lowerSpace, max + upperSpace);
-}
-
-void PlotUtils::set_ylimits(TF1* f1, TF1* f2,
- Double_t topspace, Double_t lowspace)
-{
- Double_t xlo1 = f1->GetXaxis()->GetXmin();
- Double_t xhi1 = f1->GetXaxis()->GetXmax();
- Double_t xlo2 = f2->GetXaxis()->GetXmin();
- Double_t xhi2 = f2->GetXaxis()->GetXmax();
-
- Double_t min1 = f1->GetMinimum(xlo1, xhi1);
- Double_t max1 = f1->GetMaximum(xlo1, xhi1);
-
- Double_t min2 = f2->GetMinimum(xlo2, xhi2);
- Double_t max2 = f2->GetMaximum(xlo2, xhi2);
-
- Double_t max = max1 > max2 ? max1 : max2;
- Double_t min = min1 < min2 ? min1 : min2;
-
- Double_t upperSpace = (max - min)*topspace;
- Double_t lowerSpace = (max - min)*lowspace;
-
- f1->GetYaxis()->SetRangeUser(min - lowerSpace, max + upperSpace);
- f2->GetYaxis()->SetRangeUser(min - lowerSpace, max + upperSpace);
-}
-
-void PlotUtils::multiplot(TCanvas* c, TObjArray *hArray,
- Int_t nx, Int_t ny, TString opt, Int_t iStart)
-{
- Int_t ipad = 1;
-
- c->cd();
- TList* prims = c->GetListOfPrimitives();
- if (prims->GetEntries() > 0) {
- opt.Append("same");
- }
- else padsetup(c, nx, ny);
-
- for (Int_t i = iStart; i < iStart+nx*ny; i++) {
- c->cd(ipad);
-
- TObject* obj = hArray->At(i);
-
- if (!obj)
- cout << __FILE__ << " " << __LINE__
- << ": Object in hArray is null." << endl;
- else if (obj->InheritsFrom("TH1F"))
- ((TH1F*)hArray->At(i))->Draw(opt.Data());
- else if (obj->InheritsFrom("TH1D"))
- ((TH1D*)hArray->At(i))->Draw(opt.Data());
- else if (obj->InheritsFrom("TF1"))
- ((TF1*)hArray->At(i))->Draw(opt.Data());
- else
- (hArray->At(i))->Draw(opt.Data());
-
- ipad++;
- }
-}
-
-void PlotUtils::scale_axis_labels(TObjArray* a, Double_t c)
-{
- for (Int_t i=0; i<a->GetEntries(); i++) {
- TObject* obj = a->At(i);
- if (obj->InheritsFrom("TH1F")) {
- TH1F* h = (TH1F*)obj;
- scale_axis_labels(h, c);
- }
- }
-}
-
-void PlotUtils::scale_axis_labels(TH1F* h, Double_t c)
-{
- TAxis* x = h->GetXaxis();
- TAxis* y = h->GetYaxis();
-
- x->SetLabelSize(c*x->GetLabelSize());
- y->SetLabelSize(c*y->GetLabelSize());
-
- x->SetLabelOffset(c*x->GetLabelOffset());
- y->SetLabelOffset(c*y->GetLabelOffset());
-
- x->SetTitleSize(c*x->GetTitleSize());
- y->SetTitleSize(c*y->GetTitleSize());
- return;
-}
-
-void PlotUtils::make_nice_axes(TCanvas* can, TH1F* h, Double_t c)
-{
- TAxis* x = h->GetXaxis();
- TAxis* y = h->GetYaxis();
-
- x->SetLabelSize(c*x->GetLabelSize());
- y->SetLabelSize(c*y->GetLabelSize());
-
- x->SetLabelOffset(c*x->GetLabelOffset());
- y->SetLabelOffset(c*y->GetLabelOffset());
-
- x->SetTitleSize(c*x->GetTitleSize());
- y->SetTitleSize(c*y->GetTitleSize());
-
- // new below
- x->SetTitleOffset(0.9*c*x->GetTitleOffset());
- y->SetTitleOffset(c*y->GetTitleOffset());
-
- // default left and bottom margins: 0.10
- // assume defaults, and add extra space
- can->SetLeftMargin(c*0.12);
- can->SetBottomMargin(c*0.12);
-
- x->CenterTitle();
- y->CenterTitle();
-}
-
-// I should really learn how to write templates...
-void PlotUtils::make_nice_axes(TCanvas* can, TH2F* h, Double_t c)
-{
- TAxis* x = h->GetXaxis();
- TAxis* y = h->GetYaxis();
-
- x->SetLabelSize(c*x->GetLabelSize());
- y->SetLabelSize(c*y->GetLabelSize());
-
- x->SetLabelOffset(c*x->GetLabelOffset());
- y->SetLabelOffset(c*y->GetLabelOffset());
-
- x->SetTitleSize(c*x->GetTitleSize());
- y->SetTitleSize(c*y->GetTitleSize());
-
- // new below
- x->SetTitleOffset(c*x->GetTitleOffset());
- y->SetTitleOffset(c*y->GetTitleOffset());
-
- // default left and bottom margins: 0.10
- // assume defaults, and add extra space
- can->SetLeftMargin(c*0.12);
- can->SetBottomMargin(c*0.12);
-
- x->CenterTitle();
- y->CenterTitle();
-}
+++ /dev/null
-// $Id$
-
-#ifndef Utils_h
-#define Utils_h
-
-#include <string>
-#include <vector>
-#include <TH1.h>
-#include <TH2.h>
-#include <TF1.h>
-#include <TCanvas.h>
-#include <TLegend.h>
-#include <TGraph.h>
-#include <TGraphErrors.h>
-#include <TGraphAsymmErrors.h>
-#include <TString.h>
-#include <TObjArray.h>
-
-class PlotUtils : public TObject
-{
-public:
- static void set_hist_props(TObjArray* arr,
- Int_t linecolor,
- Int_t fillcolor,
- Int_t markercolor,
- Int_t markerstyle,
- Double_t markersize = 0.75);
-
- static void set_hist_props(TH1F* h,
- Int_t linecolor,
- Int_t fillcolor,
- Int_t markercolor,
- Int_t markerstyle,
- Double_t markersize = 0.75);
-
- static void set_tgraph_props(TGraphErrors* g,
- Int_t linecolor,
- Int_t markercolor,
- Int_t markerstyle,
- Double_t markersize = 0.75);
-
- static void set_tgraph_props(TGraphAsymmErrors* g,
- Int_t linecolor,
- Int_t markercolor,
- Int_t markerstyle,
- Double_t markersize = 0.75);
-
- static void shade(TF1* f_hi, TF1* f_lo);
- static void draw_errorbox(TH1F* mid,
- TH1F* hi,
- TH1F* lo,
- Int_t color = 46, Double_t dx = 0.06);
- static void draw_errorbox(TGraph* mid,
- TGraph* hi,
- TGraph* lo,
- Int_t color = 46, Double_t dx = 0.06);
- static void draw_errorbox(TGraphErrors* mid,
- TGraphErrors* hi,
- TGraphErrors* lo,
- Int_t color = 46, Double_t dx = 0.06);
- static void draw_errorbox(TGraphAsymmErrors* mid,
- TGraphAsymmErrors* hi,
- TGraphAsymmErrors* lo,
- Int_t color = 46, Double_t dx = 0.06);
-
- static void draw_box_at_point(Double_t x, Double_t dx, Double_t y,
- Double_t dy_plus, Double_t dy_minus,
- Int_t color = 46);
-
- static Double_t maximum(Double_t, Double_t, Double_t);
- static Double_t minimum(Double_t, Double_t, Double_t);
-
- // Return position of mid, max, or mid argument (0, 1, or 2)
- // Return -1 if something is wrong.
- static Int_t midpos(Double_t x, Double_t y, Double_t z);
- static Int_t maxpos(Double_t x, Double_t y, Double_t z);
- static Int_t minpos(Double_t x, Double_t y, Double_t z);
-
- static void offset_x(TGraphErrors* g, Double_t xoff);
- static void offset_x(TGraphAsymmErrors* g, Double_t xoff);
-
- static void padsetup(TCanvas *c, Int_t nx=4, Int_t ny=3, std::string opt = "",
- Double_t lmarg = 0.12, Double_t rmarg = 0.01,
- Double_t topmarg = 0.01, Double_t botmarg = 0.12);
-
- static void multiplot(TCanvas* c, TObjArray* hArray,
- Int_t nx, Int_t ny, TString opt = "", Int_t iStart=0);
-
- static void set_ylimits(TObjArray* a1, TObjArray* a2,
- Double_t hispace=0, Double_t lospace=0);
-
- static void set_ylimits(TH1F* h1, TH1F* h2,
- Double_t hispace=0, Double_t lospace=0);
-
- static void set_ylimits(TF1* f1, TF1* f2,
- Double_t topspace = 0.1, Double_t lowspace = 0.1);
-
- static void set_012(TH1F* h0, TH1F* h1, TH1F* h2);
-
- static void set_012(TGraphErrors* gr0,
- TGraphErrors* gr1,
- TGraphErrors* gr2);
-
- static void set_012(TGraphAsymmErrors* gr0,
- TGraphAsymmErrors* gr1,
- TGraphAsymmErrors* gr2);
-
- static void scale_axis_labels(TH1F* h, Double_t c);
- static void scale_axis_labels(TObjArray* a, Double_t c);
-
- static void make_nice_axes(TCanvas* can, TH1F* h, Double_t c);
- static void make_nice_axes(TCanvas* can, TH2F* h, Double_t c);
-
- private:
- // None...allows use as a static class
-
- ClassDef(PlotUtils,0) // Plot utilities class
-};
-#endif
+++ /dev/null
-// $Id$
-
-#if !defined(__CINT__) || defined(__MAKECINT__)
-#include <Riostream.h>
-#include <TFile.h>
-#include <TTree.h>
-#include <TMath.h>
-#include <TRandom3.h>
-#include <TChain.h>
-#include <TH2F.h>
-#include <TTimeStamp.h>
-#include "TreeClasses.h"
-#include "EventPoolManager.h"
-#endif
-
-void anaCorr(const char *inFileNames,
- Int_t nEvents=-1,
- Double_t vzmin=-10,
- Double_t vzmax=+10,
- Int_t vzbins=40,
- Int_t vcmin=10,
- Int_t vcmax=500,
- Int_t psize=5,
- Double_t ptmin=0.4,
- Double_t ptmax=10,
- Double_t etamin=-1.4,
- Double_t etamax=+1.4);
-
-inline Double_t DeltaPhi(const MyPart &t1, const MyPart &t2,
- Double_t rangeMin = -TMath::Pi()/2,
- Double_t rangeMax = 3*TMath::Pi()/2);
-inline Double_t DeltaEta(const MyPart &t1, const MyPart &t2);
-inline Bool_t InBounds(Double_t val, Double_t min, Double_t max);
-inline Double_t InvMass(const MyPart &p1, const MyPart &p2);
-
-//-----------------------------------------------------------------------------------------------------
-
-void anaCorr(const char *inFileNames,
- Int_t nEvents,
- Double_t vzmin,
- Double_t vzmax,
- Int_t vzbins,
- Int_t vcmin,
- Int_t vcmax,
- Int_t psize,
- Double_t ptmin,
- Double_t ptmax,
- Double_t etamin,
- Double_t etamax)
-{
- Bool_t doExclCut = 0;
- Double_t dpmax = 0.03;
- Double_t demax = 0.01;
-
- Noti *nme = new Noti;
- TChain *tree = new TChain("MyTree");
- tree->Add(inFileNames);
- tree->SetNotify(nme);
-
- Int_t nents = tree->GetEntries();
- cout << "Found " << nents << " entries!" << endl;
- if (nents<=0)
- return;
- if (nEvents>0&&nEvents<nents) {
- nents=nEvents;
- cout << "Using " << nents << " entries!" << endl;
- }
-
- if (TClass::GetClass("MyPart"))
- TClass::GetClass("MyPart")->IgnoreTObjectStreamer();
- if (TClass::GetClass("MyTracklet"))
- TClass::GetClass("MyTracklet")->IgnoreTObjectStreamer();
-
- MyHeader *header = 0;
- TClonesArray *tracks = 0;
- TBranch *hBranch = 0;
- TBranch *tBranch = 0;
-
- Int_t multbinsa[] = {1,2,3,5,10,20,30,40,50,60,70,80,90,100,120,150,200,500};
- Int_t multbins=17;
- Double_t *multbinsa2 = new Double_t[multbins];
- for(Int_t i=0;i<multbins;++i) {
- multbinsa2[i] = multbinsa[i];
- }
-
- TH1F* hMultBins = new TH1F("hMultBins", "Event multiplicity binning",
- multbins-1, multbinsa2);
- Double_t *vzbinsa = new Double_t[vzbins+1];
- Double_t vzbinwidth = Double_t(vzmax-vzmin)/vzbins;
- for(Int_t i=0;i<vzbins+1;++i) {
- vzbinsa[i] = vzmin+i*vzbinwidth;
- }
- TH1F* hZvtxBins = new TH1F("hZvtxBins", "Event Z-vertex binning",
- vzbins, vzbinsa);
-
- EventPoolManager *pm = new EventPoolManager;
- pm->InitEventPools(psize, multbins, multbinsa, vzbins, vzbinsa);
-
-
- TH2F **hSig = new TH2F*[multbins];
- TH2F **hBkg = new TH2F*[multbins];
- TH1F **hMul = new TH1F*[multbins];
- for (Int_t i = 0; i<multbins; ++i) {
- hSig[i] = new TH2F(Form("hSig%d",i),";#Delta#eta;#Delta#phi",60,-3,3,64,-TMath::Pi()/2,3*TMath::Pi()/2);
- hBkg[i] = new TH2F(Form("hBgk%d",i),"#;Delta#eta;#Delta#phi",60,-3,3,64,-TMath::Pi()/2,3*TMath::Pi()/2);
- hMul[i] = new TH1F(Form("hMul%d",i),"#;sel tracks",500,0,500);
- }
-
- for (Int_t i=0;i<nents;++i) {
- Int_t li = tree->LoadTree(i);
- if (nme->Notified()) {
- hBranch = tree->GetBranch("header");
- hBranch->SetAddress(&header);
- tBranch = tree->GetBranch("parts");
- tBranch->SetAddress(&tracks);
- nme->Reset();
- }
-
- hBranch->GetEntry(li);
- if (header->fIsPileupSPD)
- continue;
- if ((header->fVz<vzmin) ||
- (header->fVz>vzmax) )
- continue;
- if ((header->fVc<vcmin) ||
- (header->fVc>vcmax) )
- continue;
-
- tBranch->GetEntry(li);
- Int_t ntracks = tracks->GetEntries();
- TClonesArray *seltracks = new TClonesArray("MyPart",ntracks);
- Int_t nseltracks = 0;
- for (Int_t t=0;t<ntracks;++t) {
- MyPart *part = (MyPart*)tracks->At(t);
- if (part->IsITSRefit() && !part->IsTPCRefit())
- continue;
- if (!InBounds(part->fPt,ptmin,ptmax))
- continue;
- if (!InBounds(part->fEta,etamin,etamax))
- continue;
- if (TMath::Abs(part->fD)>0.3)
- continue;
- if (TMath::Abs(part->fZ)>0.3)
- continue;
- new((*seltracks)[nseltracks]) MyPart(*part);
- ++nseltracks;
- }
-
- Int_t mBin = hMultBins->FindBin(nseltracks) - 1;
- Int_t zBin = hZvtxBins->FindBin(header->fVz) - 1;
- GenericEventPool *pool = pm->GetEventPool(mBin, zBin);
- if (!pool) {
- continue;
- }
-
- if (pool->IsReady()) {
- if (i%100==0)
- cout << "Working on event " << i << " from " << header->fRun << endl;
-
- hMul[mBin]->Fill(nseltracks);
-
- Int_t sigpairs = 0; // count pairs
- for (Int_t t1=0;t1<nseltracks;++t1) {
- MyPart *part1 = (MyPart*)seltracks->At(t1);
- for (Int_t t2=t1+1;t2<nseltracks;++t2) {
- MyPart *part2 = (MyPart*)seltracks->At(t2);
- Double_t dphi = DeltaPhi(*part1,*part2);
- Double_t deta = DeltaEta(*part1,*part2);
- Double_t dr = 1e10;
- if (doExclCut)
- dr = dphi*dphi/(dpmax*dpmax) + deta*deta/(demax*demax);
- if(dr>1) {
- ++sigpairs;
- }
- }
- }
- Double_t weight = 1./sigpairs;
- for (Int_t t1=0;t1<nseltracks;++t1) {
- MyPart *part1 = (MyPart*)seltracks->At(t1);
- for (Int_t t2=t1+1;t2<nseltracks;++t2) {
- MyPart *part2 = (MyPart*)seltracks->At(t2);
- Double_t dphi = DeltaPhi(*part1,*part2);
- Double_t deta = DeltaEta(*part1,*part2);
- Double_t dr = 1e10;
- if (doExclCut)
- dr = dphi*dphi/(dpmax*dpmax) + deta*deta/(demax*demax);
- if(dr>1) {
- ++sigpairs;
- hSig[mBin]->Fill(deta,dphi,weight);
- }
- }
- }
-
- TObjArray *mseltracks = pool->GetRandomEvent();
- Int_t nmseltracks = mseltracks->GetEntries();
- Int_t mixpairs = 0; // count mix pairs
- for (Int_t t1=0;t1<nseltracks;++t1) {
- MyPart *part1 = (MyPart*)seltracks->At(t1);
- for (Int_t t2=1;t2<nmseltracks;++t2) {
- MyPart *part2 = (MyPart*)mseltracks->At(t2);
- Double_t dphi = DeltaPhi(*part1,*part2);
- Double_t deta = DeltaEta(*part1,*part2);
- Double_t dr = 1e10;
- if (doExclCut)
- dr = dphi*dphi/(dpmax*dpmax) + deta*deta/(demax*demax);
- if(dr>1) {
- ++mixpairs;
- }
- }
- }
- Double_t weight2 = 1./mixpairs;
- for (Int_t t1=0;t1<nseltracks;++t1) {
- MyPart *part1 = (MyPart*)seltracks->At(t1);
- for (Int_t t2=1;t2<nmseltracks;++t2) {
- MyPart *part2 = (MyPart*)mseltracks->At(t2);
- Double_t dphi = DeltaPhi(*part1,*part2);
- Double_t deta = DeltaEta(*part1,*part2);
- Double_t dr = 1e10;
- if (doExclCut)
- dr = dphi*dphi/(dpmax*dpmax) + deta*deta/(demax*demax);
- if(dr>1) {
- ++sigpairs;
- hBkg[mBin]->Fill(deta,dphi,weight2);
- }
- }
- }
- }
- pool->UpdatePool(i,header,seltracks);
- if (!pool->IsReady())
- pool->PrintInfo();
- }
-
- TTimeStamp t;
- TString fname(Form("histout-%d.root",(Int_t)t.GetSec()));
- TFile outfile(fname,"recreate");
- for (Int_t i = 0; i<multbins; ++i) {
- hSig[i]->Write();
- hBkg[i]->Write();
- hMul[i]->Write();
- TH2F *hMix=(TH2F*)hSig[i]->Clone(Form("hMix%d",i));
- hMix->Divide(hSig[i],hBkg[i],1.,1.);
- hMix->Write();
- TH2F *hMix2=(TH2F*)hSig[i]->Clone(Form("hMix2%d",i));
- for(Int_t x=1;x<=hMix->GetNbinsX();++x) {
- for(Int_t y=1;y<=hMix->GetNbinsY();++y) {
- Int_t bin = hMix->GetBin(x,y);
- Double_t val = hMix->GetBinContent(bin);
- val = hMul[i]->GetMean()*(val-1);
- hMix2->SetBinContent(bin,val);
- }
- }
- hMix2->Write();
- }
-}
-
-//-----------------------------------------------------------------------------------------------------
-
-Bool_t InBounds(Double_t val, Double_t min, Double_t max)
-{
- if (val<min)
- return 0;
- if (val>max)
- return 0;
- return 1;
-}
-
-Double_t DeltaPhi(const MyPart &t1, const MyPart &t2,
- Double_t rangeMin, Double_t rangeMax)
-{
- Double_t dphi = -999;
- Double_t pi = TMath::Pi();
- Double_t phia = t1.fPhi;
- Double_t phib = t2.fPhi;
-
- if (phia < 0) phia += 2*pi;
- else if (phia > 2*pi) phia -= 2*pi;
- if (phib < 0) phib += 2*pi;
- else if (phib > 2*pi) phib -= 2*pi;
- dphi = phib - phia;
- if (dphi < rangeMin) dphi += 2*pi;
- else if (dphi > rangeMax) dphi -= 2*pi;
-
- return dphi;
-}
-
-Double_t DeltaEta(const MyPart &t1, const MyPart &t2)
-{
- return t1.fEta - t2.fEta;
-}
-
-Double_t InvMass(const MyPart &p1, const MyPart &p2)
-{
- Double_t px1 = p1.Px();
- Double_t py1 = p1.Py();
- Double_t pz1 = p1.Pz();
- Double_t px2 = p2.Px();
- Double_t py2 = p2.Py();
- Double_t pz2 = p2.Pz();
-
- Double_t pm1 = TMath::Sqrt(px1*px1+py1*py1+pz1*pz1);
- Double_t pm2 = TMath::Sqrt(px2*px2+py2*py2+pz1*pz2);
- Double_t p12 = px1*px2+py1*py2+pz1*pz2;
-
- Double_t m = TMath::Sqrt(pm1*pm2-p12);
- return m;
-}
-
+++ /dev/null
-void batch()
-{
- ifstream ifs("goodruns.txt");
- int run = 0;
- const char* dir = "../rootfiles/res_LHC10e_09122010/mergedruns";
- char* inFileName = "";
- char* outFileName = "";
-
- while (ifs >> run){
- inFileName = Form("%s/merged_run%i.root", dir, run);
-
- outFileName = Form("output/histos_%i.root", run);
-
- char* cmd = Form("root -b -q 'runAutoCorr.C+g(\"%s\", \"%s\")'",
- inFileName, outFileName);
-
- cout << cmd << endl;
-
- gSystem->Exec(cmd);
-
- }
- return;
-}
+++ /dev/null
-// sum.root made by doing
-// cd output/
-// hadd sum.root histos_*.root
-
-void draw(const char* histFile = "output/sum.root")
-{
- TFile* inFile = new TFile(histFile, "read");
- int nMultBinsRead = hMultBins->GetNbinsX();
- const int nMultBins = nMultBinsRead;
- int nPtBinsRead = hPtBins->GetNbinsX();
- const int nPtBins = nPtBinsRead;
-
- TH2F* hSig[nMultBins];
- TH2F* hBkg[nMultBins];
- TH2F* hSB[nMultBins];
-
- TH2F* hSigPt[nMultBins][nPtBins];
- TH2F* hBkgPt[nMultBins][nPtBins];
- TH2F* hSBPt[nMultBins][nPtBins];
-
-
- for (int iM=0; iM<nMultBins; iM++) {
- hSig[iM] = (TH2F*)inFile->Get(Form("hSig_%i", iM));
- hBkg[iM] = (TH2F*)inFile->Get(Form("hBkg_%i", iM));
- hSB[iM] = (TH2F*)inFile->Get(Form("hSB_%i", iM));
-
- for (int iPt=0; iPt<nPtBins; iPt++) {
- hSigPt[iM][iPt] = (TH2F*)inFile->Get(Form("hSigPt_%i_%i",iM,iPt));
- hBkgPt[iM][iPt] = (TH2F*)inFile->Get(Form("hBkgPt_%i_%i",iM,iPt));
- hSBPt[iM][iPt] = (TH2F*)inFile->Get(Form("hSBPt_%i_%i", iM,iPt));
- }
-
- }
-
- TCanvas* c1 = new TCanvas("c1", "c1", 1);
- c1->cd();
- hDEta->Draw();
- // hEta->Draw("same");
-
- TCanvas* c2 = new TCanvas("c2", "c2", 1);
- c2->cd();
- hDPhi->Draw();
-
- TCanvas* c3 = new TCanvas("c3", "c3", 1);
- c3->cd();
- gPad->SetLogy();
- hPtAll->Draw("ep");
- hPtSel->SetLineColor(kRed);
- hPtSel->Draw("epsame");
-
- TCanvas* c4 = new TCanvas("c4", "c4", 1);
- c4->Divide(2, 1, 0.001, 0.001);
- c4->cd(1);
- hDR->Draw("surf2");
- c4->cd(2);
- hDRcut->Draw("colz");
-
- TCanvas* cEff = new TCanvas("cEff", "cEff", 1);
- cEff->cd();
- hMixEff->Draw("colz");
-
- TCanvas* cSig = new TCanvas("cSig", "cSig", 1);
- cSig->Divide(nMultBins, 1, 0.001, 0.001);
- for (int iM=0; iM<nMultBins; iM++) {
- cSig->cd(iM+1);
- hSig[iM]->Draw("surf1");
- }
- TCanvas* cBkg = new TCanvas("cBkg", "cBkg", 1);
- cBkg->Divide(nMultBins, 1, 0.001, 0.001);
- for (int iM=0; iM<nMultBins; iM++) {
- cBkg->cd(iM+1);
- hBkg[iM]->Draw("surf1");
- }
-
- TCanvas* cSB = new TCanvas("cSB", "cSB", 1);
- cSB->Divide(nMultBins, 1, 0.001, 0.001);
- for (int iM=0; iM<nMultBins; iM++) {
- cSB->cd(iM+1);
- hSB[iM]->Draw("surf1");
- }
-
- TCanvas* cSBpt[nPtBins];
- for (int iPt=0; iPt<nPtBins; iPt++) {
- cSBpt[iPt]= new TCanvas(Form("cSBpt_%i",iPt), Form("cSBpt_%i",iPt), 1);
- cSBpt[iPt]->Divide(nMultBins, 1, 0.001, 0.001);
- for (int iM=0; iM<nMultBins; iM++) {
- cSBpt[iPt]->cd(iM+1);
- hSBPt[iM][iPt]->Draw("surf1");
- }
- }
-
- TCanvas* cMult = new TCanvas("cMult", "cMult", 1);
- // cMult->Divide(nMultBins, 1, 0.001, 0.001);
- cMult->cd();
- hMultAll->Draw();
- hMultSel->SetLineColor(kRed);
- hMultSel->Draw("same");
-
- return;
-}
+++ /dev/null
-130517
-130519
-130520
-130524
-130526
-130601
-130608
-130620
-130621
-130623
-130627
-130628
-130640
-130696
-130704
-130793
-130795
-130798
-130799
-130833
-130834
-130840
-130842
-130844
-130847
-130848
+++ /dev/null
-// $Id:$
-
-#if !defined(__CINT__) || defined(__MAKECINT__)
-#include <vector>
-#include <map>
-#include <Riostream.h>
-#include <TChain.h>
-#include <TClonesArray.h>
-#include <TError.h>
-#include <TFile.h>
-#include <TDirectory.h>
-#include <TMath.h>
-#include <TRandom3.h>
-#include <TSystem.h>
-#include <TObjectTable.h>
-#include "TreeClasses.C"
-#endif
-
-#define DUPLICATECHECK
-#ifdef DUPLICATECHECK
-#define DUPLICATEVERBOSE
-#endif
-
-struct EvInfo
-{
- Int_t fRun;
- Short_t fArrId;
- UInt_t fEntry;
-};
-
-
-void mergeCorr(Int_t run_from=-1,
- Int_t run_to=-1,
- const char *outPrefix = "merged",
- const char *inFileNames = "res/output_*.root");
-
-void mergeCorr(Int_t nEvents,
- const char *outFileName,
- const char *inFileNames);
-
-//-----------------------------------------------------------------------------------------------------
-
-void mergeCorr(Int_t run_from,
- Int_t run_to,
- const char *outPrefix,
- const char *inFileNames)
-{
- TChain *c = new TChain("MyTree");
- c->Add(inFileNames);
-
- std::map<Int_t,Int_t> goodruns;
-
- TObjArray *filenames = c->GetListOfFiles();
- Int_t nfiles = filenames->GetEntries();
- for(Int_t i=0; i<nfiles; ++i) {
- TString fname(filenames->At(i)->GetTitle());
- TFile *file = TFile::Open(fname);
- if (!file)
- continue;
- TList *output = dynamic_cast<TList*>(file->Get("output"));
- output->SetOwner(1);
- TTree *tree = dynamic_cast<TTree*>(output->FindObject("MyTree"));
- Int_t nent = tree->GetEntries();
- if (nent<1) {
- delete output;
- file->Close();
- delete file;
- continue;
- }
-
- MyHeader *header = 0;
- TBranch *br = tree->GetBranch("header");
- br->SetAddress(&header);
-
- std::vector<Int_t> multruns;
-
- Int_t runno = -1;
- for (Int_t ev=0;ev<nent;++ev) {
- br->GetEntry(ev);
- if (header->fRun==0) {
- continue;
- }
- if (runno<0)
- runno = header->fRun;
- else if (runno!=header->fRun) {
- Bool_t found = 0;
- for (UInt_t j=0;j<multruns.size();++j) {
- if (multruns.at(j)==header->fRun) {
- found = 1;
- break;
- }
- }
- if (!found)
- multruns.push_back(header->fRun);
- }
- }
- if (multruns.size()<=0) {
- if ((runno>run_from) && (run_to<0 || runno<run_to)) {
- cout << "Found run " << runno << " associated to " << fname << endl;
- goodruns.insert(pair<Int_t,Int_t>(runno,0));
- TString dir(gSystem->DirName(inFileNames));
- dir+="/runs";
- TString link(Form("%s/%s",gSystem->pwd(),fname.Data()));
- TString rundir(Form("%s/%d", dir.Data(), runno));
- TString cmd(Form("mkdir -p %s && cd %s && ln -sf %s",
- rundir.Data(), rundir.Data(), link.Data()));
- gSystem->Exec(cmd);
- }
- } else {
- cout << "Multiple runs (" << runno;
- for (UInt_t j=0;j<multruns.size();++j)
- cout << ", " << multruns.at(j);
- cout << ") found in " << fname << endl;
- cout << " Cannot deal with this case yet!!!" << endl;
- }
- delete output;
- file->Close();
- delete file;
- }
-
- for (map<Int_t, Int_t>::iterator it = goodruns.begin();it!=goodruns.end();++it) {
- TString infiles(Form("%s/runs/%d/*.root",gSystem->DirName(inFileNames),it->first));
- TString outdir(Form("%s/mergedruns",gSystem->DirName(inFileNames)));
- TString outFileName(Form("%s/%s_run%d.root",outdir.Data(),outPrefix,it->first));
- gSystem->Exec(Form("mkdir -p %s",outdir.Data()));
- gSystem->Exec(Form("root -b -q mergeCorr.C+'(-1,\"%s\",\"%s\")'",
- outFileName.Data(),infiles.Data()));
- }
-}
-
-//-----------------------------------------------------------------------------------------------------
-
-void mergeCorr(Int_t nEvents,
- const char *outFileName,
- const char *inFileNames)
-{
- TChain *c = new TChain("MyTree");
- c->Add(inFileNames);
-
- map<ULong64_t, EvInfo*> einfos;
-#ifdef DUPLICATEVERBOSE
- map<ULong64_t, MyHeader*> eheaders;
-#endif
-
- Int_t totnent = 0;
- Int_t totnsus = 0;
-#ifdef DUPLICATECHECK
- Int_t totncan = 0;
-#ifdef DUPLICATEVERBOSE
- Int_t totndup = 0;
-#endif
-#endif
-
- TObjArray *maps = new TObjArray;
- maps->SetOwner(1);
- Int_t lrun = 0;
- TObjArray *filenames = c->GetListOfFiles();
- Int_t nfiles = filenames->GetEntries();
- for(Int_t i=0; i<nfiles; ++i) {
- TString fname(filenames->At(i)->GetTitle());
- TFile *file = TFile::Open(fname);
- if (!file)
- continue;
- TList *output = dynamic_cast<TList*>(file->Get("output"));
- output->SetOwner(1);
- TTree *tree = dynamic_cast<TTree*>(output->FindObject("MyTree"));
- Int_t nent = tree->GetEntries();
- if (nent<1) {
- delete output;
- file->Close();
- delete file;
- continue;
- }
-
- totnent += nent;
- cout << "Found "<< nent << " entries in " << fname << endl;
- TObjArray *filesnames = dynamic_cast<TObjArray*>(output->FindObject("filesnames"));
-
- MyHeader *header = 0;
- TBranch *br = tree->GetBranch("header");
- br->SetAddress(&header);
-
- for (Int_t ev=0;ev<nent;++ev) {
- br->GetEntry(ev);
- if (header->fRun==0) {
-#ifdef DUPLICATEVERBOSE
- cout << "--> Suspicious found for " << header->fRun << " "
- << header->fTime << " " << header->fPeriod << " " << header->fBx << endl;
-#endif
- ++totnsus;
- continue;
- }
-
- //ULong64_t evtid = header->fTime;
- ULong64_t evtid = header->GetEventId();
-#ifdef DUPLICATECHECK
- if (einfos.find(evtid)!=einfos.end()) {
-#ifdef DUPLICATEVERBOSE
- cout << "--> Potential duplicate found for " << header->fRun << " "
- << header->fTime << " " << header->fPeriod << " " << header->fBx << endl;
-#endif
- ++totncan;
-#ifdef DUPLICATEVERBOSE
- map<ULong64_t, MyHeader*>::iterator it = eheaders.find(evtid);
- if ( (it->second->fRun==header->fRun) &&
- (it->second->fNTracks==header->fNTracks) &&
- (it->second->fNSelTracks==header->fNSelTracks) &&
- (it->second->fNTracklets==header->fNTracklets) &&
- (it->second->fVz==header->fVz)) {
- cout << "--> Duplicate found:" << endl;
- cout << " Run: " << it->second->fRun << " vs " << header->fRun << endl;
- cout << " Ntracks: " << it->second->fNTracks << " vs " << header->fNTracks << endl;
- cout << " Nseltracks: " << it->second->fNSelTracks << " vs " << header->fNSelTracks << endl;
- cout << " Ntracklets: " << it->second->fNTracklets << " vs " << header->fNTracklets << endl;
- cout << " Vz: " << it->second->fVz << " vs " << header->fVz << endl;
- cout << " Duplicate event " << header->fEvNumberInFile << " in file "
- << filesnames->At(header->fFileId)->GetName() << endl;
- ++totndup;
- }
-#endif
- continue;
- }
-#ifdef DUPLICATEVERBOSE
- MyHeader *cheader = new MyHeader(*header);
- eheaders.insert( pair<ULong64_t, MyHeader*>(evtid,cheader) );
-#endif
-#endif
- EvInfo *cinfo = new EvInfo;
- cinfo->fRun = header->fRun;
- cinfo->fArrId = i;
- cinfo->fEntry = ev;
- einfos.insert( pair<ULong64_t, EvInfo*>(evtid,cinfo) );
- if (lrun != header->fRun) {
- for (Int_t lx=0;lx<4;++lx) {
- lrun = header->fRun;
- TString cname(Form("L%d_run%d_map",lx,lrun));
- if (lx==3)
- cname=Form("TrigClass_run%d_map",lrun);
- TMap *map = dynamic_cast<TMap*>(output->FindObject(cname));
- if (map) {
- maps->Add(map);
- output->Remove(map);
- }
- }
- }
- }
- delete output;
- file->Close();
- delete file;
- if ((nEvents>0) && (einfos.size()>=(UInt_t)nEvents))
- break;
- }
- cout << "Found total: " << totnent << endl;
- cout << "Found suspicious: " << totnsus << endl;
-#ifdef DUPLICATECHECK
- cout << "Found potential duplicates: " << totncan << endl;
-#ifdef DUPLICATEVERBOSE
- cout << "Found duplicates: " << totndup << endl;
-#endif
-#endif
-
- // sort output tree
- std::vector<TFile*> files(nfiles);
- std::vector<TTree*> trees(nfiles);
- std::vector<MyHeader*> headers(nfiles);
- std::vector<TClonesArray*> partas(nfiles);
- std::vector<TClonesArray*> trklas(nfiles);
-
- MyHeader *newHeader = new MyHeader;
- if (TClass::GetClass("MyPart"))
- TClass::GetClass("MyPart")->IgnoreTObjectStreamer();
- TClonesArray *newParts = new TClonesArray("MyPart",10000);
- if (TClass::GetClass("MyTracklet"))
- TClass::GetClass("MyTracklet")->IgnoreTObjectStreamer();
- TClonesArray *newTracklets = new TClonesArray("MyTracklet",10000);
-
- TFile *newFile = TFile::Open(outFileName,"recreate","Merged/Sorted MyTree created by mergeCorr.C",5);
- TDirectory::TContext cd(newFile);
- for (Int_t m=0;m<maps->GetEntries();++m) {
- TMap *map = (TMap*)maps->At(m);
- //map->Print();
- map->Write(map->GetName(),TObject::kSingleKey);
- }
- TTree *newTree = new TTree("MyTree", "MyTree created by MyAliTask.cxx");
- newTree->Branch("header",&newHeader, 32*1024, 99);
- newTree->Branch("parts",&newParts, 256*1024, 99);
- newTree->Branch("tracklets",&newTracklets, 256*1024, 99);
- newTree->SetDirectory(newFile);
- Int_t newEntries = 0;
- for (map<ULong64_t, EvInfo*>::iterator it = einfos.begin();it!=einfos.end();++it) {
- Short_t id = it->second->fArrId;
- UInt_t entry = it->second->fEntry;
- if (!trees.at(id)) { // connect new tree
- TString fname(filenames->At(id)->GetTitle());
- TFile *file = TFile::Open(fname);
- TDirectory::TContext cd2(newFile,file);
- files[id] = file;
- TList *output = dynamic_cast<TList*>(file->Get("output"));
- TTree *tree = dynamic_cast<TTree*>(output->FindObject("MyTree"));
- trees[id] = tree;
- output->Remove(tree);
- output->SetOwner(1);
- delete output;
- tree->SetBranchAddress("header",&headers[id]);
- tree->SetBranchAddress("parts",&partas[id]);
- tree->SetBranchAddress("tracklets",&trklas[id]);
- }
- files[id]->cd();
- trees.at(id)->GetEntry(entry);
- *newHeader = *(headers[id]);
- newParts->Clear();
- Int_t nparts = partas[id]->GetEntries();
- for (Int_t p = 0; p < nparts; ++p) {
- MyPart *orig = dynamic_cast<MyPart*>(partas[id]->At(p));
- new((*newParts)[p]) MyPart(*orig);
- }
- Int_t ntrkls = trklas[id]->GetEntries();
- for (Int_t p = 0; p < ntrkls; ++p) {
- MyTracklet *orig = dynamic_cast<MyTracklet*>(trklas[id]->At(p));
- new((*newTracklets)[p]) MyTracklet(*orig);
- }
- newFile->cd();
- newTree->Fill();
- ++newEntries;
- delete it->second;
- if ((nEvents>0) && (newEntries>=nEvents))
- break;
- }
-
- newTree->Write();
- newFile->Close();
- delete newFile;
-
- for(Int_t i=0;i<nfiles;++i) {
- delete trees.at(i);
- delete files.at(i);
- }
-}
+++ /dev/null
-{
- if (1) {
- cout << "Loading " << gSystem->Getenv("PWD")
- << "/rootlogon.C" << endl;
- cout << "Using ROOT version "
- << gROOT->GetVersion() << endl;
- }
-
- gROOT->SetStyle("Plain");
- gStyle->SetPalette(1);
-
- gSystem->Load("libTree.so");
- if (gSystem->Getenv("TMPDIR"))
- gSystem->SetBuildDir(gSystem->Getenv("TMPDIR"));
- if (1) {
- delete gRandom;
- gRandom = new TRandom3(0);
- }
- gROOT->LoadMacro("TreeClasses.C+g");
- gROOT->LoadMacro("EventPool.C+g");
- gROOT->LoadMacro("AutoCorr.C+g");
- gROOT->LoadMacro("Utils.C+g");
- if (0) {
- gROOT->LoadMacro("EventPoolManager.C+g");
- gROOT->LoadMacro("anaCorr.C+g");
- }
-}
+++ /dev/null
-// $Id$
-
-// To run this:
-// .x rootlogon.C
-// .x runAutoCorr.C+
-
-#if !defined(__CINT__) || defined(__MAKECINT__)
-#include "TreeClasses.h"
-#include "EventPool.h"
-#include "AutoCorr.h"
-#include "TStopwatch.h"
-#include "TString.h"
-#include "TROOT.h"
-#endif
-
-bool debug = true;
-
-void runAutoCorr(const char* dataFile =
- "../rootfiles/res_LHC10c_09212010/merged_run120824.root",
- // "../rootfiles/res_LHC10e_09122010/mergedruns/merged_run130601.root",
- const char* outFileName = "output/120824.root")
-{
- TFile* outFile = new TFile(outFileName, "recreate");
- const double PI = TMath::Pi();
- Int_t poolsize = 50;
- Int_t nMix = 5;
- const Int_t nMultBins = 1;
- Double_t multbin[nMultBins+1] = {50, 400};
- const Int_t nZvtxBins = 3;
- Double_t zvtxbin[nZvtxBins+1] = {-10, -3, 4, 11};
- const Int_t nPtBins = 1;
- Double_t ptbin[nPtBins+1] = {1.0, 2.0};
- Double_t etaMin = -1.4; // for track cuts
- Double_t etaMax = 1.4;
- Double_t ptMin = ptbin[0]; // for track cuts
- Double_t ptMax = ptbin[nPtBins];
- Double_t dEtaMax = 3.0; // for histogram limits
- TString sDataSet = "res_LHC10e_09122010";
-
- // Event cuts
- Int_t minVc = 25;
- Int_t maxNTracklets = 200;
- Double_t zMin = zvtxbin[0];
- Double_t zMax = zvtxbin[nZvtxBins];
- Int_t runNumber = -1;
- Int_t firstRunNumber = -1;
-
- TString sMultBins(""), sZvtxBins(""), sPtBins("");
- for (int i=0; i<nMultBins+1; i++) sMultBins.Append(Form("%.0f ", multbin[i]));
- for (int i=0; i<nZvtxBins+1; i++) sZvtxBins.Append(Form("%.1f ", zvtxbin[i]));
- for (int i=0; i<nPtBins+1; i++) sPtBins.Append(Form("%.2f ", ptbin[i]));
-
- // Document the applied cuts, get quick look with cuts->Print()
- TList* cuts = new TList();
- cuts->Add(new TNamed("data_sample", sDataSet.Data()));
- cuts->Add(new TNamed("mult_bins", sMultBins.Data()));
- cuts->Add(new TNamed("zvtx_bins", sZvtxBins.Data()));
- cuts->Add(new TNamed("pt1_pt2_bins", sPtBins.Data()));
- cuts->Add(new TNamed("mass_cut", "No mass cut"));
- cuts->Add(new TNamed("eta_range", Form("%.1f < #eta < %.1f", etaMin, etaMax)));
- cuts->Add(new TNamed("z_vtx_range", Form("%.1f < z-vtx < %.1f", zMin, zMax)));
- cuts->Add(new TNamed("min_vc_ntracks", Form("No. vtx. contributors > %i", minVc)));
- cuts->Add(new TNamed("max_ntracklets", Form("No. trackets > %i", maxNTracklets)));
- cuts->Add(new TNamed("pool_size", Form("Event pool size = %i", poolsize )));
- cuts->Add(new TNamed("ntracks_mix", Form("Mixed tracks per real track = %i", nMix )));
- cuts->Add(new TNamed("comments", "Using dphi x deta = 0.03 x 0.03 cut"));
-
- // Histos to hold binning info
- TH1F* hMultBins = new TH1F("hMultBins", "Event multiplicity binning",
- nMultBins, multbin);
- TH1F* hZvtxBins = new TH1F("hZvtxBins", "Event Z-vertex binning",
- nZvtxBins, zvtxbin);
- TH1F* hPtBins = new TH1F("hPtBins", "p_{T} binning",
- nPtBins, ptbin);
-
- TChain *tree = new TChain("MyTree");
- tree->Add(dataFile);
- cout << "Found " << tree->GetEntries() << " entries in tree!" << endl;
-
- MyHeader* ev = 0;
- TClonesArray* trk = 0;
-
- TBranch* evBranch = tree->GetBranch("header");
- evBranch->SetAddress(&ev);
-
- TBranch* trBranch = tree->GetBranch("parts");
- trBranch->SetAddress(&trk);
-
- AutoCorr* ac = new AutoCorr();
- // ac->InitEventPools(poolsize, nMultBins, multbin, nZvtxBins, zvtxbin);
- ac->InitEventPools(poolsize, nMultBins, multbin, nZvtxBins, zvtxbin, ptMin, ptMax);
-
- if (debug) {
- cout << nMultBins << " x " << nZvtxBins
- << " event pool(s) to initialize." << endl;
- cout << "Mult. binning: ";
- for (int j=0; j<=nMultBins; j++) cout << multbin[j] << " ";
- cout << endl;
- cout << "Zvtx. binning: ";
- for (int j=0; j<=nZvtxBins; j++) cout << zvtxbin[j] << " ";
- cout << endl;
- }
-
- TH2F* hSig[nMultBins];
- TH2F* hBkg[nMultBins];
- TH2F* hSB[nMultBins];
- TH2F* hSigPt[nMultBins][nPtBins];
- TH2F* hBkgPt[nMultBins][nPtBins];
- TH2F* hSBPt[nMultBins][nPtBins];
-
- TH1F* hMult[nMultBins];
- TH1F* hEta = new TH1F("hEta", "hEta", 600, -3, 3);
- TH1F* hDEta = new TH1F("hDEta", "hDEta", 600, -3, 3);
- TH1F* hDPhi = new TH1F("hDPhi", "hDPhi", 60, -PI/2, 3*PI/2);
- TH1F* hMultAll = new TH1F("hMultAll", "hMultAll", 1000, 0, 1000);
- TH1F* hMultSel = new TH1F("hMultSel", "hMultSel", 1000, 0, 1000);
- TH1F* hPtAll = new TH1F("hPtAll", "hPtAll", 100, 0., 10.);
- TH1F* hPtSel = new TH1F("hPtSel", "hPtSel", 100, 0., 10.);
- TH1F* hVc = new TH1F("hVc", "Vertex contributor tracks", 1000, 0, 1000);
- TH1F* hNTracklets = new TH1F("hNTracklets", "Tracklets", 1000, 0, 1000);
- TH1F* hZvtx = new TH1F("hZvtx", "Event Z-vertex", 60, -30, 30);
- TH1F* hMass = new TH1F("hMass", "hMass", 1000, 0, 10);
- TH1F* hMassBg = new TH1F("hMassBg", "hMassBg", 1000, 0, 10);
- TH1F* hNevts = new TH1F("hNevts", "# events: available, sampled, accepted", 3, 0, 3);
-
- //----------------------
- TString sMultBin[nMultBins];
- for (int i=0; i<nMultBins; i++) {
- double m1 = hMultBins->GetBinLowEdge(i+1);
- double m2 = m1 + hMultBins->GetBinWidth(i+1);
- sMultBin[i] = Form("%.0f - %.0f tracks", m1, m2);
- }
- TString sPtBin[nPtBins];
- for (int i=0; i<nPtBins; i++) {
- double pt1 = hPtBins->GetBinLowEdge(i+1);
- double pt2 = pt1 + hPtBins->GetBinWidth(i+1);
- sPtBin[i] = Form("%.2f - %.2f GeV/c", pt1, pt2);
- }
- //---------------------
-
- for (int iM=0; iM<nMultBins; iM++) {
- char* name;
- char* title;
- name = Form("hMult_%i", iM);
- title = Form("%s, %.2f - %.2f GeV/c", sMultBin[iM].Data(), ptMin, ptMax);
-
- hMult[iM] = new TH1F(name, title, 1000, 0, 1000);
-
- name = Form("hSig_%i", iM);
- hSig[iM] = new TH2F(name, title, 60, -PI/2, 3*PI/2, 60, -dEtaMax, dEtaMax);
- hSig[iM]->Sumw2();
- name = Form("hBkg_%i", iM);
- hBkg[iM] = new TH2F(name, title, 60, -PI/2, 3*PI/2, 60, -dEtaMax, dEtaMax);
- hBkg[iM]->Sumw2();
- name = Form("hSB_%i", iM);
- hSB[iM] = new TH2F(name, title, 60, -PI/2, 3*PI/2, 60, -dEtaMax, dEtaMax);
- hSB[iM]->Sumw2();
-
- for (int iPt=0; iPt<nPtBins; iPt++) {
- name = Form("hSigPt_%i_%i", iM, iPt);
- title = Form("%s, %s", sMultBin[iM].Data(), sPtBin[iPt].Data());
- hSigPt[iM][iPt] = new TH2F(name, title, 60, -PI/2, 3*PI/2, 60, -dEtaMax, dEtaMax);
- hSigPt[iM][iPt]->Sumw2();
- name = Form("hBkgPt_%i_%i", iM, iPt);
- hBkgPt[iM][iPt] = new TH2F(name, title, 60, -PI/2, 3*PI/2, 60, -dEtaMax, dEtaMax);
- hBkgPt[iM][iPt]->Sumw2();
- name = Form("hSBPt_%i_%i", iM, iPt);
- hSBPt[iM][iPt] = new TH2F(name, title, 60, -PI/2, 3*PI/2, 60, -dEtaMax, dEtaMax);
- hSBPt[iM][iPt]->Sumw2();
- }
-
- }
-
- TH2F* hDR = new TH2F("hDR", "hDR", 100, -0.01, 0.01, 100, -0.005, 0.005);
- hDR->SetTitle("#Delta#phi-#Delta#eta Distribution (no cut);#Delta#phi;#Delta#eta");
- TH2F* hDRcut = new TH2F("hDRcut", "hDRcut", 100, -0.2, 0.2, 100, -0.1, 0.1);
- hDRcut->SetTitle("#Delta#phi-#Delta#eta Distribution;#Delta#phi;#Delta#eta");
- hDEta->SetLineColor(kRed);
- hDPhi->SetLineColor(kBlue);
- TH2F* hMixEff = new TH2F("hMixEff", "", 100, 0., 500., 11, -0.05, 1.05);
- hMixEff->SetTitle("Mixing efficiency vs. multiplicity;"
- "multiplicity;"
- "accepted/sampled tracks");
- // ---------------------------------
-
- Int_t nEvents = tree->GetEntries();
- Int_t nAcceptedEvents = 0;
-
- TStopwatch* watch = new TStopwatch();
- watch->Start();
- for (int n=0; n<nEvents; n++) {
- evBranch->GetEntry(n);
-
- // need to add trigger mask here!!
- if (!ac->IsEventOk(*ev, minVc, maxNTracklets, zMin, zMax))
- continue;
-
- trBranch->GetEntry(n);
- if (n % 1 == 0) cout << n << endl;
-
- hMultAll->Fill(ev->fNSelTracks); // No cuts
-
- Int_t multBin = hMultBins->FindBin(ev->fNSelTracks) - 1;
- Int_t zvtxBin = hZvtxBins->FindBin(ev->fVz) - 1;
- if (!ac->InBounds(multBin, 0, nMultBins-1))
- continue;
- if (!ac->InBounds(zvtxBin, 0, nZvtxBins-1))
- continue;
- if (runNumber == -1) {
- runNumber = firstRunNumber = ev->fRun;
- }
- runNumber = ev->fRun;
- if (runNumber != firstRunNumber) {
- cout << "Warning: Run number has changed!" << endl;
- firstRunNumber = runNumber;
- }
- nAcceptedEvents++;
-
- Int_t ntracks = trk->GetEntries();
- hMultSel->Fill(ntracks); // After event selection
- hMult[multBin]->Fill(ntracks);
- hVc->Fill(ev->fVc);
- hZvtx->Fill(ev->fVz);
- hNTracklets->Fill(ev->fNTracklets);
- EventPool* pool = ac->GetEventPool(multBin, zvtxBin);
-
- if (pool->IsPoolReady()) {
- for (int i=0; i<ntracks; i++) {
- MyPart* t1 = (MyPart*)trk->At(i);
- hPtAll->Fill(t1->fPt);
- Int_t ptBin1 = hPtBins->FindBin(t1->fPt) - 1;
- if (ptBin1 < 0 || ptBin1 >= nPtBins)
- continue;
- if (!ac->IsTrackOk(*t1, etaMin, etaMax, ptMin, ptMax))
- continue;
-
- hEta->Fill(t1->fEta);
- hPtSel->Fill(t1->fPt);
-
- // --- Same-event correlation loop ---
- for (int j=i+1; j<ntracks; j++) {
- MyPart* t2 = (MyPart*)trk->At(j);
- if (!ac->IsTrackOk(*t2, etaMin, etaMax, ptMin, ptMax))
- continue;
-
- Double_t deta = ac->DeltaEta(*t1, *t2);
- Double_t dphi = ac->DeltaPhi(*t1, *t2);
- hDR->Fill(dphi, deta); // No cut - contamination pk.
-
- if (ac->IsPairOk(*t1, *t2)) {
-
- hDEta->Fill(deta);
- hDPhi->Fill(dphi);
- hSig[multBin]->Fill(dphi, deta);
- hDRcut->Fill(dphi, deta); // Show elliptical hole from cut
- hMass->Fill(ac->InvMass(*t1, *t2));
-
- Int_t ptBin2 = hPtBins->FindBin(t2->fPt) - 1;
- if (ptBin2 < 0 || ptBin2 >= nPtBins)
- continue;
- if (ptBin1==ptBin2)
- hSigPt[multBin][ptBin1]->Fill(dphi, deta);
- }
-
- Int_t jMix = 0, jWatcher = 0, jLimit = 50*nMix;
- for (jMix=0; jMix<nMix;) {
- if (++jWatcher > jLimit) {
- cout << "Warning: mix loop count > " << jLimit
- << ". Mixed with " << jMix << " tracks."
- << endl;
- break;
- }
- MyPart* tbg = pool->GetRandomTrack();
- if (!ac->IsTrackOk(*tbg, etaMin, etaMax, ptMin, ptMax))
- continue;
-
- if (ac->IsMixedPairOk(*t1, *tbg)) {
- hMassBg->Fill(ac->InvMass(*t1, *tbg));
- Double_t deta_bg = ac->DeltaEta(*t1, *tbg);
- Double_t dphi_bg = ac->DeltaPhi(*t1, *tbg);
- hBkg[multBin]->Fill(dphi_bg, deta_bg);
-
- Int_t ptBin3 = hPtBins->FindBin(tbg->fPt) - 1;
- if (ptBin1==ptBin3)
- hBkgPt[multBin][ptBin3]->Fill(dphi_bg, deta_bg);
-
- jMix++;
- }
- }
- hMixEff->Fill(ntracks, (double)jMix/jWatcher);
- }
- }
- }
- else {
- ac->UpdatePools(n, ev, trk);
- }
- }
- watch->Stop();
- cout << "\nFinished event loop.\n" << endl;
- cout << nEvents << " events sampled from " << dataFile << endl;
- cout << nAcceptedEvents << " events selected." << endl;
- cout << "Used ROOT " << gROOT->GetVersion() << endl;
- cout << " CPU time: " << watch->CpuTime() << endl;
- cout << " Real time: " << watch->RealTime() << endl;
-
- cuts->Add(new TNamed("run_number", Form("Run number = %i", runNumber )));
- cuts->Add(new TNamed("n_events_file",Form("Events in file = %i", tree->GetEntries())));
- cuts->Add(new TNamed("n_events_loop", Form("Events looped = %i", nEvents )));
- cuts->Add(new TNamed("n_events_kept", Form("Events kept = %i", nAcceptedEvents )));
- cuts->Add(new TNamed("cpu_time", Form("CPU time = %f", watch->RealTime() )));
- cuts->Add(new TNamed("real_time", Form("Real time = %f", watch->RealTime() )));
- cuts->Add(new TNamed("root_version", Form("ROOT version = %s", gROOT->GetVersion() )));
- cuts->Add(new TNamed("hostname", Form("Hostname = %s", gSystem->HostName() )));
-
- hMassBg->Scale(1./nMix);
- for (int iM=0; iM<nMultBins; iM++) {
- hBkg[iM]->Scale(1./nMix);
- hSB[iM]->Divide(hSig[iM], hBkg[iM]);
-
- for (int iPt=0; iPt<nPtBins; iPt++) {
- hBkgPt[iM][iPt]->Scale(1./nMix);
- hSBPt[iM][iPt]->Divide(hSigPt[iM][iPt], hBkgPt[iM][iPt]);
- }
- }
- hNevts->SetBinContent(1, tree->GetEntries());
- hNevts->SetBinContent(2, nEvents);
- hNevts->SetBinContent(3, nAcceptedEvents);
-
- outFile->Write();
- cuts->Write("cuts", TObject::kSingleKey);
- outFile->Close();
- return;
-}
-
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff