During simulation: fill STU region w/ non null time sums
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibBase.cxx
index fa33ab7..e263bbd 100644 (file)
 #include "TGraph.h"
 #include "TGraphErrors.h"
 #include "TF1.h"
+#include "TH1.h"
+#include "THnSparse.h"
 #include "TH1D.h"
 #include "TH2D.h"
-#include "THnSparse.h"
+#include "TAxis.h"
+#include "AliRecoParam.h"
 
 
 #include "AliLog.h"
@@ -72,6 +75,10 @@ AliTPCcalibBase::AliTPCcalibBase():
     fHasLaser(kFALSE),                    //flag the laser is overlayed with given event 
     fRejectLaser(kTRUE),                 //flag- reject laser
     fTriggerClass(),
+    fCurrentEvent(0),           //! current event
+    fCurrentTrack(0),           //! current esd track
+    fCurrentFriendTrack(0),           //! current esd track
+    fCurrentSeed(0),            //! current seed
     fDebugLevel(0)
 {
   //
@@ -93,6 +100,10 @@ AliTPCcalibBase::AliTPCcalibBase(const char * name, const char * title):
   fHasLaser(kFALSE),                    //flag the laser is overlayed with given event 
   fRejectLaser(kTRUE),                 //flag- reject laser
   fTriggerClass(),
+  fCurrentEvent(0),           //! current event
+  fCurrentTrack(0),           //! current esd track
+  fCurrentFriendTrack(0),           //! current esd track
+  fCurrentSeed(0),            //! current seed
   fDebugLevel(0)
 {
   //
@@ -114,6 +125,10 @@ AliTPCcalibBase::AliTPCcalibBase(const AliTPCcalibBase&calib):
   fHasLaser(calib.fHasLaser),                    //flag the laser is overlayed with given event
   fRejectLaser(calib.fRejectLaser),                 //flag- reject laser
   fTriggerClass(calib.fTriggerClass),
+  fCurrentEvent(0),           //! current event
+  fCurrentTrack(0),           //! current esd track
+  fCurrentFriendTrack(0),           //! current esd track
+  fCurrentSeed(0),            //! current seed
   fDebugLevel(calib.fDebugLevel)
 {
   //
@@ -123,8 +138,10 @@ AliTPCcalibBase::AliTPCcalibBase(const AliTPCcalibBase&calib):
 
 AliTPCcalibBase &AliTPCcalibBase::operator=(const AliTPCcalibBase&calib){
   //
+  // operator=
   //
-  //
+  if (this == &calib) return (*this);
+
   ((TNamed *)this)->operator=(calib);
   fDebugStreamer=0;
   fStreamLevel=calib.fStreamLevel;
@@ -187,26 +204,13 @@ Bool_t AliTPCcalibBase::HasLaser(AliESDEvent *event){
   //
   //
   //
-  // Thresholds more than 8 tracks with small dip angle
-  
-  const Int_t kMinLaserTracks = 8;
-  const Float_t kThrLaser       = 0.3;
-  const Float_t kLaserTgl       = 0.01;
-
-  Int_t ntracks = event->GetNumberOfTracks();
-  if (ntracks<kMinLaserTracks) return kFALSE;
-  Float_t nlaser=0;
-  Float_t nall=0;
-  for (Int_t i=0;i<ntracks;++i) {
-    AliESDtrack *track=event->GetTrack(i);
-    if (!track) continue;
-    if (track->GetTPCNcls()<=0) continue; 
-    nall++;
-    if (TMath::Abs(track->GetTgl())<kLaserTgl) nlaser++;
+  // Use event specie
+  Bool_t isLaser=kFALSE;
+  UInt_t specie = event->GetEventSpecie();  // select only cosmic events
+  if (specie==AliRecoParam::kCalib) {
+    isLaser = kTRUE;
   }
-  if (nlaser>kMinLaserTracks) return kTRUE;
-  if (nall>0 && nlaser/nall>kThrLaser) return kTRUE;
-  return kFALSE;
+  return isLaser;
 }
 
 
@@ -249,10 +253,11 @@ void AliTPCcalibBase::RegisterDebugOutput(const char *path){
 
 
 
-TGraphErrors * AliTPCcalibBase::FitSlices(THnSparse *h, Int_t axisDim1, Int_t axisDim2, Int_t minEntries, Int_t nmaxBin, Float_t fracLow, Float_t fracUp){
+TGraphErrors * AliTPCcalibBase::FitSlices(THnSparse *h, Int_t axisDim1, Int_t axisDim2, Int_t minEntries, Int_t nmaxBin, Float_t fracLow, Float_t fracUp, Bool_t useMedian, TTreeSRedirector *cstream, Int_t ival){
   //
   // Fitting slices of the projection(axisDim1,axisDim2) of a sparse histogram
   // 
+  TF1 funcGaus("funcGaus","gaus");
   TH2D * hist = h->Projection(axisDim1, axisDim2);
   Double_t *xvec = new Double_t[hist->GetNbinsX()];
   Double_t *yvec = new Double_t[hist->GetNbinsX()];
@@ -262,7 +267,7 @@ TGraphErrors * AliTPCcalibBase::FitSlices(THnSparse *h, Int_t axisDim1, Int_t ax
   TH1D * projectionHist =0;
   //
 
-  for(Int_t i=1; i < hist->GetNbinsX(); i++) {
+  for(Int_t i=1; i <= hist->GetNbinsX(); i++) {
     Int_t nsum=0;
     Int_t imin   =  i;
     Int_t imax   =  i;
@@ -283,15 +288,19 @@ TGraphErrors * AliTPCcalibBase::FitSlices(THnSparse *h, Int_t axisDim1, Int_t ax
     Float_t xMin = projectionHist->GetXaxis()->GetXmin();
     Float_t xMax = projectionHist->GetXaxis()->GetXmax();
     Float_t xMedian = (xMin+xMax)*0.5;
-    Float_t integral = projectionHist->GetSum();
+    Float_t integral = 0;
+    for(Int_t jbin=1; jbin<projectionHist->GetNbinsX()-1; jbin++) {
+      integral+=projectionHist->GetBinContent(jbin);
+    }
+    //printf("Integral %f\t%f\n",integral, projectionHist->GetSum());
     //
     //
     Float_t currentSum=0;
     for(Int_t jbin=1; jbin<projectionHist->GetNbinsX()-1; jbin++) {
       currentSum += projectionHist->GetBinContent(jbin);
-      if (currentSum/integral<fracLow) xMin = projectionHist->GetBinCenter(jbin);
-      if (currentSum/integral<fracUp)  xMax = projectionHist->GetBinCenter(jbin+1);      
-      if (currentSum/integral<0.5 && projectionHist->GetBinContent(jbin)>0){
+      if (currentSum<fracLow*integral) xMin = projectionHist->GetBinCenter(jbin);
+      if (currentSum<fracUp*integral)  xMax = projectionHist->GetBinCenter(jbin+1);      
+      if (currentSum<0.5*integral && projectionHist->GetBinContent(jbin)>0){
        xMedian = (projectionHist->GetBinCenter(jbin)*projectionHist->GetBinContent(jbin)+
                   projectionHist->GetBinCenter(jbin+1)*projectionHist->GetBinContent(jbin+1))/
          (projectionHist->GetBinContent(jbin)+projectionHist->GetBinContent(jbin+1));
@@ -305,14 +314,26 @@ TGraphErrors * AliTPCcalibBase::FitSlices(THnSparse *h, Int_t axisDim1, Int_t ax
     Double_t xrms    =  hist->GetRMS()+hist->GetXaxis()->GetBinWidth(1)/TMath::Sqrt(12.); 
     Double_t binWidth = projectionHist->GetXaxis()->GetBinWidth(1);
     if (rms>0){
-      TF1 funcGaus("funcGaus","gaus");
       // cut on +- 4 RMS
       projectionHist->Fit(&funcGaus,"QN","",xMin, xMax);
+      Double_t chi2 = funcGaus.GetChisquare();
       //  
       xvec[counter] = xcenter;
-      yvec[counter] = funcGaus.GetParameter(1);
+      yvec[counter] = funcGaus.GetParameter(ival);
       xerr[counter] = xrms;
-      yerr[counter] = funcGaus.GetParError(1); 
+      yerr[counter] = funcGaus.GetParError(ival); 
+      if (useMedian) yvec[counter] = xMedian;
+      if (cstream){
+       (*cstream)<<"fitDebug"<<
+         "xcenter="<<xcenter<<
+         "xMin="<<xMin<<
+         "xMax="<<xMax<<
+         "xMedian="<<xMedian<<
+         "xFitM"<<yvec[counter]<<
+         "xFitS"<<yerr[counter]<<
+         "chi2="<<chi2<<         
+         "\n";
+      }
       counter++;
     }else{
       xvec[counter] = xcenter;
@@ -332,3 +353,47 @@ TGraphErrors * AliTPCcalibBase::FitSlices(THnSparse *h, Int_t axisDim1, Int_t ax
   delete hist;
   return graphErrors;
 }
+
+
+void AliTPCcalibBase::BinLogX(THnSparse *h, Int_t axisDim) {
+
+  // Method for the correct logarithmic binning of histograms
+
+  TAxis *axis = h->GetAxis(axisDim);
+  int bins = axis->GetNbins();
+
+  Double_t from = axis->GetXmin();
+  Double_t to = axis->GetXmax();
+  Double_t *new_bins = new Double_t[bins + 1];
+
+  new_bins[0] = from;
+  Double_t factor = pow(to/from, 1./bins);
+
+  for (int i = 1; i <= bins; i++) {
+   new_bins[i] = factor * new_bins[i-1];
+  }
+  axis->Set(bins, new_bins);
+  delete [] new_bins;
+
+}
+void AliTPCcalibBase::BinLogX(TH1 *h) {
+
+  // Method for the correct logarithmic binning of histograms
+
+  TAxis *axis = h->GetXaxis();
+  int bins = axis->GetNbins();
+
+  Double_t from = axis->GetXmin();
+  Double_t to = axis->GetXmax();
+  Double_t *new_bins = new Double_t[bins + 1];
+
+  new_bins[0] = from;
+  Double_t factor = pow(to/from, 1./bins);
+
+  for (int i = 1; i <= bins; i++) {
+   new_bins[i] = factor * new_bins[i-1];
+  }
+  axis->Set(bins, new_bins);
+  delete [] new_bins;
+
+}