[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, 13000,  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, 13000, 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 = event->GetNumberOfTracks();
  //
  // 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};
012c4694	131	Double_t xmaxPadEqual[6] = { 1.5, 1.5, 3.5, 13000, 250, 3};
f72219cb	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};
012c4694	145	Double_t xmaxGainMult[6] = {300., 300., 250, 3.5, 13000, 159.5};
f72219cb	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();
012c4694	194	Int_t nContributors = event->GetNumberOfTracks();
f72219cb	195	//
	196	// track loop
	197	//
	198	for (Int_t i=0;i<ntracks;++i) {
	199	//
	200	AliESDtrack *track = event->GetTrack(i);
	201	if (!track) continue;
	202	//
	203	AliExternalTrackParam * trackIn = (AliExternalTrackParam *)track->GetInnerParam();
	204	if (!trackIn) continue;
	205
	206	// calculate necessary track parameters
	207	Double_t meanP = trackIn->GetP();
	208	Int_t ncls = track->GetTPCNcls();
	209
	210	if (ncls < 80) continue;
	211
	212	// exclude tracks which do not look like primaries or are simply too short or on wrong sectors
	213
	214	if (TMath::Abs(trackIn->Eta()) > 0.8) continue;
	215	UInt_t status = track->GetStatus();
	216	if ((status&AliESDtrack::kTPCrefit)==0) continue;
	217	//if (track->GetNcls(0) < 3) continue; // ITS clusters
	218	Float_t dca[2], cov[3];
	219	track->GetImpactParameters(dca,cov);
	220	Float_t primVtxDCA = TMath::Sqrt(dca[0]*dca[0]);
	221	if (primVtxDCA > 10 \|\| primVtxDCA < 0.00001) continue;
	222	if (TMath::Abs(dca[1]) > 5) continue;
	223	//
	224	// require that the track does not cross any dead area
	225	//
	226	//if (track->GetTPCNclsF() < 158) continue;
	227	//
	228	//if (seed->CookShape(1) > 1) continue;
	229	//if (TMath::Abs(trackIn->GetY()) > 20) continue;
	230	//if (TMath::Abs(d)>20) continue; // distance to the 0,0; select only tracks which cross chambers under proper angle
	231	//if (TMath::Abs(trackIn->GetSnp()) > 0.6) continue;
	232	if (primVtxDCA < 3 && track->GetNcls(0) > 3 && track->GetKinkIndex(0) == 0 && ncls > 100) fHistQA->Fill(meanP, track->GetTPCsignal(), 5);
	233
	234	// Get seeds
	235	AliESDfriendTrack *friendTrack = esdFriend->GetTrack(i);
	236	if (!friendTrack) continue;
	237	TObject *calibObject;
	238	AliTPCseed *seed = 0;
	239	for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
	240	if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
	241	}
	242
	243	if (seed) {
	244	//
	245	const AliExternalTrackParam * trackOut = friendTrack->GetTPCOut();
	246	if (!trackIn) continue;
	247	if (!trackOut) continue;
	248	Double_t meanDrift = 250 - 0.5*TMath::Abs(trackIn->GetZ() + trackOut->GetZ());
	249	//
	250	for (Int_t irow =0; irow<160;irow++) {
	251	AliTPCTrackerPoint * point = seed->GetTrackPoint(irow);
	252	if (point==0) continue;
	253	AliTPCclusterMI * cl = seed->GetClusterPointer(irow);
	254	if (cl==0) continue;
	255	//
	256	Float_t rsigmay = TMath::Sqrt(point->GetSigmaY());
	257	fHistClusterShape->Fill(rsigmay);
	258	}
259	//
260	Int_t row0 = 0;
261	Int_t row1 = 160;
262	//
263	Double_t signalShortMax = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,1,0,62);
264	Double_t signalMedMax = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,1,63,126);
265	Double_t signalLongMax = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,1,127,159);
266	Double_t signalMax = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,1,row0,row1);
267	Double_t signalArrayMax[4] = {signalShortMax, signalMedMax, signalLongMax, signalMax};
268	//
269	Double_t signalShortTot = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,0,0,62);
270	Double_t signalMedTot = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,0,63,126);
271	Double_t signalLongTot = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,0,127,159);
272	Double_t signalTot = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,0,row0,row1);
273	Double_t signalArrayTot[4] = {signalShortTot, signalMedTot, signalLongTot, signalTot};
274	//
275	Double_t mipSignalShort = fUseMax ? signalShortMax : signalShortTot;
276	Double_t mipSignalMed = fUseMax ? signalMedMax : signalMedTot;
277	Double_t mipSignalLong = fUseMax ? signalLongMax : signalLongTot;
278	Double_t mipSignalOroc = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,fUseMax,63,159);
279	Double_t signal = fUseMax ? signalMax : signalTot;
280	//
281	fHistQA->Fill(meanP, mipSignalShort, 0);
282	fHistQA->Fill(meanP, mipSignalMed, 1);
283	fHistQA->Fill(meanP, mipSignalLong, 2);
284	fHistQA->Fill(meanP, signal, 3);
285	fHistQA->Fill(meanP, mipSignalOroc, 4);
286	//
287	// "dEdxRatioMax","dEdxRatioTot","padType","mult","driftlength", "1_pt"
288	Float_t meanMax = (1/3.)*(signalArrayMax[0] + signalArrayMax[1] + signalArrayMax[2]);
b4ffdb90	289	Float_t meanTot = (1/3.)*(signalArrayTot[0] + signalArrayTot[1] + signalArrayTot[2]);
4a4c6be2	290	if (meanMax < 1e-5 \|\| meanTot < 1e-5) continue;
f72219cb	291	for(Int_t ipad = 0; ipad < 4; ipad ++) {
	292	Double_t vecPadEqual[6] = {signalArrayMax[ipad]/meanMax, signalArrayTot[ipad]/meanTot, ipad, nContributors, meanDrift, track->OneOverPt()};
	293	fHistPadEqual->Fill(vecPadEqual);
	294	}
	295	//
	296	if (meanP > 0.4 && meanP < 0.55) {
	297	Double_t vecMult[6] = {seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,1,row0,row1),
	298	seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,0,row0,row1),
	299	meanDrift,
	300	3,
	301	nContributors,
	302	ncls};
	303	//
	304	fHistGainMult->Fill(vecMult);
	305	vecMult[0]=mipSignalShort; vecMult[1]=mipSignalShort; vecMult[3]=0;
	306	fHistGainMult->Fill(vecMult);
	307	vecMult[0]=mipSignalMed; vecMult[1]=mipSignalMed; vecMult[3]=1;
	308	fHistGainMult->Fill(vecMult);
	309	vecMult[0]=mipSignalLong; vecMult[1]=mipSignalLong; vecMult[3]=2;
	310	fHistGainMult->Fill(vecMult);
	311	//
	312	}
	313	//
	314	//
	315	if (meanP > 0.5 \|\| meanP < 0.4) continue; // only MIP pions
	316	//
	317	// 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
	318	//
	319	Bool_t isNotSplit[3] = {kTRUE, kTRUE, kTRUE}; // short, medium, long (true if the track is not split between two chambers)
	320	//
	321	Double_t sector[4] = {-1, -1, -1, -1}; // sector number short, medium, long, all
	322	Int_t ncl[3] = {0,0,0};
	323	//
	324	for (Int_t irow=0; irow < 159; irow++){
	325	Int_t padRegion = 0;
	326	if (irow > 62) padRegion = 1;
	327	if (irow > 126) padRegion = 2;
	328	//
	329	AliTPCclusterMI* cluster = seed->GetClusterPointer(irow);
	330	if (!cluster) continue;
	331	if (sector[padRegion] == -1) {
	332	sector[padRegion] = cluster->GetDetector();
	333	continue;
	334	}
	335	if (sector[padRegion] != -1 && sector[padRegion] != cluster->GetDetector()) isNotSplit[padRegion] = kFALSE;
	336	ncl[padRegion]++;
	337	}
	338	//
	339	// MIP, sect, pad, run
	340	//
	341	Double_t vecMip[5] = {mipSignalShort, mipSignalMed, mipSignalLong, signal, mipSignalOroc};
	342	//
	343	for(Int_t ipad = 0; ipad < 3; ipad++) {
	344	//
	345	Double_t vecGainSec[5] = {vecMip[ipad], sector[ipad], ipad, runNumber, ncl[ipad]};
	346	if (isNotSplit[ipad]) fHistGainSector->Fill(vecGainSec);
	347	}
	348	}
	349
	350	}
	351	}
	352
	353
a8ef8a9c	354	void AliTPCcalibGainMult::MakeLookup(THnSparse * /hist/, Char_t * /outputFile/) {
	355	//
	356	// Not yet implemented
	357	//
4a4c6be2	358	}
f72219cb	359
	360
	361	void AliTPCcalibGainMult::Analyze() {
	362
	363
	364	return;
	365
	366	}
	367
	368
	369	Long64_t AliTPCcalibGainMult::Merge(TCollection *li) {
	370
	371	TIterator* iter = li->MakeIterator();
	372	AliTPCcalibGainMult* cal = 0;
	373
	374	while ((cal = (AliTPCcalibGainMult*)iter->Next())) {
	375	if (!cal->InheritsFrom(AliTPCcalibGainMult::Class())) {
	376	Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
	377	return -1;
	378	}
	379
	380	if (cal->GetHistNTracks()) fHistNTracks->Add(cal->GetHistNTracks());
	381	if (cal->GetHistClusterShape()) fHistClusterShape->Add(cal->GetHistClusterShape());
	382	if (cal->GetHistQA()) fHistQA->Add(cal->GetHistQA());
	383	if (cal->GetHistGainSector()) fHistGainSector->Add(cal->GetHistGainSector());
	384	if (cal->GetHistPadEqual()) fHistPadEqual->Add(cal->GetHistPadEqual());
	385	if (cal->GetHistGainMult()) fHistGainMult->Add(cal->GetHistGainMult());
	386
	387	}
	388
	389	return 0;
	390
	391	}
	392
	393
	394
	395	void AliTPCcalibGainMult::BinLogX(const TH1 *h) {
	396
	397	// Method for the correct logarithmic binning of histograms
	398
	399	TAxis *axis = h->GetXaxis();
	400	int bins = axis->GetNbins();
	401
	402	Double_t from = axis->GetXmin();
	403	Double_t to = axis->GetXmax();
	404	Double_t *newBins = new Double_t[bins + 1];
	405
	406	newBins[0] = from;
	407	Double_t factor = pow(to/from, 1./bins);
	408
	409	for (int i = 1; i <= bins; i++) {
	410	newBins[i] = factor * newBins[i-1];
	411	}
	412	axis->Set(bins, newBins);
	413	delete [] newBins;
	414
	415
	416	}
	417
	418
	419	void AliTPCcalibGainMult::UpdateGainMap() {
	420	//
	421	// read in the old gain map and scale it appropriately...
	422	//
423	/*
424	gSystem->Load("libANALYSIS");
425	gSystem->Load("libTPCcalib");
426	//
427	TFile jj("Run0_999999999_v1_s0.root");
428	AliTPCCalPad * pad = AliCDBEntry->GetObject()->Clone();
429	TFile hh("output.root");
430	AliTPCcalibGainMult * gain = calibTracksGain;
431	TH2D * histGainSec = gain->GetHistGainSector()->Projection(0,1);
432	//
433	TObjArray arr;
434	histGainSec->FitSlicesY(0, 0, -1, 0, "QNR", &arr);
435	TH1D * meanGainSec = arr->At(1);
436	Double_t gainsIROC[36];
437	Double_t gainsOROC[36];
438	Double_t gains[72];
439	//
440	for(Int_t isec = 1; isec < meanGainSec->GetNbinsX() + 1; isec++) {
441	cout << isec << " " << meanGainSec->GetXaxis()->GetBinCenter(isec) << " " <<meanGainSec->GetBinContent(isec) << endl;
442	gains[isec-1] = meanGainSec->GetBinContent(isec);
443	if (isec < 37) {
444	gainsIROC[isec-1] = meanGainSec->GetBinContent(isec);
445	} else {
446	gainsOROC[isec - 36 -1] = meanGainSec->GetBinContent(isec);
447	}
448	}
449	Double_t meanIroc = TMath::Mean(36, gainsIROC);
450	Double_t meanOroc = TMath::Mean(36, gainsIROC);
451	for(Int_t i = 0; i < 36; i++) gains[i] /= meanIroc;
452	for(Int_t i = 36; i < 72; i++) gains[i] /= meanOroc;
453	//
454	for(Int_t i = 0; i < 72; i++) {
455	AliTPCCalROC * chamber = pad->GetCalROC(i);
456	chamber->Multiply(gains[i]);
457	cout << i << " "<< chamber->GetMean() << endl;
458	}
459	//
460	// update the OCDB
461	//
462	AliCDBMetaData *metaData= new AliCDBMetaData();
463	metaData->SetObjectClassName("AliTPCCalPad");
464	metaData->SetResponsible("Alexander Kalweit");
465	metaData->SetBeamPeriod(1);
466	metaData->SetAliRootVersion("04-19-05"); //root version
467	metaData->SetComment("New gain map for 1600V OROC gain increase and equalization. Valid for runs starting after technical stop beginning of September.");
468	AliCDBId id1("TPC/Calib/GainFactorDedx", 131541, AliCDBRunRange::Infinity()); // important: new gain runs here..
469	AliCDBStorage * gStorage = AliCDBManager::Instance()->GetStorage("local:///d/alice05/akalweit/projects/OCDBupdate/HighGain_2010-09-03/OCDB/");
470	gStorage->Put(pad, id1, metaData);
471	*/
472
473	}
474
475	void AliTPCcalibGainMult::UpdateClusterParam() {
476	//
477	//
478	//
479	/*
480	gSystem->Load("libANALYSIS");
481	gSystem->Load("libTPCcalib");
482	//
483	TFile ff("OldClsParam.root");
484	AliTPCClusterParam * param = AliCDBEntry->GetObject()->Clone();
485
486	TFile hh("output.root");
487	AliTPCcalibGainMult * gain = calibGainMult;
488	TH2D * histGainSec = gain->GetHistGainSector()->Projection(0,2);
489	TObjArray arr;
490	histGainSec->FitSlicesY(0, 0, -1, 0, "QNR", &arr);
491	histGainSec->Draw("colz");
492	TH1D * fitVal = arr.At(1);
493	fitVal->Draw("same");
494	param->GetQnormCorrMatrix()->Print();
495	param->GetQnormCorrMatrix()(0,5) *= fitVal->GetBinContent(1)/fitVal->GetBinContent(1); // short pads Qtot
496	param->GetQnormCorrMatrix()(1,5) *= fitVal->GetBinContent(2)/fitVal->GetBinContent(1); // med pads Qtot
497	param->GetQnormCorrMatrix()(2,5) *= fitVal->GetBinContent(3)/fitVal->GetBinContent(1); // long pads Qtot
498	//
499	param->GetQnormCorrMatrix()(3,5) *= fitVal->GetBinContent(1)/fitVal->GetBinContent(1); // short pads Qmax -> scaling assumed
500	param->GetQnormCorrMatrix()(4,5) *= fitVal->GetBinContent(2)/fitVal->GetBinContent(1); // med pads Qmax -> scaling assumed
501	param->GetQnormCorrMatrix()(5,5) *= fitVal->GetBinContent(3)/fitVal->GetBinContent(1); // long pads Qmax -> scaling assumed
502	//
503	TFile jj("newClusterParam.root","RECREATE");
504	param->Write();
505	param->GetQnormCorrMatrix()->Print();
506	//
507	// update the OCDB
508	//
509	AliCDBMetaData *metaData= new AliCDBMetaData();
510	metaData->SetObjectClassName("AliTPCClusterParam");
511	metaData->SetResponsible("Alexander Kalweit");
512	metaData->SetBeamPeriod(1);
513	metaData->SetAliRootVersion("04-19-04"); //root version
514	metaData->SetComment("1600V OROC / hard thres. / new algorithm");
515	AliCDBId id1("TPC/Calib/ClusterParam", 0, AliCDBRunRange::Infinity()); // important: new gain runs here..
516	AliCDBStorage * gStorage = AliCDBManager::Instance()->GetStorage("local:///lustre/alice/akalweit/baseline/CalibrationEntries/OldThres_NewAlgo_PP");
517	gStorage->Put(param, id1, metaData);
518	*/
519
520
521	}
522