[u/mrichter/AliRoot.git] / TPC / AliTPCcalibGainMult.cxx


/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/*


Basic calibration and QA class for the TPC gain calibration based on tracks from BEAM events.


Send comments etc. to: A.Kalweit@gsi.de, marian.ivanov@cern.ch
*/


#include "Riostream.h"
#include "TChain.h"
#include "TTree.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TList.h"
#include "TMath.h"
#include "TCanvas.h"
#include "TFile.h"
#include "TF1.h"
#include "TVectorD.h"
#include "TProfile.h"

#include "AliTPCcalibDB.h"
#include "AliTPCclusterMI.h"
#include "AliTPCClusterParam.h"
#include "AliTPCseed.h"
#include "AliESDVertex.h"
#include "AliESDEvent.h"
#include "AliESDfriend.h"
#include "AliESDInputHandler.h"
#include "AliAnalysisManager.h"
#include "AliTPCParam.h"

#include "AliComplexCluster.h"
#include "AliTPCclusterMI.h"

#include "AliLog.h"

#include "AliTPCcalibGainMult.h"

#include "TTreeStream.h"


ClassImp(AliTPCcalibGainMult)


AliTPCcalibGainMult::AliTPCcalibGainMult() 
  :AliTPCcalibBase(),
   fMIP(0),
   fLowerTrunc(0),
   fUpperTrunc(0),
   fUseMax(kFALSE),
   fHistNTracks(0),
   fHistClusterShape(0),
   fHistQA(0),
   fHistGainSector(0),
   fHistPadEqual(0),
   fHistGainMult(0)
{  
  //
  // Empty default cosntructor
  //
  AliInfo("Default Constructor");  
}


AliTPCcalibGainMult::AliTPCcalibGainMult(const Text_t *name, const Text_t *title) 
  :AliTPCcalibBase(),
   fMIP(0),
   fLowerTrunc(0),
   fUpperTrunc(0),
   fUseMax(kFALSE),
   fHistNTracks(0),
   fHistClusterShape(0),
   fHistQA(0),
   fHistGainSector(0),
   fHistPadEqual(0),
   fHistGainMult(0)
{
  //
  //
  //  
  SetName(name);
  SetTitle(title);
  //
  fMIP = 50.;
  fLowerTrunc = 0.02; // IMPORTANT CHANGE --> REMOVE HARDWIRED TRUNCATION FROM TRACKER
  fUpperTrunc = 0.6;
  fUseMax = kTRUE; // IMPORTANT CHANGE FOR PbPb; standard: kFALSE;
  //
  fHistNTracks = new TH1F("ntracks","Number of Tracks per Event; number of tracks per event; number of tracks",1001,-0.5,1000.5);
  fHistClusterShape = new TH1F("fHistClusterShape","cluster shape; rms meas. / rms exp.;",300,0,3);
  fHistQA = new TH3F("fHistQA","dEdx; momentum p (GeV); TPC signal (a.u.); pad",500,0.1,20.,500,0.,500,6,-0.5,5.5);
  BinLogX(fHistQA);
  //
  //
  //                          MIP, sect,  pad (short,med,long,full,oroc),   run,      ncl
  Int_t binsGainSec[5]    = { 145,   72,    4,  10000000,   65};
  Double_t xminGainSec[5] = { 10., -0.5, -0.5,      -0.5, -0.5}; 
  Double_t xmaxGainSec[5] = {300., 71.5,  3.5, 9999999.5, 64.5};
  TString axisNameSec[5]={"Q","sector","pad type","run", "ncl"};
  TString axisTitleSec[5]={"Q (a.u)","sector","pad type","run","ncl"};
  //
  fHistGainSector = new THnSparseF("fHistGainSector","0:MIP, 1:sect, 2:pad, 3:run, 4:ncl", 5, binsGainSec, xminGainSec, xmaxGainSec);
  for (Int_t iaxis=0; iaxis<5;iaxis++){
    fHistGainSector->GetAxis(iaxis)->SetName(axisNameSec[iaxis]);
    fHistGainSector->GetAxis(iaxis)->SetTitle(axisTitleSec[iaxis]);
  }
  //
  //
  //
  Int_t binsPadEqual[6]    = { 200, 200,    4,   20,   50, 100};
  Double_t xminPadEqual[6] = { 0.5, 0.5, -0.5,    0, -250,   0}; 
  Double_t xmaxPadEqual[6] = { 1.5, 1.5,  3.5, 4000,  250,   3};
  TString axisNamePadEqual[6]   = {"dEdxRatioMax","dEdxRatioTot","padType","mult","driftlength", "1_pt"};
  TString axisTitlePadEqual[6]  = {"dEdx_padRegion/mean_dEdx Qmax", "dEdx_padRegion/mean_dEdx Qtot","padType","mult","driftlength", "1/pt"};
  //
  fHistPadEqual = new THnSparseF("fHistPadEqual","0:dEdx_pad/dEdx_mean, 1:pad, 2:mult, 3:drift, 4:1/pt", 6, binsPadEqual, xminPadEqual, xmaxPadEqual);
  for (Int_t iaxis=0; iaxis<6;iaxis++){
    fHistPadEqual->GetAxis(iaxis)->SetName(axisNamePadEqual[iaxis]);
    fHistPadEqual->GetAxis(iaxis)->SetTitle(axisTitlePadEqual[iaxis]);
  }
  //
  //
  //                    MIP Qmax, MIP Qtot,  z,  pad, vtx. contribut., ncl
  Int_t binsGainMult[6]    = { 145,  145,   25,    4,  100,  80};
  Double_t xminGainMult[6] = { 10.,  10.,    0, -0.5,    0, -0.5}; 
  Double_t xmaxGainMult[6] = {300., 300.,  250,  3.5, 5000, 159.5};
  TString axisNameMult[6]={"Qmax","Qtot","drift","padtype""multiplicity","ncl"};
  TString axisTitleMult[6]={"Qmax (a.u)","Qtot (a.u.)","driftlenght l (cm)","Pad Type","multiplicity","ncl"};
  //
  fHistGainMult = new THnSparseF("fHistGainMult","MIP Qmax, MIP Qtot, z, type, vtx. contribut.", 6, binsGainMult, xminGainMult, xmaxGainMult); 
  for (Int_t iaxis=0; iaxis<6;iaxis++){
    fHistGainMult->GetAxis(iaxis)->SetName(axisNameMult[iaxis]);
    fHistGainMult->GetAxis(iaxis)->SetTitle(axisTitleMult[iaxis]);
  }
  //
  AliInfo("Non Default Constructor");  
  //
}


AliTPCcalibGainMult::~AliTPCcalibGainMult(){
  //
  //
  //
  delete fHistNTracks;            //  histogram showing number of ESD tracks per event
  delete fHistClusterShape;       //  histogram to check the cluster shape
  delete fHistQA;                 //  dE/dx histogram showing the final spectrum
  //
  delete fHistGainSector;   //  histogram which shows MIP peak for each of the 3x36 sectors (pad region)
  delete fHistPadEqual;     //  histogram for the equalization of the gain in the different pad regions -> pass0
  delete fHistGainMult;     //  histogram which shows decrease of MIP signal as a function


}


void AliTPCcalibGainMult::Process(AliESDEvent *event) {
  //
  //
  //
  if (!event) {
    Printf("ERROR: ESD not available");
    return;
  }  
  Int_t ntracks=event->GetNumberOfTracks(); 
  fHistNTracks->Fill(ntracks);
  
  AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
  if (!esdFriend) {
   Printf("ERROR: esdFriend not available");
   return;
  }
  UInt_t runNumber = event->GetRunNumber();
  Int_t nContributors = 0;
  if (event->GetPrimaryVertexTPC()) nContributors = event->GetPrimaryVertexTPC()->GetNContributors();
  //
  // track loop
  //
  for (Int_t i=0;i<ntracks;++i) {
    //
    AliESDtrack *track = event->GetTrack(i);
    if (!track) continue;
    //   
    AliExternalTrackParam * trackIn  = (AliExternalTrackParam *)track->GetInnerParam();
    if (!trackIn) continue;
  
    // calculate necessary track parameters
    Double_t meanP = trackIn->GetP();
    Int_t ncls = track->GetTPCNcls();

    if (ncls < 80) continue;     
    
    // exclude tracks which do not look like primaries or are simply too short or on wrong sectors

    if (TMath::Abs(trackIn->Eta()) > 0.8) 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);
    Float_t primVtxDCA = TMath::Sqrt(dca[0]*dca[0]);
    if (primVtxDCA > 10 || primVtxDCA < 0.00001) continue;
    if (TMath::Abs(dca[1]) > 5) continue;
    //
    // require that the track does not cross any dead area
    //
    //if (track->GetTPCNclsF() < 158) continue;
    //
    //if (seed->CookShape(1) > 1) continue;
    //if (TMath::Abs(trackIn->GetY()) > 20) continue;
    //if (TMath::Abs(d)>20) continue;   // distance to the 0,0; select only tracks which cross chambers under proper angle
    //if (TMath::Abs(trackIn->GetSnp()) > 0.6) continue;
    if (primVtxDCA < 3 && track->GetNcls(0) > 3 && track->GetKinkIndex(0) == 0 && ncls > 100) fHistQA->Fill(meanP, track->GetTPCsignal(), 5);

    // 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) {
      //
      const AliExternalTrackParam * trackOut = friendTrack->GetTPCOut();
      if (!trackIn) continue;
      if (!trackOut) continue;
      Double_t meanDrift = 250 - 0.5*TMath::Abs(trackIn->GetZ() + trackOut->GetZ());
      //
      for (Int_t irow =0; irow<160;irow++)    {
	AliTPCTrackerPoint * point = seed->GetTrackPoint(irow);
	if (point==0) continue;
	AliTPCclusterMI * cl = seed->GetClusterPointer(irow);
	if (cl==0) continue;	
	//
	Float_t rsigmay =  TMath::Sqrt(point->GetSigmaY());
	fHistClusterShape->Fill(rsigmay);
      }
      //
      Int_t row0 = 0;
      Int_t row1 = 160;
      //
      Double_t signalShortMax = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,1,0,62);
      Double_t signalMedMax   = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,1,63,126);
      Double_t signalLongMax  = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,1,127,159);
      Double_t signalMax      = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,1,row0,row1);
      Double_t signalArrayMax[4] = {signalShortMax, signalMedMax, signalLongMax, signalMax};
      //
      Double_t signalShortTot = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,0,0,62);
      Double_t signalMedTot   = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,0,63,126);
      Double_t signalLongTot  = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,0,127,159);
      Double_t signalTot      = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,0,row0,row1);
      Double_t signalArrayTot[4] = {signalShortTot, signalMedTot, signalLongTot, signalTot};
      //
      Double_t mipSignalShort = fUseMax ? signalShortMax : signalShortTot;
      Double_t mipSignalMed   = fUseMax ? signalMedMax   : signalMedTot;
      Double_t mipSignalLong  = fUseMax ? signalLongMax  : signalLongTot;
      Double_t mipSignalOroc  = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,fUseMax,63,159);
      Double_t signal =  fUseMax ? signalMax  : signalTot;
      //
      fHistQA->Fill(meanP, mipSignalShort, 0);
      fHistQA->Fill(meanP, mipSignalMed, 1);
      fHistQA->Fill(meanP, mipSignalLong, 2);
      fHistQA->Fill(meanP, signal, 3);
      fHistQA->Fill(meanP, mipSignalOroc, 4);
      //
      // "dEdxRatioMax","dEdxRatioTot","padType","mult","driftlength", "1_pt"
      Float_t meanMax = (1/3.)*(signalArrayMax[0] + signalArrayMax[1] + signalArrayMax[2]);
      Float_t meanTot = (1/3.)*(signalArrayTot[0] + signalArrayTot[1] + signalArrayTot[2]); 
      if (meanMax < 1e-5 || meanTot < 1e-5) continue;
      for(Int_t ipad = 0; ipad < 4; ipad ++) {
	Double_t vecPadEqual[6] = {signalArrayMax[ipad]/meanMax, signalArrayTot[ipad]/meanTot, ipad, nContributors, meanDrift, track->OneOverPt()};
	fHistPadEqual->Fill(vecPadEqual);
      }
      //
      if (meanP > 0.4 && meanP < 0.55) {
	Double_t vecMult[6] = {seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,1,row0,row1),
			       seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,0,row0,row1),
			       meanDrift,
			       3,
			       nContributors,
			       ncls};
	//
	fHistGainMult->Fill(vecMult);
	vecMult[0]=mipSignalShort; vecMult[1]=mipSignalShort; vecMult[3]=0;
	fHistGainMult->Fill(vecMult);
	vecMult[0]=mipSignalMed; vecMult[1]=mipSignalMed; vecMult[3]=1;
	fHistGainMult->Fill(vecMult);
	vecMult[0]=mipSignalLong; vecMult[1]=mipSignalLong; vecMult[3]=2;
	fHistGainMult->Fill(vecMult);
	//
      }
      //
      //
      if (meanP > 0.5 || meanP < 0.4) continue; // only MIP pions
      //
      // for each track, we look at the three different pad regions, split it into tracklets, check if the sector does not change and fill the histogram
      //
      Bool_t isNotSplit[3] = {kTRUE, kTRUE, kTRUE}; //  short, medium, long (true if the track is not split between two chambers)
      //
      Double_t sector[4] = {-1, -1, -1, -1}; // sector number short, medium, long, all
      Int_t ncl[3] = {0,0,0};
      //
      for (Int_t irow=0; irow < 159; irow++){
	Int_t padRegion = 0;
	if (irow > 62) padRegion = 1;
	if (irow > 126) padRegion = 2;
	//
	AliTPCclusterMI* cluster = seed->GetClusterPointer(irow);
	if (!cluster) continue;
	if (sector[padRegion] == -1) {
	  sector[padRegion] = cluster->GetDetector();
	  continue;
	}
	if (sector[padRegion] != -1 && sector[padRegion] != cluster->GetDetector()) isNotSplit[padRegion] = kFALSE;
	ncl[padRegion]++;
      }
      //
      //                        MIP, sect,  pad,   run
      //
      Double_t vecMip[5] = {mipSignalShort, mipSignalMed, mipSignalLong, signal, mipSignalOroc};
      //
      for(Int_t ipad = 0; ipad < 3; ipad++) {
	//
	Double_t vecGainSec[5] = {vecMip[ipad], sector[ipad], ipad, runNumber, ncl[ipad]};
	if (isNotSplit[ipad]) fHistGainSector->Fill(vecGainSec);
      }
    }
   
  }    
}  


void AliTPCcalibGainMult::MakeLookup(THnSparse * /*hist*/, Char_t * /*outputFile*/) {
  //
  // Not  yet implemented
  //
}


void AliTPCcalibGainMult::Analyze() {


  return;

}


Long64_t AliTPCcalibGainMult::Merge(TCollection *li) {

  TIterator* iter = li->MakeIterator();
  AliTPCcalibGainMult* cal = 0;

  while ((cal = (AliTPCcalibGainMult*)iter->Next())) {
    if (!cal->InheritsFrom(AliTPCcalibGainMult::Class())) {
      Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
      return -1;
    }
    
    if (cal->GetHistNTracks()) fHistNTracks->Add(cal->GetHistNTracks());
    if (cal->GetHistClusterShape()) fHistClusterShape->Add(cal->GetHistClusterShape());
    if (cal->GetHistQA()) fHistQA->Add(cal->GetHistQA());
    if (cal->GetHistGainSector()) fHistGainSector->Add(cal->GetHistGainSector());
    if (cal->GetHistPadEqual()) fHistPadEqual->Add(cal->GetHistPadEqual());
    if (cal->GetHistGainMult()) fHistGainMult->Add(cal->GetHistGainMult());
 
  }
  
  return 0;
  
}


void AliTPCcalibGainMult::BinLogX(const 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 *newBins = new Double_t[bins + 1];
   
  newBins[0] = from;
  Double_t factor = pow(to/from, 1./bins);
  
  for (int i = 1; i <= bins; i++) {
   newBins[i] = factor * newBins[i-1];
  }
  axis->Set(bins, newBins);
  delete [] newBins;
  
  
}


void AliTPCcalibGainMult::UpdateGainMap() {
  //
  // read in the old gain map and scale it appropriately...
  //
  /*
  gSystem->Load("libANALYSIS");
  gSystem->Load("libTPCcalib");
  //
  TFile jj("Run0_999999999_v1_s0.root");
  AliTPCCalPad * pad = AliCDBEntry->GetObject()->Clone();
  TFile hh("output.root");
  AliTPCcalibGainMult * gain = calibTracksGain;
  TH2D * histGainSec = gain->GetHistGainSector()->Projection(0,1);
  //
  TObjArray arr;
  histGainSec->FitSlicesY(0, 0, -1, 0, "QNR", &arr);
  TH1D * meanGainSec = arr->At(1);
  Double_t gainsIROC[36];
  Double_t gainsOROC[36];
  Double_t gains[72];
  //
  for(Int_t isec = 1; isec < meanGainSec->GetNbinsX() + 1; isec++) {
    cout << isec << " " << meanGainSec->GetXaxis()->GetBinCenter(isec) << " " <<meanGainSec->GetBinContent(isec) << endl;
    gains[isec-1] = meanGainSec->GetBinContent(isec);
    if (isec < 37) {
      gainsIROC[isec-1] = meanGainSec->GetBinContent(isec);
    } else {
      gainsOROC[isec - 36 -1] = meanGainSec->GetBinContent(isec);
    }
  }
  Double_t meanIroc = TMath::Mean(36, gainsIROC);
  Double_t meanOroc = TMath::Mean(36, gainsIROC);
  for(Int_t i = 0; i < 36; i++) gains[i] /= meanIroc;
  for(Int_t i = 36; i < 72; i++) gains[i] /= meanOroc;
  //
  for(Int_t i = 0; i < 72; i++) {
    AliTPCCalROC * chamber = pad->GetCalROC(i);
    chamber->Multiply(gains[i]);
    cout << i << " "<< chamber->GetMean() << endl;
  }
  //
  // update the OCDB
  //
  AliCDBMetaData *metaData= new AliCDBMetaData();
  metaData->SetObjectClassName("AliTPCCalPad");
  metaData->SetResponsible("Alexander Kalweit");
  metaData->SetBeamPeriod(1);
  metaData->SetAliRootVersion("04-19-05"); //root version
  metaData->SetComment("New gain map for 1600V OROC gain increase and equalization. Valid for runs starting after technical stop beginning of September.");
  AliCDBId id1("TPC/Calib/GainFactorDedx", 131541, AliCDBRunRange::Infinity()); // important: new gain runs here..
  AliCDBStorage * gStorage = AliCDBManager::Instance()->GetStorage("local:///d/alice05/akalweit/projects/OCDBupdate/HighGain_2010-09-03/OCDB/");
  gStorage->Put(pad, id1, metaData);
  */
  
}

void AliTPCcalibGainMult::UpdateClusterParam() {
  //
  //
  //
  /*
  gSystem->Load("libANALYSIS");
  gSystem->Load("libTPCcalib");
  //
  TFile ff("OldClsParam.root");
  AliTPCClusterParam * param = AliCDBEntry->GetObject()->Clone();
 
  TFile hh("output.root");
  AliTPCcalibGainMult * gain = calibGainMult;
  TH2D * histGainSec = gain->GetHistGainSector()->Projection(0,2);
  TObjArray arr;
  histGainSec->FitSlicesY(0, 0, -1, 0, "QNR", &arr);
  histGainSec->Draw("colz");
  TH1D * fitVal = arr.At(1);
  fitVal->Draw("same");
  param->GetQnormCorrMatrix()->Print();
  param->GetQnormCorrMatrix()(0,5) *= fitVal->GetBinContent(1)/fitVal->GetBinContent(1); // short pads Qtot
  param->GetQnormCorrMatrix()(1,5) *= fitVal->GetBinContent(2)/fitVal->GetBinContent(1); // med pads Qtot
  param->GetQnormCorrMatrix()(2,5) *= fitVal->GetBinContent(3)/fitVal->GetBinContent(1); // long pads Qtot
  //
  param->GetQnormCorrMatrix()(3,5) *= fitVal->GetBinContent(1)/fitVal->GetBinContent(1); // short pads Qmax -> scaling assumed
  param->GetQnormCorrMatrix()(4,5) *= fitVal->GetBinContent(2)/fitVal->GetBinContent(1); // med pads Qmax -> scaling assumed
  param->GetQnormCorrMatrix()(5,5) *= fitVal->GetBinContent(3)/fitVal->GetBinContent(1); // long pads Qmax -> scaling assumed
  //
  TFile jj("newClusterParam.root","RECREATE");
  param->Write();
  param->GetQnormCorrMatrix()->Print();
  //
  // update the OCDB
  // 
  AliCDBMetaData *metaData= new AliCDBMetaData();
  metaData->SetObjectClassName("AliTPCClusterParam");
  metaData->SetResponsible("Alexander Kalweit");
  metaData->SetBeamPeriod(1);
  metaData->SetAliRootVersion("04-19-04"); //root version
  metaData->SetComment("1600V OROC / hard thres. / new algorithm");
  AliCDBId id1("TPC/Calib/ClusterParam", 0, AliCDBRunRange::Infinity()); // important: new gain runs here..
  AliCDBStorage * gStorage = AliCDBManager::Instance()->GetStorage("local:///lustre/alice/akalweit/baseline/CalibrationEntries/OldThres_NewAlgo_PP");
  gStorage->Put(param, id1, metaData);
  */
  

}
Commit	Line	Data
f72219cb	1
	2	/**************************************************************************
	3	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	4	* *
	5	* Author: The ALICE Off-line Project. *
	6	* Contributors are mentioned in the code where appropriate. *
	7	* *
	8	* Permission to use, copy, modify and distribute this software and its *
	9	* documentation strictly for non-commercial purposes is hereby granted *
	10	* without fee, provided that the above copyright notice appears in all *
	11	* copies and that both the copyright notice and this permission notice *
	12	* appear in the supporting documentation. The authors make no claims *
	13	* about the suitability of this software for any purpose. It is *
	14	* provided "as is" without express or implied warranty. *
	15	**************************************************************************/
	16
	17	/*
	18
	19
	20	Basic calibration and QA class for the TPC gain calibration based on tracks from BEAM events.
	21
	22
	23	Send comments etc. to: A.Kalweit@gsi.de, marian.ivanov@cern.ch
	24	*/
	25
	26
	27	#include "Riostream.h"
	28	#include "TChain.h"
	29	#include "TTree.h"
	30	#include "TH1F.h"
	31	#include "TH2F.h"
	32	#include "TList.h"
	33	#include "TMath.h"
	34	#include "TCanvas.h"
	35	#include "TFile.h"
	36	#include "TF1.h"
	37	#include "TVectorD.h"
	38	#include "TProfile.h"
	39
	40	#include "AliTPCcalibDB.h"
	41	#include "AliTPCclusterMI.h"
	42	#include "AliTPCClusterParam.h"
	43	#include "AliTPCseed.h"
	44	#include "AliESDVertex.h"
	45	#include "AliESDEvent.h"
	46	#include "AliESDfriend.h"
	47	#include "AliESDInputHandler.h"
	48	#include "AliAnalysisManager.h"
	49	#include "AliTPCParam.h"
	50
	51	#include "AliComplexCluster.h"
	52	#include "AliTPCclusterMI.h"
	53
	54	#include "AliLog.h"
	55
	56	#include "AliTPCcalibGainMult.h"
	57
	58	#include "TTreeStream.h"
	59
	60
	61	ClassImp(AliTPCcalibGainMult)
	62
	63
	64	AliTPCcalibGainMult::AliTPCcalibGainMult()
65	:AliTPCcalibBase(),
66	fMIP(0),
67	fLowerTrunc(0),
68	fUpperTrunc(0),
69	fUseMax(kFALSE),
70	fHistNTracks(0),
71	fHistClusterShape(0),
72	fHistQA(0),
73	fHistGainSector(0),
74	fHistPadEqual(0),
75	fHistGainMult(0)
76	{
77	//
78	// Empty default cosntructor
79	//
80	AliInfo("Default Constructor");
81	}
82
83
84	AliTPCcalibGainMult::AliTPCcalibGainMult(const Text_t name, const Text_t title)
85	:AliTPCcalibBase(),
86	fMIP(0),
87	fLowerTrunc(0),
88	fUpperTrunc(0),
89	fUseMax(kFALSE),
90	fHistNTracks(0),
91	fHistClusterShape(0),
92	fHistQA(0),
93	fHistGainSector(0),
94	fHistPadEqual(0),
95	fHistGainMult(0)
96	{
97	//
98	//
99	//
100	SetName(name);
101	SetTitle(title);
102	//
103	fMIP = 50.;
104	fLowerTrunc = 0.02; // IMPORTANT CHANGE --> REMOVE HARDWIRED TRUNCATION FROM TRACKER
105	fUpperTrunc = 0.6;
106	fUseMax = kTRUE; // IMPORTANT CHANGE FOR PbPb; standard: kFALSE;
107	//
108	fHistNTracks = new TH1F("ntracks","Number of Tracks per Event; number of tracks per event; number of tracks",1001,-0.5,1000.5);
109	fHistClusterShape = new TH1F("fHistClusterShape","cluster shape; rms meas. / rms exp.;",300,0,3);
110	fHistQA = new TH3F("fHistQA","dEdx; momentum p (GeV); TPC signal (a.u.); pad",500,0.1,20.,500,0.,500,6,-0.5,5.5);
111	BinLogX(fHistQA);
112	//
113	//
114	// MIP, sect, pad (short,med,long,full,oroc), run, ncl
115	Int_t binsGainSec[5] = { 145, 72, 4, 10000000, 65};
116	Double_t xminGainSec[5] = { 10., -0.5, -0.5, -0.5, -0.5};
117	Double_t xmaxGainSec[5] = {300., 71.5, 3.5, 9999999.5, 64.5};
118	TString axisNameSec[5]={"Q","sector","pad type","run", "ncl"};
119	TString axisTitleSec[5]={"Q (a.u)","sector","pad type","run","ncl"};
120	//
121	fHistGainSector = new THnSparseF("fHistGainSector","0:MIP, 1:sect, 2:pad, 3:run, 4:ncl", 5, binsGainSec, xminGainSec, xmaxGainSec);
122	for (Int_t iaxis=0; iaxis<5;iaxis++){
123	fHistGainSector->GetAxis(iaxis)->SetName(axisNameSec[iaxis]);
124	fHistGainSector->GetAxis(iaxis)->SetTitle(axisTitleSec[iaxis]);
125	}
126	//
127	//
128	//
129	Int_t binsPadEqual[6] = { 200, 200, 4, 20, 50, 100};
130	Double_t xminPadEqual[6] = { 0.5, 0.5, -0.5, 0, -250, 0};
131	Double_t xmaxPadEqual[6] = { 1.5, 1.5, 3.5, 4000, 250, 3};
132	TString axisNamePadEqual[6] = {"dEdxRatioMax","dEdxRatioTot","padType","mult","driftlength", "1_pt"};
133	TString axisTitlePadEqual[6] = {"dEdx_padRegion/mean_dEdx Qmax", "dEdx_padRegion/mean_dEdx Qtot","padType","mult","driftlength", "1/pt"};
134	//
135	fHistPadEqual = new THnSparseF("fHistPadEqual","0:dEdx_pad/dEdx_mean, 1:pad, 2:mult, 3:drift, 4:1/pt", 6, binsPadEqual, xminPadEqual, xmaxPadEqual);
136	for (Int_t iaxis=0; iaxis<6;iaxis++){
137	fHistPadEqual->GetAxis(iaxis)->SetName(axisNamePadEqual[iaxis]);
138	fHistPadEqual->GetAxis(iaxis)->SetTitle(axisTitlePadEqual[iaxis]);
139	}
140	//
141	//
142	// MIP Qmax, MIP Qtot, z, pad, vtx. contribut., ncl
143	Int_t binsGainMult[6] = { 145, 145, 25, 4, 100, 80};
144	Double_t xminGainMult[6] = { 10., 10., 0, -0.5, 0, -0.5};
145	Double_t xmaxGainMult[6] = {300., 300., 250, 3.5, 5000, 159.5};
146	TString axisNameMult[6]={"Qmax","Qtot","drift","padtype""multiplicity","ncl"};
147	TString axisTitleMult[6]={"Qmax (a.u)","Qtot (a.u.)","driftlenght l (cm)","Pad Type","multiplicity","ncl"};
148	//
149	fHistGainMult = new THnSparseF("fHistGainMult","MIP Qmax, MIP Qtot, z, type, vtx. contribut.", 6, binsGainMult, xminGainMult, xmaxGainMult);
150	for (Int_t iaxis=0; iaxis<6;iaxis++){
151	fHistGainMult->GetAxis(iaxis)->SetName(axisNameMult[iaxis]);
152	fHistGainMult->GetAxis(iaxis)->SetTitle(axisTitleMult[iaxis]);
153	}
154	//
155	AliInfo("Non Default Constructor");
156	//
157	}
158
159
160	AliTPCcalibGainMult::~AliTPCcalibGainMult(){
161	//
162	//
163	//
b9ab8e40	164	delete fHistNTracks; // histogram showing number of ESD tracks per event
	165	delete fHistClusterShape; // histogram to check the cluster shape
	166	delete fHistQA; // dE/dx histogram showing the final spectrum
	167	//
	168	delete fHistGainSector; // histogram which shows MIP peak for each of the 3x36 sectors (pad region)
	169	delete fHistPadEqual; // histogram for the equalization of the gain in the different pad regions -> pass0
	170	delete fHistGainMult; // histogram which shows decrease of MIP signal as a function
	171
f72219cb	172
	173	}
	174
	175
	176
	177	void AliTPCcalibGainMult::Process(AliESDEvent *event) {
	178	//
	179	//
	180	//
	181	if (!event) {
	182	Printf("ERROR: ESD not available");
	183	return;
	184	}
	185	Int_t ntracks=event->GetNumberOfTracks();
	186	fHistNTracks->Fill(ntracks);
	187
	188	AliESDfriend esdFriend=static_cast<AliESDfriend>(event->FindListObject("AliESDfriend"));
	189	if (!esdFriend) {
	190	Printf("ERROR: esdFriend not available");
	191	return;
	192	}
	193	UInt_t runNumber = event->GetRunNumber();
	194	Int_t nContributors = 0;
	195	if (event->GetPrimaryVertexTPC()) nContributors = event->GetPrimaryVertexTPC()->GetNContributors();
	196	//
	197	// track loop
	198	//
	199	for (Int_t i=0;i<ntracks;++i) {
	200	//
	201	AliESDtrack *track = event->GetTrack(i);
	202	if (!track) continue;
	203	//
	204	AliExternalTrackParam * trackIn = (AliExternalTrackParam *)track->GetInnerParam();
	205	if (!trackIn) continue;
	206
	207	// calculate necessary track parameters
	208	Double_t meanP = trackIn->GetP();
	209	Int_t ncls = track->GetTPCNcls();
	210
	211	if (ncls < 80) continue;
	212
	213	// exclude tracks which do not look like primaries or are simply too short or on wrong sectors
	214
	215	if (TMath::Abs(trackIn->Eta()) > 0.8) continue;
	216	UInt_t status = track->GetStatus();
	217	if ((status&AliESDtrack::kTPCrefit)==0) continue;
	218	//if (track->GetNcls(0) < 3) continue; // ITS clusters
	219	Float_t dca[2], cov[3];
	220	track->GetImpactParameters(dca,cov);
	221	Float_t primVtxDCA = TMath::Sqrt(dca[0]*dca[0]);
	222	if (primVtxDCA > 10 \|\| primVtxDCA < 0.00001) continue;
	223	if (TMath::Abs(dca[1]) > 5) continue;
	224	//
	225	// require that the track does not cross any dead area
	226	//
	227	//if (track->GetTPCNclsF() < 158) continue;
	228	//
	229	//if (seed->CookShape(1) > 1) continue;
	230	//if (TMath::Abs(trackIn->GetY()) > 20) continue;
	231	//if (TMath::Abs(d)>20) continue; // distance to the 0,0; select only tracks which cross chambers under proper angle
	232	//if (TMath::Abs(trackIn->GetSnp()) > 0.6) continue;
	233	if (primVtxDCA < 3 && track->GetNcls(0) > 3 && track->GetKinkIndex(0) == 0 && ncls > 100) fHistQA->Fill(meanP, track->GetTPCsignal(), 5);
	234
	235	// Get seeds
236	AliESDfriendTrack *friendTrack = esdFriend->GetTrack(i);
237	if (!friendTrack) continue;
238	TObject *calibObject;
239	AliTPCseed *seed = 0;
240	for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
241	if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
242	}
243
244	if (seed) {
245	//
246	const AliExternalTrackParam * trackOut = friendTrack->GetTPCOut();
247	if (!trackIn) continue;
248	if (!trackOut) continue;
249	Double_t meanDrift = 250 - 0.5*TMath::Abs(trackIn->GetZ() + trackOut->GetZ());
250	//
251	for (Int_t irow =0; irow<160;irow++) {
252	AliTPCTrackerPoint * point = seed->GetTrackPoint(irow);
253	if (point==0) continue;
254	AliTPCclusterMI * cl = seed->GetClusterPointer(irow);
255	if (cl==0) continue;
256	//
257	Float_t rsigmay = TMath::Sqrt(point->GetSigmaY());
258	fHistClusterShape->Fill(rsigmay);
259	}
260	//
261	Int_t row0 = 0;
262	Int_t row1 = 160;
263	//
264	Double_t signalShortMax = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,1,0,62);
265	Double_t signalMedMax = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,1,63,126);
266	Double_t signalLongMax = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,1,127,159);
267	Double_t signalMax = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,1,row0,row1);
268	Double_t signalArrayMax[4] = {signalShortMax, signalMedMax, signalLongMax, signalMax};
269	//
270	Double_t signalShortTot = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,0,0,62);
271	Double_t signalMedTot = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,0,63,126);
272	Double_t signalLongTot = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,0,127,159);
273	Double_t signalTot = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,0,row0,row1);
274	Double_t signalArrayTot[4] = {signalShortTot, signalMedTot, signalLongTot, signalTot};
275	//
276	Double_t mipSignalShort = fUseMax ? signalShortMax : signalShortTot;
277	Double_t mipSignalMed = fUseMax ? signalMedMax : signalMedTot;
278	Double_t mipSignalLong = fUseMax ? signalLongMax : signalLongTot;
279	Double_t mipSignalOroc = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,fUseMax,63,159);
280	Double_t signal = fUseMax ? signalMax : signalTot;
281	//
282	fHistQA->Fill(meanP, mipSignalShort, 0);
283	fHistQA->Fill(meanP, mipSignalMed, 1);
284	fHistQA->Fill(meanP, mipSignalLong, 2);
285	fHistQA->Fill(meanP, signal, 3);
286	fHistQA->Fill(meanP, mipSignalOroc, 4);
287	//
288	// "dEdxRatioMax","dEdxRatioTot","padType","mult","driftlength", "1_pt"
289	Float_t meanMax = (1/3.)*(signalArrayMax[0] + signalArrayMax[1] + signalArrayMax[2]);
b4ffdb90	290	Float_t meanTot = (1/3.)*(signalArrayTot[0] + signalArrayTot[1] + signalArrayTot[2]);
4a4c6be2	291	if (meanMax < 1e-5 \|\| meanTot < 1e-5) continue;
f72219cb	292	for(Int_t ipad = 0; ipad < 4; ipad ++) {
	293	Double_t vecPadEqual[6] = {signalArrayMax[ipad]/meanMax, signalArrayTot[ipad]/meanTot, ipad, nContributors, meanDrift, track->OneOverPt()};
	294	fHistPadEqual->Fill(vecPadEqual);
	295	}
	296	//
	297	if (meanP > 0.4 && meanP < 0.55) {
	298	Double_t vecMult[6] = {seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,1,row0,row1),
	299	seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,0,row0,row1),
	300	meanDrift,
	301	3,
	302	nContributors,
	303	ncls};
	304	//
	305	fHistGainMult->Fill(vecMult);
	306	vecMult[0]=mipSignalShort; vecMult[1]=mipSignalShort; vecMult[3]=0;
	307	fHistGainMult->Fill(vecMult);
	308	vecMult[0]=mipSignalMed; vecMult[1]=mipSignalMed; vecMult[3]=1;
	309	fHistGainMult->Fill(vecMult);
	310	vecMult[0]=mipSignalLong; vecMult[1]=mipSignalLong; vecMult[3]=2;
	311	fHistGainMult->Fill(vecMult);
	312	//
	313	}
	314	//
	315	//
	316	if (meanP > 0.5 \|\| meanP < 0.4) continue; // only MIP pions
	317	//
	318	// for each track, we look at the three different pad regions, split it into tracklets, check if the sector does not change and fill the histogram
	319	//
	320	Bool_t isNotSplit[3] = {kTRUE, kTRUE, kTRUE}; // short, medium, long (true if the track is not split between two chambers)
	321	//
	322	Double_t sector[4] = {-1, -1, -1, -1}; // sector number short, medium, long, all
	323	Int_t ncl[3] = {0,0,0};
	324	//
	325	for (Int_t irow=0; irow < 159; irow++){
	326	Int_t padRegion = 0;
	327	if (irow > 62) padRegion = 1;
	328	if (irow > 126) padRegion = 2;
	329	//
	330	AliTPCclusterMI* cluster = seed->GetClusterPointer(irow);
	331	if (!cluster) continue;
	332	if (sector[padRegion] == -1) {
	333	sector[padRegion] = cluster->GetDetector();
	334	continue;
	335	}
	336	if (sector[padRegion] != -1 && sector[padRegion] != cluster->GetDetector()) isNotSplit[padRegion] = kFALSE;
	337	ncl[padRegion]++;
	338	}
	339	//
	340	// MIP, sect, pad, run
	341	//
	342	Double_t vecMip[5] = {mipSignalShort, mipSignalMed, mipSignalLong, signal, mipSignalOroc};
	343	//
	344	for(Int_t ipad = 0; ipad < 3; ipad++) {
	345	//
	346	Double_t vecGainSec[5] = {vecMip[ipad], sector[ipad], ipad, runNumber, ncl[ipad]};
	347	if (isNotSplit[ipad]) fHistGainSector->Fill(vecGainSec);
	348	}
	349	}
	350
	351	}
	352	}
	353
	354
a8ef8a9c	355	void AliTPCcalibGainMult::MakeLookup(THnSparse * /hist/, Char_t * /outputFile/) {
	356	//
	357	// Not yet implemented
	358	//
4a4c6be2	359	}
f72219cb	360
	361
	362	void AliTPCcalibGainMult::Analyze() {
	363
	364
	365	return;
	366
	367	}
	368
	369
	370	Long64_t AliTPCcalibGainMult::Merge(TCollection *li) {
	371
	372	TIterator* iter = li->MakeIterator();
	373	AliTPCcalibGainMult* cal = 0;
	374
	375	while ((cal = (AliTPCcalibGainMult*)iter->Next())) {
	376	if (!cal->InheritsFrom(AliTPCcalibGainMult::Class())) {
	377	Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
	378	return -1;
	379	}
	380
	381	if (cal->GetHistNTracks()) fHistNTracks->Add(cal->GetHistNTracks());
	382	if (cal->GetHistClusterShape()) fHistClusterShape->Add(cal->GetHistClusterShape());
	383	if (cal->GetHistQA()) fHistQA->Add(cal->GetHistQA());
	384	if (cal->GetHistGainSector()) fHistGainSector->Add(cal->GetHistGainSector());
	385	if (cal->GetHistPadEqual()) fHistPadEqual->Add(cal->GetHistPadEqual());
	386	if (cal->GetHistGainMult()) fHistGainMult->Add(cal->GetHistGainMult());
	387
	388	}
	389
	390	return 0;
	391
	392	}
	393
	394
	395
	396	void AliTPCcalibGainMult::BinLogX(const TH1 *h) {
	397
	398	// Method for the correct logarithmic binning of histograms
	399
	400	TAxis *axis = h->GetXaxis();
	401	int bins = axis->GetNbins();
	402
	403	Double_t from = axis->GetXmin();
	404	Double_t to = axis->GetXmax();
	405	Double_t *newBins = new Double_t[bins + 1];
	406
	407	newBins[0] = from;
	408	Double_t factor = pow(to/from, 1./bins);
	409
	410	for (int i = 1; i <= bins; i++) {
	411	newBins[i] = factor * newBins[i-1];
	412	}
	413	axis->Set(bins, newBins);
	414	delete [] newBins;
	415
	416
	417	}
	418
	419
	420	void AliTPCcalibGainMult::UpdateGainMap() {
	421	//
	422	// read in the old gain map and scale it appropriately...
	423	//
424	/*
425	gSystem->Load("libANALYSIS");
426	gSystem->Load("libTPCcalib");
427	//
428	TFile jj("Run0_999999999_v1_s0.root");
429	AliTPCCalPad * pad = AliCDBEntry->GetObject()->Clone();
430	TFile hh("output.root");
431	AliTPCcalibGainMult * gain = calibTracksGain;
432	TH2D * histGainSec = gain->GetHistGainSector()->Projection(0,1);
433	//
434	TObjArray arr;
435	histGainSec->FitSlicesY(0, 0, -1, 0, "QNR", &arr);
436	TH1D * meanGainSec = arr->At(1);
437	Double_t gainsIROC[36];
438	Double_t gainsOROC[36];
439	Double_t gains[72];
440	//
441	for(Int_t isec = 1; isec < meanGainSec->GetNbinsX() + 1; isec++) {
442	cout << isec << " " << meanGainSec->GetXaxis()->GetBinCenter(isec) << " " <<meanGainSec->GetBinContent(isec) << endl;
443	gains[isec-1] = meanGainSec->GetBinContent(isec);
444	if (isec < 37) {
445	gainsIROC[isec-1] = meanGainSec->GetBinContent(isec);
446	} else {
447	gainsOROC[isec - 36 -1] = meanGainSec->GetBinContent(isec);
448	}
449	}
450	Double_t meanIroc = TMath::Mean(36, gainsIROC);
451	Double_t meanOroc = TMath::Mean(36, gainsIROC);
452	for(Int_t i = 0; i < 36; i++) gains[i] /= meanIroc;
453	for(Int_t i = 36; i < 72; i++) gains[i] /= meanOroc;
454	//
455	for(Int_t i = 0; i < 72; i++) {
456	AliTPCCalROC * chamber = pad->GetCalROC(i);
457	chamber->Multiply(gains[i]);
458	cout << i << " "<< chamber->GetMean() << endl;
459	}
460	//
461	// update the OCDB
462	//
463	AliCDBMetaData *metaData= new AliCDBMetaData();
464	metaData->SetObjectClassName("AliTPCCalPad");
465	metaData->SetResponsible("Alexander Kalweit");
466	metaData->SetBeamPeriod(1);
467	metaData->SetAliRootVersion("04-19-05"); //root version
468	metaData->SetComment("New gain map for 1600V OROC gain increase and equalization. Valid for runs starting after technical stop beginning of September.");
469	AliCDBId id1("TPC/Calib/GainFactorDedx", 131541, AliCDBRunRange::Infinity()); // important: new gain runs here..
470	AliCDBStorage * gStorage = AliCDBManager::Instance()->GetStorage("local:///d/alice05/akalweit/projects/OCDBupdate/HighGain_2010-09-03/OCDB/");
471	gStorage->Put(pad, id1, metaData);
472	*/
473
474	}
475
476	void AliTPCcalibGainMult::UpdateClusterParam() {
477	//
478	//
479	//
480	/*
481	gSystem->Load("libANALYSIS");
482	gSystem->Load("libTPCcalib");
483	//
484	TFile ff("OldClsParam.root");
485	AliTPCClusterParam * param = AliCDBEntry->GetObject()->Clone();
486
487	TFile hh("output.root");
488	AliTPCcalibGainMult * gain = calibGainMult;
489	TH2D * histGainSec = gain->GetHistGainSector()->Projection(0,2);
490	TObjArray arr;
491	histGainSec->FitSlicesY(0, 0, -1, 0, "QNR", &arr);
492	histGainSec->Draw("colz");
493	TH1D * fitVal = arr.At(1);
494	fitVal->Draw("same");
495	param->GetQnormCorrMatrix()->Print();
496	param->GetQnormCorrMatrix()(0,5) *= fitVal->GetBinContent(1)/fitVal->GetBinContent(1); // short pads Qtot
497	param->GetQnormCorrMatrix()(1,5) *= fitVal->GetBinContent(2)/fitVal->GetBinContent(1); // med pads Qtot
498	param->GetQnormCorrMatrix()(2,5) *= fitVal->GetBinContent(3)/fitVal->GetBinContent(1); // long pads Qtot
499	//
500	param->GetQnormCorrMatrix()(3,5) *= fitVal->GetBinContent(1)/fitVal->GetBinContent(1); // short pads Qmax -> scaling assumed
501	param->GetQnormCorrMatrix()(4,5) *= fitVal->GetBinContent(2)/fitVal->GetBinContent(1); // med pads Qmax -> scaling assumed
502	param->GetQnormCorrMatrix()(5,5) *= fitVal->GetBinContent(3)/fitVal->GetBinContent(1); // long pads Qmax -> scaling assumed
503	//
504	TFile jj("newClusterParam.root","RECREATE");
505	param->Write();
506	param->GetQnormCorrMatrix()->Print();
507	//
508	// update the OCDB
509	//
510	AliCDBMetaData *metaData= new AliCDBMetaData();
511	metaData->SetObjectClassName("AliTPCClusterParam");
512	metaData->SetResponsible("Alexander Kalweit");
513	metaData->SetBeamPeriod(1);
514	metaData->SetAliRootVersion("04-19-04"); //root version
515	metaData->SetComment("1600V OROC / hard thres. / new algorithm");
516	AliCDBId id1("TPC/Calib/ClusterParam", 0, AliCDBRunRange::Infinity()); // important: new gain runs here..
517	AliCDBStorage * gStorage = AliCDBManager::Instance()->GetStorage("local:///lustre/alice/akalweit/baseline/CalibrationEntries/OldThres_NewAlgo_PP");
518	gStorage->Put(param, id1, metaData);
519	*/
520
521
522	}
523