[u/mrichter/AliRoot.git] / TRD / AliTRDdEdxCalibUtils.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.                  *
 **************************************************************************/

//
// class to calculate TRD dEdx
// xx
// xx
// xx
// xx
//
//  Xianguo Lu 
//  lu@physi.uni-heidelberg.de
//  Xianguo.Lu@cern.ch
//  
//


#include "TF1.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "THnSparse.h"
#include "TMath.h"
#include "TMatrixD.h"
#include "TMinuit.h"
#include "TObjArray.h"
#include "TRandom3.h"
#include "TStopwatch.h"
#include "TVectorD.h"

#include "TTreeStream.h"

#include "AliCDBId.h"
#include "AliCDBMetaData.h"
#include "AliCDBStorage.h"
#include "AliESDEvent.h"
#include "AliESDfriendTrack.h"
#include "AliESDtrack.h"
#include "AliTRDcalibDB.h"
#include "AliTRDCalROC.h"
#include "AliTRDtrackV1.h"

#include "AliTRDdEdxBaseUtils.h"
#include "AliTRDdEdxReconUtils.h"
#include "AliTRDdEdxCalibHistArray.h"
#include "AliTRDdEdxCalibUtils.h"

#define EPSILON 1e-8

AliTRDdEdxCalibHistArray * AliTRDdEdxCalibUtils::fgHistArray = 0x0;
TObjArray * AliTRDdEdxCalibUtils::fgObjArray = 0x0;

//============================================================
//                           CalibObj
//============================================================
TObjArray * AliTRDdEdxCalibUtils::GetObjArray()
{
  //
  //return fgObjArray, initialized if null
  //

  if(!fgObjArray){
    fgObjArray = new TObjArray(8);
  }

  return fgObjArray;
}

TObjArray * AliTRDdEdxCalibUtils::GetObj(const Bool_t kinvq, const Double_t mag, const Int_t charge)
{
  //
  //return calib obj
  //
  if(!fgObjArray){
    printf("AliTRDdEdxCalibUtils::GetObjArray(kinvq, mag, charge) error fgObjArray null!!\n"); exit(1);
  }

  return (TObjArray*) fgObjArray->At(AliTRDdEdxCalibHistArray::GetIterator(kinvq, mag, charge));
}

void AliTRDdEdxCalibUtils::DeleteObjArray()
{
  //
  //delete calib obj
  //
  if(fgObjArray){
    fgObjArray->SetOwner();
    delete fgObjArray;
    fgObjArray = 0x0;
  }
}

Bool_t AliTRDdEdxCalibUtils::GenerateOCDB(const Int_t run, const TString path)
{
  //
  //generate OCDB object PHQ, do like
  //AliTRDdEdxCalibUtils::GenerateOCDB(run, "local://./")
  //if fgObjArray==0x0, generate default one
  //else generate according to fgObjArray
  //

  TObjArray * arr8 = 0x0;
  if(fgObjArray){
    arr8 = fgObjArray;
  }
  else{
    arr8 = new TObjArray(8);
    arr8->SetOwner();

    for(Int_t ii=0; ii<8; ii++){
      TObjArray * arr1 = new TObjArray(1);
      arr1->SetOwner();
      TString objn(AliTRDdEdxCalibHistArray::GetNameAt(ii));
      objn.ReplaceAll("Hist","Obj");
      arr1->SetName(objn);
      
      const Int_t nbins = AliTRDdEdxBaseUtils::NTRDtimebin();
      TVectorD * vec = new TVectorD(nbins);
      for(Int_t jj=0; jj<nbins; jj++){
        (*vec)[jj] = 1;
      }
      arr1->Add(vec);
      arr8->Add(arr1);
    }
  }

  AliCDBMetaData *metaData= new AliCDBMetaData();
  metaData->SetObjectClassName("TObjArray");
  metaData->SetResponsible("Raphaelle Bailhache and Xianguo Lu");

  AliCDBId id1("TRD/Calib/PHQ", run<0? 0 : run , run<0 ? 999999999 : run);
  AliCDBStorage * gStorage = AliCDBManager::Instance()->GetStorage(path);
  gStorage->Put(arr8, id1, metaData);

  delete metaData;
  delete arr8;

  return kTRUE;
}

//============================================================
//                           CalibHist
//============================================================
void AliTRDdEdxCalibUtils::DeleteHistArray()
{
  //
  //delete calib hist
  //
  delete fgHistArray;
  fgHistArray = 0x0;
}

THnBase * AliTRDdEdxCalibUtils::GetHistAt(const Int_t iter)
{
  //
  //
  //
  if(iter<fgHistArray->GetSize())
    return (THnBase*)fgHistArray->At(iter);
  else
    return 0x0;
}

void AliTRDdEdxCalibUtils::IniHistArray(TList *list, const Bool_t kNoInv)
{
  //
  //initialize calib hist, list should not own the hist, or list->Clear/delete hist should not be called
  //

  delete fgHistArray;
  fgHistArray = new AliTRDdEdxCalibHistArray(kNoInv);
  list->Add(fgHistArray);
}

Bool_t AliTRDdEdxCalibUtils::ReadHistArray(const TString filename, const TString listname)
{
  //
  //used in AliTRDPreprocessorOffline
  //read in calib hist from file, only for PHQ
  //

  //maybe already open by others... don't close
  TFile fcalib(filename);
  TObjArray * array = (TObjArray*)fcalib.Get(listname);
  const TString histname = AliTRDdEdxCalibHistArray::GetArrayName();
 
  AliTRDdEdxCalibHistArray * tmph=0x0;
  if(array){
    tmph = (AliTRDdEdxCalibHistArray *) array->FindObject(histname);
  }
  else{
    tmph = (AliTRDdEdxCalibHistArray *) fcalib.Get(histname);
  }
  if(!tmph){
    printf("AliTRDdEdxCalibUtils::ReadCalibHist warning calib hist not found! %s %s %s\n", filename.Data(), listname.Data(), histname.Data());
    fcalib.ls();
    if(array){
      array->ls();
    }
    return kFALSE;
  }
  
  delete fgHistArray; 
  fgHistArray = new AliTRDdEdxCalibHistArray(*tmph);

  return kTRUE;
}

void AliTRDdEdxCalibUtils::FillHist(const Int_t ncls, const TVectorD *arrayQ, const TVectorD *arrayX, THnBase * hcalib, const Double_t scale)
{
  //
  //fill calibration hist
  //
  if(!hcalib){printf("AliTRDdEdxCalibUtils::FillCalibHist errro hcalib null!!\n"); exit(1);}

  for(Int_t ii=0; ii<ncls; ii++){
    Double_t dq = (*arrayQ)[ii]/scale;
    const Double_t qmin = hcalib->GetAxis(1)->GetXmin() +0.5 * hcalib->GetAxis(1)->GetBinWidth(1);
    const Double_t qmax = hcalib->GetAxis(1)->GetXmax() -0.5 * hcalib->GetAxis(1)->GetBinWidth(1);

    if(dq<qmin)
      dq = qmin;
    if(dq>qmax)
      dq = qmax;

    const Double_t xx = (*arrayX)[ii];
    const Double_t xmin = hcalib->GetAxis(0)->GetXmin();
    const Double_t xmax = hcalib->GetAxis(0)->GetXmax();

    if(xx>=xmax || xx<xmin){
      printf("AliTRDdEdxCalibUtils::FillCalibHist error x overflow or underflow! %s %15f %15f %15f\n", hcalib->GetName(),  xx, xmin, xmax); exit(1);
    }

    const Double_t var[]={xx, dq};
    hcalib->Fill(var);
  }
}

void AliTRDdEdxCalibUtils::FillHist(const AliTRDtrackV1 *trdv1, const Bool_t kinvq, const Double_t mag, const Int_t charge, const Double_t scale) 
{
  //
  //get cluster Q and fill calib hist, if kinvq = kTRUE, 1/Q is filled
  //
  if(!fgHistArray){
    printf("AliTRDdEdxCalibUtils::FillHist fgHistArray not initialized!!\n"); exit(1);
  }

  THnBase * hcalib = fgHistArray->GetHist(kinvq, mag, charge);

  TVectorD arrayQ(200), arrayX(200);
  const Int_t ncls = AliTRDdEdxReconUtils::GetArrayClusterQ(kinvq, &arrayQ, &arrayX, trdv1);
  FillHist(ncls, &arrayQ, &arrayX, hcalib, kinvq ? 1/scale : scale);

  static Int_t kprint = 100;
  if(kprint<0){
    printf("\nAliTRDdEdxCalibUtils::FillHist summary: \n");
    printf("\nkinvq= %d;\n", kinvq);
    for(Int_t iq=0; iq<ncls; iq++){
      printf("arrayX[%3d] = %15.0f; arrayQ[%3d] =  %15f;\n", iq, arrayX[iq], iq, arrayQ[iq]);
    }
    printf("\n");
  }
  kprint++;
}

//============================================================
//
//============================================================

Double_t AliTRDdEdxCalibUtils::GetCalibTPCscale(const Int_t tpcncls, const Double_t tpcsig)
{
  //
  //the scale used in calibration
  //

  if(tpcncls < AliTRDdEdxBaseUtils::CalibTPCnclsCut())
    return -999;

  if(tpcsig<EPSILON)
    return -999;

  return tpcsig/45;

}

void AliTRDdEdxCalibUtils::GetPHCountMeanRMS(const TH1D *hnor, TH1D *&hmean)
{
  //
  //calculate from the ph calib hist the (mean-3sigma) ph-count in the chamber, save in the TH1D output
  //
  const Int_t ndet = 540;
  TObjArray *obj=new TObjArray(ndet);
  obj->SetOwner();
  for(Int_t ii=0; ii<ndet; ii++){
    obj->Add(new TVectorD(AliTRDseedV1::kNtb));
  }

  //ibin = binlowedge of bin(ibin+1) = the number fills this bin
  for(Int_t ibin=0; ibin<hnor->GetNbinsX(); ibin++){
    const Double_t stat = hnor->GetBinContent(ibin+1);
    if(stat<EPSILON){
      continue;
    }

    const Int_t idet = AliTRDdEdxBaseUtils::ToDetector(ibin);
    const Int_t itb  = AliTRDdEdxBaseUtils::ToTimeBin(ibin);
    TVectorD *vec=(TVectorD *)obj->At(idet);
    (*vec)[itb] = stat;
  }

  hmean = new TH1D(Form("%sdetmean", hnor->GetName()), "", hnor->GetNbinsX(), hnor->GetXaxis()->GetXmin(), hnor->GetXaxis()->GetXmax());
  for(Int_t ibin=0; ibin<hnor->GetNbinsX(); ibin++){
    const Int_t idet = AliTRDdEdxBaseUtils::ToDetector(ibin);
    const TVectorD *vec=(TVectorD *)obj->At(idet);

    Int_t nonzero = 0;
    for(Int_t ii=0; ii<vec->GetNrows(); ii++){
      if((*vec)[ii]>EPSILON){
        nonzero++;
      }
    }

    Double_t mean = 0;
    const Double_t lowfrac = 0.6;
    //if there are too many 0's, reject this chamber by settig mean=rms=0
    if(nonzero> (AliTRDseedV1::kNtb*(1-lowfrac)) ){
      //only highest (1-lowfrac)*31 timebins are used to estimate the mean and rms! important! otherwise the 0' will make rms very large!
      mean = AliTRDdEdxBaseUtils::TruncatedMean(AliTRDseedV1::kNtb, vec->GetMatrixArray(), lowfrac, 1);
    }

    hmean->SetBinContent(ibin+1, mean);
  }

  delete obj;
}

  /*
    http://root.cern.ch/root/html/TH1.html
    When an histogram is created, a reference to it is automatically added to the list of in-memory objects for the current file or directory. This default behaviour can be changed by:
    
    h->SetDirectory(0);          for the current histogram h
    TH1::AddDirectory(kFALSE);   sets a global switch disabling the reference
    
    When the histogram is deleted, the reference to it is removed from the list of objects in memory. When a file is closed, all histograms in memory associated with this file are automatically deleted. 
  */

TObjArray* AliTRDdEdxCalibUtils::HistToObj(const THnBase *hh, Int_t run, TList *lout, TTreeSRedirector *calibStream)
{
  //
  //produce calibration objects
  //

  const TString hname = hh->GetName();
  const Bool_t kinvq = TString(hname(hname.First('Q')+1,1)).Atoi()&4;

  //----------------------------------------
  //               Define nbin, tag, cobj0
  //----------------------------------------
  const Int_t nbin = AliTRDdEdxBaseUtils::NTRDtimebin();

    
  TString tag(hname);
  tag.ReplaceAll("Hist","Obj");
  
  TObjArray * cobj0 = new TObjArray(1);
  cobj0->SetOwner();
  cobj0->SetName(tag);
  cobj0->Add(new TVectorD(nbin));
  
  //----------------------------------------
  //               Define lowFrac, highFrac
  //----------------------------------------
  Double_t lowFrac = -999, highFrac = -999;
  if(!kinvq) {
    lowFrac = 0.01; highFrac = AliTRDdEdxBaseUtils::Q0Frac();
  }
  else{
    lowFrac = AliTRDdEdxBaseUtils::Q1Frac(); highFrac = 0.99;
  }
  
  //----------------------------------------
  //               Get analysis result
  //----------------------------------------
  TH1::AddDirectory(kFALSE);//important!
  TH1D *hnor=0x0, *hmpv=0x0, *hres=0x0, *hwid=0x0, *htrdphmean = 0x0;//if(!lout), these have to be deleted
  TH2D *hpj = hh->Projection(1,0);
  AliTRDdEdxBaseUtils::FitSlicesY(hpj, hnor, hmpv, hwid, hres, 0, lowFrac, highFrac);
  GetPHCountMeanRMS(hnor, htrdphmean);
  if(lout){
    lout->Add(htrdphmean);
  }
  delete hpj;
  
  if(lout){
    lout->Add(hnor);
    lout->Add(hmpv);
    lout->Add(hwid);
    lout->Add(hres);
  }
  
  //----------------------------------------
  //               Define Counter
  //----------------------------------------
  TVectorD *countDet=0x0;
  TObjArray *countSSL=0x0;

  if(!kinvq){
    countDet=new TVectorD(540);
    countSSL=new TObjArray(90);//SectorStackLayer
    countSSL->SetOwner();
    for(Int_t ii=0; ii<90; ii++){
      countSSL->Add(new TVectorD(6));
    }
  }

  //----------------------------------------
  //               Fill cobj0
  //----------------------------------------

  //ibin = binlowedge of bin(ibin+1) = the number fills this bin
  for(Int_t ibin=0; ibin<nbin; ibin++){
    Double_t gnor = hnor->GetBinContent(ibin+1);
    Double_t gmpv = hmpv->GetBinContent(ibin+1);
    Double_t gwid = hwid->GetBinContent(ibin+1);
    Double_t gres = hres->GetBinContent(ibin+1);

    //--- set additional cut by kpass
    Bool_t kpass = kTRUE;
    Double_t gtrdphmean = -999;
    if(htrdphmean){
      gtrdphmean = htrdphmean->GetBinContent(ibin+1);
      //chamber no statistics (e.g. too many 0's), not usual, not seen in run 143237
      if(gtrdphmean<EPSILON){
        kpass = kFALSE;
      }
      if(gnor<AliTRDdEdxBaseUtils::TimeBinCountCut()*gtrdphmean){
        kpass = kFALSE;
      }
    }
    
    //--- set calibration constant p0
    Double_t p0= 0;
    
    //reason for gmpv=0:
    //1)gnor<=3; truncation in hist: (0, 0.6*ntot=1.8 with ntot=3]={1}, in hist entries can pile up so that ntot=2, or 3, and (ntot>nstart && ntot<=nstop) is skipped;
    //2)TruncatedMean(hist) out of range (only for Q0, not Q1).
    
    if(gmpv>EPSILON && kpass){ 
      if(tag.Contains("T0")){
        p0 = gmpv;
      }
      else{
        p0 = 1/gmpv;
      }
      //printf("outcalibobj%s %d %15.6e\n", tag.Data(), ibin, p0);
    }

    (*( (TVectorD*)cobj0->At(0) ))[ibin] = p0;

    //--- save optional record
    if(p0>EPSILON && countDet && countSSL){
      const Int_t idet = AliTRDdEdxBaseUtils::ToDetector(ibin);
      (*countDet)[idet]=1;
      
      const Int_t isector = AliTRDdEdxBaseUtils::ToSector(ibin);
      const Int_t istack = AliTRDdEdxBaseUtils::ToStack(ibin);
      const Int_t ilayer = AliTRDdEdxBaseUtils::ToLayer(ibin);
      TVectorD * vecsectorstack = (TVectorD*)countSSL->At(istack*18+isector);
      (*vecsectorstack)[ilayer]=1;
    }
    
    if(calibStream){
      (*calibStream)<<tag<<
        "run="<<run<<
        "p0="<<p0<<
        
        "nor="<<gnor<<
        "mpv="<<gmpv<<
        "wid="<<gwid<<
        "res="<<gres<<
        "gtrdphmean="<<gtrdphmean<<
        
        "ibin="<<ibin<<
        "\n";
    }
  }
  
  //----------------------------------------
  //               Status Report
  //----------------------------------------
  if(countDet && countSSL){
    TVectorD count2Dstack(90);
    for(Int_t ii=0; ii<90; ii++){
      TVectorD * vecsectorstack = (TVectorD*)countSSL->At(ii);
      const Int_t nlayer = (Int_t)vecsectorstack->Sum();
      if(nlayer==6){
        count2Dstack[ii]=1;
      }
    }

    printf("\nAliTRDdEdxCalibUtils::GetCalibObj Summary run: %d name: %s entries: %.0f ndetector: %03.0f n2dstack %02.0f\n\n", run, hname.Data(), hh->GetEntries(), countDet->Sum(), count2Dstack.Sum());
  }

  //----------------------------------------
  //               Clean Up
  //----------------------------------------
  
  TH1D **hhs[]={&hnor, &hmpv, &hwid, &hres, &htrdphmean};
  const Int_t nhh=sizeof(hhs)/sizeof(TH1D**);
  for(Int_t ihh=0; ihh<nhh; ihh++){
    if(!lout){
      delete (*hhs[ihh]);
    }
  }
  
  delete countDet;
  delete countSSL;

  //----------------------------------------

  return cobj0;
}
Commit	Line	Data
6951a056	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	//
	17	// class to calculate TRD dEdx
	18	// xx
	19	// xx
	20	// xx
	21	// xx
	22	//
	23	// Xianguo Lu
	24	// lu@physi.uni-heidelberg.de
	25	// Xianguo.Lu@cern.ch
	26	//
	27	//
	28
	29
	30	#include "TF1.h"
	31	#include "TFile.h"
	32	#include "TH1D.h"
	33	#include "TH2D.h"
	34	#include "THnSparse.h"
	35	#include "TMath.h"
	36	#include "TMatrixD.h"
	37	#include "TMinuit.h"
	38	#include "TObjArray.h"
	39	#include "TRandom3.h"
	40	#include "TStopwatch.h"
	41	#include "TVectorD.h"
	42
	43	#include "TTreeStream.h"
	44
	45	#include "AliCDBId.h"
	46	#include "AliCDBMetaData.h"
	47	#include "AliCDBStorage.h"
	48	#include "AliESDEvent.h"
	49	#include "AliESDfriendTrack.h"
	50	#include "AliESDtrack.h"
	51	#include "AliTRDcalibDB.h"
	52	#include "AliTRDCalROC.h"
	53	#include "AliTRDtrackV1.h"
	54
	55	#include "AliTRDdEdxBaseUtils.h"
	56	#include "AliTRDdEdxReconUtils.h"
	57	#include "AliTRDdEdxCalibHistArray.h"
	58	#include "AliTRDdEdxCalibUtils.h"
	59
6fa94550	60	#define EPSILON 1e-8
6951a056	61
	62	AliTRDdEdxCalibHistArray * AliTRDdEdxCalibUtils::fgHistArray = 0x0;
	63	TObjArray * AliTRDdEdxCalibUtils::fgObjArray = 0x0;
	64
	65	//============================================================
	66	// CalibObj
	67	//============================================================
	68	TObjArray * AliTRDdEdxCalibUtils::GetObjArray()
	69	{
	70	//
	71	//return fgObjArray, initialized if null
	72	//
	73
	74	if(!fgObjArray){
	75	fgObjArray = new TObjArray(8);
	76	}
	77
	78	return fgObjArray;
	79	}
	80
	81	TObjArray * AliTRDdEdxCalibUtils::GetObj(const Bool_t kinvq, const Double_t mag, const Int_t charge)
	82	{
	83	//
	84	//return calib obj
	85	//
	86	if(!fgObjArray){
	87	printf("AliTRDdEdxCalibUtils::GetObjArray(kinvq, mag, charge) error fgObjArray null!!\n"); exit(1);
	88	}
	89
	90	return (TObjArray*) fgObjArray->At(AliTRDdEdxCalibHistArray::GetIterator(kinvq, mag, charge));
	91	}
	92
	93	void AliTRDdEdxCalibUtils::DeleteObjArray()
	94	{
	95	//
	96	//delete calib obj
	97	//
	98	if(fgObjArray){
	99	fgObjArray->SetOwner();
	100	delete fgObjArray;
	101	fgObjArray = 0x0;
	102	}
	103	}
	104
687aa844	105	Bool_t AliTRDdEdxCalibUtils::GenerateOCDB(const Int_t run, const TString path)
6951a056	106	{
6951a056	107	//
687aa844	108	//generate OCDB object PHQ, do like
	109	//AliTRDdEdxCalibUtils::GenerateOCDB(run, "local://./")
	110	//if fgObjArray==0x0, generate default one
	111	//else generate according to fgObjArray
6951a056	112	//
6951a056	113
687aa844	114	TObjArray * arr8 = 0x0;
	115	if(fgObjArray){
	116	arr8 = fgObjArray;
	117	}
	118	else{
	119	arr8 = new TObjArray(8);
	120	arr8->SetOwner();
	121
	122	for(Int_t ii=0; ii<8; ii++){
	123	TObjArray * arr1 = new TObjArray(1);
	124	arr1->SetOwner();
	125	TString objn(AliTRDdEdxCalibHistArray::GetNameAt(ii));
	126	objn.ReplaceAll("Hist","Obj");
	127	arr1->SetName(objn);
	128
	129	const Int_t nbins = AliTRDdEdxBaseUtils::NTRDtimebin();
	130	TVectorD * vec = new TVectorD(nbins);
	131	for(Int_t jj=0; jj<nbins; jj++){
	132	(*vec)[jj] = 1;
	133	}
	134	arr1->Add(vec);
	135	arr8->Add(arr1);
6951a056	136	}
6951a056	137	}
	138
	139	AliCDBMetaData *metaData= new AliCDBMetaData();
	140	metaData->SetObjectClassName("TObjArray");
	141	metaData->SetResponsible("Raphaelle Bailhache and Xianguo Lu");
	142
fe829aa0	143	AliCDBId id1("TRD/Calib/PHQ", run<0? 0 : run , run<0 ? 999999999 : run);
6951a056	144	AliCDBStorage * gStorage = AliCDBManager::Instance()->GetStorage(path);
	145	gStorage->Put(arr8, id1, metaData);
	146
	147	delete metaData;
	148	delete arr8;
	149
	150	return kTRUE;
	151	}
	152
	153	//============================================================
	154	// CalibHist
	155	//============================================================
	156	void AliTRDdEdxCalibUtils::DeleteHistArray()
	157	{
	158	//
	159	//delete calib hist
	160	//
	161	delete fgHistArray;
	162	fgHistArray = 0x0;
	163	}
	164
	165	THnBase * AliTRDdEdxCalibUtils::GetHistAt(const Int_t iter)
	166	{
	167	//
	168	//
	169	//
	170	if(iter<fgHistArray->GetSize())
	171	return (THnBase*)fgHistArray->At(iter);
	172	else
	173	return 0x0;
	174	}
	175
	176	void AliTRDdEdxCalibUtils::IniHistArray(TList *list, const Bool_t kNoInv)
	177	{
	178	//
	179	//initialize calib hist, list should not own the hist, or list->Clear/delete hist should not be called
	180	//
	181
	182	delete fgHistArray;
	183	fgHistArray = new AliTRDdEdxCalibHistArray(kNoInv);
	184	list->Add(fgHistArray);
	185	}
	186
	187	Bool_t AliTRDdEdxCalibUtils::ReadHistArray(const TString filename, const TString listname)
	188	{
	189	//
	190	//used in AliTRDPreprocessorOffline
	191	//read in calib hist from file, only for PHQ
	192	//
	193
	194	//maybe already open by others... don't close
	195	TFile fcalib(filename);
	196	TObjArray * array = (TObjArray*)fcalib.Get(listname);
	197	const TString histname = AliTRDdEdxCalibHistArray::GetArrayName();
	198
	199	AliTRDdEdxCalibHistArray * tmph=0x0;
	200	if(array){
	201	tmph = (AliTRDdEdxCalibHistArray *) array->FindObject(histname);
	202	}
	203	else{
	204	tmph = (AliTRDdEdxCalibHistArray *) fcalib.Get(histname);
	205	}
	206	if(!tmph){
	207	printf("AliTRDdEdxCalibUtils::ReadCalibHist warning calib hist not found! %s %s %s\n", filename.Data(), listname.Data(), histname.Data());
208	fcalib.ls();
209	if(array){
210	array->ls();
211	}
212	return kFALSE;
213	}
214
215	delete fgHistArray;
216	fgHistArray = new AliTRDdEdxCalibHistArray(*tmph);
217
218	return kTRUE;
219	}
220
221	void AliTRDdEdxCalibUtils::FillHist(const Int_t ncls, const TVectorD arrayQ, const TVectorD arrayX, THnBase * hcalib, const Double_t scale)
222	{
223	//
224	//fill calibration hist
225	//
226	if(!hcalib){printf("AliTRDdEdxCalibUtils::FillCalibHist errro hcalib null!!\n"); exit(1);}
227
228	for(Int_t ii=0; ii<ncls; ii++){
6856b6a8	229	Double_t dq = (*arrayQ)[ii]/scale;
6856b6a8	230	const Double_t qmin = hcalib->GetAxis(1)->GetXmin() +0.5 * hcalib->GetAxis(1)->GetBinWidth(1);
6951a056	231	const Double_t qmax = hcalib->GetAxis(1)->GetXmax() -0.5 * hcalib->GetAxis(1)->GetBinWidth(1);
6856b6a8	232
	233	if(dq<qmin)
	234	dq = qmin;
	235	if(dq>qmax)
	236	dq = qmax;
	237
	238	const Double_t xx = (*arrayX)[ii];
6951a056	239	const Double_t xmin = hcalib->GetAxis(0)->GetXmin();
	240	const Double_t xmax = hcalib->GetAxis(0)->GetXmax();
	241
	242	if(xx>=xmax \|\| xx<xmin){
	243	printf("AliTRDdEdxCalibUtils::FillCalibHist error x overflow or underflow! %s %15f %15f %15f\n", hcalib->GetName(), xx, xmin, xmax); exit(1);
	244	}
	245
6856b6a8	246	const Double_t var[]={xx, dq};
6951a056	247	hcalib->Fill(var);
	248	}
	249	}
	250
	251	void AliTRDdEdxCalibUtils::FillHist(const AliTRDtrackV1 *trdv1, const Bool_t kinvq, const Double_t mag, const Int_t charge, const Double_t scale)
	252	{
	253	//
	254	//get cluster Q and fill calib hist, if kinvq = kTRUE, 1/Q is filled
	255	//
	256	if(!fgHistArray){
	257	printf("AliTRDdEdxCalibUtils::FillHist fgHistArray not initialized!!\n"); exit(1);
	258	}
	259
	260	THnBase * hcalib = fgHistArray->GetHist(kinvq, mag, charge);
	261
	262	TVectorD arrayQ(200), arrayX(200);
	263	const Int_t ncls = AliTRDdEdxReconUtils::GetArrayClusterQ(kinvq, &arrayQ, &arrayX, trdv1);
	264	FillHist(ncls, &arrayQ, &arrayX, hcalib, kinvq ? 1/scale : scale);
	265
	266	static Int_t kprint = 100;
	267	if(kprint<0){
	268	printf("\nAliTRDdEdxCalibUtils::FillHist summary: \n");
	269	printf("\nkinvq= %d;\n", kinvq);
	270	for(Int_t iq=0; iq<ncls; iq++){
	271	printf("arrayX[%3d] = %15.0f; arrayQ[%3d] = %15f;\n", iq, arrayX[iq], iq, arrayQ[iq]);
	272	}
	273	printf("\n");
	274	}
	275	kprint++;
	276	}
	277
	278	//============================================================
	279	//
	280	//============================================================
	281
	282	Double_t AliTRDdEdxCalibUtils::GetCalibTPCscale(const Int_t tpcncls, const Double_t tpcsig)
	283	{
	284	//
	285	//the scale used in calibration
	286	//
	287
	288	if(tpcncls < AliTRDdEdxBaseUtils::CalibTPCnclsCut())
	289	return -999;
	290
	291	if(tpcsig<EPSILON)
	292	return -999;
	293
687aa844	294	return tpcsig/45;
6951a056	295
	296	}
	297
	298	void AliTRDdEdxCalibUtils::GetPHCountMeanRMS(const TH1D hnor, TH1D &hmean)
	299	{
	300	//
	301	//calculate from the ph calib hist the (mean-3sigma) ph-count in the chamber, save in the TH1D output
	302	//
	303	const Int_t ndet = 540;
	304	TObjArray *obj=new TObjArray(ndet);
	305	obj->SetOwner();
	306	for(Int_t ii=0; ii<ndet; ii++){
	307	obj->Add(new TVectorD(AliTRDseedV1::kNtb));
	308	}
	309
	310	//ibin = binlowedge of bin(ibin+1) = the number fills this bin
	311	for(Int_t ibin=0; ibin<hnor->GetNbinsX(); ibin++){
	312	const Double_t stat = hnor->GetBinContent(ibin+1);
	313	if(stat<EPSILON){
	314	continue;
	315	}
	316
	317	const Int_t idet = AliTRDdEdxBaseUtils::ToDetector(ibin);
	318	const Int_t itb = AliTRDdEdxBaseUtils::ToTimeBin(ibin);
	319	TVectorD vec=(TVectorD )obj->At(idet);
	320	(*vec)[itb] = stat;
	321	}
	322
	323	hmean = new TH1D(Form("%sdetmean", hnor->GetName()), "", hnor->GetNbinsX(), hnor->GetXaxis()->GetXmin(), hnor->GetXaxis()->GetXmax());
	324	for(Int_t ibin=0; ibin<hnor->GetNbinsX(); ibin++){
	325	const Int_t idet = AliTRDdEdxBaseUtils::ToDetector(ibin);
	326	const TVectorD vec=(TVectorD )obj->At(idet);
	327
	328	Int_t nonzero = 0;
	329	for(Int_t ii=0; ii<vec->GetNrows(); ii++){
	330	if((*vec)[ii]>EPSILON){
	331	nonzero++;
	332	}
	333	}
	334
	335	Double_t mean = 0;
	336	const Double_t lowfrac = 0.6;
	337	//if there are too many 0's, reject this chamber by settig mean=rms=0
	338	if(nonzero> (AliTRDseedV1::kNtb*(1-lowfrac)) ){
	339	//only highest (1-lowfrac)*31 timebins are used to estimate the mean and rms! important! otherwise the 0' will make rms very large!
	340	mean = AliTRDdEdxBaseUtils::TruncatedMean(AliTRDseedV1::kNtb, vec->GetMatrixArray(), lowfrac, 1);
	341	}
	342
	343	hmean->SetBinContent(ibin+1, mean);
	344	}
	345
	346	delete obj;
	347	}
	348
6951a056	349	/*
	350	http://root.cern.ch/root/html/TH1.html
	351	When an histogram is created, a reference to it is automatically added to the list of in-memory objects for the current file or directory. This default behaviour can be changed by:
	352
	353	h->SetDirectory(0); for the current histogram h
	354	TH1::AddDirectory(kFALSE); sets a global switch disabling the reference
	355
	356	When the histogram is deleted, the reference to it is removed from the list of objects in memory. When a file is closed, all histograms in memory associated with this file are automatically deleted.
	357	*/
6951a056	358
	359	TObjArray* AliTRDdEdxCalibUtils::HistToObj(const THnBase hh, Int_t run, TList lout, TTreeSRedirector *calibStream)
	360	{
	361	//
	362	//produce calibration objects
	363	//
	364
	365	const TString hname = hh->GetName();
	366	const Bool_t kinvq = TString(hname(hname.First('Q')+1,1)).Atoi()&4;
	367
	368	//----------------------------------------
	369	// Define nbin, tag, cobj0
	370	//----------------------------------------
	371	const Int_t nbin = AliTRDdEdxBaseUtils::NTRDtimebin();
	372
	373
	374	TString tag(hname);
	375	tag.ReplaceAll("Hist","Obj");
	376
	377	TObjArray * cobj0 = new TObjArray(1);
	378	cobj0->SetOwner();
	379	cobj0->SetName(tag);
	380	cobj0->Add(new TVectorD(nbin));
	381
	382	//----------------------------------------
	383	// Define lowFrac, highFrac
	384	//----------------------------------------
	385	Double_t lowFrac = -999, highFrac = -999;
	386	if(!kinvq) {
	387	lowFrac = 0.01; highFrac = AliTRDdEdxBaseUtils::Q0Frac();
	388	}
	389	else{
	390	lowFrac = AliTRDdEdxBaseUtils::Q1Frac(); highFrac = 0.99;
	391	}
	392
	393	//----------------------------------------
	394	// Get analysis result
	395	//----------------------------------------
	396	TH1::AddDirectory(kFALSE);//important!
	397	TH1D hnor=0x0, hmpv=0x0, hres=0x0, hwid=0x0, *htrdphmean = 0x0;//if(!lout), these have to be deleted
	398	TH2D *hpj = hh->Projection(1,0);
	399	AliTRDdEdxBaseUtils::FitSlicesY(hpj, hnor, hmpv, hwid, hres, 0, lowFrac, highFrac);
	400	GetPHCountMeanRMS(hnor, htrdphmean);
	401	if(lout){
	402	lout->Add(htrdphmean);
	403	}
	404	delete hpj;
	405
	406	if(lout){
	407	lout->Add(hnor);
	408	lout->Add(hmpv);
	409	lout->Add(hwid);
	410	lout->Add(hres);
	411	}
	412
	413	//----------------------------------------
	414	// Define Counter
	415	//----------------------------------------
	416	TVectorD *countDet=0x0;
	417	TObjArray *countSSL=0x0;
	418
	419	if(!kinvq){
	420	countDet=new TVectorD(540);
	421	countSSL=new TObjArray(90);//SectorStackLayer
422	countSSL->SetOwner();
423	for(Int_t ii=0; ii<90; ii++){
424	countSSL->Add(new TVectorD(6));
425	}
426	}
427
428	//----------------------------------------
429	// Fill cobj0
430	//----------------------------------------
431
432	//ibin = binlowedge of bin(ibin+1) = the number fills this bin
433	for(Int_t ibin=0; ibin<nbin; ibin++){
434	Double_t gnor = hnor->GetBinContent(ibin+1);
435	Double_t gmpv = hmpv->GetBinContent(ibin+1);
436	Double_t gwid = hwid->GetBinContent(ibin+1);
437	Double_t gres = hres->GetBinContent(ibin+1);
438
439	//--- set additional cut by kpass
440	Bool_t kpass = kTRUE;
441	Double_t gtrdphmean = -999;
442	if(htrdphmean){
443	gtrdphmean = htrdphmean->GetBinContent(ibin+1);
444	//chamber no statistics (e.g. too many 0's), not usual, not seen in run 143237
445	if(gtrdphmean<EPSILON){
446	kpass = kFALSE;
447	}
448	if(gnor<AliTRDdEdxBaseUtils::TimeBinCountCut()*gtrdphmean){
449	kpass = kFALSE;
450	}
451	}
452
453	//--- set calibration constant p0
454	Double_t p0= 0;
455
456	//reason for gmpv=0:
457	//1)gnor<=3; truncation in hist: (0, 0.6*ntot=1.8 with ntot=3]={1}, in hist entries can pile up so that ntot=2, or 3, and (ntot>nstart && ntot<=nstop) is skipped;
458	//2)TruncatedMean(hist) out of range (only for Q0, not Q1).
459
460	if(gmpv>EPSILON && kpass){
461	if(tag.Contains("T0")){
462	p0 = gmpv;
463	}
464	else{
465	p0 = 1/gmpv;
466	}
467	//printf("outcalibobj%s %d %15.6e\n", tag.Data(), ibin, p0);
468	}
469
470	(( (TVectorD)cobj0->At(0) ))[ibin] = p0;
471
472	//--- save optional record
473	if(p0>EPSILON && countDet && countSSL){
474	const Int_t idet = AliTRDdEdxBaseUtils::ToDetector(ibin);
475	(*countDet)[idet]=1;
476
477	const Int_t isector = AliTRDdEdxBaseUtils::ToSector(ibin);
478	const Int_t istack = AliTRDdEdxBaseUtils::ToStack(ibin);
479	const Int_t ilayer = AliTRDdEdxBaseUtils::ToLayer(ibin);
480	TVectorD * vecsectorstack = (TVectorD)countSSL->At(istack18+isector);
481	(*vecsectorstack)[ilayer]=1;
482	}
483
484	if(calibStream){
485	(*calibStream)<<tag<<
486	"run="<<run<<
487	"p0="<<p0<<
488
489	"nor="<<gnor<<
490	"mpv="<<gmpv<<
491	"wid="<<gwid<<
492	"res="<<gres<<
493	"gtrdphmean="<<gtrdphmean<<
494
495	"ibin="<<ibin<<
496	"\n";
497	}
498	}
499
500	//----------------------------------------
501	// Status Report
502	//----------------------------------------
503	if(countDet && countSSL){
504	TVectorD count2Dstack(90);
505	for(Int_t ii=0; ii<90; ii++){
506	TVectorD * vecsectorstack = (TVectorD*)countSSL->At(ii);
507	const Int_t nlayer = (Int_t)vecsectorstack->Sum();
508	if(nlayer==6){
509	count2Dstack[ii]=1;
510	}
511	}
512
513	printf("\nAliTRDdEdxCalibUtils::GetCalibObj Summary run: %d name: %s entries: %.0f ndetector: %03.0f n2dstack %02.0f\n\n", run, hname.Data(), hh->GetEntries(), countDet->Sum(), count2Dstack.Sum());
514	}
515
516	//----------------------------------------
517	// Clean Up
518	//----------------------------------------
519
520	TH1D **hhs[]={&hnor, &hmpv, &hwid, &hres, &htrdphmean};
521	const Int_t nhh=sizeof(hhs)/sizeof(TH1D**);
522	for(Int_t ihh=0; ihh<nhh; ihh++){
523	if(!lout){
524	delete (*hhs[ihh]);
525	}
526	}
527
528	delete countDet;
529	delete countSSL;
530
531	//----------------------------------------
532
533	return cobj0;
534	}
535