]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG0/dNdEta/dNdEtaAnalysis.cxx
Fixing bugs in FMD reconstruction. Everything should work now
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / dNdEtaAnalysis.cxx
index e0956e3e3c2a7e129410409f85341a447ceb9477..40d02f346e6167d8d084fdfdadae9e3415933f75 100644 (file)
@@ -43,7 +43,7 @@ dNdEtaAnalysis::dNdEtaAnalysis() :
 }
 
 //____________________________________________________________________
-dNdEtaAnalysis::dNdEtaAnalysis(Char_t* name, Char_t* title, AliPWG0Helper::AnalysisMode analysisMode) :
+dNdEtaAnalysis::dNdEtaAnalysis(const Char_t* name, const Char_t* title, AliPWG0Helper::AnalysisMode analysisMode) :
   TNamed(name, title),
   fData(0),
   fMult(0),
@@ -201,7 +201,7 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
   // the measured result is not used up to a multiplicity of multCut (the bin at multCut is the first that is used!)
   //
 
-  fTag.Form("Correcting dN/deta spectrum >>> %s <<<. Correction type: %d, pt cut: %.2f.", tag, (Int_t) correctionType, ptCut);
+  fTag.Form("Correcting dN/deta spectrum (data: %d) >>> %s <<<. Correction type: %d, pt cut: %.2f.", (Int_t) fAnalysisMode, tag, (Int_t) correctionType, ptCut);
   Printf("\n\n%s", fTag.Data());
 
   // set corrections to 1
@@ -209,8 +209,8 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
 
   if (correction && correctionType != AlidNdEtaCorrection::kNone)
   {
-    TH3F* trackCorr = fData->GetTrackCorrection()->GetCorrectionHistogram();
-    TH2F* eventCorr = fData->GetEventCorrection()->GetCorrectionHistogram();
+    TH3* trackCorr = fData->GetTrackCorrection()->GetCorrectionHistogram();
+    TH2* eventCorr = fData->GetEventCorrection()->GetCorrectionHistogram();
 
     if (correctionType >= AlidNdEtaCorrection::kTrack2Particle)
       trackCorr->Multiply(correction->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetCorrectionHistogram());
@@ -268,6 +268,7 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
 
     TH2* eTrig =    correction->GetVertexRecoCorrection()->GetEventCorrection()->GetGeneratedHistogram();
     TH2* eTrigVtx = correction->GetVertexRecoCorrection()->GetEventCorrection()->GetMeasuredHistogram();
+    //TH1* eTrigVtx_projx = eTrigVtx->ProjectionX("eTrigVtx_projx", 2, rawMeasured->GetNbinsY()+1);
     TH1* eTrigVtx_projx = eTrigVtx->ProjectionX("eTrigVtx_projx", 2, rawMeasured->GetNbinsY()+1);
 
     //new TCanvas; correctedEvents->DrawCopy("TEXT");
@@ -281,6 +282,9 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
 
     Printf("%d triggered events with 0 mult. -- %d triggered events with vertex", triggeredEventsWith0Mult, allEventsWithVertex);
 
+    TH1* kineBias = (TH1*) vertexDist->Clone("kineBias");
+    kineBias->Reset();
+
     for (Int_t i = 1; i <= rawMeasured->GetNbinsX(); i++)
     {
       Double_t alpha = (Double_t) vertexDist->GetBinContent(i) / allEventsWithVertex;
@@ -291,6 +295,7 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
 
       Double_t fZ = eTrigVtx_projx->Integral(0, eTrigVtx_projx->GetNbinsX()+1) / eTrigVtx_projx->GetBinContent(i) *
         eTrig->GetBinContent(i, 1) / eTrig->Integral(0, eTrig->GetNbinsX()+1, 1, 1);
+      kineBias->SetBinContent(i, fZ);
 
       events *= fZ;
 
@@ -303,18 +308,20 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
     }
 
     //new TCanvas; correctedEvents->DrawCopy("TEXT");
+    //new TCanvas; kineBias->DrawCopy();
   }
 
   fData->PrintInfo(ptCut);
 
-  TH3F* dataHist = fData->GetTrackCorrection()->GetGeneratedHistogram();
+  TH3* dataHist = fData->GetTrackCorrection()->GetGeneratedHistogram();
 
   // integrate multiplicity axis out (include under/overflow bins!!!)
-  TH2F* tmp = fData->GetEventCorrection()->GetGeneratedHistogram();
+  TH2* tmp = fData->GetEventCorrection()->GetGeneratedHistogram();
 
   TH1D* vertexHist = (TH1D*) tmp->ProjectionX("_px", 0, tmp->GetNbinsY() + 1, "e");
 
   // create pt hist
+  if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS)
   {
     // reset all ranges
     dataHist->GetXaxis()->SetRange(0, 0);
@@ -353,14 +360,19 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
   dataHist->GetYaxis()->SetRange(0, 0);
   dataHist->GetZaxis()->SetRange(0, 0);
 
-  // integrate over pt (with pt cut)
+  // integrate over pt (with pt cut) (TPC, TPCITS) or multiplicity (SPD)
   Int_t ptLowBin = 1;
-  if (ptCut > 0)
+  if (ptCut > 0 && fAnalysisMode != AliPWG0Helper::kSPD)
     ptLowBin = dataHist->GetZaxis()->FindBin(ptCut);
+    
+  //new TCanvas; dataHist->DrawCopy();
 
+  //dataHist->Sumw2();
   dataHist->GetZaxis()->SetRange(ptLowBin, dataHist->GetZaxis()->GetNbins()+1);
-  printf("pt range %d %d\n", ptLowBin, dataHist->GetZaxis()->GetNbins()+1);
+  printf("pt/multiplicity range %d %d\n", ptLowBin, dataHist->GetZaxis()->GetNbins()+1);
   TH2D* vtxVsEta = dynamic_cast<TH2D*> (dataHist->Project3D("yx2e"));
+  
+  //new TCanvas; vtxVsEta->Draw("COLZ");
 
   dataHist->GetZaxis()->SetRange(0, 0);
   vtxVsEta->GetXaxis()->SetTitle(dataHist->GetXaxis()->GetTitle());
@@ -368,10 +380,12 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
 
   if (vtxVsEta == 0)
   {
-    printf("ERROR: pt integration failed\n");
+    printf("ERROR: pt/multiplicity integration failed\n");
     return;
   }
 
+  //new TCanvas(tag, tag, 800, 600);  vtxVsEta->DrawCopy("COLZ");
+
   // clear result histograms
   for (Int_t vertexPos=0; vertexPos<kVertexBinning; ++vertexPos)
   {
@@ -391,19 +405,27 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
       Int_t vertexBinEnd   = vertexHist->GetXaxis()->FindBin(vertexRangeEnd[vertexPos]);
 
       const Int_t *binBegin = 0;
-      // adjust acceptance range for SPD
+      const Int_t maxBins = 60;
+
+      // adjust acceptance range
+      // produce with drawPlots.C: DetermineAcceptance(...)
       if (fAnalysisMode == AliPWG0Helper::kSPD)
       {
-        //const Int_t binBegin[30] = { 18, 16, 15, 14, 13, 13, 12, 11, 9, 7, 6, 5, 5, 5, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 1 };
-        const Int_t binBeginSPD[30] = { -1, -1, -1, -1, 16, 14, 12, 10, 9, 7, 6, 5, 5, 4, 4, 4, 4, 3, 3, 3, 3, 2, 2, 2, 2, 2, -1, -1, -1, -1}; // limit in correction map is 5
+        //const Int_t binBeginSPD[30] = { 18, 16, 15, 14, 13, 13, 12, 11, 9, 7, 6, 5, 5, 5, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 1 }; // by eye
+        //const Int_t binBeginSPD[30] = { -1, -1, -1, -1, 16, 14, 12, 10, 9, 7, 6, 5, 5, 4, 4, 4, 4, 3, 3, 3, 3, 2, 2, 2, 2, 2, -1, -1, -1, -1}; // limit in correction map is 5
+
         //const Int_t binBegin[30] = { -1, -1, -1, 17, 15, 14, 12, 10, 8, 7, 6, 5, 4, 4, 4, 4, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, -1, -1, -1};  // limit in correction map is 10
-        binBegin = binBeginSPD;
 
+        //const Int_t binBeginSPD[30] = { -1, -1, -1, -1, 16, 15, 13, 11, 9, 8, 7, 6, 5, 5, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 3, -1, -1, -1, -1}; // limit 2
+        const Int_t binBeginSPD[maxBins] = {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 15, 14, 13, 12, 11, 10, 9, 9, 8, 7, 7, 6, 6, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}; //limit 5
+
+        binBegin = binBeginSPD;
       }
       else if (fAnalysisMode == AliPWG0Helper::kTPC)
       {
-//        const Int_t binBegin[30] = { -1, -1, -1, -1, 16, 14, 12, 10, 9, 7, 6, 5, 5, 4, 4, 4, 4, 3, 3, 3, 3, 2, 2, 2, 2, 2, -1, -1, -1, -1}; // limit in correction map is 5
-        const Int_t binBeginTPC[30] = {-1, -1, -1, -1, -1, -1, -1, -1, 9, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, -1, -1, -1, -1, -1, -1, -1, -1};
+        //const Int_t binBeginTPC[30] = {-1, -1, -1, -1, -1, -1, -1, -1, -1, 4, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1}; // limit 5, pt cut off 0.2 mev/c
+        const Int_t binBeginTPC[maxBins] = {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 9, 5, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}; // limit 5
+
         binBegin = binBeginTPC;
       }
       else if (fAnalysisMode == AliPWG0Helper::kTPCITS)
@@ -411,21 +433,30 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
         // TODO create map
       }
 
-      Int_t vtxBegin = binBegin[iEta - 1];
-      Int_t vtxEnd   = 18 + 1 - binBegin[30 - iEta];
+      Int_t vtxBegin = 1;
+      Int_t vtxEnd   = maxBins;
+
+      if (binBegin)
+      {
+        vtxBegin = binBegin[iEta - 1];
+        vtxEnd = 18 + 1 - binBegin[maxBins - iEta];
+      }
+      else
+        Printf("WARNING: No acceptance applied!");
       
       // eta range not accessible
       if (vtxBegin == -1)
         continue;
       
-      //Float_t eta = vtxVsEta->GetYaxis()->GetBinCenter(iEta);
       //Printf("%d %d | %d %d", vtxBegin, vertexHist->GetXaxis()->FindBin(GetVtxMin(eta)), vtxEnd, vertexHist->GetXaxis()->FindBin(-GetVtxMin(-eta)));
       //vtxBegin = vertexHist->GetXaxis()->FindBin(GetVtxMin(eta));
       //vtxEnd = vertexHist->GetXaxis()->FindBin(-GetVtxMin(-eta));
-      //printf("Eta bin: %d (%f) Vertex range: %d; before: %d %d ", iEta, eta, vertexPos, vertexBinBegin, vertexBinEnd);
+
+      //Float_t eta = vtxVsEta->GetYaxis()->GetBinCenter(iEta);
+      //printf("Eta bin: %d (%f) Vertex range: %d; before: %d %d (range) %d %d (acceptance)", iEta, eta, vertexPos, vertexBinBegin, vertexBinEnd, vtxBegin, vtxEnd);
       vertexBinBegin = TMath::Max(vertexBinBegin, vtxBegin);
       vertexBinEnd =   TMath::Min(vertexBinEnd, vtxEnd);
-      //Printf("after:  %d %d", vertexBinBegin, vertexBinEnd);
+      //Printf(" after:  %d %d", vertexBinBegin, vertexBinEnd);
 
       // no data for this bin
       if (vertexBinBegin > vertexBinEnd)
@@ -434,39 +465,55 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
         continue;
       }
 
-      Float_t totalEvents = vertexHist->Integral(vertexBinBegin, vertexBinEnd);
-      if (totalEvents == 0)
-      {
-        printf("WARNING: No events for hist %d %d %d\n", vertexPos, vertexBinBegin, vertexBinEnd);
-        continue;
-      }
-
+      Float_t totalEvents = 0;
       Float_t sum = 0;
       Float_t sumError2 = 0;
-      for (Int_t iVtx = vertexBinBegin; iVtx <= vertexBinEnd; iVtx++)
+      Float_t unusedTracks = 0;
+      Float_t unusedEvents = 0;
+      for (Int_t iVtx = 1; iVtx <= vtxVsEta->GetNbinsX(); iVtx++)
       {
-        if (vtxVsEta->GetBinContent(iVtx, iEta) != 0)
+        if (iVtx >= vertexBinBegin && iVtx <= vertexBinEnd)
         {
-          sum = sum + vtxVsEta->GetBinContent(iVtx, iEta);
+          if (vtxVsEta->GetBinContent(iVtx, iEta) != 0)
+          {
+            sum += vtxVsEta->GetBinContent(iVtx, iEta);
 
-          if (sumError2 > 10e30)
-            Printf("WARNING: sum of error2 is dangerously large - be prepared for crash... ");
+            if (sumError2 > 10e30)
+              Printf("WARNING: sum of error2 is dangerously large - be prepared for crash... ");
 
-          sumError2 = sumError2 + TMath::Power(vtxVsEta->GetBinError(iVtx, iEta),2);
+            sumError2 = sumError2 + TMath::Power(vtxVsEta->GetBinError(iVtx, iEta),2);
+          }
+          totalEvents += vertexHist->GetBinContent(iVtx);
         }
+        else
+        {
+          unusedTracks += vtxVsEta->GetBinContent(iVtx, iEta);
+          unusedEvents += vertexHist->GetBinContent(iVtx);
+        }
+      }
+
+      if (totalEvents == 0)
+      {
+        printf("WARNING: No events for hist %d %d %d\n", vertexPos, vertexBinBegin, vertexBinEnd);
+        continue;
       }
 
       Float_t ptCutOffCorrection = 1;
-      if (correction && ptCut > 0)
-        ptCutOffCorrection = correction->GetMeasuredFraction(correctionType, ptCut, vtxVsEta->GetYaxis()->GetBinCenter(iEta));
 
-      if (ptCutOffCorrection <= 0)
+      // find pt cut off correction factor
+      if (fAnalysisMode != AliPWG0Helper::kSPD)
       {
-        printf("UNEXPECTED: ptCutOffCorrection is %f for hist %d %d %d\n", ptCutOffCorrection, vertexPos, vertexBinBegin, vertexBinEnd);
-        continue;
+        if (correction && ptCut > 0)
+            ptCutOffCorrection = correction->GetMeasuredFraction(correctionType, ptCut, vtxVsEta->GetYaxis()->GetBinCenter(iEta), vertexBinBegin, vertexBinEnd);
+
+        if (ptCutOffCorrection <= 0)
+        {
+            printf("UNEXPECTED: ptCutOffCorrection is %f for hist %d %d %d\n", ptCutOffCorrection, vertexPos, vertexBinBegin, vertexBinEnd);
+            continue;
+        }
       }
 
-      //printf("Eta: %d Vertex Range: %d %d, Event Count %f, Track Sum: %f, Track Sum corrected: %f\n", iEta, vertexBinBegin, vertexBinEnd, totalEvents, sum, sum / ptCutOffCorrection);
+      //printf("Eta: %d (%f) Vertex Range: %d %d, Event Count %f, Track Sum: %f, Track Sum corrected: %f \n", iEta, vtxVsEta->GetYaxis()->GetBinCenter(iEta), vertexBinBegin, vertexBinEnd, totalEvents, sum, sum / ptCutOffCorrection);
 
       Int_t bin = fdNdEta[vertexPos]->FindBin(vtxVsEta->GetYaxis()->GetBinCenter(iEta));
       if (bin > 0 && bin <= fdNdEta[vertexPos]->GetNbinsX())
@@ -486,7 +533,7 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
         fdNdEtaPtCutOffCorrected[vertexPos]->SetBinContent(bin, dndeta);
         fdNdEtaPtCutOffCorrected[vertexPos]->SetBinError(bin, error);
 
-        //Printf("Bin %d has dN/deta = %f +- %f; %f +- %f", bin, dndeta, error, fdNdEta[vertexPos]->GetBinContent(bin), fdNdEta[vertexPos]->GetBinError(bin));
+        //Printf("Bin %d has dN/deta = %f +- %f; %.2f tracks %.2f events (outside acceptance: %.2f tracks, %.2f events)", bin, dndeta, error, sum, totalEvents, unusedTracks, unusedEvents);
       }
     }
   }
@@ -592,11 +639,11 @@ void dNdEtaAnalysis::LoadHistograms(const Char_t* dir)
 
   fPtDist = dynamic_cast<TH1F*> (gDirectory->Get(fPtDist->GetName()));
 
-  if (dynamic_cast<TNamed*> (gFile->Get("fTag")))
-    fTag = (dynamic_cast<TNamed*> (gFile->Get("fTag")))->GetTitle();
+  if (dynamic_cast<TNamed*> (gDirectory->Get("fTag")))
+    fTag = (dynamic_cast<TNamed*> (gDirectory->Get("fTag")))->GetTitle();
 
-  if (dynamic_cast<TParameter<Int_t>*> (gFile->Get("fAnalysisMode")))
-    fAnalysisMode = (AliPWG0Helper::AnalysisMode) (dynamic_cast<TParameter<Int_t>*> (gFile->Get("fAnalysisMode")))->GetVal();
+  if (dynamic_cast<TParameter<Int_t>*> (gDirectory->Get("fAnalysisMode")))
+    fAnalysisMode = (AliPWG0Helper::AnalysisMode) (dynamic_cast<TParameter<Int_t>*> (gDirectory->Get("fAnalysisMode")))->GetVal();
 
   gDirectory->cd("../");
 }