#include <TF1.h>
#include <TSpline.h>
#include <TFile.h>
+#include <TArrayF.h>
#include <AliVEvent.h>
#include <AliVTrack.h>
#include <AliLog.h>
#include <AliPID.h>
#include <AliOADBContainer.h>
-#include <AliTRDPIDReference.h>
+#include <AliTRDPIDResponseObject.h>
+#include <AliTOFPIDParams.h>
#include "AliPIDResponse.h"
+#include "AliCentrality.h"
+
ClassImp(AliPIDResponse);
AliPIDResponse::AliPIDResponse(Bool_t isMC/*=kFALSE*/) :
fOldRun(0),
fArrPidResponseMaster(0x0),
fResolutionCorrection(0x0),
-fTRDPIDParams(0x0),
-fTRDPIDReference(0x0),
-fTOFTimeZeroType(kBest_T0),
-fTOFres(100.),
-fCurrentEvent(0x0)
+fTRDPIDResponseObject(0x0),
+fTOFtail(1.1),
+fTOFPIDParams(0x0),
+fEMCALPIDParams(0x0),
+fCurrentEvent(0x0),
+fCurrCentrality(0.0)
{
//
// default ctor
//
- AliLog::SetClassDebugLevel("AliPIDResponse",10);
- AliLog::SetClassDebugLevel("AliESDpid",10);
- AliLog::SetClassDebugLevel("AliAODpidUtil",10);
+ AliLog::SetClassDebugLevel("AliPIDResponse",0);
+ AliLog::SetClassDebugLevel("AliESDpid",0);
+ AliLog::SetClassDebugLevel("AliAODpidUtil",0);
memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
}
//
// dtor
//
- delete fArrPidResponseMaster;
- delete fTRDPIDParams;
- delete fTRDPIDReference;
+ delete fArrPidResponseMaster;
+ delete fTRDPIDResponseObject;
+ if (fTOFPIDParams) delete fTOFPIDParams;
}
//______________________________________________________________________________
fOldRun(0),
fArrPidResponseMaster(0x0),
fResolutionCorrection(0x0),
-fTRDPIDParams(0x0),
-fTRDPIDReference(0x0),
-fTOFTimeZeroType(AliPIDResponse::kBest_T0),
-fTOFres(100.),
-fCurrentEvent(0x0)
+fTRDPIDResponseObject(0x0),
+fTOFtail(1.1),
+fTOFPIDParams(0x0),
+fEMCALPIDParams(0x0),
+fCurrentEvent(0x0),
+fCurrCentrality(0.0)
{
//
// copy ctor
fOldRun=0;
fArrPidResponseMaster=0x0;
fResolutionCorrection=0x0;
- fTRDPIDParams=0x0;
- fTRDPIDReference=0x0;
+ fTRDPIDResponseObject=0x0;
+ fEMCALPIDParams=0x0;
memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
- fTOFTimeZeroType=AliPIDResponse::kBest_T0;
- fTOFres=100.;
+ fTOFtail=1.1;
+ fTOFPIDParams=0x0;
fCurrentEvent=other.fCurrentEvent;
}
return *this;
// Track matching
nMatchClus = track->GetEMCALcluster();
if(nMatchClus > -1){
-
+
mom = track->P();
pt = track->Pt();
charge = track->Charge();
if(matchedClus){
- // matched cluster is EMCAL
- if(matchedClus->IsEMCAL()){
-
- fClsE = matchedClus->E();
- EovP = fClsE/mom;
-
-
- // NSigma value really meaningful only for electrons!
- return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge);
+ // matched cluster is EMCAL
+ if(matchedClus->IsEMCAL()){
+
+ fClsE = matchedClus->E();
+ EovP = fClsE/mom;
+
+
+ // NSigma value really meaningful only for electrons!
+ return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge);
+ }
}
}
- }
-
+
return -999;
}
+//______________________________________________________________________________
+Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVTrack *track, AliPID::EParticleType type, Double_t &eop, Double_t showershape[4]) const {
+
+ AliVCluster *matchedClus = NULL;
+
+ Double_t mom = -1.;
+ Double_t pt = -1.;
+ Double_t EovP = -1.;
+ Double_t fClsE = -1.;
+
+ Int_t nMatchClus = -1;
+ Int_t charge = 0;
+
+ // Track matching
+ nMatchClus = track->GetEMCALcluster();
+ if(nMatchClus > -1){
+
+ mom = track->P();
+ pt = track->Pt();
+ charge = track->Charge();
+
+ matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
+
+ if(matchedClus){
+
+ // matched cluster is EMCAL
+ if(matchedClus->IsEMCAL()){
+
+ fClsE = matchedClus->E();
+ EovP = fClsE/mom;
+
+ // fill used EMCAL variables here
+ eop = EovP; // E/p
+ showershape[0] = matchedClus->GetNCells(); // number of cells in cluster
+ showershape[1] = matchedClus->GetM02(); // long axis
+ showershape[2] = matchedClus->GetM20(); // short axis
+ showershape[3] = matchedClus->GetDispersion(); // dispersion
+
+ // NSigma value really meaningful only for electrons!
+ return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge);
+ }
+ }
+ }
+ return -999;
+}
+
+
//______________________________________________________________________________
AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability (EDetCode detCode, const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
{
// Compute PID response for the
//
+ Double_t meanCorrFactor = 0.11/fTOFtail; // Correction factor on the mean because of the tail (should be ~ 0.1 with tail = 1.1)
+
// set flat distribution (no decision)
for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
Double_t time[AliPID::kSPECIESN];
track->GetIntegratedTimes(time);
- Double_t sigma[AliPID::kSPECIES];
- for (Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++) {
- sigma[iPart] = fTOFResponse.GetExpectedSigma(track->P(),time[iPart],AliPID::ParticleMass(iPart));
- }
+ // Double_t sigma[AliPID::kSPECIES];
+ // for (Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++) {
+ // sigma[iPart] = fTOFResponse.GetExpectedSigma(track->P(),time[iPart],AliPID::ParticleMass(iPart));
+ // }
Bool_t mismatch = kTRUE/*, heavy = kTRUE*/;
for (Int_t j=0; j<AliPID::kSPECIES; j++) {
AliPID::EParticleType type=AliPID::EParticleType(j);
- Double_t nsigmas=NumberOfSigmasTOF(track,type);
-
- Double_t sig = sigma[j];
+ Double_t nsigmas=NumberOfSigmasTOF(track,type) + meanCorrFactor;
+
+ // Double_t sig = sigma[j];
+ Double_t expTime = fTOFResponse.GetExpectedSignal(track,type);
+ Double_t sig = fTOFResponse.GetExpectedSigma(track->P(),expTime,AliPID::ParticleMassZ(type));
+ if (TMath::Abs(nsigmas) > (fRange+2)) {
+ if(nsigmas < fTOFtail)
+ p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
+ else
+ p[j] = TMath::Exp(-(fRange+2 - fTOFtail*0.5)*fTOFtail)/sig;
+ } else{
+ if(nsigmas < fTOFtail)
+ p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
+ else
+ p[j] = TMath::Exp(-(nsigmas - fTOFtail*0.5)*fTOFtail)/sig;
+ }
+
+ /* OLD Gaussian shape
if (TMath::Abs(nsigmas) > (fRange+2)) {
p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
} else
p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
+ */
if (TMath::Abs(nsigmas)<5.){
Double_t nsigmasTPC=NumberOfSigmasTPC(track,type);
}
//TOF resolution
- SetTOFResponse(event, (AliPIDResponse::EStartTimeType_t)fTOFTimeZeroType);
-
+ SetTOFResponse(event, (AliPIDResponse::EStartTimeType_t)fTOFPIDParams->GetStartTimeMethod());
+
+
+ // Get and set centrality
+ AliCentrality *centrality = event->GetCentrality();
+ if(centrality){
+ fCurrCentrality = centrality->GetCentralityPercentile("V0M");
+ }
+ else{
+ fCurrCentrality = -1;
+ }
}
//______________________________________________________________________________
SetTRDPidResponseMaster();
InitializeTRDResponse();
+
+ SetEMCALPidResponseMaster();
+ InitializeEMCALResponse();
- fTOFResponse.SetTimeResolution(fTOFres);
+ SetTOFPidResponseMaster();
+ InitializeTOFResponse();
}
//_____________________________________________________
fBeamType="PP";
- TPRegexp reg(".*(LHC11[a-z]+[0-9]+[a-z_]*)/.*");
+
+ TPRegexp reg(".*(LHC1[1-2][a-z]+[0-9]+[a-z_]*)/.*");
+ TPRegexp reg12a17(".*(LHC12a17[a-z]+)/.*");
+
//find the period by run number (UGLY, but not stored in ESD and AOD... )
if (fRun>=114737&&fRun<=117223) { fLHCperiod="LHC10B"; fMCperiodTPC="LHC10D1"; }
else if (fRun>=118503&&fRun<=121040) { fLHCperiod="LHC10C"; fMCperiodTPC="LHC10D1"; }
else if (fRun>=122195&&fRun<=126437) { fLHCperiod="LHC10D"; fMCperiodTPC="LHC10F6A"; }
- else if (fRun>=127719&&fRun<=130850) { fLHCperiod="LHC10E"; fMCperiodTPC="LHC10F6A"; }
+ else if (fRun>=127710&&fRun<=130850) { fLHCperiod="LHC10E"; fMCperiodTPC="LHC10F6A"; }
else if (fRun>=133004&&fRun<=135029) { fLHCperiod="LHC10F"; fMCperiodTPC="LHC10F6A"; }
else if (fRun>=135654&&fRun<=136377) { fLHCperiod="LHC10G"; fMCperiodTPC="LHC10F6A"; }
else if (fRun>=136851&&fRun<=139517) {
if (reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11A10";
fBeamType="PBPB";
}
- else if (fRun>=139699) { fLHCperiod="LHC11A"; fMCperiodTPC="LHC10F6A"; }
+ else if (fRun>=139699&&fRun<=146860) { fLHCperiod="LHC11A"; fMCperiodTPC="LHC10F6A"; }
+ //TODO: periods 11B, 11C are not yet treated assume 11d for the moment
+ else if (fRun>=148531&&fRun<=155384) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
+ else if (fRun>=156477&&fRun<=159635) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
+ else if (fRun>=166529) {
+ fLHCperiod="LHC11H";
+ fMCperiodTPC="LHC11A10";
+ fBeamType="PBPB";
+ if (reg12a17.MatchB(fCurrentFile)) fMCperiodTPC="LHC12A17";
+ }
+
//exception new pp MC productions from 2011
if (fBeamType=="PP" && reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11B2";
+ // exception for 11f1
+ if (fCurrentFile.Contains("LHC11f1/")) fMCperiodTPC="LHC11F1";
}
//______________________________________________________________________________
for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,0x0);
}
-
+
+ // period
+ TString period=fLHCperiod;
+ if (fIsMC) period=fMCperiodTPC;
+
+ AliInfo(Form("Searching splines for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
+ Bool_t found=kFALSE;
//
//set the new PID splines
//
- TString period=fLHCperiod;
if (fArrPidResponseMaster){
TObject *grAll=0x0;
//for MC don't use period information
// if (fIsMC) period="[A-Z0-9]*";
//for MC use MC period information
- if (fIsMC) period=fMCperiodTPC;
//pattern for the default entry (valid for all particles)
TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,responseFunction);
fTPCResponse.SetUseDatabase(kTRUE);
AliInfo(Form("Adding graph: %d - %s",ispec,responseFunction->GetName()));
+ found=kTRUE;
break;
}
}
if (!fTPCResponse.GetResponseFunction((AliPID::EParticleType)ispec)){
fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,grAll);
AliInfo(Form("Adding graph: %d - %s",ispec,grAll->GetName()));
+ found=kTRUE;
}
}
}
}
-
+
+ if (!found){
+ AliError(Form("No splines found for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
+ }
//
// Setup resolution parametrisation
//
//
// Load the TRD pid params and references from the OADB
//
- if(fTRDPIDParams) return;
+ if(fTRDPIDResponseObject) return;
AliOADBContainer contParams("contParams");
- contParams.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDParams.root", fOADBPath.Data()), "AliTRDPIDParams");
- fTRDPIDParams = (TObjArray *)contParams.GetObject(fRun);
-
+ Int_t statusResponse = contParams.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()), "AliTRDPIDResponseObject");
+ if(statusResponse){
+ AliError("Failed initializing PID Response Object from OADB");
+ } else {
+ AliInfo(Form("Loading TRD Response from %s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()));
+ fTRDPIDResponseObject = dynamic_cast<AliTRDPIDResponseObject *>(contParams.GetObject(fRun));
+ if(!fTRDPIDResponseObject){
+ AliError(Form("TRD Response not found in run %d", fRun));
+ }
+ }
+ /*
AliOADBContainer contRefs("contRefs");
- contRefs.InitFromFile(Form("%s/COMMON/PID/dReferencesLQ1D.root", fOADBPath.Data()), "AliTRDPIDReference");
- fTRDPIDReference = (AliTRDPIDReference *)contRefs.GetObject(fRun);
+ Int_t statusRefs = contRefs.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDReferenceLQ1D.root", fOADBPath.Data()), "AliTRDPIDReference");
+ if(statusRefs){
+ AliInfo("Failed Loading References for TRD");
+ } else {
+ AliInfo(Form("Loading TRD References from %s/COMMON/PID/data/TRDPIDReferenceLQ1D.root", fOADBPath.Data()));
+ fTRDPIDReference = dynamic_cast<AliTRDPIDReference *>(contRefs.GetObject(fRun));
+ if(!fTRDPIDReference){
+ AliError(Form("TRD References not found in OADB Container for run %d", fRun));
+ }
+ }
+ */
}
//______________________________________________________________________________
//
// Set PID Params and references to the TRD PID response
//
- fTRDResponse.SetPIDParams(fTRDPIDParams);
- fTRDResponse.Load(fTRDPIDReference);
- if(fLHCperiod == "LHC10b" || fLHCperiod == "LHC10c" || fLHCperiod == "LHC10d" || fLHCperiod == "LHC10e"){
- fTRDslicesForPID[0] = 0;
- fTRDslicesForPID[1] = 7;
+ fTRDResponse.SetPIDResponseObject(fTRDPIDResponseObject);
+ if(fLHCperiod == "LHC10b" || fLHCperiod == "LHC10c" || fLHCperiod == "LHC10d" || fLHCperiod == "LHC10e"){
+ fTRDslicesForPID[0] = 0;
+ fTRDslicesForPID[1] = 7;
+ }
+}
+
+//______________________________________________________________________________
+void AliPIDResponse::SetTOFPidResponseMaster()
+{
+ //
+ // Load the TOF pid params from the OADB
+ //
+ TFile *oadbf = new TFile(Form("%s/COMMON/PID/data/TOFPIDParams.root",fOADBPath.Data()));
+ if (oadbf->IsOpen()) {
+ AliInfo(Form("Loading TOF Params from %s/COMMON/PID/data/TOFPIDParams.root", fOADBPath.Data()));
+ AliOADBContainer *oadbc = (AliOADBContainer *)oadbf->Get("TOFoadb");
+ if (fTOFPIDParams) delete fTOFPIDParams;
+ fTOFPIDParams = dynamic_cast<AliTOFPIDParams *>(oadbc->GetObject(fRun,"TOFparams"));
+ oadbf->Close();
+ delete oadbc;
+ } else {
+ AliError(Form("TOFPIDParams.root not found in %s/COMMON/PID/data !!",fOADBPath.Data()));
+ }
+ delete oadbf;
+
+ }
+
+//______________________________________________________________________________
+void AliPIDResponse::InitializeTOFResponse(){
+ //
+ // Set PID Params to the TOF PID response
+ //
+ for (Int_t i=0;i<4;i++) {
+ fTOFResponse.SetTrackParameter(i,fTOFPIDParams->GetSigParams(i));
}
+ fTOFResponse.SetTimeResolution(fTOFPIDParams->GetTOFresolution());
+
}
+
//_________________________________________________________________________
Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack, Double_t efficiencyLevel) const {
//
Double_t probs[AliPID::kSPECIES];
ComputeTRDProbability(vtrack, AliPID::kSPECIES, probs);
- Int_t ntracklets=0;
- Double_t p = 0;
+ Int_t ntracklets = vtrack->GetTRDntrackletsPID();
+ // Take mean of the TRD momenta in the given tracklets
+ Float_t p = 0, trdmomenta[AliVTrack::kTRDnPlanes];
+ Int_t nmomenta = 0;
for(Int_t iPl=0;iPl<AliVTrack::kTRDnPlanes;iPl++){
if(vtrack->GetTRDmomentum(iPl) > 0.){
- ntracklets++;
- p = vtrack->GetTRDmomentum(iPl);
+ trdmomenta[nmomenta++] = vtrack->GetTRDmomentum(iPl);
}
}
+ p = TMath::Mean(nmomenta, trdmomenta);
return fTRDResponse.IdentifiedAsElectron(ntracklets, probs, p, efficiencyLevel);
}
+//______________________________________________________________________________
+void AliPIDResponse::SetEMCALPidResponseMaster()
+{
+ //
+ // Load the EMCAL pid response functions from the OADB
+ //
+ TObjArray* fEMCALPIDParamsRun = NULL;
+ TObjArray* fEMCALPIDParamsPass = NULL;
+
+ if(fEMCALPIDParams) return;
+ AliOADBContainer contParams("contParams");
+
+ Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()), "AliEMCALPIDParams");
+ if(statusPars){
+ AliError("Failed initializing PID Params from OADB");
+ }
+ else {
+ AliInfo(Form("Loading EMCAL Params from %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
+
+ fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(fRun));
+ if(fEMCALPIDParamsRun) fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",fRecoPass)));
+ if(fEMCALPIDParamsPass) fEMCALPIDParams = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
+
+ if(!fEMCALPIDParams){
+ AliInfo(Form("EMCAL Params not found in run %d pass %d", fRun, fRecoPass));
+ AliInfo("Will take the standard LHC11d instead ...");
+
+ fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(156477));
+ if(fEMCALPIDParamsRun) fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",1)));
+ if(fEMCALPIDParamsPass) fEMCALPIDParams = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
+
+ if(!fEMCALPIDParams){
+ AliError(Form("DEFAULT EMCAL Params (LHC11d) not found in file %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
+ }
+ }
+ }
+}
+
+//______________________________________________________________________________
+void AliPIDResponse::InitializeEMCALResponse(){
+ //
+ // Set PID Params to the EMCAL PID response
+ //
+ fEMCALResponse.SetPIDParams(fEMCALPIDParams);
+
+}
+//_________________________________________________________________________
+void AliPIDResponse::SetTOFResponse(AliVEvent *vevent,EStartTimeType_t option){
+ //
+ // Set TOF response function
+ // Input option for event_time used
+ //
+
+ Float_t t0spread = 0.; //vevent->GetEventTimeSpread();
+ if(t0spread < 10) t0spread = 80;
+
+ // T0 from TOF algorithm
+
+ Bool_t flagT0TOF=kFALSE;
+ Bool_t flagT0T0=kFALSE;
+ Float_t *startTime = new Float_t[fTOFResponse.GetNmomBins()];
+ Float_t *startTimeRes = new Float_t[fTOFResponse.GetNmomBins()];
+ Int_t *startTimeMask = new Int_t[fTOFResponse.GetNmomBins()];
+
+ // T0-TOF arrays
+ Float_t *estimatedT0event = new Float_t[fTOFResponse.GetNmomBins()];
+ Float_t *estimatedT0resolution = new Float_t[fTOFResponse.GetNmomBins()];
+ for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
+ estimatedT0event[i]=0.0;
+ estimatedT0resolution[i]=0.0;
+ startTimeMask[i] = 0;
+ }
+
+ Float_t resT0A=75,resT0C=65,resT0AC=55;
+ if(vevent->GetT0TOF()){ // check if T0 detector information is available
+ flagT0T0=kTRUE;
+ }
+
+
+ AliTOFHeader *tofHeader = (AliTOFHeader*)vevent->GetTOFHeader();
+
+ if (tofHeader) { // read global info and T0-TOF
+ fTOFResponse.SetTimeResolution(tofHeader->GetTOFResolution());
+ t0spread = tofHeader->GetT0spread(); // read t0 sprad
+ if(t0spread < 10) t0spread = 80;
+
+ flagT0TOF=kTRUE;
+ for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){ // read T0-TOF default value
+ startTime[i]=tofHeader->GetDefaultEventTimeVal();
+ startTimeRes[i]=tofHeader->GetDefaultEventTimeRes();
+ if(startTimeRes[i] < 1.e-5) startTimeRes[i] = t0spread;
+ }
+
+ TArrayI *ibin=(TArrayI*)tofHeader->GetNvalues();
+ TArrayF *t0Bin=(TArrayF*)tofHeader->GetEventTimeValues();
+ TArrayF *t0ResBin=(TArrayF*)tofHeader->GetEventTimeRes();
+ for(Int_t j=0;j < tofHeader->GetNbins();j++){ // fill T0-TOF in p-bins
+ Int_t icurrent = (Int_t)ibin->GetAt(j);
+ startTime[icurrent]=t0Bin->GetAt(j);
+ startTimeRes[icurrent]=t0ResBin->GetAt(j);
+ if(startTimeRes[icurrent] < 1.e-5) startTimeRes[icurrent] = t0spread;
+ }
+ }
+
+ // for cut of 3 sigma on t0 spread
+ Float_t t0cut = 3 * t0spread;
+ if(t0cut < 500) t0cut = 500;
+
+ if(option == kFILL_T0){ // T0-FILL is used
+ for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
+ estimatedT0event[i]=0.0;
+ estimatedT0resolution[i]=t0spread;
+ }
+ fTOFResponse.SetT0event(estimatedT0event);
+ fTOFResponse.SetT0resolution(estimatedT0resolution);
+ }
+
+ if(option == kTOF_T0){ // T0-TOF is used when available (T0-FILL otherwise) from ESD
+ if(flagT0TOF){
+ fTOFResponse.SetT0event(startTime);
+ fTOFResponse.SetT0resolution(startTimeRes);
+ for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
+ if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
+ fTOFResponse.SetT0binMask(i,startTimeMask[i]);
+ }
+ }
+ else{
+ for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
+ estimatedT0event[i]=0.0;
+ estimatedT0resolution[i]=t0spread;
+ fTOFResponse.SetT0binMask(i,startTimeMask[i]);
+ }
+ fTOFResponse.SetT0event(estimatedT0event);
+ fTOFResponse.SetT0resolution(estimatedT0resolution);
+ }
+ }
+ else if(option == kBest_T0){ // T0-T0 or T0-TOF are used when available (T0-FILL otherwise) from ESD
+ Float_t t0AC=-10000;
+ Float_t t0A=-10000;
+ Float_t t0C=-10000;
+ if(flagT0T0){
+ t0AC= vevent->GetT0TOF()[0];
+ t0A= vevent->GetT0TOF()[1];
+ t0C= vevent->GetT0TOF()[2];
+ }
+
+ Float_t t0t0Best = 0;
+ Float_t t0t0BestRes = 9999;
+ Int_t t0used=0;
+ if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
+ t0t0Best = t0AC;
+ t0t0BestRes = resT0AC;
+ t0used=6;
+ }
+ else if(TMath::Abs(t0C) < t0cut){
+ t0t0Best = t0C;
+ t0t0BestRes = resT0C;
+ t0used=4;
+ }
+ else if(TMath::Abs(t0A) < t0cut){
+ t0t0Best = t0A;
+ t0t0BestRes = resT0A;
+ t0used=2;
+ }
+
+ if(flagT0TOF){ // if T0-TOF info is available
+ for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
+ if(t0t0BestRes < 999){
+ if(startTimeRes[i] < t0spread){
+ Double_t wtot = 1./startTimeRes[i]/startTimeRes[i] + 1./t0t0BestRes/t0t0BestRes;
+ Double_t t0best = startTime[i]/startTimeRes[i]/startTimeRes[i] + t0t0Best/t0t0BestRes/t0t0BestRes;
+ estimatedT0event[i]=t0best / wtot;
+ estimatedT0resolution[i]=1./TMath::Sqrt(wtot);
+ startTimeMask[i] = t0used+1;
+ }
+ else {
+ estimatedT0event[i]=t0t0Best;
+ estimatedT0resolution[i]=t0t0BestRes;
+ startTimeMask[i] = t0used;
+ }
+ }
+ else{
+ estimatedT0event[i]=startTime[i];
+ estimatedT0resolution[i]=startTimeRes[i];
+ if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
+ }
+ fTOFResponse.SetT0binMask(i,startTimeMask[i]);
+ }
+ fTOFResponse.SetT0event(estimatedT0event);
+ fTOFResponse.SetT0resolution(estimatedT0resolution);
+ }
+ else{ // if no T0-TOF info is available
+ for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
+ fTOFResponse.SetT0binMask(i,t0used);
+ if(t0t0BestRes < 999){
+ estimatedT0event[i]=t0t0Best;
+ estimatedT0resolution[i]=t0t0BestRes;
+ }
+ else{
+ estimatedT0event[i]=0.0;
+ estimatedT0resolution[i]=t0spread;
+ }
+ }
+ fTOFResponse.SetT0event(estimatedT0event);
+ fTOFResponse.SetT0resolution(estimatedT0resolution);
+ }
+ }
+
+ else if(option == kT0_T0){ // T0-T0 is used when available (T0-FILL otherwise)
+ Float_t t0AC=-10000;
+ Float_t t0A=-10000;
+ Float_t t0C=-10000;
+ if(flagT0T0){
+ t0AC= vevent->GetT0TOF()[0];
+ t0A= vevent->GetT0TOF()[1];
+ t0C= vevent->GetT0TOF()[2];
+ }
+
+ if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
+ for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
+ estimatedT0event[i]=t0AC;
+ estimatedT0resolution[i]=resT0AC;
+ fTOFResponse.SetT0binMask(i,6);
+ }
+ }
+ else if(TMath::Abs(t0C) < t0cut){
+ for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
+ estimatedT0event[i]=t0C;
+ estimatedT0resolution[i]=resT0C;
+ fTOFResponse.SetT0binMask(i,4);
+ }
+ }
+ else if(TMath::Abs(t0A) < t0cut){
+ for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
+ estimatedT0event[i]=t0A;
+ estimatedT0resolution[i]=resT0A;
+ fTOFResponse.SetT0binMask(i,2);
+ }
+ }
+ else{
+ for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
+ estimatedT0event[i]=0.0;
+ estimatedT0resolution[i]=t0spread;
+ fTOFResponse.SetT0binMask(i,0);
+ }
+ }
+ fTOFResponse.SetT0event(estimatedT0event);
+ fTOFResponse.SetT0resolution(estimatedT0resolution);
+ }
+ delete [] startTime;
+ delete [] startTimeRes;
+ delete [] startTimeMask;
+ delete [] estimatedT0event;
+ delete [] estimatedT0resolution;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff