[u/mrichter/AliRoot.git] / TOF / AliTOFT0maker.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.                  *
 **************************************************************************/
/* $Id: AliTOFT0maker.cxx,v 1.8 2010/01/19 16:32:20 noferini Exp $ */

/////////////////////////////////////////////////////////////////////////////
//                                                                         //
//  This class contains the basic functions for the time zero              //
//  evaluation with TOF detector informations.                             //
// Use case in an analysis task:                                           //
//                                                                         //
// Create the object in the task constructor (fTOFmaker is a private var)  //
// AliESDpid *extPID=new AliESDpid();                                      //
// fTOFmaker = new AliTOFT0maker(extPID);                                  //
// fTOFmaker->SetTimeResolution(100.0); // if you want set the TOF res     //
// 115 ps is the TOF default resolution value                              //
//                                                                         //
// Use the RemakePID method in the task::Exec                              //
// Double_t* calcolot0;                                                    //
// calcolot0=fTOFmaker->RemakePID(fESD);                                   //
// //calcolot0[0] = calculated event time                                  // 
// //calcolot0[1] = event time time resolution                             //
// //calcolot0[2] = average event time for the current fill                //
// //calcolot0[3] = tracks at TOF                                          // 
// //calcolot0[4] = calculated event time (only TOF)                       //
// //calcolot0[5] = event time time resolution (only TOF)                  //
// //calcolot0[6] = sigma t0 fill                                          //
// //calcolot0[7] = tracks at TOF really used in tht algorithm             // 
//                                                                         //
// Let consider that:                                                      //
// - the PIF is automatically recalculated with the event time subtrction  //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////

#include "AliTOFT0v1.h"
#include "AliTOFT0maker.h"
#include "AliPID.h"
#include "AliLog.h"
#include "AliESDpid.h"
#include "AliESDEvent.h"
#include "TFile.h"
#include "TH1F.h"
#include "AliTOFcalib.h"
#include "AliTOFRunParams.h"
#include "TRandom.h"
#include "AliTOFHeader.h"

ClassImp(AliTOFT0maker)
           
//____________________________________________________________________________ 
AliTOFT0maker::AliTOFT0maker():
  TObject(),
  fT0TOF(NULL),
  fPIDesd(NULL),
  fExternalPIDFlag(kFALSE),
  fTOFcalib(NULL),
  fNoTOFT0(0),
  fNmomBins(0),
  fTimeResolution(100),
  fT0sigma(1000),
  fHmapChannel(0),
  fKmask(0),
  fT0width(150.),
  fT0spreadExt(-1.),
  fT0fillExt(0)
{
  // ctr
  fCalculated[0] = 0;
  fCalculated[1] = 0;
  fCalculated[2] = 0;
  fCalculated[3] = 0;

  fT0cur[0]=0.;
  fT0cur[1]=0.;

  if(AliPID::ParticleMass(0) == 0) new AliPID();

  fPIDesd = new AliESDpid();

  fNmomBins = fPIDesd->GetTOFResponse().GetNmomBins();
  SetTOFResponse();

  fT0TOF = new AliTOFT0v1(fPIDesd);

}
//____________________________________________________________________________ 
AliTOFT0maker::AliTOFT0maker(AliESDpid *externalPID, AliTOFcalib *tofCalib):
    TObject(),
    fT0TOF(NULL),
    fPIDesd(externalPID),
    fExternalPIDFlag(kTRUE),
    fTOFcalib(tofCalib),
    fNoTOFT0(0),
    fNmomBins(0),
    fTimeResolution(100),
    fT0sigma(1000),
    fHmapChannel(0),
    fKmask(0),
    fT0width(150.),
    fT0spreadExt(-1.),
    fT0fillExt(0)
{
  // ctr
  fCalculated[0] = 0;
  fCalculated[1] = 0;
  fCalculated[2] = 0;
  fCalculated[3] = 0;

  fT0cur[0]=0.;
  fT0cur[1]=0.;

  if(AliPID::ParticleMass(0) == 0) new AliPID();

  if(!fPIDesd){
    fPIDesd = new AliESDpid();
    fExternalPIDFlag = kFALSE;
  }

  fNmomBins = fPIDesd->GetTOFResponse().GetNmomBins();
  SetTOFResponse();

  fT0TOF = new AliTOFT0v1(fPIDesd);

}

//____________________________________________________________________________ 
AliTOFT0maker::~AliTOFT0maker()
{
  // dtor

  delete fT0TOF;
  if (!fExternalPIDFlag) delete fPIDesd;
}
//____________________________________________________________________________ 
Double_t* AliTOFT0maker::ComputeT0TOF(AliESDEvent *esd,Double_t t0time,Double_t t0sigma){
  //
  // Remake TOF PID probabilities
  //
  Double_t t0tof[6];

  if(fKmask) ApplyMask(esd);

  Double_t t0fill = 0.;

  fPIDesd->GetTOFResponse().ResetT0info();

  /* get T0 spread from TOFcalib if available otherwise use default value */
  if (fTOFcalib && esd) {
    AliTOFRunParams *runParams = fTOFcalib->GetRunParams();
    if (runParams && runParams->GetTimestamp(0) != 0) {
      Float_t t0spread = runParams->EvalT0Spread(esd->GetTimeStamp());
      if(fT0spreadExt > 0) SetT0FillWidth(fT0spreadExt);
      else{
	SetT0FillWidth(t0spread);
	t0fill = fT0fillExt;
      }
    }
    else{
      if(fT0spreadExt > 0) SetT0FillWidth(fT0spreadExt);
      t0fill = fT0fillExt;
    }
  }
  else if(esd){
      Float_t t0spread = esd->GetSigma2DiamondZ(); // vertex pread ^2
      if(t0spread > 0) t0spread = TMath::Sqrt(t0spread)/0.0299792458;

      if(fT0spreadExt > 0) SetT0FillWidth(fT0spreadExt);
      else{
	  SetT0FillWidth(t0spread);
	  t0fill = fT0fillExt;
      }
  }

  Float_t thrGood = TMath::Max(Float_t(500.),fT0width*3);

  fT0TOF->Init(esd);
  AliTOFT0v1* t0maker= fT0TOF;

  t0maker->DefineT0("all",1.5,3.0);
  t0tof[0] = t0maker->GetResult(0);
  t0tof[1] = t0maker->GetResult(1);
  t0tof[2] = t0maker->GetResult(2);
  t0tof[3] = t0maker->GetResult(3);
  t0tof[4] = t0maker->GetResult(4);
  t0tof[5] = t0maker->GetResult(5);

  Float_t lT0Current=0.;
  fT0sigma=1000;

//   Int_t nrun = esd->GetRunNumber();

  t0time += t0fill;

  Float_t sigmaFill = fT0width;

  if(sigmaFill < 20) sigmaFill = 140;

  fCalculated[0]=-1000*t0tof[0]; // best t0
  fCalculated[1]=1000*t0tof[1]; // sigma best t0
  fCalculated[2] = t0fill;    //t0 fill
  fCalculated[3] = t0tof[2];  // n TOF tracks
  fCalculated[4]=-1000*t0tof[0]; // TOF t0
  fCalculated[5]=1000*t0tof[1]; // TOF t0 sigma
  fCalculated[6]=sigmaFill; // sigma t0 fill
  fCalculated[7] = t0tof[3];  // n TOF tracks used for T0

  if(fCalculated[7] > 30) thrGood = 10000000;

  //statistics
  fCalculated[8] = t0tof[4]; // real time in s
  fCalculated[9] = t0tof[5]; // cpu time in s

  if(fCalculated[1] < sigmaFill && TMath::Abs(fCalculated[0] - t0fill) < thrGood && fCalculated[1] < fTimeResolution*1.2){
    fT0sigma=fCalculated[1];
    lT0Current=fCalculated[0];
  }
  else{
    fCalculated[4] = t0fill;
    fCalculated[5] = sigmaFill;
  }

  if(fCalculated[1] < 1 || fT0sigma > sigmaFill || fCalculated[1] > fTimeResolution* 1.2){
    fT0sigma =1000;
    fCalculated[4] = t0fill;
    fCalculated[5] = sigmaFill;
  }

  if(t0sigma < 1000){
    if(fT0sigma < 1000){
      Double_t w1 = 1./t0sigma/t0sigma;
      Double_t w2 = 1./fCalculated[1]/fCalculated[1];

      Double_t wtot = w1+w2;

      lT0Current = (w1*t0time + w2*fCalculated[0]) / wtot;
      fT0sigma = TMath::Sqrt(1./wtot);
    }
    else{
      lT0Current=t0time;
      fT0sigma=t0sigma;
    }
  }

  if(fT0sigma < sigmaFill && TMath::Abs(lT0Current - t0fill) < thrGood){
    fCalculated[1]=fT0sigma;
    fCalculated[0]=lT0Current;
  }

  if(fT0sigma >= 1000 || fNoTOFT0){
    lT0Current = t0fill;
    fT0sigma = sigmaFill;

    fCalculated[0] = t0fill;
    fCalculated[1] = sigmaFill;
  }

  // T0 pt bin
  Float_t *t0values = new Float_t[fNmomBins];
  Float_t *t0resolution = new Float_t[fNmomBins];
  if(fCalculated[7] < 100){
    for(Int_t i=0;i<fNmomBins;i++){
      t0maker->DefineT0("all",fPIDesd->GetTOFResponse().GetMinMom(i),fPIDesd->GetTOFResponse().GetMaxMom(i));
      t0tof[0] = t0maker->GetResult(0);
      t0tof[1] = t0maker->GetResult(1);
      t0tof[2] = t0maker->GetResult(2);
      t0tof[3] = t0maker->GetResult(3);

      Float_t t0bin =-1000*t0tof[0]; // best t0
      Float_t t0binRes =1000*t0tof[1]; // sigma best t0
      
      if(t0binRes < sigmaFill  && t0binRes < fTimeResolution * 1.2 && TMath::Abs(t0bin - t0fill) < thrGood){
	// Ok T0
	if(t0sigma < 1000){
	  Double_t w1 = 1./t0sigma/t0sigma;
	  Double_t w2 = 1./t0binRes/t0binRes;
	  
	  Double_t wtot = w1+w2;
	  
	  t0bin = (w1*t0time + w2*t0bin) / wtot;
	  t0binRes = TMath::Sqrt(1./wtot);
	}
      }
      else{
	t0bin = t0fill;
	t0binRes = sigmaFill;
	if(t0sigma < 1000){
	  t0bin = t0time;
	   t0binRes= t0sigma;	  
	}
      }
      t0values[i] = t0bin;
      t0resolution[i] = t0binRes;
    }
  }
  else{
    for(Int_t i=0;i<fNmomBins;i++){
      t0values[i] = lT0Current;
      t0resolution[i] = fT0sigma;
    }
  }
  for(Int_t i=0;i<fNmomBins;i++){
    fPIDesd->GetTOFResponse().SetT0bin(i,t0values[i]);
    fPIDesd->GetTOFResponse().SetT0binRes(i,t0resolution[i]);
  }
  
  delete[] t0values;
  delete[] t0resolution;

  return fCalculated;
}
//____________________________________________________________________________ 
Double_t  *AliTOFT0maker::GetT0p(Float_t p){// [0]=to -- [1] = sigma T0
  Int_t i=fPIDesd->GetTOFResponse().GetMomBin(p);
  
  fT0cur[0] = fPIDesd->GetTOFResponse().GetT0bin(i);
  fT0cur[1] = fPIDesd->GetTOFResponse().GetT0binRes(i);

  return fT0cur;
}
//____________________________________________________________________________ 
void AliTOFT0maker::SetTOFResponse(){
  fPIDesd->GetTOFResponse().SetTimeResolution(fTimeResolution);
}
//____________________________________________________________________________ 
Float_t AliTOFT0maker::GetExpectedSigma(Float_t mom, Float_t tof, Float_t mass){
  Float_t sigma = fPIDesd->GetTOFResponse().GetExpectedSigma(mom,tof,mass);

  return sigma;
}
//____________________________________________________________________________ 
void AliTOFT0maker::ApplyT0TOF(AliESDEvent *esd){
  //
  // Recalculate TOF PID probabilities
  //

// subtruct t0 for each track
   Int_t ntracks = esd->GetNumberOfTracks();
  
   while (ntracks--) {
       AliESDtrack *t=esd->GetTrack(ntracks);
    
       if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
    
       Double_t time=t->GetTOFsignal();
       Float_t p = t->GetP();

       Double_t *t0=GetT0p(p);
       time -= t0[0];
       t->SetTOFsignal(time);
   }

   for(Int_t i=0;i<fNmomBins;i++){
       fPIDesd->GetTOFResponse().SetT0bin(i,0.0);
   }
   
  //
}
//____________________________________________________________________________ 
void  AliTOFT0maker::LoadChannelMap(char *filename){
  // Load the histo with the channel off map
  TFile *f= new TFile(filename);
  if(!f){
    printf("Cannot open the channel map file (%s)\n",filename);
    return;
  }
  
  fHmapChannel = (TH1F *) f->Get("hChEnabled");
  
  if(!fHmapChannel){
    printf("Cannot laod the channel map histo (from %s)\n",filename);
    return;
  }
    
}
//____________________________________________________________________________ 
void AliTOFT0maker::ApplyMask(AliESDEvent * const esd){
  // Switch off the disable channel
  if(!fHmapChannel){
    printf("Channel Map is not available\n");
    return;
  }
  
  Int_t ntracks = esd->GetNumberOfTracks();
  
  while (ntracks--) {
    AliESDtrack *t=esd->GetTrack(ntracks);    

    if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;

    Int_t chan = t->GetTOFCalChannel();
 
    if(fHmapChannel->GetBinContent(chan) < 0.01){
      t->ResetStatus(AliESDtrack::kTOFout);
    }
  }
}

Float_t  
AliTOFT0maker::TuneForMC(AliESDEvent *esd){ // return true T0 event
  //
  // tune for MC data
  //

  Float_t TOFtimeResolutionDefault=80;

  Float_t t0 = gRandom->Gaus(0.,fT0width); 

  Float_t extraSmearing = 0;

  if(fTimeResolution > TOFtimeResolutionDefault){
    extraSmearing = TMath::Sqrt(fTimeResolution*fTimeResolution - TOFtimeResolutionDefault*TOFtimeResolutionDefault);
  }

  // subtruct t0 for each track
  Int_t ntracks = esd->GetNumberOfTracks();
  
  while (ntracks--) {
    AliESDtrack *t=esd->GetTrack(ntracks);
    
    if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
    
    /* check if channel is enabled */
    if (fTOFcalib){
	if(!fTOFcalib->IsChannelEnabled(t->GetTOFCalChannel())) {
	    /* reset TOF status */
	    t->ResetStatus(AliESDtrack::kTOFin);
	    t->ResetStatus(AliESDtrack::kTOFout);
	    t->ResetStatus(AliESDtrack::kTOFmismatch);
	    t->ResetStatus(AliESDtrack::kTOFpid);
	}
    }

    Double_t time=t->GetTOFsignal();

    time += t0;

    if(extraSmearing>0){
      Float_t smearing = gRandom->Gaus(0.,extraSmearing);
      time += smearing;
    }

    t->SetTOFsignal(time);
  }
  //
  return t0;
}
//_________________________________________________________________________
void     AliTOFT0maker::WriteInESD(AliESDEvent *esd){
  //
  // Write t0Gen, t0ResGen, nt0;
  //       t0resESD[0:nt0], it0ESD[0:nt0]
  // in the AliESDEvent
  //
    Int_t nMomBins = fPIDesd->GetTOFResponse().GetNmomBins();

    Int_t nt0=0;
    Float_t *t0 = new Float_t[nMomBins];
    Float_t *t0res = new Float_t[nMomBins];
    Int_t *multT0 = new Int_t[nMomBins];

    for(Int_t i=0;i<nMomBins;i++){
//	printf("START %i) %f %f\n",i,fT0event[i],fT0resolution[i]);
	multT0[i]=0;
	Bool_t kNewT0=kTRUE;
	for(Int_t j=0;j < nt0;j++){
	    if(TMath::Abs(fPIDesd->GetTOFResponse().GetT0bin(i) - t0[j])<0.1){
		kNewT0=kFALSE;
		multT0[j]++;
		j=nMomBins*10;
	    }
	}
	if(kNewT0){
	    t0[nt0]=fPIDesd->GetTOFResponse().GetT0bin(i);
	    t0res[nt0]=fPIDesd->GetTOFResponse().GetT0binRes(i);
	    nt0++;
	}
    }

    Int_t iMultT0=0,nmult=0;
    for(Int_t j=0;j < nt0;j++){ // find the most frequent
	if(multT0[j] > nmult){
	    iMultT0 = j;
	    nmult = multT0[j];
	}
    }

    Float_t *t0ESD = new Float_t[nMomBins];
    Float_t *t0resESD = new Float_t[nMomBins];
    Int_t *it0ESD = new Int_t[nMomBins];
    
    Float_t t0Gen,t0ResGen;
    t0Gen = t0[iMultT0];
    t0ResGen = t0res[iMultT0];
    nt0=0;
    //  printf("T0 to ESD\n%f %f\n",t0Gen,t0ResGen);
    for(Int_t i=0;i<nMomBins;i++){
	if(TMath::Abs(fPIDesd->GetTOFResponse().GetT0bin(i) - t0Gen)>0.1){
	    t0ESD[nt0]=fPIDesd->GetTOFResponse().GetT0bin(i);
	    t0resESD[nt0]=fPIDesd->GetTOFResponse().GetT0binRes(i);
	    it0ESD[nt0]=i;
//	    printf("%i) %f %f %i\n",nt0,t0ESD[nt0],t0resESD[nt0],it0ESD[nt0]);
	    nt0++;
	}
    }

    // Write  t0Gen,t0ResGen; nt0; t0resESD[0:nt0],it0ESD[0:nt0] in the AliESDEvent

   AliTOFHeader *tofHeader =
	new AliTOFHeader(t0Gen,t0ResGen,nt0,
			 t0ESD,t0resESD,it0ESD,fTimeResolution,fT0width);
  
    esd->SetTOFHeader(tofHeader);

    AliDebug(1,Form("resTOF=%f T0spread=%f t0Gen=%f t0resGen=%f",fTimeResolution,fT0width,t0Gen,t0ResGen));
    AliDebug(1,Form("%d ",nt0));
    for (Int_t ii=0; ii<nt0; ii++)
      AliDebug(1,Form("pBin=%d t0val=%f t0res=%f",it0ESD[ii],t0ESD[ii],t0resESD[ii]));
    
    delete[] t0;
    t0 = NULL;
    delete[] t0res;
    t0res = NULL;
    delete[] multT0;
    multT0 = NULL;
    delete[] t0ESD;
    t0ESD = NULL;
    delete[] t0resESD;
    t0resESD = NULL;
    delete[] it0ESD;
    it0ESD = NULL;

}
Commit	Line	Data
8f589502	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	**************************************************************************/
6a4e212e	15	/* $Id: AliTOFT0maker.cxx,v 1.8 2010/01/19 16:32:20 noferini Exp $ */
8f589502	16
	17	/////////////////////////////////////////////////////////////////////////////
	18	// //
	19	// This class contains the basic functions for the time zero //
6a4e212e	20	// evaluation with TOF detector informations. //
	21	// Use case in an analysis task: //
	22	// //
	23	// Create the object in the task constructor (fTOFmaker is a private var) //
5b4ed716	24	// AliESDpid *extPID=new AliESDpid(); //
	25	// fTOFmaker = new AliTOFT0maker(extPID); //
	26	// fTOFmaker->SetTimeResolution(100.0); // if you want set the TOF res //
6a4e212e	27	// 115 ps is the TOF default resolution value //
	28	// //
	29	// Use the RemakePID method in the task::Exec //
	30	// Double_t* calcolot0; //
	31	// calcolot0=fTOFmaker->RemakePID(fESD); //
	32	// //calcolot0[0] = calculated event time //
	33	// //calcolot0[1] = event time time resolution //
	34	// //calcolot0[2] = average event time for the current fill //
03bd764d	35	// //calcolot0[3] = tracks at TOF //
	36	// //calcolot0[4] = calculated event time (only TOF) //
	37	// //calcolot0[5] = event time time resolution (only TOF) //
	38	// //calcolot0[6] = sigma t0 fill //
	39	// //calcolot0[7] = tracks at TOF really used in tht algorithm //
6a4e212e	40	// //
	41	// Let consider that: //
	42	// - the PIF is automatically recalculated with the event time subtrction //
8f589502	43	// //
	44	/////////////////////////////////////////////////////////////////////////////
	45
536031f2	46	#include "AliTOFT0v1.h"
536031f2	47	#include "AliTOFT0maker.h"
536031f2	48	#include "AliPID.h"
f858b00e	49	#include "AliLog.h"
536031f2	50	#include "AliESDpid.h"
03bd764d	51	#include "AliESDEvent.h"
	52	#include "TFile.h"
	53	#include "TH1F.h"
5b4ed716	54	#include "AliTOFcalib.h"
	55	#include "AliTOFRunParams.h"
	56	#include "TRandom.h"
f858b00e	57	#include "AliTOFHeader.h"
536031f2	58
	59	ClassImp(AliTOFT0maker)
	60
	61	//____________________________________________________________________________
5b4ed716	62	AliTOFT0maker::AliTOFT0maker():
	63	TObject(),
	64	fT0TOF(NULL),
	65	fPIDesd(NULL),
	66	fExternalPIDFlag(kFALSE),
	67	fTOFcalib(NULL),
	68	fNoTOFT0(0),
2a258f40	69	fNmomBins(0),
5b4ed716	70	fTimeResolution(100),
	71	fT0sigma(1000),
	72	fHmapChannel(0),
	73	fKmask(0),
2a258f40	74	fT0width(150.),
	75	fT0spreadExt(-1.),
	76	fT0fillExt(0)
5b4ed716	77	{
	78	// ctr
	79	fCalculated[0] = 0;
	80	fCalculated[1] = 0;
	81	fCalculated[2] = 0;
	82	fCalculated[3] = 0;
	83
020dbfc2	84	fT0cur[0]=0.;
	85	fT0cur[1]=0.;
	86
5b4ed716	87	if(AliPID::ParticleMass(0) == 0) new AliPID();
	88
	89	fPIDesd = new AliESDpid();
	90
2a258f40	91	fNmomBins = fPIDesd->GetTOFResponse().GetNmomBins();
	92	SetTOFResponse();
	93
	94	fT0TOF = new AliTOFT0v1(fPIDesd);
	95
5b4ed716	96	}
	97	//____________________________________________________________________________
	98	AliTOFT0maker::AliTOFT0maker(AliESDpid externalPID, AliTOFcalib tofCalib):
	99	TObject(),
	100	fT0TOF(NULL),
	101	fPIDesd(externalPID),
	102	fExternalPIDFlag(kTRUE),
	103	fTOFcalib(tofCalib),
03bd764d	104	fNoTOFT0(0),
2a258f40	105	fNmomBins(0),
5b4ed716	106	fTimeResolution(100),
03bd764d	107	fT0sigma(1000),
03bd764d	108	fHmapChannel(0),
5b4ed716	109	fKmask(0),
2a258f40	110	fT0width(150.),
	111	fT0spreadExt(-1.),
	112	fT0fillExt(0)
536031f2	113	{
6a4e212e	114	// ctr
6a4e212e	115	fCalculated[0] = 0;
	116	fCalculated[1] = 0;
	117	fCalculated[2] = 0;
03bd764d	118	fCalculated[3] = 0;
6a4e212e	119
020dbfc2	120	fT0cur[0]=0.;
	121	fT0cur[1]=0.;
	122
536031f2	123	if(AliPID::ParticleMass(0) == 0) new AliPID();
03bd764d	124
5b4ed716	125	if(!fPIDesd){
	126	fPIDesd = new AliESDpid();
	127	fExternalPIDFlag = kFALSE;
	128	}
	129
2a258f40	130	fNmomBins = fPIDesd->GetTOFResponse().GetNmomBins();
2a258f40	131	SetTOFResponse();
5b4ed716	132
2a258f40	133	fT0TOF = new AliTOFT0v1(fPIDesd);
8f589502	134
8f589502	135	}
5b4ed716	136
8f589502	137	//____________________________________________________________________________
536031f2	138	AliTOFT0maker::~AliTOFT0maker()
	139	{
	140	// dtor
5b4ed716	141
	142	delete fT0TOF;
	143	if (!fExternalPIDFlag) delete fPIDesd;
536031f2	144	}
536031f2	145	//____________________________________________________________________________
5b4ed716	146	Double_t* AliTOFT0maker::ComputeT0TOF(AliESDEvent *esd,Double_t t0time,Double_t t0sigma){
8f589502	147	//
	148	// Remake TOF PID probabilities
	149	//
2a258f40	150	Double_t t0tof[6];
536031f2	151
03bd764d	152	if(fKmask) ApplyMask(esd);
03bd764d	153
2a258f40	154	Double_t t0fill = 0.;
	155
	156	fPIDesd->GetTOFResponse().ResetT0info();
	157
5b4ed716	158	/* get T0 spread from TOFcalib if available otherwise use default value */
	159	if (fTOFcalib && esd) {
	160	AliTOFRunParams *runParams = fTOFcalib->GetRunParams();
	161	if (runParams && runParams->GetTimestamp(0) != 0) {
	162	Float_t t0spread = runParams->EvalT0Spread(esd->GetTimeStamp());
2a258f40	163	if(fT0spreadExt > 0) SetT0FillWidth(fT0spreadExt);
	164	else{
	165	SetT0FillWidth(t0spread);
	166	t0fill = fT0fillExt;
	167	}
	168	}
	169	else{
90beec49	170	if(fT0spreadExt > 0) SetT0FillWidth(fT0spreadExt);
2a258f40	171	t0fill = fT0fillExt;
5b4ed716	172	}
536031f2	173	}
bc4fdef6	174	else if(esd){
	175	Float_t t0spread = esd->GetSigma2DiamondZ(); // vertex pread ^2
	176	if(t0spread > 0) t0spread = TMath::Sqrt(t0spread)/0.0299792458;
	177
	178	if(fT0spreadExt > 0) SetT0FillWidth(fT0spreadExt);
	179	else{
	180	SetT0FillWidth(t0spread);
	181	t0fill = fT0fillExt;
	182	}
	183	}
03bd764d	184
bf40b7fa	185	Float_t thrGood = TMath::Max(Float_t(500.),fT0width*3);
bf40b7fa	186
5b4ed716	187	fT0TOF->Init(esd);
5b4ed716	188	AliTOFT0v1* t0maker= fT0TOF;
5b4ed716	189
	190	t0maker->DefineT0("all",1.5,3.0);
	191	t0tof[0] = t0maker->GetResult(0);
	192	t0tof[1] = t0maker->GetResult(1);
	193	t0tof[2] = t0maker->GetResult(2);
	194	t0tof[3] = t0maker->GetResult(3);
2a258f40	195	t0tof[4] = t0maker->GetResult(4);
2a258f40	196	t0tof[5] = t0maker->GetResult(5);
536031f2	197
8f589502	198	Float_t lT0Current=0.;
536031f2	199	fT0sigma=1000;
536031f2	200
03bd764d	201	// Int_t nrun = esd->GetRunNumber();
5b4ed716	202
03bd764d	203	t0time += t0fill;
8f589502	204
5b4ed716	205	Float_t sigmaFill = fT0width;
	206
	207	if(sigmaFill < 20) sigmaFill = 140;
8f589502	208
03bd764d	209	fCalculated[0]=-1000*t0tof[0]; // best t0
	210	fCalculated[1]=1000*t0tof[1]; // sigma best t0
	211	fCalculated[2] = t0fill; //t0 fill
	212	fCalculated[3] = t0tof[2]; // n TOF tracks
	213	fCalculated[4]=-1000*t0tof[0]; // TOF t0
	214	fCalculated[5]=1000*t0tof[1]; // TOF t0 sigma
	215	fCalculated[6]=sigmaFill; // sigma t0 fill
	216	fCalculated[7] = t0tof[3]; // n TOF tracks used for T0
	217
bf40b7fa	218	if(fCalculated[7] > 30) thrGood = 10000000;
bf40b7fa	219
2a258f40	220	//statistics
	221	fCalculated[8] = t0tof[4]; // real time in s
	222	fCalculated[9] = t0tof[5]; // cpu time in s
	223
bf40b7fa	224	if(fCalculated[1] < sigmaFill && TMath::Abs(fCalculated[0] - t0fill) < thrGood && fCalculated[1] < fTimeResolution*1.2){
6a4e212e	225	fT0sigma=fCalculated[1];
6a4e212e	226	lT0Current=fCalculated[0];
536031f2	227	}
03bd764d	228	else{
	229	fCalculated[4] = t0fill;
	230	fCalculated[5] = sigmaFill;
	231	}
	232
5b4ed716	233	if(fCalculated[1] < 1 \|\| fT0sigma > sigmaFill \|\| fCalculated[1] > fTimeResolution* 1.2){
03bd764d	234	fT0sigma =1000;
	235	fCalculated[4] = t0fill;
	236	fCalculated[5] = sigmaFill;
	237	}
536031f2	238
	239	if(t0sigma < 1000){
	240	if(fT0sigma < 1000){
	241	Double_t w1 = 1./t0sigma/t0sigma;
6a4e212e	242	Double_t w2 = 1./fCalculated[1]/fCalculated[1];
536031f2	243
	244	Double_t wtot = w1+w2;
	245
6a4e212e	246	lT0Current = (w1t0time + w2fCalculated[0]) / wtot;
536031f2	247	fT0sigma = TMath::Sqrt(1./wtot);
	248	}
	249	else{
8f589502	250	lT0Current=t0time;
536031f2	251	fT0sigma=t0sigma;
	252	}
	253	}
	254
bf40b7fa	255	if(fT0sigma < sigmaFill && TMath::Abs(lT0Current - t0fill) < thrGood){
03bd764d	256	fCalculated[1]=fT0sigma;
	257	fCalculated[0]=lT0Current;
	258	}
	259
	260	if(fT0sigma >= 1000 \|\| fNoTOFT0){
8f589502	261	lT0Current = t0fill;
03bd764d	262	fT0sigma = sigmaFill;
8f589502	263
6a4e212e	264	fCalculated[0] = t0fill;
03bd764d	265	fCalculated[1] = sigmaFill;
536031f2	266	}
536031f2	267
5b4ed716	268	// T0 pt bin
721ed687	269	Float_t *t0values = new Float_t[fNmomBins];
721ed687	270	Float_t *t0resolution = new Float_t[fNmomBins];
2a258f40	271	if(fCalculated[7] < 100){
	272	for(Int_t i=0;i<fNmomBins;i++){
	273	t0maker->DefineT0("all",fPIDesd->GetTOFResponse().GetMinMom(i),fPIDesd->GetTOFResponse().GetMaxMom(i));
	274	t0tof[0] = t0maker->GetResult(0);
	275	t0tof[1] = t0maker->GetResult(1);
	276	t0tof[2] = t0maker->GetResult(2);
	277	t0tof[3] = t0maker->GetResult(3);
2a258f40	278
	279	Float_t t0bin =-1000*t0tof[0]; // best t0
	280	Float_t t0binRes =1000*t0tof[1]; // sigma best t0
	281
bf40b7fa	282	if(t0binRes < sigmaFill && t0binRes < fTimeResolution * 1.2 && TMath::Abs(t0bin - t0fill) < thrGood){
2a258f40	283	// Ok T0
	284	if(t0sigma < 1000){
	285	Double_t w1 = 1./t0sigma/t0sigma;
	286	Double_t w2 = 1./t0binRes/t0binRes;
	287
	288	Double_t wtot = w1+w2;
	289
	290	t0bin = (w1t0time + w2t0bin) / wtot;
	291	t0binRes = TMath::Sqrt(1./wtot);
	292	}
	293	}
	294	else{
	295	t0bin = t0fill;
	296	t0binRes = sigmaFill;
	297	if(t0sigma < 1000){
	298	t0bin = t0time;
	299	t0binRes= t0sigma;
	300	}
	301	}
721ed687	302	t0values[i] = t0bin;
721ed687	303	t0resolution[i] = t0binRes;
5b4ed716	304	}
2a258f40	305	}
	306	else{
	307	for(Int_t i=0;i<fNmomBins;i++){
721ed687	308	t0values[i] = lT0Current;
721ed687	309	t0resolution[i] = fT0sigma;
5b4ed716	310	}
5b4ed716	311	}
721ed687	312	for(Int_t i=0;i<fNmomBins;i++){
	313	fPIDesd->GetTOFResponse().SetT0bin(i,t0values[i]);
	314	fPIDesd->GetTOFResponse().SetT0binRes(i,t0resolution[i]);
	315	}
	316
	317	delete[] t0values;
	318	delete[] t0resolution;
536031f2	319
6a4e212e	320	return fCalculated;
536031f2	321	}
536031f2	322	//____________________________________________________________________________
5b4ed716	323	Double_t *AliTOFT0maker::GetT0p(Float_t p){// [0]=to -- [1] = sigma T0
2a258f40	324	Int_t i=fPIDesd->GetTOFResponse().GetMomBin(p);
03bd764d	325
2a258f40	326	fT0cur[0] = fPIDesd->GetTOFResponse().GetT0bin(i);
	327	fT0cur[1] = fPIDesd->GetTOFResponse().GetT0binRes(i);
	328
5b4ed716	329	return fT0cur;
536031f2	330	}
536031f2	331	//____________________________________________________________________________
5b4ed716	332	void AliTOFT0maker::SetTOFResponse(){
2a258f40	333	fPIDesd->GetTOFResponse().SetTimeResolution(fTimeResolution);
5b4ed716	334	}
	335	//____________________________________________________________________________
	336	Float_t AliTOFT0maker::GetExpectedSigma(Float_t mom, Float_t tof, Float_t mass){
5b4ed716	337	Float_t sigma = fPIDesd->GetTOFResponse().GetExpectedSigma(mom,tof,mass);
8f589502	338
5b4ed716	339	return sigma;
536031f2	340	}
536031f2	341	//____________________________________________________________________________
5b4ed716	342	void AliTOFT0maker::ApplyT0TOF(AliESDEvent *esd){
8f589502	343	//
5b4ed716	344	// Recalculate TOF PID probabilities
8f589502	345	//
8f589502	346
90beec49	347	// subtruct t0 for each track
90beec49	348	Int_t ntracks = esd->GetNumberOfTracks();
536031f2	349
90beec49	350	while (ntracks--) {
90beec49	351	AliESDtrack *t=esd->GetTrack(ntracks);
03bd764d	352
90beec49	353	if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
5b4ed716	354
90beec49	355	Double_t time=t->GetTOFsignal();
	356	Float_t p = t->GetP();
	357
	358	Double_t *t0=GetT0p(p);
	359	time -= t0[0];
	360	t->SetTOFsignal(time);
	361	}
	362
	363	for(Int_t i=0;i<fNmomBins;i++){
	364	fPIDesd->GetTOFResponse().SetT0bin(i,0.0);
	365	}
	366
5b4ed716	367	//
536031f2	368	}
03bd764d	369	//____________________________________________________________________________
	370	void AliTOFT0maker::LoadChannelMap(char *filename){
	371	// Load the histo with the channel off map
	372	TFile *f= new TFile(filename);
	373	if(!f){
	374	printf("Cannot open the channel map file (%s)\n",filename);
	375	return;
	376	}
	377
	378	fHmapChannel = (TH1F *) f->Get("hChEnabled");
	379
	380	if(!fHmapChannel){
	381	printf("Cannot laod the channel map histo (from %s)\n",filename);
	382	return;
	383	}
	384
	385	}
	386	//____________________________________________________________________________
	387	void AliTOFT0maker::ApplyMask(AliESDEvent * const esd){
	388	// Switch off the disable channel
	389	if(!fHmapChannel){
	390	printf("Channel Map is not available\n");
	391	return;
	392	}
	393
	394	Int_t ntracks = esd->GetNumberOfTracks();
	395
	396	while (ntracks--) {
	397	AliESDtrack *t=esd->GetTrack(ntracks);
	398
	399	if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
	400
	401	Int_t chan = t->GetTOFCalChannel();
	402
	403	if(fHmapChannel->GetBinContent(chan) < 0.01){
	404	t->ResetStatus(AliESDtrack::kTOFout);
	405	}
	406	}
	407	}
	408
5b4ed716	409	Float_t
	410	AliTOFT0maker::TuneForMC(AliESDEvent *esd){ // return true T0 event
	411	//
	412	// tune for MC data
	413	//
	414
	415	Float_t TOFtimeResolutionDefault=80;
	416
	417	Float_t t0 = gRandom->Gaus(0.,fT0width);
	418
	419	Float_t extraSmearing = 0;
	420
	421	if(fTimeResolution > TOFtimeResolutionDefault){
	422	extraSmearing = TMath::Sqrt(fTimeResolutionfTimeResolution - TOFtimeResolutionDefaultTOFtimeResolutionDefault);
	423	}
	424
	425	// subtruct t0 for each track
	426	Int_t ntracks = esd->GetNumberOfTracks();
	427
	428	while (ntracks--) {
	429	AliESDtrack *t=esd->GetTrack(ntracks);
	430
	431	if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
	432
	433	/* check if channel is enabled */
bc4fdef6	434	if (fTOFcalib){
	435	if(!fTOFcalib->IsChannelEnabled(t->GetTOFCalChannel())) {
	436	/* reset TOF status */
	437	t->ResetStatus(AliESDtrack::kTOFin);
	438	t->ResetStatus(AliESDtrack::kTOFout);
	439	t->ResetStatus(AliESDtrack::kTOFmismatch);
	440	t->ResetStatus(AliESDtrack::kTOFpid);
	441	}
5b4ed716	442	}
	443
	444	Double_t time=t->GetTOFsignal();
	445
	446	time += t0;
	447
	448	if(extraSmearing>0){
	449	Float_t smearing = gRandom->Gaus(0.,extraSmearing);
	450	time += smearing;
	451	}
	452
	453	t->SetTOFsignal(time);
	454	}
	455	//
	456	return t0;
	457	}
f858b00e	458	//_________________________________________________________________________
	459	void AliTOFT0maker::WriteInESD(AliESDEvent *esd){
	460	//
	461	// Write t0Gen, t0ResGen, nt0;
	462	// t0resESD[0:nt0], it0ESD[0:nt0]
	463	// in the AliESDEvent
	464	//
	465	Int_t nMomBins = fPIDesd->GetTOFResponse().GetNmomBins();
	466
	467	Int_t nt0=0;
	468	Float_t *t0 = new Float_t[nMomBins];
	469	Float_t *t0res = new Float_t[nMomBins];
	470	Int_t *multT0 = new Int_t[nMomBins];
	471
	472	for(Int_t i=0;i<nMomBins;i++){
	473	// printf("START %i) %f %f\n",i,fT0event[i],fT0resolution[i]);
	474	multT0[i]=0;
	475	Bool_t kNewT0=kTRUE;
	476	for(Int_t j=0;j < nt0;j++){
	477	if(TMath::Abs(fPIDesd->GetTOFResponse().GetT0bin(i) - t0[j])<0.1){
	478	kNewT0=kFALSE;
	479	multT0[j]++;
	480	j=nMomBins*10;
	481	}
	482	}
	483	if(kNewT0){
	484	t0[nt0]=fPIDesd->GetTOFResponse().GetT0bin(i);
	485	t0res[nt0]=fPIDesd->GetTOFResponse().GetT0binRes(i);
	486	nt0++;
	487	}
	488	}
	489
	490	Int_t iMultT0=0,nmult=0;
	491	for(Int_t j=0;j < nt0;j++){ // find the most frequent
	492	if(multT0[j] > nmult){
	493	iMultT0 = j;
	494	nmult = multT0[j];
	495	}
	496	}
	497
	498	Float_t *t0ESD = new Float_t[nMomBins];
	499	Float_t *t0resESD = new Float_t[nMomBins];
	500	Int_t *it0ESD = new Int_t[nMomBins];
	501
	502	Float_t t0Gen,t0ResGen;
	503	t0Gen = t0[iMultT0];
	504	t0ResGen = t0res[iMultT0];
	505	nt0=0;
	506	// printf("T0 to ESD\n%f %f\n",t0Gen,t0ResGen);
	507	for(Int_t i=0;i<nMomBins;i++){
	508	if(TMath::Abs(fPIDesd->GetTOFResponse().GetT0bin(i) - t0Gen)>0.1){
	509	t0ESD[nt0]=fPIDesd->GetTOFResponse().GetT0bin(i);
	510	t0resESD[nt0]=fPIDesd->GetTOFResponse().GetT0binRes(i);
	511	it0ESD[nt0]=i;
	512	// printf("%i) %f %f %i\n",nt0,t0ESD[nt0],t0resESD[nt0],it0ESD[nt0]);
	513	nt0++;
	514	}
	515	}
	516
	517	// Write t0Gen,t0ResGen; nt0; t0resESD[0:nt0],it0ESD[0:nt0] in the AliESDEvent
	518
	519	AliTOFHeader *tofHeader =
	520	new AliTOFHeader(t0Gen,t0ResGen,nt0,
c85a1b24	521	t0ESD,t0resESD,it0ESD,fTimeResolution,fT0width);
f858b00e	522
	523	esd->SetTOFHeader(tofHeader);
	524
c85a1b24	525	AliDebug(1,Form("resTOF=%f T0spread=%f t0Gen=%f t0resGen=%f",fTimeResolution,fT0width,t0Gen,t0ResGen));
c85a1b24	526	AliDebug(1,Form("%d ",nt0));
f858b00e	527	for (Int_t ii=0; ii<nt0; ii++)
c85a1b24	528	AliDebug(1,Form("pBin=%d t0val=%f t0res=%f",it0ESD[ii],t0ESD[ii],t0resESD[ii]));
f858b00e	529
	530	delete[] t0;
	531	t0 = NULL;
	532	delete[] t0res;
	533	t0res = NULL;
	534	delete[] multT0;
	535	multT0 = NULL;
	536	delete[] t0ESD;
	537	t0ESD = NULL;
	538	delete[] t0resESD;
	539	t0resESD = NULL;
	540	delete[] it0ESD;
	541	it0ESD = NULL;
	542
	543	}