apply corrections PbPb - Markus Koehler

diff --git a/PWG0/dNdPt/macros/ExtractSpectra.C b/PWG0/dNdPt/macros/ExtractSpectra.C
new file mode 100644 (file)
index 0000000..5f4de8b
--- /dev/null
+++ b/PWG0/dNdPt/macros/ExtractSpectra.C
@@ -0,0 +1,187 @@
+
+
+void ExtractSpectra() {
+
+
+       //macro to extract and plot 1d and 2d histograms
+       //for efficiency and efficiency corrected spectra
+
+       gROOT->SetStyle("Plain");
+       gStyle->SetPalette(1);
+
+
+       TFile *f = TFile::Open("rootfile");
+       // rootfile has to be generated by the AlidNdPtAnalysisPbPb class
+
+
+       AlidNdPtAnalysisPbPb *obj =(AlidNdPtAnalysisPbPb*)dNdPtPbPb->FindObject("dNdPtAnalysisPbPb");
+       //dNdPtPbPb is container
+
+       ////////////////
+       // still to do : calculate errors properly
+       ////////////////
+
+
+////////////////////////////////////
+// Event level
+////////////////////////////////////
+
+       //The projections are 0->z Vertex ; 1->Multiplicity
+
+       THnSparse *fSparseEventsRec = obj->GetRecEventMatrix(); //all reconstructed events
+       THnSparse *fSparseEventsTrigger = obj->GetTriggerEventMatrix(); //all triggered events
+       THnSparse *fSparseGenEvents = obj->GetGenEventMatrix(); //all triggered events
+       
+       TH2D * fHist2EventTriggerEff = (TH2D*) fSparseEventsTrigger->Projection(0,1)->Clone();
+       fHist2EventTriggerEff->Divide(obj->GetGenEventMatrix()->Projection(0,1));
+
+       TH2D * fHist2Event = (TH2D*) fSparseEventsRec->Projection(0,1)->Clone();
+       fHist2Event->Divide(fHist2EventTriggerEff);
+
+       Double_t NumberOfEvents = (Double_t) fHist2Event->Integral();
+
+//     printf("========================\n");
+//     printf("Number of events is %lf \n",NumberOfEvents);
+//     printf("========================\n");
+
+////////////////////////////////////
+// Track level
+////////////////////////////////////
+
+       /////////////
+       // 1. extract primaries
+       // 2. extract efficiency for trigger, vertex and tracking (1d and 2d)
+       // 3. correct for efficiency 
+       /////////////
+
+       THnSparse *fSparseTrackRecPrim = obj->GetRecTrackMatrix();//all reconstructed track
+       fSparseTrackRecPrim->Add(obj->GetRecSecTrackMatrix(),-1);//primary tracks
+
+       THnSparse *fSparseGenPrimTrack = obj->GetGenPrimTrackMatrix();//generated primary tracks
+       THnSparse *fSparseGenSecTrack = obj->GetRecSecTrackMatrix();//secondary tracks
+
+       THnSparse *fSparseTriggerTrackEvent = obj->GetTriggerTrackEventMatrix();//Tracks from triggered events
+       THnSparse *fSparseVtxTrackEvent = obj->GetRecTrackEventMatrix();//Tracks from events with rec. vtx
+       THnSparse *fSparseGenTrackEvent = obj->GetGenTrackEventMatrix();//generated TrackEvent matrix
+
+       //////////////
+       // efficiencies and normalization for 1d
+       //////////////
+
+       //The projections are 0->z Vertex ; 1->pT ; 2->Eta
+
+       //tracking
+       TH1D * fHistTrackEffzV = (TH1D*) fSparseTrackRecPrim->Projection(0)->Clone();
+       fHistTrackEffzV->Divide(fHistTrackEffzV,fSparseGenPrimTrack->Projection(0),1,1,"B");
+       fHistTrackEffzV->Divide(fSparseGenPrimTrack->Projection(0));
+       TH1D * fHistTrackEffPt = (TH1D*) fSparseTrackRecPrim->Projection(1)->Clone();
+       fHistTrackEffPt->Divide(fHistTrackEffPt,fSparseGenPrimTrack->Projection(1) , 1 ,1 , "B");
+       TH1D * fHistTrackEffEta = (TH1D*) fSparseTrackRecPrim->Projection(2)->Clone();
+       fHistTrackEffEta->Divide(fHistTrackEffEta,fSparseGenPrimTrack->Projection(2),1,1,"B");
+
+       //trigger
+       TH1D * fHistTriggerEffzV = (TH1D*) fSparseTriggerTrackEvent->Projection(0)->Clone();
+       fHistTriggerEffzV->Divide(fHistTriggerEffzV,fSparseGenTrackEvent->Projection(0),1,1,"B");
+       TH1D * fHistTriggerEffPt = (TH1D*) fSparseTriggerTrackEvent->Projection(1)->Clone();
+       fHistTriggerEffPt->Divide(fHistTriggerEffPt,fSparseGenTrackEvent->Projection(1),1,1,"B");
+       TH1D * fHistTriggerEffEta = (TH1D*) fSparseTriggerTrackEvent->Projection(2)->Clone();
+       fHistTriggerEffEta->Divide(fHistTriggerEffEta,fSparseGenTrackEvent->Projection(2),1,1,"B");
+
+
+       //rec Vertex
+       TH1D * fHistVtxEffzV = (TH1D*) fSparseVtxTrackEvent->Projection(0)->Clone();
+       fHistVtxEffzV->Divide(fHistVtxEffzV,fSparseGenTrackEvent->Projection(0),1,1,"B");
+       TH1D * fHistVtxEffPt = (TH1D*) fSparseVtxTrackEvent->Projection(1)->Clone();
+       fHistVtxEffPt->Divide(fHistVtxEffPt,fSparseGenTrackEvent->Projection(1),1,1,"B");
+       TH1D * fHistVtxEffEta = (TH1D*) fSparseVtxTrackEvent->Projection(2)->Clone();
+       fHistVtxEffEta->Divide(fHistVtxEffEta,fSparseGenTrackEvent->Projection(2),1,1,"B");
+
+       TH1D *fHistTrackzV = (TH1D*)fSparseTrackRecPrim->Projection(0)->Clone();
+       fHistTrackzV->Divide(fHistTrackEffzV);
+       fHistTrackzV->Divide(fHistTriggerEffzV);
+       fHistTrackzV->Divide(fHistVtxEffzV);
+       fHistTrackzV->Scale(1./NumberOfEvents);
+       fHistTrackzV->Scale(1,"width");
+
+
+       TH1D *fHistTrackPt = (TH1D*)fSparseTrackRecPrim->Projection(1)->Clone();
+       fHistTrackPt->Divide(fHistTrackEffPt);
+       fHistTrackPt->Divide(fHistTriggerEffPt);
+       fHistTrackPt->Divide(fHistVtxEffPt);
+       fHistTrackPt->Scale(1./NumberOfEvents);
+       fHistTrackPt->Scale(1,"width");
+
+
+       TH1D *fHistTrackEta = (TH1D*)fSparseTrackRecPrim->Projection(2)->Clone();
+       fHistTrackEta->Divide(fHistTrackEffEta);
+       fHistTrackEta->Divide(fHistTriggerEffEta);
+       fHistTrackEta->Divide(fHistVtxEffEta);
+       fHistTrackEta->Scale(1./NumberOfEvents);
+       fHistTrackEta->Scale(1,"width");
+
+       //////////////
+       // 2 dim efficiencies 
+       //////////////
+
+       //The projections are 0->z Vertex ; 1->pT ; 2->Eta
+
+       //tracking
+       TH2D * fHist2TrackEffzV_Pt = (TH2D*) fSparseTrackRecPrim->Projection(0,1)->Clone();
+       fHist2TrackEffzV_Pt->Divide(fSparseGenPrimTrack->Projection(0,1));
+       TH2D * fHist2TrackEffzV_Eta = (TH2D*) fSparseTrackRecPrim->Projection(0,2)->Clone();
+       fHist2TrackEffzV_Eta->Divide(fSparseGenPrimTrack->Projection(0,2));
+       TH2D * fHist2TrackEffPt_Eta = (TH2D*) fSparseTrackRecPrim->Projection(1,2)->Clone();
+       fHist2TrackEffPt_Eta->Divide(fSparseGenPrimTrack->Projection(1,2));
+
+       //trigger
+       TH2D * fHist2TriggerEffzV_Pt = (TH2D*) fSparseTriggerTrackEvent->Projection(0,1)->Clone();
+       fHist2TriggerEffzV_Pt->Divide(fSparseGenTrackEvent->Projection(0,1));
+       TH2D * fHist2TriggerEffzV_Eta = (TH2D*) fSparseTriggerTrackEvent->Projection(0,2)->Clone();
+       fHist2TriggerEffzV_Eta->Divide(fSparseGenTrackEvent->Projection(0,2));
+       TH2D * fHist2TriggerEffPt_Eta = (TH2D*) fSparseTriggerTrackEvent->Projection(1,2)->Clone();
+       fHist2TriggerEffPt_Eta->Divide(fSparseGenTrackEvent->Projection(1,2));
+
+       //rec Vertex
+       TH2D * fHist2VtxEffzV_Pt = (TH2D*) fSparseVtxTrackEvent->Projection(0,1)->Clone();
+       fHist2VtxEffzV_Pt->Divide(fSparseGenTrackEvent->Projection(0,1));
+       TH2D * fHist2VtxEffzV_Eta = (TH2D*) fSparseVtxTrackEvent->Projection(0,2)->Clone();
+       fHist2VtxEffzV_Eta->Divide(fSparseGenTrackEvent->Projection(0,2));
+       TH2D * fHist2VtxEffPt_Eta = (TH2D*) fSparseVtxTrackEvent->Projection(1,2)->Clone();
+       fHist2VtxEffPt_Eta->Divide(fSparseGenTrackEvent->Projection(1,2));
+
+       //////////////
+       // correct for efficiencies (2d)
+       //////////////
+
+       TH2D *fHist2TrackzV_Pt = (TH2D*)fSparseTrackRecPrim->Projection(0,1)->Clone();
+       fHist2TrackzV_Pt->Divide(fHist2TrackEffzV_Pt);
+       fHist2TrackzV_Pt->Divide(fHist2TriggerEffzV_Pt);
+       fHist2TrackzV_Pt->Divide(fHist2VtxEffzV_Pt);
+
+       TH2D *fHist2TrackPt_Eta = (TH2D*)fSparseTrackRecPrim->Projection(1,2)->Clone();
+       fHist2TrackPt_Eta->Divide(fHist2TrackEffPt_Eta);
+       fHist2TrackPt_Eta->Divide(fHist2TriggerEffPt_Eta);
+       fHist2TrackPt_Eta->Divide(fHist2VtxEffPt_Eta);
+
+       TH2D *fHist2TrackzV_Eta = (TH2D*)fSparseTrackRecPrim->Projection(0,2)->Clone();
+       fHist2TrackzV_Eta->Divide(fHist2TrackEffzV_Eta);
+       fHist2TrackzV_Eta->Divide(fHist2TriggerEffzV_Eta);
+       fHist2TrackzV_Eta->Divide(fHist2VtxEffzV_Eta);
+
+/*
+
+       //for example
+       TCanvas *cPtEta = new TCanvas("cPtEta","cPtEta",900,700);
+       fHist2TrackEffPt_Eta->SetStats(0);
+       fHist2TrackEffPt_Eta->Draw("col");
+
+*/
+
+
+}
+
+
+
+
+
+