// //
// Create the object in the task constructor (fTOFmaker is a private var) //
// fTOFmaker = new AliTOFT0maker(); //
-// fTOFmaker->SetTimeResolution(115.0e-12); // if you want set the TOF res //
+// fTOFmaker->SetTimeResolution(130.0); // if you want set the TOF res //
// 115 ps is the TOF default resolution value //
// //
// Use the RemakePID method in the task::Exec //
// //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 <Riostream.h>
-#include <stdlib.h>
-
#include "AliTOFT0v1.h"
#include "AliTOFT0maker.h"
#include "AliTOFcalibHisto.h"
#include "AliPID.h"
#include "AliESDpid.h"
+#include "AliESDEvent.h"
+#include "TFile.h"
+#include "TH1F.h"
ClassImp(AliTOFT0maker)
//____________________________________________________________________________
-AliTOFT0maker::AliTOFT0maker() :
-TObject(),
- fCalib(new AliTOFcalibHisto()),
- fESDswitch(0),
- fTimeResolution(115),
- fT0sigma(1000)
+ AliTOFT0maker::AliTOFT0maker():
+ fCalib(new AliTOFcalibHisto()),
+ fnT0(0),
+ fiT0(0),
+ fNoTOFT0(0),
+ fESDswitch(0),
+ fTimeResolution(115),
+ fT0sigma(1000),
+ fHmapChannel(0),
+ fKmask(0)
{
- //
// ctr
- //
-
fCalculated[0] = 0;
fCalculated[1] = 0;
fCalculated[2] = 0;
+ fCalculated[3] = 0;
+ // fCalib->SetCalibParFileName("./AliTOFcalibPar.LHC10b.7000GeV.20100405.root");
fCalib->LoadCalibPar();
if(AliPID::ParticleMass(0) == 0) new AliPID();
+
+ SetESDdata();
}
//____________________________________________________________________________
AliTOFT0maker::AliTOFT0maker(const AliTOFT0maker & t) :
-TObject(),
+ TObject(),
fCalib(t.fCalib),
+ fnT0(t.fnT0),
+ fiT0(t.fiT0),
+ fNoTOFT0(t.fNoTOFT0),
fESDswitch(t.fESDswitch),
fTimeResolution(t.fTimeResolution),
- fT0sigma(t.fT0sigma)
+ fT0sigma(t.fT0sigma),
+ fHmapChannel(t.fHmapChannel),
+ fKmask(t.fKmask)
{
+ // copy ctr
}
//____________________________________________________________________________
Double_t *t0tof;
+ if(fKmask) ApplyMask(esd);
+
AliTOFT0v1* t0maker=new AliTOFT0v1(esd);
- t0maker->SetCalib(fCalib);
+// t0maker->SetCalib(fCalib);
t0maker->SetTimeResolution(fTimeResolution*1e-12);
if(! fESDswitch){
- t0tof=t0maker->DefineT0RawCorrection("all");
TakeTimeRawCorrection(esd);
}
- else t0tof=t0maker->DefineT0("all");
+
+ t0tof=t0maker->DefineT0("all");
Float_t lT0Current=0.;
fT0sigma=1000;
- Int_t nrun = esd->GetRunNumber();
- Double_t t0fill = GetT0Fill(nrun);
+// Int_t nrun = esd->GetRunNumber();
+ Double_t t0fill = GetT0Fill();
+ t0time += t0fill;
- fCalculated[0]=-1000*t0tof[0];
- fCalculated[1]=1000*t0tof[1];
- fCalculated[2] = t0fill;
+ Float_t sigmaFill = (t0fill - Int_t(t0fill))*1000;
+ if(sigmaFill < 0) sigmaFill += 1000;
- if(fCalculated[1] < 150 && TMath::Abs(fCalculated[0] - t0fill) < 500){
+ 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[1] < sigmaFill){
+ if(fnT0 < 10){
+ fT0fill[fiT0] = fCalculated[0];
+ fT0sigmaTOF[fiT0] = fCalculated[1];
+ fiT0++;
+ fnT0++;
+ }
+ else if(TMath::Abs(fCalculated[0] - t0fill) < 500){
+ fT0fill[fiT0] = fCalculated[0];
+ fT0sigmaTOF[fiT0] = fCalculated[1];
+ fiT0++;
+ fnT0++;
+ }
+
+ // printf("%i - %i) %f\n",fiT0,fnT0,t0fill);
+ }
+ if(fnT0==10) fiT0=0;
+
+ if(fiT0 > fgkNmaxT0step-1) fiT0=0;
+
+ if(fnT0 < 100){
+ t0time -= t0fill;
+ sigmaFill=200;
+ t0fill=0;
+ fCalculated[2] = t0fill; //t0 fill
+ }
+
+ if(fCalculated[1] < sigmaFill && TMath::Abs(fCalculated[0] - t0fill) < 500){
fT0sigma=fCalculated[1];
lT0Current=fCalculated[0];
}
+ else{
+ fCalculated[4] = t0fill;
+ fCalculated[5] = sigmaFill;
+ }
+
+ if(fCalculated[1] < 1 || fT0sigma > sigmaFill){
+ fT0sigma =1000;
+ fCalculated[4] = t0fill;
+ fCalculated[5] = sigmaFill;
+ }
if(t0sigma < 1000){
if(fT0sigma < 1000){
}
}
- if(fT0sigma >= 1000){
+ if(fT0sigma < sigmaFill && TMath::Abs(lT0Current - t0fill) < 500){
+ fCalculated[1]=fT0sigma;
+ fCalculated[0]=lT0Current;
+ }
+
+ if(fT0sigma >= 1000 || fNoTOFT0){
lT0Current = t0fill;
- fT0sigma = 135;
+ fT0sigma = sigmaFill;
fCalculated[0] = t0fill;
- fCalculated[1] = 150;
+ fCalculated[1] = sigmaFill;
}
- RemakeTOFpid(esd,lT0Current);
+
+ RemakeTOFpid(esd,lT0Current);
+
return fCalculated;
}
//____________________________________________________________________________
//
// Take raw corrections for time measurements
//
-
+
Int_t ntracks = esd->GetNumberOfTracks();
-
+
while (ntracks--) {
AliESDtrack *t=esd->GetTrack(ntracks);
if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
Double_t time=t->GetTOFsignalRaw();
+
Double_t tot = t->GetTOFsignalToT();
Int_t chan = t->GetTOFCalChannel();
Double_t corr = fCalib->GetFullCorrection(chan,tot) - fCalib->GetCorrection(AliTOFcalibHisto::kTimeSlewingCorr,chan,0);
time -= corr*1000;
- Int_t crate = Int_t(fCalib->GetCalibMap(AliTOFcalibHisto::kDDL,chan));
-
- if(crate == 63 || crate == 62){
- time += 9200;
- }
-
t->SetTOFsignal(time);
+
}
}
//____________________________________________________________________________
}
//____________________________________________________________________________
-Double_t AliTOFT0maker::GetT0Fill(Int_t nrun) const {
+Double_t AliTOFT0maker::GetT0Fill() const {
//
// Return T0 of filling
//
+
+ Double_t t0=0.200;
+
+ Int_t n=fnT0;
- Double_t t0;
- if(nrun==104065) t0= 1771614;
- else if(nrun==104068) t0= 1771603;
- else if(nrun==104070) t0= 1771594;
- else if(nrun==104073) t0= 1771610;
- else if(nrun==104080) t0= 1771305;
- else if(nrun==104083) t0= 1771613;
- else if(nrun==104157) t0= 1771665;
- else if(nrun==104159) t0= 1771679;
- else if(nrun==104160) t0= 1771633;
- else if(nrun==104316) t0= 1764344;
- else if(nrun==104320) t0= 1764342;
- else if(nrun==104321) t0= 1764371;
- else if(nrun==104439) t0= 1771750;
- else if(nrun==104792) t0= 1771755;
- else if(nrun==104793) t0= 1771762;
- else if(nrun==104799) t0= 1771828;
- else if(nrun==104800) t0= 1771788;
- else if(nrun==104801) t0= 1771796;
- else if(nrun==104802) t0= 1771775;
- else if(nrun==104803) t0= 1771795;
- else if(nrun==104824) t0= 1771751;
- else if(nrun==104825) t0= 1771763;
- else if(nrun==104845) t0= 1771792;
- else if(nrun==104852) t0= 1771817;
- else if(nrun==104864) t0= 1771825;
- else if(nrun==104865) t0= 1771827;
- else if(nrun==104867) t0= 1771841;
- else if(nrun==104876) t0= 1771856;
- else if(nrun==104878) t0= 1771847;
- else if(nrun==104879) t0= 1771830;
- else if(nrun==104892) t0= 1771837;
- else t0= 487;
-
- if(fESDswitch) t0 -= 487;
+ if(n >10 && n <= 20) n = 10;
+ else if(n > 20){
+ n -= 10;
+ }
+ if(n > fgkNmaxT0step) n = fgkNmaxT0step;
+
+ if(n>1){
+ Double_t lT0av=0;
+ Double_t lT0sigmaav=0;
+ Double_t lT0avErr=0;
+ for(Int_t i=0;i<n;i++){
+ lT0av+=fT0fill[i];
+ lT0sigmaav += fT0sigmaTOF[fiT0];
+ lT0avErr+=fT0fill[i]*fT0fill[i];
+ }
+ lT0avErr -= lT0av*lT0av/n;
+ lT0av /= n;
+ lT0sigmaav /= n;
+ lT0avErr = TMath::Sqrt(TMath::Max(lT0avErr/(n-1) - lT0sigmaav*lT0sigmaav,0.00001));
+
+
+ if(lT0avErr > 300) lT0avErr = 300;
+
+ lT0av = Int_t(lT0av) + lT0avErr/1000.;
+
+ return lT0av;
+ }
+
+
return t0;
}
+//____________________________________________________________________________
+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);
+ }
+ }
+}
+
+//____________________________________________________________________________
//-- Mod By Silvia Arcelli, Francesco Noferini, Barbara Guerzoni
//////////////////////////////////////////////////////////////////////////////
-#include "Riostream.h"
-
-#include "AliTOFT0v1.h"
#include "AliESDtrack.h"
-#include "AliTOFcalibHisto.h"
#include "AliESDEvent.h"
+#include "AliTOFT0v1.h"
ClassImp(AliTOFT0v1)
//____________________________________________________________________________
AliTOFT0v1::AliTOFT0v1():
- fLowerMomBound(0.4),
- fUpperMomBound(2.0),
+ fLowerMomBound(0.5),
+ fUpperMomBound(1.5),
fTimeResolution(0.80e-10),
fTimeCorr(0.),
- fEvent(0x0),
- fCalib(0x0)
+ fEvent(0x0)
+// fCalib(0x0)
{
//
// default constructor
//____________________________________________________________________________
AliTOFT0v1::AliTOFT0v1(AliESDEvent* event):
- fLowerMomBound(0.4),
- fUpperMomBound(2.0),
+ fLowerMomBound(0.5),
+ fUpperMomBound(1.5),
fTimeResolution(0.80e-10),
fTimeCorr(0.),
- fEvent(event),
- fCalib(0x0)
+ fEvent(event)
+// fCalib(0x0)
{
//
// real constructor
fUpperMomBound(tzero.fUpperMomBound),
fTimeResolution(tzero.fTimeResolution),
fTimeCorr(tzero.fTimeCorr),
- fEvent(tzero.fEvent),
- fCalib(tzero.fCalib)
+ fEvent(tzero.fEvent)
+// fCalib(tzero.fCalib)
{
//
// copy constructor
fTimeResolution=tzero.fTimeResolution;
fTimeCorr=tzero.fTimeCorr;
fEvent=tzero.fEvent;
- fCalib=tzero.fCalib;
+// fCalib=tzero.fCalib;
fT0SigmaT0def[0]=tzero.fT0SigmaT0def[0];
fT0SigmaT0def[1]=tzero.fT0SigmaT0def[1];
fT0SigmaT0def[2]=tzero.fT0SigmaT0def[2];
AliTOFT0v1::~AliTOFT0v1()
{
// dtor
- fCalib=NULL;
+// fCalib=NULL;
fEvent=NULL;
}
{
// Caluclate the Event Time using the ESD TOF time
- Float_t timeresolutioninns=fTimeResolution*(1.e+9) * TMath::Sqrt(2.); // convert in [ns]
+ Float_t timeresolutioninns=fTimeResolution*(1.e+9); // convert in [ns]
const Int_t nmaxtracksinset=10;
// if(strstr(option,"all")){
Double_t momOld=t->GetP();
Double_t mom=momOld-0.0036*momOld;
if ((t->GetStatus()&AliESDtrack::kTIME)==0) continue;
+ if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
Double_t time=t->GetTOFsignal();
time*=1.E-3; // tof given in nanoseconds
}
-// cout << " N. of ESD tracks : " << ntrk << endl;
-// cout << " N. of preselected tracks : " << ngoodtrk << endl;
-// cout << " Minimum tof time in set (in ns) : " << mintime << endl;
+// cout << " N. of ESD tracks : " << ntrk << endl;
+// cout << " N. of preselected tracks : " << ngoodtrk << endl;
+// cout << " Minimum tof time in set (in ns) : " << mintime << endl;
AliESDtrack **gtracks=new AliESDtrack*[ngoodtrk];
}
if(ngoodtrkt0>=2){
// Decide how many tracks in set
- Int_t ntracksinset = min(ngoodtrkt0,nmaxtracksinsetCurrent);
+ Int_t ntracksinset = std::min(ngoodtrkt0,nmaxtracksinsetCurrent);
Int_t nset=1;
if(ngoodtrkt0>nmaxtracksinsetCurrent) {nset= (Int_t)(ngoodtrkt0/ntracksinset)+1;}
Double_t chi2singlecut = chi2cut[ntracksinsetmy-1]/ntracksinsetmy + TMath::Abs(chisquarebest-chi2cut[ntracksinsetmy-1])/ntracksinsetmy;
-// printf("tracks removed with a chi2 > %f (chi2total = %f w.r.t. the limit of %f)\n",chi2singlecut,chisquarebest,chi2cut[ntracksinsetmy-1]);
+// printf("tracks removed with a chi2 > %f (chi2total = %f w.r.t. the limit of %f)\n",chi2singlecut,chisquarebest,chi2cut[ntracksinsetmy-1]);
Bool_t kRedoT0 = kFALSE;
ntracksinsetmyCut = ntracksinsetmy;
if(chisquarebest<999.){
Double_t dblechisquare=(Double_t)chisquarebest;
confLevel=(Float_t)TMath::Prob(dblechisquare,ntracksinsetmyCut-1);
-// cout << " Set Number " << nsets << endl;
-// cout << "Best Assignment, selection " << assparticle[0] <<
-// assparticle[1] << assparticle[2] <<
-// assparticle[3] << assparticle[4] <<
-// assparticle[5] << endl;
-// cout << " Chisquare of the set "<< chisquarebest <<endl;
-// cout << " C.L. of the set "<< confLevel <<endl;
-// cout << " T0 for this set (in ns) " << t0best << endl;
+// cout << " Set Number " << nsets << endl;
+// cout << "Best Assignment, selection " << assparticle[0] <<
+// assparticle[1] << assparticle[2] <<
+// assparticle[3] << assparticle[4] <<
+// assparticle[5] << endl;
+// cout << " Chisquare of the set "<< chisquarebest <<endl;
+// cout << " C.L. of the set "<< confLevel <<endl;
+// cout << " T0 for this set (in ns) " << t0best << endl;
for(Int_t icsq=0; icsq<ntracksinsetmy;icsq++){
if(! usetrack[icsq]) continue;
-// cout << "Track # " << icsq << " T0 offsets = "
-// << besttimezero[icsq]-t0best <<
-// " track error = " << bestsqTrackError[icsq]
-// << " Chisquare = " << bestchisquare[icsq]
-// << " Momentum = " << bestmomentum[icsq]
-// << " TOF = " << besttimeofflight[icsq]
-// << " TOF tracking = " << besttexp[icsq]
-// << " is used = " << usetrack[icsq] << endl;
+// cout << "Track # " << icsq << " T0 offsets = "
+// << besttimezero[icsq]-t0best <<
+// " track error = " << bestsqTrackError[icsq]
+// << " Chisquare = " << bestchisquare[icsq]
+// << " Momentum = " << bestmomentum[icsq]
+// << " TOF = " << besttimeofflight[icsq]
+// << " TOF tracking = " << besttexp[icsq]
+// << " is used = " << usetrack[icsq] << endl;
}
// Pick up only those with C.L. >1%
}
fT0SigmaT0def[0]=t0def;
- fT0SigmaT0def[1]=TMath::Sqrt(deltat0def*deltat0def*(ngoodtrktrulyused/(ngoodtrktrulyused-1)));
+ fT0SigmaT0def[1]=TMath::Sqrt(deltat0def*deltat0def);//*ngoodtrktrulyused/(ngoodtrktrulyused-1));
fT0SigmaT0def[2]=ngoodtrkt0;
fT0SigmaT0def[3]=ngoodtrktrulyused;
}
return fT0SigmaT0def;
}
//__________________________________________________________________
-Double_t * AliTOFT0v1::DefineT0RawCorrection(Option_t *option)
-{
- // Caluclate the Event Time using the RAW+correction TOF time
-
- Float_t timeresolutioninns=fTimeResolution*(1.e+9); // convert in [ns]
-
- const Int_t nmaxtracksinset=10;
-// if(strstr(option,"all")){
-// cout << "Selecting primary tracks with momentum between " << fLowerMomBound << " GeV/c and " << fUpperMomBound << " GeV/c" << endl;
-// cout << "Memorandum: 0 means PION | 1 means KAON | 2 means PROTON" << endl;
-// }
-
- Float_t stripmean = 0;
-
- Int_t nsets=0;
- Int_t nUsedTracks=0;
- Int_t ngoodsetsSel= 0;
- Float_t t0bestSel[300];
- Float_t eT0bestSel[300];
- Float_t chiSquarebestSel[300];
- Float_t confLevelbestSel[300];
- Float_t t0bestallSel=0.;
- Float_t eT0bestallSel=0.;
- Float_t sumWt0bestallSel=0.;
- Float_t eMeanTzeroPi=0.;
- Float_t meantzeropi=0.;
- Float_t sumAllweightspi=0.;
- Double_t t0def=-999;
- Double_t deltat0def=999;
- Int_t ngoodtrktrulyused=0;
- Int_t ntracksinsetmyCut = 0;
-
- Int_t ntrk=fEvent->GetNumberOfTracks();
-
- AliESDtrack **tracks=new AliESDtrack*[ntrk];
- Int_t ngoodtrk=0;
- Int_t ngoodtrkt0 =0;
- Float_t mintime =1E6;
-
- // First Track loop, Selection of good tracks
-
- for (Int_t itrk=0; itrk<ntrk; itrk++) {
- AliESDtrack *t=fEvent->GetTrack(itrk);
- Double_t momOld=t->GetP();
- Double_t mom=momOld-0.0036*momOld;
- if ((t->GetStatus()&AliESDtrack::kTIME)==0) continue;
- Double_t tot = t->GetTOFsignalToT();
- Double_t time=t->GetTOFsignalRaw();
- Int_t chan = t->GetTOFCalChannel();
- Double_t corr = fCalib->GetFullCorrection(chan,tot) - fCalib->GetCorrection(AliTOFcalibHisto::kTimeSlewingCorr,chan,0);
- time -= corr*1000.;
-
- Int_t crate = (Int_t) fCalib->GetCalibMap(AliTOFcalibHisto::kDDL,chan);
- if(crate == 63 || crate == 62){
- time += 9200;
- }
-
- Int_t strip = (Int_t) fCalib->GetCalibMap(AliTOFcalibHisto::kSectorStrip,chan);
-
- time*=1.E-3; // tof given in nanoseconds
- if (!(mom<=fUpperMomBound && mom>=fLowerMomBound))continue;
-
- if (!AcceptTrack(t)) continue;
-
- if(t->GetP() < fLowerMomBound || t->GetIntegratedLength() < 350 || t->GetTOFsignalToT() < 0.000000001)continue; //skip decays
- if(time <= mintime) mintime=time;
- tracks[ngoodtrk]=t;
- ngoodtrk++;
- stripmean += strip;
- }
- if(ngoodtrk) stripmean /= ngoodtrk;
-
-// cout << " N. of ESD tracks : " << ntrk << endl;
-// cout << " N. of preselected tracks : " << ngoodtrk << endl;
-// cout << " Minimum tof time in set (in ns) : " << mintime << endl;
-
- AliESDtrack **gtracks=new AliESDtrack*[ngoodtrk];
-
- for (Int_t jtrk=0; jtrk< ngoodtrk; jtrk++) {
- AliESDtrack *t=tracks[jtrk];
- Double_t tot = t->GetTOFsignalToT();
- Double_t time=t->GetTOFsignalRaw();
- Int_t chan = t->GetTOFCalChannel();
- Double_t corr = fCalib->GetFullCorrection(chan,tot) - fCalib->GetCorrection(AliTOFcalibHisto::kTimeSlewingCorr,chan,0);
- time -= corr*1000.;
-
- Int_t crate = (Int_t) fCalib->GetCalibMap(AliTOFcalibHisto::kDDL,chan);
- if(crate == 63 || crate == 62){
- time += 9200;
- }
-
- if((time-mintime*1.E3)<50.E3){ // For pp and per
- gtracks[ngoodtrkt0]=t;
- ngoodtrkt0++;
- }
- }
-
- Int_t nseteq = (ngoodtrkt0-1)/nmaxtracksinset + 1;
- Int_t nmaxtracksinsetCurrent=ngoodtrkt0/nseteq;
- if(nmaxtracksinsetCurrent*nseteq < ngoodtrkt0) nmaxtracksinsetCurrent++;
-
- if(ngoodtrkt0 > 1){
- Int_t nlastset = (ngoodtrkt0 % nmaxtracksinsetCurrent);
-
- while(nlastset-nseteq+1 > 2 ){
- nmaxtracksinsetCurrent++;
- nlastset -= nseteq-1;
- }
- if(nmaxtracksinsetCurrent > nmaxtracksinset) nmaxtracksinsetCurrent = nmaxtracksinset;
- }
-
- if(ngoodtrkt0<2){
-// cout << "less than 2 tracks, skip event " << endl;
- t0def=-999;
- deltat0def=0.600;
- fT0SigmaT0def[0]=t0def;
- fT0SigmaT0def[1]=deltat0def;
- fT0SigmaT0def[2]=ngoodtrkt0;
- fT0SigmaT0def[3]=ngoodtrkt0;
- //goto finish;
- }
- if(ngoodtrkt0>=2){
- // Decide how many tracks in set
- Int_t ntracksinset = min(ngoodtrkt0,nmaxtracksinsetCurrent);
- Int_t nset=1;
- if(ngoodtrkt0>nmaxtracksinset) {nset= (Int_t)(ngoodtrkt0/ntracksinset)+1;}
-
- // Loop over selected sets
-
- if(nset>=1){
- for (Int_t i=0; i< nset; i++) {
-
- Float_t t0best=999.;
- Float_t eT0best=999.;
- Float_t chisquarebest=99999.;
- Int_t npionbest=0;
-
- Int_t ntracksinsetmy=0;
- AliESDtrack **tracksT0=new AliESDtrack*[ntracksinset];
- for (Int_t itrk=0; itrk<ntracksinset; itrk++) {
- Int_t index = itrk+i*ntracksinset;
- if(index < ngoodtrkt0){
- AliESDtrack *t=gtracks[index];
- tracksT0[itrk]=t;
- ntracksinsetmy++;
- }
- }
-
- // Analyse it
-
- Int_t assparticle[nmaxtracksinset];
- Float_t exptof[nmaxtracksinset][3];
- Float_t timeofflight[nmaxtracksinset];
- Float_t momentum[nmaxtracksinset];
- Float_t timezero[nmaxtracksinset];
- Float_t weightedtimezero[nmaxtracksinset];
- Float_t beta[nmaxtracksinset];
- Float_t texp[nmaxtracksinset];
- Float_t dtexp[nmaxtracksinset];
- Float_t sqMomError[nmaxtracksinset];
- Float_t sqTrackError[nmaxtracksinset];
- Float_t massarray[3]={0.13957,0.493677,0.9382723};
- Float_t tracktoflen[nmaxtracksinset];
- Float_t besttimezero[nmaxtracksinset];
- Float_t besttexp[nmaxtracksinset];
- Float_t besttimeofflight[nmaxtracksinset];
- Float_t bestmomentum[nmaxtracksinset];
- Float_t bestchisquare[nmaxtracksinset];
- Float_t bestweightedtimezero[nmaxtracksinset];
- Float_t bestsqTrackError[nmaxtracksinset];
- Int_t imass[nmaxtracksinset];
-
- for (Int_t j=0; j<ntracksinset; j++) {
- assparticle[j] = 3;
- timeofflight[j] = 0;
- momentum[j] = 0;
- timezero[j] = 0;
- weightedtimezero[j] = 0;
- beta[j] = 0;
- texp[j] = 0;
- dtexp[j] = 0;
- sqMomError[j] = 0;
- sqTrackError[j] = 0;
- tracktoflen[j] = 0;
- besttimezero[j] = 0;
- besttexp[j] = 0;
- besttimeofflight[j] = 0;
- bestmomentum[j] = 0;
- bestchisquare[j] = 0;
- bestweightedtimezero[j] = 0;
- bestsqTrackError[j] = 0;
- imass[j] = 1;
- }
-
- for (Int_t j=0; j<ntracksinsetmy; j++) {
- AliESDtrack *t=tracksT0[j];
- Double_t momOld=t->GetP();
- Double_t mom=momOld-0.0036*momOld;
- Double_t tot = t->GetTOFsignalToT();
- Double_t time=t->GetTOFsignalRaw();
- Int_t chan = t->GetTOFCalChannel();
- Double_t corr = fCalib->GetFullCorrection(chan,tot) - fCalib->GetCorrection(AliTOFcalibHisto::kTimeSlewingCorr,chan,0);
- time -= corr*1000.;
-
- Int_t crate = (Int_t) fCalib->GetCalibMap(AliTOFcalibHisto::kDDL,chan);
- if(crate == 63 || crate == 62){
- time += 9200;
- }
-
- time*=1.E-3; // tof given in nanoseconds
- Double_t exptime[10]; t->GetIntegratedTimes(exptime);
- Double_t toflen=t->GetIntegratedLength();
- toflen=toflen/100.; // toflen given in m
-
- timeofflight[j]=time;
- tracktoflen[j]=toflen;
- exptof[j][0]=exptime[2]*1.E-3+fTimeCorr;// in ns
- exptof[j][1]=exptime[3]*1.E-3+fTimeCorr;
- exptof[j][2]=exptime[4]*1.E-3+fTimeCorr;
- momentum[j]=mom;
- assparticle[j]=3;
-
- } //end for (Int_t j=0; j<ntracksinsetmy; j++) {
-
- for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
- beta[itz]=momentum[itz]/sqrt(massarray[0]*massarray[0]
- +momentum[itz]*momentum[itz]);
- sqMomError[itz]= ((1.-beta[itz]*beta[itz])*0.01)*((1.-beta[itz]*beta[itz])*0.01)*(tracktoflen[itz]/(0.299792*beta[itz]))*(tracktoflen[itz]/(0.299792*beta[itz]));
- sqTrackError[itz]=(timeresolutioninns*timeresolutioninns+sqMomError[itz]); //in ns
- timezero[itz]=exptof[itz][0]-timeofflight[itz];// in ns
- weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
- sumAllweightspi+=1./sqTrackError[itz];
- meantzeropi+=weightedtimezero[itz];
- } // end loop for (Int_t itz=0; itz< ntracksinset;itz++)
-
-
- // Then, Combinatorial Algorithm
-
- if(ntracksinsetmy<2 )break;
-
- for (Int_t j=0; j<ntracksinsetmy; j++) {
- imass[j] = 3;
- }
-
- Int_t ncombinatorial = Int_t(TMath::Power(3,ntracksinsetmy));
-
- // Loop on mass hypotheses
- for (Int_t k=0; k < ncombinatorial;k++) {
- for (Int_t j=0; j<ntracksinsetmy; j++) {
- imass[j] = (k % Int_t(TMath::Power(3,ntracksinsetmy-j)))/Int_t(TMath::Power(3,ntracksinsetmy-j-1));
- texp[j]=exptof[j][imass[j]];
- dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
- }
- Float_t sumAllweights=0.;
- Float_t meantzero=0.;
- Float_t eMeanTzero=0.;
- Double_t sumAllSquare=0.;
-
- for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
- sqTrackError[itz]=
- (timeresolutioninns*
- timeresolutioninns
- +dtexp[itz]*dtexp[itz]*1E-6); //in ns2
-
-// printf("pt=%f -- TOF res=%f -- track res=%f -- all res=%f\n",momentum[itz],timeresolutioninns,dtexp[itz],TMath::Sqrt(sqTrackError[itz]));
-// getchar();
- timezero[itz]=texp[itz]-timeofflight[itz];// in ns
-
- weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
- sumAllSquare += timezero[itz]*timezero[itz];
- sumAllweights+=1./sqTrackError[itz];
- meantzero+=weightedtimezero[itz];
-
- } // end loop for (Int_t itz=0; itz<15;itz++)
-
- meantzero=meantzero/sumAllweights; // it is given in [ns]
- eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns]
-
- //changed
- for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
- sumAllSquare+= (timezero[itz] - meantzero)*(timezero[itz] - meantzero);
- }
- // eMeanTzero = TMath::Sqrt(sumAllSquare/ntracksinsetmy);
-
- // calculate chisquare
-
- Float_t chisquare=0.;
- for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
- chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
-
- } // end loop for (Int_t icsq=0; icsq<15;icsq++)
-
- if(chisquare<=chisquarebest){
- for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) {
-
- bestsqTrackError[iqsq]=sqTrackError[iqsq];
- besttimezero[iqsq]=timezero[iqsq];
- bestmomentum[iqsq]=momentum[iqsq];
- besttimeofflight[iqsq]=timeofflight[iqsq];
- besttexp[iqsq]=texp[iqsq];
- bestweightedtimezero[iqsq]=weightedtimezero[iqsq];
- bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq];
- }
-
- Int_t npion=0;
- for (Int_t j=0; j<ntracksinsetmy; j++) {
- assparticle[j]=imass[j];
- if(imass[j] == 0) npion++;
- }
- npionbest=npion;
- chisquarebest=chisquare;
- t0best=meantzero;
- eT0best=eMeanTzero;
- } // close if(dummychisquare<=chisquare)
-
- }
-
- Double_t chi2cut[nmaxtracksinset];
- chi2cut[0] = 0;
- chi2cut[1] = 6.6; // corresponding to a C.L. of 0.01
- for (Int_t j=2; j<ntracksinset; j++) {
- chi2cut[j] = chi2cut[1] * TMath::Sqrt(j*1.);
- }
-
- Double_t chi2singlecut = chi2cut[ntracksinsetmy-1]/ntracksinsetmy + TMath::Abs(chisquarebest-chi2cut[ntracksinsetmy-1])/ntracksinsetmy;
-
-// printf("tracks removed with a chi2 > %f (chi2total = %f w.r.t. the limit of %f)\n",chi2singlecut,chisquarebest,chi2cut[ntracksinsetmy-1]);
-
- Bool_t kRedoT0 = kFALSE;
- ntracksinsetmyCut = ntracksinsetmy;
- Bool_t usetrack[nmaxtracksinset];
- for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
- usetrack[icsq] = kTRUE;
- if((bestchisquare[icsq] > chisquarebest*0.5 && ntracksinsetmy > 2) || (bestchisquare[icsq] > chi2singlecut)){
- kRedoT0 = kTRUE;
- ntracksinsetmyCut--;
- usetrack[icsq] = kFALSE;
- }
- } // end loop for (Int_t icsq=0; icsq<15;icsq++)
-
- // printf("ntrackinsetmy = %i - %i\n",ntracksinsetmy,ntracksinsetmyCut);
-
- // Loop on mass hypotheses Redo
- if(kRedoT0 && ntracksinsetmyCut > 1){
- // printf("Redo T0\n");
- for (Int_t k=0; k < ncombinatorial;k++) {
- for (Int_t j=0; j<ntracksinsetmy; j++) {
- imass[j] = (k % Int_t(TMath::Power(3,ntracksinsetmy-j))) / Int_t(TMath::Power(3,ntracksinsetmy-j-1));
- texp[j]=exptof[j][imass[j]];
- dtexp[j]=GetMomError(imass[j], momentum[j], texp[j]);
- }
-
- Float_t sumAllweights=0.;
- Float_t meantzero=0.;
- Float_t eMeanTzero=0.;
- Double_t sumAllSquare=0;
-
- for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
- if(! usetrack[itz]) continue;
- sqTrackError[itz]=
- (timeresolutioninns*
- timeresolutioninns
- +dtexp[itz]*dtexp[itz]*1E-6); //in ns2
-
- timezero[itz]=texp[itz]-timeofflight[itz];// in ns
-
- weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
- sumAllweights+=1./sqTrackError[itz];
- meantzero+=weightedtimezero[itz];
- } // end loop for (Int_t itz=0; itz<15;itz++)
-
- meantzero=meantzero/sumAllweights; // it is given in [ns]
- eMeanTzero=sqrt(1./sumAllweights); // it is given in [ns]
-
- //changed
- for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
- if(! usetrack[itz]) continue;
- sumAllSquare+= (timezero[itz] - meantzero)*(timezero[itz] - meantzero);
- }
- // eMeanTzero = TMath::Sqrt(sumAllSquare/ntracksinsetmyCut);
-
- // calculate chisquare
-
- Float_t chisquare=0.;
- for (Int_t icsq=0; icsq<ntracksinsetmy;icsq++) {
- if(! usetrack[icsq]) continue;
- chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
-
- } // end loop for (Int_t icsq=0; icsq<15;icsq++)
-
- Int_t npion=0;
- for (Int_t j=0; j<ntracksinsetmy; j++) {
- assparticle[j]=imass[j];
- if(imass[j] == 0) npion++;
- }
-
- if(chisquare<=chisquarebest){
- for(Int_t iqsq = 0; iqsq<ntracksinsetmy; iqsq++) {
- if(! usetrack[iqsq]) continue;
- bestsqTrackError[iqsq]=sqTrackError[iqsq];
- besttimezero[iqsq]=timezero[iqsq];
- bestmomentum[iqsq]=momentum[iqsq];
- besttimeofflight[iqsq]=timeofflight[iqsq];
- besttexp[iqsq]=texp[iqsq];
- bestweightedtimezero[iqsq]=weightedtimezero[iqsq];
- bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq];
- }
-
- npionbest=npion;
- chisquarebest=chisquare;
- t0best=meantzero;
- eT0best=eMeanTzero;
- } // close if(dummychisquare<=chisquare)
-
- }
- }
-
- if(chisquarebest >= 999){
- printf("How is it possible (chi2 = %f)? T0best = %f\n",chisquarebest,t0best);
-
-// for(Int_t icsq=0; icsq<ntracksinsetmy;icsq++){
-// cout << "Track # " << icsq << " T0 offsets = "
-// << besttimezero[icsq]-t0best <<
-// " track error = " << bestsqTrackError[icsq]
-// << " Chisquare = " << bestchisquare[icsq]
-// << " Momentum = " << bestmomentum[icsq]
-// << " TOF = " << besttimeofflight[icsq]
-// << " TOF tracking = " << besttexp[icsq]
-// << " is used = " << usetrack[icsq] << endl;
-// }
- }
-
- // filling histos
- Float_t confLevel=999;
-
- // Sets with decent chisquares
-
- if(chisquarebest<999.){
- Double_t dblechisquare=(Double_t)chisquarebest;
- confLevel=(Float_t)TMath::Prob(dblechisquare,ntracksinsetmyCut-1);
-// cout << " Set Number " << nsets << endl;
-// cout << "Best Assignment, selection " << assparticle[0] <<
-// assparticle[1] << assparticle[2] <<
-// assparticle[3] << assparticle[4] <<
-// assparticle[5] << endl;
-// cout << " Chisquare of the set "<< chisquarebest <<endl;
-// cout << " C.L. of the set "<< confLevel <<endl;
-// cout << " T0 for this set (in ns) " << t0best << endl;
-
- for(Int_t icsq=0; icsq<ntracksinsetmy;icsq++){
-
- if(! usetrack[icsq]) continue;
-
-// cout << "Track # " << icsq << " T0 offsets = "
-// << besttimezero[icsq]-t0best <<
-// " track error = " << bestsqTrackError[icsq]
-// << " Chisquare = " << bestchisquare[icsq]
-// << " Momentum = " << bestmomentum[icsq]
-// << " TOF = " << besttimeofflight[icsq]
-// << " TOF tracking = " << besttexp[icsq]
-// << " is used = " << usetrack[icsq] << endl;
- }
-
- // Pick up only those with C.L. >1%
- // if(confLevel>0.01 && ngoodsetsSel<200){
- if(confLevel>0.01 && ngoodsetsSel<200){
- chiSquarebestSel[ngoodsetsSel]=chisquarebest;
- confLevelbestSel[ngoodsetsSel]=confLevel;
- t0bestSel[ngoodsetsSel]=t0best/eT0best/eT0best;
- eT0bestSel[ngoodsetsSel]=1./eT0best/eT0best;
- t0bestallSel += t0best/eT0best/eT0best;
- sumWt0bestallSel += 1./eT0best/eT0best;
-
- ngoodsetsSel++;
- ngoodtrktrulyused+=ntracksinsetmyCut;
- }
- else{
- // printf("conflevel = %f -- ngoodsetsSel = %i -- ntrackset = %i\n",confLevel,ngoodsetsSel,ntracksinsetmy);
- }
- }
- delete[] tracksT0;
- nsets++;
-
- } // end for the current set
-
- nUsedTracks = ngoodtrkt0;
- if(strstr(option,"all")){
- if(sumAllweightspi>0.){
- meantzeropi=meantzeropi/sumAllweightspi; // it is given in [ns]
- eMeanTzeroPi=sqrt(1./sumAllweightspi); // it is given in [ns]
- }
-
- if(sumWt0bestallSel>0){
- t0bestallSel = t0bestallSel/sumWt0bestallSel;
- eT0bestallSel = sqrt(1./sumWt0bestallSel);
- }// end of if(sumWt0bestallSel>0){
-
- }
-
- t0def=t0bestallSel;
- deltat0def=eT0bestallSel;
- if ((TMath::Abs(t0bestallSel)<0.001)&&(TMath::Abs(eT0bestallSel)<0.001)){
- t0def=-999; deltat0def=0.600;
- }
-
- fT0SigmaT0def[0]=t0def;
- fT0SigmaT0def[1]=TMath::Sqrt(deltat0def*deltat0def*(ngoodtrktrulyused/(ngoodtrktrulyused-1)));
- fT0SigmaT0def[2]=ngoodtrkt0;
- fT0SigmaT0def[3]=ngoodtrktrulyused;
- }
- }
-
- return fT0SigmaT0def;
- }
-
-//__________________________________________________________________
Float_t AliTOFT0v1::GetMomError(Int_t index, Float_t mom, Float_t texp) const
{
// Take the error extimate for the TOF time in the track reconstruction
Double_t mass=kMasses[index+2];
Double_t dpp=0.01; //mean relative pt resolution;
+ if(mom > 1) dpp = 0.01*mom;
Double_t sigma=dpp*texp*1E3/(1.+ mom*mom/(mass*mass));
sigma =TMath::Sqrt(sigma*sigma);
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff