[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),
  fTOFT0algorithm(1)
{
  // 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),
    fTOFT0algorithm(1)
{
  // 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;
  if (fTOFT0algorithm==2) t0maker->SetOptimization(kTRUE);
  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);

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