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