]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TPC/AliTPCcalibTimeGain.cxx
removes OCDB access during merging by splitting out the functionality
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibTimeGain.cxx
index 21b19cb944e2b50d3813b7809e696d266e1f34ec..bfd1bd739fe0afe3aa0c4aaa25d515d3ee52b162 100644 (file)
@@ -164,6 +164,7 @@ TGaxis *axis = new TGaxis(xmax,ymin,xmax,ymax,ymin,ymax,50510,"+L");
 #include "AliTracker.h"
 #include "AliMagF.h"
 #include "AliTPCCalROC.h"
+#include "AliESDv0.h"
 
 #include "AliLog.h"
 
@@ -197,7 +198,10 @@ AliTPCcalibTimeGain::AliTPCcalibTimeGain()
    fIsCosmic(0),
    fLowMemoryConsumption(0)
 {  
-  AliInfo("Default Constructor");  
+  //
+  // Default constructor
+  //
+  AliDebug(5,"Default Constructor");  
 }
 
 
@@ -218,23 +222,25 @@ AliTPCcalibTimeGain::AliTPCcalibTimeGain(const Text_t *name, const Text_t *title
    fIsCosmic(0),
    fLowMemoryConsumption(0)
 {
-  
+  //
+  // No default constructor
+  //
   SetName(name);
   SetTitle(title);
   
-  AliInfo("Non Default Constructor");
+  AliDebug(5,"Non Default Constructor");
   
   fIntegrationTimeDeDx = deltaIntegrationTimeGain;
   
   Double_t deltaTime = EndTime - StartTime;
   
   
-  // main histogram for time dependence: dE/dx, time, type (1-muon cosmic,2-pion beam data), meanDriftlength, momenta (only filled if enough space is available), run number
+  // main histogram for time dependence: dE/dx, time, type (1-muon cosmic,2-pion beam data, 3&4 - proton points at higher dE/dx), meanDriftlength, momenta (only filled if enough space is available), run number, eta
   Int_t timeBins = TMath::Nint(deltaTime/deltaIntegrationTimeGain);
-  Int_t binsGainTime[6]    = {150,  timeBins,    2,  25, 200, 10000000};
-  Double_t xminGainTime[6] = {0.5, StartTime,  0.5,   0, 0.1,    -0.5};
-  Double_t xmaxGainTime[6] = {  8,   EndTime,  2.5, 250,  50, 9999999.5};
-  fHistGainTime = new THnSparseF("HistGainTime","dEdx time dep.;dEdx,time,type,driftlength,momenta,run number;dEdx",6,binsGainTime,xminGainTime,xmaxGainTime);
+  Int_t binsGainTime[7]    = {150,  timeBins,    4,  25, 200, 10000000, 20};
+  Double_t xminGainTime[7] = {0.5, StartTime,  0.5,   0, 0.1,    -0.5,  -1};
+  Double_t xmaxGainTime[7] = {  8,   EndTime,  4.5, 250,  50, 9999999.5, 1};
+  fHistGainTime = new THnSparseF("HistGainTime","dEdx time dep.;dEdx,time,type,driftlength,momenta,run number, eta;dEdx",7,binsGainTime,xminGainTime,xmaxGainTime);
   BinLogX(fHistGainTime, 4);
   //
   fHistDeDxTotal = new TH2F("DeDx","dEdx; momentum p (GeV); TPC signal (a.u.)",250,0.01,100.,1000,0.,8);
@@ -243,8 +249,8 @@ AliTPCcalibTimeGain::AliTPCcalibTimeGain(const Text_t *name, const Text_t *title
   // default values for dE/dx
   fMIP = 50.;
   fUseMax = kTRUE;
-  fLowerTrunc = 0.0;
-  fUpperTrunc = 0.7;
+  fLowerTrunc = 0.02;
+  fUpperTrunc = 0.6;
   fUseShapeNorm = kTRUE;
   fUsePosNorm = kFALSE;
   fUsePadNorm = kFALSE;
@@ -260,8 +266,11 @@ AliTPCcalibTimeGain::AliTPCcalibTimeGain(const Text_t *name, const Text_t *title
 
 AliTPCcalibTimeGain::~AliTPCcalibTimeGain(){
   //
+  // Destructor
   //
-  //
+  delete fHistGainTime;
+  delete fGainVsTime;
+  delete fHistDeDxTotal;
 }
 
 
@@ -273,7 +282,12 @@ void AliTPCcalibTimeGain::Process(AliESDEvent *event) {
     Printf("ERROR: ESD not available");
     return;
   }
-  
+  AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
+  if (!ESDfriend) {
+   return;
+  }
+  if (ESDfriend->TestSkipBit()) return;
+
   if (fIsCosmic) { // this should be removed at some point based on trigger mask !?
     ProcessCosmicEvent(event);
   } else {
@@ -287,9 +301,11 @@ void AliTPCcalibTimeGain::Process(AliESDEvent *event) {
 
 
 void AliTPCcalibTimeGain::ProcessCosmicEvent(AliESDEvent *event) {
-
-  AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
-  if (!ESDfriend) {
+  //
+  // Process in case of cosmic event
+  //
+  AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
+  if (!esdFriend) {
    Printf("ERROR: ESDfriend not available");
    return;
   }
@@ -304,25 +320,24 @@ void AliTPCcalibTimeGain::ProcessCosmicEvent(AliESDEvent *event) {
 
     AliESDtrack *track = event->GetTrack(i);
     if (!track) continue;
-        
+    AliESDfriendTrack *friendTrack = esdFriend->GetTrack(i);
+    if (!friendTrack) continue;        
     const AliExternalTrackParam * trackIn = track->GetInnerParam();
-    const AliExternalTrackParam * trackOut = track->GetOuterParam();
+    const AliExternalTrackParam * trackOut = friendTrack->GetTPCOut();
     if (!trackIn) continue;
     if (!trackOut) continue;
 
     // calculate necessary track parameters
     Double_t meanP = trackIn->GetP();
     Double_t meanDrift = 250 - 0.5*TMath::Abs(trackIn->GetZ() + trackOut->GetZ());
-    Int_t NclsDeDx = track->GetTPCNcls();
+    Int_t nclsDeDx = track->GetTPCNcls();
 
     // exclude tracks which do not look like primaries or are simply too short or on wrong sectors
-    if (NclsDeDx < 60) continue;     
+    if (nclsDeDx < 60) continue;     
     if (TMath::Abs(trackIn->GetTgl()) > 1) continue;
     if (TMath::Abs(trackIn->GetSnp()) > 0.6) continue;
     
     // Get seeds
-    AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
-    if (!friendTrack) continue;
     TObject *calibObject;
     AliTPCseed *seed = 0;
     for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
@@ -330,15 +345,15 @@ void AliTPCcalibTimeGain::ProcessCosmicEvent(AliESDEvent *event) {
     }    
 
     if (seed) { 
-      Double_t TPCsignal = GetTPCdEdx(seed);
-      if (NclsDeDx > 100) fHistDeDxTotal->Fill(meanP, TPCsignal);
+      Double_t tpcSignal = GetTPCdEdx(seed);
+      if (nclsDeDx > 100) fHistDeDxTotal->Fill(meanP, tpcSignal);
       //
       if (fLowMemoryConsumption) {
        if (meanP < 20) continue;
        meanP = 30; // set all momenta to one in order to save memory
       }
       //dE/dx, time, type (1-muon cosmic,2-pion beam data), momenta
-      Double_t vec[6] = {TPCsignal,time,1,meanDrift,meanP,runNumber};
+      Double_t vec[6] = {tpcSignal,time,1,meanDrift,meanP,runNumber};
       fHistGainTime->Fill(vec);
     }
     
@@ -349,9 +364,11 @@ void AliTPCcalibTimeGain::ProcessCosmicEvent(AliESDEvent *event) {
 
 
 void AliTPCcalibTimeGain::ProcessBeamEvent(AliESDEvent *event) {
-
-  AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
-  if (!ESDfriend) {
+  //
+  // Process in case of beam event
+  //
+  AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
+  if (!esdFriend) {
    Printf("ERROR: ESDfriend not available");
    return;
   }
@@ -362,55 +379,116 @@ void AliTPCcalibTimeGain::ProcessBeamEvent(AliESDEvent *event) {
   //
   // track loop
   //
-  for (Int_t i=0;i<ntracks;++i) {
+  for (Int_t i=0;i<ntracks;++i) { // begin track loop
 
     AliESDtrack *track = event->GetTrack(i);
     if (!track) continue;
+    AliESDfriendTrack *friendTrack = esdFriend->GetTrack(i);
+    if (!friendTrack) continue;
         
     const AliExternalTrackParam * trackIn = track->GetInnerParam();
-    const AliExternalTrackParam * trackOut = track->GetOuterParam();
+    const AliExternalTrackParam * trackOut = friendTrack->GetTPCOut();
     if (!trackIn) continue;
     if (!trackOut) continue;
 
     // calculate necessary track parameters
     Double_t meanP = trackIn->GetP();
     Double_t meanDrift = 250 - 0.5*TMath::Abs(trackIn->GetZ() + trackOut->GetZ());
-    Int_t NclsDeDx = track->GetTPCNcls();
+    Int_t nclsDeDx = track->GetTPCNcls();
 
     // exclude tracks which do not look like primaries or are simply too short or on wrong sectors
-    if (NclsDeDx < 60) continue;     
-    if (TMath::Abs(trackIn->GetTgl()) > 1) continue;
-    if (TMath::Abs(trackIn->GetSnp()) > 0.6) continue;
+    if (nclsDeDx < 60) continue;     
+    //if (TMath::Abs(trackIn->GetTgl()) > 1) continue;
+    //if (TMath::Abs(trackIn->GetSnp()) > 0.6) continue;
+    if (TMath::Abs(trackIn->Eta()) > 1) continue;
+    UInt_t status = track->GetStatus();
+    if ((status&AliESDtrack::kTPCrefit)==0) continue;
+    //if (track->GetNcls(0) < 3) continue; // ITS clusters
+    Float_t dca[2], cov[3];
+    track->GetImpactParameters(dca,cov);
+    if (TMath::Abs(dca[0]) > 7 || TMath::Abs(dca[0]) < 0.0000001 || TMath::Abs(dca[1]) > 25 ) continue; // cut in xy
+    Double_t eta = trackIn->Eta();
     
     // Get seeds
-    AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
-    if (!friendTrack) continue;
     TObject *calibObject;
     AliTPCseed *seed = 0;
     for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
       if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
     }    
 
-    if (seed) { 
-      Double_t TPCsignal = GetTPCdEdx(seed);
-      fHistDeDxTotal->Fill(meanP, TPCsignal);
+    if (seed) {
+      Int_t particleCase = 0;
+      if (meanP < 0.5  && meanP > 0.4)  particleCase = 2; // MIP pions
+      if (meanP < 0.57 && meanP > 0.56) particleCase = 3; // protons 1
+      if (meanP < 0.66 && meanP > 0.65) particleCase = 4; // protons 2
       //
-      if (fLowMemoryConsumption) {
-       if (meanP > 0.5 || meanP < 0.4) continue;
-       meanP = 0.45; // set all momenta to one in order to save memory
-      }
-      //dE/dx, time, type (1-muon cosmic,2-pion beam data), momenta
-      Double_t vec[6] = {TPCsignal,time,2,meanDrift,meanP,runNumber};
+      if (fLowMemoryConsumption && particleCase == 0) continue;
+      //
+      Double_t tpcSignal = GetTPCdEdx(seed);
+      fHistDeDxTotal->Fill(meanP, tpcSignal);
+      //
+      //dE/dx, time, type (1-muon cosmic,2-pion beam data, 3&4 protons), momenta, runNumner, eta
+      Double_t vec[7] = {tpcSignal,time,particleCase,meanDrift,meanP,runNumber, eta};
       fHistGainTime->Fill(vec);
     }
     
+  } // end track loop
+  //
+  // V0 loop -- in beam events the cosmic part of the histogram is filled with GammaConversions
+  //
+  for(Int_t iv0 = 0; iv0 < event->GetNumberOfV0s(); iv0++) {
+    AliESDv0 * v0 = event->GetV0(iv0);
+    if (!v0->GetOnFlyStatus()) continue;
+    if (v0->GetEffMass(0,0) > 0.02) continue; // select low inv. mass
+    Double_t xyz[3];
+    v0->GetXYZ(xyz[0], xyz[1], xyz[2]);
+    if (TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]) < 3 || TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]) > 30) continue;
+    //
+    // "loop over daughters" 
+    //
+    for(Int_t idaughter = 0; idaughter < 2; idaughter++) { // daughter loop
+      Int_t index = idaughter == 0 ? v0->GetPindex() : v0->GetNindex();
+      AliESDtrack * trackP = event->GetTrack(index);
+      AliESDfriendTrack *friendTrackP = esdFriend->GetTrack(index);
+      if (!friendTrackP) continue;
+      const AliExternalTrackParam * trackPIn = trackP->GetInnerParam();
+      const AliExternalTrackParam * trackPOut = friendTrackP->GetTPCOut();
+      if (!trackPIn) continue;
+      if (!trackPOut) continue;
+      // calculate necessary track parameters
+      Double_t meanP = trackPIn->GetP();
+      Double_t meanDrift = 250 - 0.5*TMath::Abs(trackPIn->GetZ() + trackPOut->GetZ());
+      Int_t nclsDeDx = trackP->GetTPCNcls();
+      // exclude tracks which do not look like primaries or are simply too short or on wrong sectors
+      if (nclsDeDx < 60) continue;     
+      if (TMath::Abs(trackPIn->GetTgl()) > 1) continue;
+      //
+      TObject *calibObject;
+      AliTPCseed *seed = 0;
+      for (Int_t l=0;(calibObject=friendTrackP->GetCalibObject(l));++l) {
+      if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
+      }    
+      if (seed) { 
+       if (fLowMemoryConsumption) {
+         if (meanP > 0.5 || meanP < 0.4) continue;
+         meanP = 0.45; // set all momenta to one in order to save memory
+      }
+       Double_t tpcSignal = GetTPCdEdx(seed);
+       //dE/dx, time, type (1-muon cosmic,2-pion beam data), momenta
+       Double_t vec[6] = {tpcSignal,time,1,meanDrift,meanP,runNumber};
+       fHistGainTime->Fill(vec);
+      }
+    }
+    
   }
 
 }
 
 
 Float_t AliTPCcalibTimeGain::GetTPCdEdx(AliTPCseed * seed) {
-
+  //
+  // calculate tpc dEdx
+  //
   Double_t signal = 0;
   //
   if (!fUseCookAnalytical) {
@@ -425,7 +503,7 @@ Float_t AliTPCcalibTimeGain::GetTPCdEdx(AliTPCseed * seed) {
 
 void AliTPCcalibTimeGain::AnalyzeRun(Int_t minEntries) {
   //
-  //
+  // Analyze results of calibration
   //
   if (fIsCosmic) {
     fHistGainTime->GetAxis(2)->SetRangeUser(0.51,1.49); // only cosmics
@@ -443,7 +521,7 @@ void AliTPCcalibTimeGain::AnalyzeRun(Int_t minEntries) {
 
 TGraphErrors * AliTPCcalibTimeGain::GetGraphGainVsTime(Int_t runNumber, Int_t minEntries) {
   //
-  //
+  // Analyze results and get the graph 
   //
   if (runNumber == 0) {
     if (!fGainVsTime) {
@@ -459,7 +537,9 @@ TGraphErrors * AliTPCcalibTimeGain::GetGraphGainVsTime(Int_t runNumber, Int_t mi
 }
 
 Long64_t AliTPCcalibTimeGain::Merge(TCollection *li) {
-
+  //
+  // merge component
+  //
   TIterator* iter = li->MakeIterator();
   AliTPCcalibTimeGain* cal = 0;
 
@@ -474,7 +554,7 @@ Long64_t AliTPCcalibTimeGain::Merge(TCollection *li) {
     if (cal->GetHistDeDxTotal()) fHistDeDxTotal->Add(cal->GetHistDeDxTotal());
 
   }
-  
+  delete iter;
   return 0;
   
 }
@@ -482,7 +562,7 @@ Long64_t AliTPCcalibTimeGain::Merge(TCollection *li) {
 
 AliSplineFit * AliTPCcalibTimeGain::MakeSplineFit(TGraphErrors * graph) {
   //
-  //
+  // make spline fit of gain
   //
   AliSplineFit *fit = new AliSplineFit();
   fit->SetGraph(graph);
@@ -495,7 +575,7 @@ AliSplineFit * AliTPCcalibTimeGain::MakeSplineFit(TGraphErrors * graph) {
 
 
 
-TGraphErrors * AliTPCcalibTimeGain::GetGraphAttachment(Int_t minEntries, Int_t nmaxBin, Float_t fracLow, Float_t fracUp) {
+TGraphErrors * AliTPCcalibTimeGain::GetGraphAttachment(Int_t minEntries, Int_t nmaxBin, Float_t /*fracLow*/, Float_t /*fracUp*/) {
   //
   // For each time bin the driftlength-dependence of the signal is fitted.
   //
@@ -551,53 +631,56 @@ TGraphErrors * AliTPCcalibTimeGain::GetGraphAttachment(Int_t minEntries, Int_t n
 
 
 void AliTPCcalibTimeGain::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];
+  Double_t *newBins = new Double_t[bins + 1];
    
-  new_bins[0] = from;
+  newBins[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];
+   newBins[i] = factor * newBins[i-1];
   }
-  axis->Set(bins, new_bins);
-  delete new_bins;
+  axis->Set(bins, newBins);
+  delete[] newBins;
   
 }
 
 
 void AliTPCcalibTimeGain::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];
+  Double_t *newBins = new Double_t[bins + 1];
    
-  new_bins[0] = from;
+  newBins[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];
+   newBins[i] = factor * newBins[i-1];
   }
-  axis->Set(bins, new_bins);
-  delete new_bins;
+  axis->Set(bins, newBins);
+  delete[] newBins;
   
 }
 
 
 
 void AliTPCcalibTimeGain::CalculateBetheAlephParams(TH2F *hist, Double_t * ini) {
+  //
+  // Fit the bethe bloch params
+  //
   //{0.0762*MIP,10.632,1.34e-05,1.863,1.948}
   const Double_t sigma = 0.06;
 
@@ -615,7 +698,7 @@ void AliTPCcalibTimeGain::CalculateBetheAlephParams(TH2F *hist, Double_t * ini)
   foKaon->SetParameters(ini);
   foProton->SetParameters(ini);
   
-  TCanvas *CanvCheck1 = new TCanvas();
+  TCanvas *canvCheck1 = new TCanvas();
   hist->Draw("colz");
   foElectron->Draw("same");
   foMuon->Draw("same");
@@ -678,7 +761,7 @@ void AliTPCcalibTimeGain::CalculateBetheAlephParams(TH2F *hist, Double_t * ini)
   TF1 * funcAlephD = new TF1("AlephParametrizationD", "AliExternalTrackParam::BetheBlochAleph(x,[0],[1],[2],[3],[4])",0.3,10000);
   funcAlephD->SetParameters(ini);
 
-  TCanvas *CanvCheck2 = new TCanvas();
+  TCanvas *canvCheck2 = new TCanvas();
   histBG->Draw();
   
   //FitSlices
@@ -701,7 +784,7 @@ void AliTPCcalibTimeGain::CalculateBetheAlephParams(TH2F *hist, Double_t * ini)
   foKaon->SetParameters(ini);
   foProton->SetParameters(ini);
   
-  TCanvas *CanvCheck3 = new TCanvas();
+  TCanvas *canvCheck3 = new TCanvas();
   hist->Draw("colz");
   foElectron->Draw("same");
   foMuon->Draw("same");
@@ -709,9 +792,9 @@ void AliTPCcalibTimeGain::CalculateBetheAlephParams(TH2F *hist, Double_t * ini)
   foKaon->Draw("same");  
   foProton->Draw("same");
   
-  CanvCheck1->Print();
-  CanvCheck2->Print();
-  CanvCheck3->Print();
+  canvCheck1->Print();
+  canvCheck2->Print();
+  canvCheck3->Print();
 
   return;