--- /dev/null
+/**************************************************************************
+ * 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: AliTOFT0v1.cxx,v 1.8 2003/10/23 16:32:20 hristov Exp $ */
+
+//_________________________________________________________________________
+// This is a TTask that made the calculation of the Time zero using TOF.
+// Description: The algorithm used to calculate the time zero of interaction
+// using TOF detector is the following.
+// We select in the MonteCarlo some primary particles - or tracks in the following -
+// that strike the TOF detector (the larger part are pions, kaons or protons).
+// We choose a set of 10 selected tracks, for each track You have the length
+// of the track when the TOF is reached (a standard TOF hit does not contain this
+// additional information, this is the reason why we implemented a new time zero
+// dedicated TOF hit class AliTOFhitT0; in order to store this type of hit You
+// have to use the AliTOFv4T0 as TOF class in Your Config.C. In AliTOFv4T0 the
+// StepManager was modified in order to fill the TOF hit branch with this type
+// of hits; in fact the AliTOF::AddT0Hit is called rather that the usual AliTOF::AddHit),
+// the momentum at generation (from TreeK) and the time of flight
+// given by the TOF detector.
+// (Observe that the ctor of the AliTOF class, when the AliTOFv4T0 class is used, is called
+// with the "tzero" option: it is in order create the fHits TClonesArray filled with
+// AliTOFhitT0 objects, rather than with normal AliTOFhit)
+// Then Momentum and time of flight for each track are smeared according to
+// known experimental resolution (all sources of error have been token into account).
+// Let consider now only one set of 10 tracks (the algorithm is the same for all sets).
+// Assuming the (mass) hypothesis that each track can be AUT a pion, AUT a kaon, AUT a proton,
+// we consider all the 3 at 10 possible cases.
+// For each track in each (mass) configuration
+// (a configuration can be e.g. pion/pion/kaon/proton/pion/proton/kaon/kaon/pion/pion)
+// we calculate the time zero (we know in fact the velocity of the track after
+// the assumption about its mass, the time of flight given by the TOF, and the
+// corresponding path travelled till the TOF detector). Then for each mass configuration we have
+// 10 time zero and we can calculate the ChiSquare for the current configuration using the
+// weighted mean over all 10 time zero.
+// We call the best assignment the mass configuration that gives the minimum value of the ChiSquare.
+// We plot the weighted mean over all 10 time zero for the best assignment,
+// the ChiSquare for the best assignment and the corresponding confidence level.
+// The strong assumption is the MC selection of primary particles. It will be introduced
+// in the future also some more realistic simulation about this point.
+// Use case:
+// root [0] AliTOFT0v1 * tzero = new AliTOFT0v1("galice.root")
+// Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
+// root [1] tzero->ExecuteTask()
+// root [2] tzero->ExecuteTask("tim")
+// // available parameters:
+// tim - print benchmarking information
+// all - print usefull informations about the number of misidentified tracks
+// and a comparison about the true configuration (known from MC) and the best
+// assignment
+// Different Selections for pp and Pb-Pb: Momentum Range, Max Time, # pions
+//-- Author: F. Pierella
+//-- Mod By SA.
+//////////////////////////////////////////////////////////////////////////////
+
+#include <Riostream.h>
+#include <stdlib.h>
+
+#include <TBenchmark.h>
+#include <TCanvas.h>
+#include <TClonesArray.h>
+#include <TFile.h>
+#include <TFolder.h>
+#include <TFrame.h>
+#include <TH1.h>
+#include <TH2.h>
+#include <TProfile.h>
+#include <TParticle.h>
+#include <TROOT.h>
+#include <TSystem.h>
+#include <TTree.h>
+#include <TVirtualMC.h>
+#include "AliTOFT0v1.h"
+#include "AliMC.h"
+#include "AliESDtrack.h"
+#include "AliESD.h"
+#include <AliStack.h>
+#include "AliESDtrackCuts.h"
+#include "AliTOFcalibHisto.h"
+
+ClassImp(AliTOFT0v1)
+
+//____________________________________________________________________________
+AliTOFT0v1::AliTOFT0v1(AliESDEvent* event)
+{
+ fLowerMomBound=0.4; // [GeV/c] default value pp
+ fUpperMomBound=2.0 ; // [GeV/c] default value pp
+ fTimeResolution = 0.80e-10; // 80 ps by default
+ fTimeCorr = 0.0; // in ns by default
+ fLOffset=0.0;
+ fT0Offset=0.0;
+ fEvent=event;
+
+ fT0SigmaT0def[0]=-999.;
+ fT0SigmaT0def[1]=999.;
+ fT0SigmaT0def[2]=-999.;
+ fT0SigmaT0def[3]=-999.;
+
+ fDeltaTfromMisallinement = 0.; // in ps
+}
+
+//____________________________________________________________________________
+AliTOFT0v1::AliTOFT0v1(const AliTOFT0v1 & tzero)
+{
+ ( (AliTOFT0v1 &)tzero ).Copy(*this);
+}
+
+//____________________________________________________________________________
+AliTOFT0v1::~AliTOFT0v1()
+{
+ // dtor
+
+}
+void AliTOFT0v1::SetTimeResolution(Double_t timeresolution){
+ fTimeResolution=timeresolution;
+}
+//____________________________________________________________________________
+//____________________________________________________________________________
+Double_t * AliTOFT0v1::DefineT0(Option_t *option)
+{
+ Float_t timeresolutioninns=fTimeResolution*(1.e+9) * TMath::Sqrt(2.); // 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;
+ }
+
+
+ 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 time=t->GetTOFsignal();
+
+ time*=1.E-3; // tof given in nanoseconds
+ if (!(mom<=fUpperMomBound && mom>=fLowerMomBound))continue;
+
+ AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts();
+ Bool_t tpcRefit = kTRUE;
+ Double_t nSigma = 4;
+ Bool_t sigmaToVertex = kTRUE;
+ esdTrackCuts->SetRequireSigmaToVertex(sigmaToVertex);
+ if (sigmaToVertex) {
+ esdTrackCuts->SetMaxNsigmaToVertex(nSigma);
+ }
+ else{
+ esdTrackCuts->SetMaxDCAToVertexZ(3.0);
+ esdTrackCuts->SetMaxDCAToVertexXY(3.0);
+ }
+ esdTrackCuts->SetRequireTPCRefit(tpcRefit);
+ esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
+ esdTrackCuts->SetMinNClustersTPC(50);
+ esdTrackCuts->SetMaxChi2PerClusterTPC(3.5);
+ Bool_t accepted;
+ accepted=esdTrackCuts->AcceptTrack(t);
+ if(!accepted) continue;
+
+ if(t->GetP() < fLowerMomBound || t->GetIntegratedLength() < 350 || t->GetTOFsignalToT() < 0.000000001)continue; //skip decays
+ if(time <= mintime) mintime=time;
+ tracks[ngoodtrk]=t;
+ ngoodtrk++;
+ }
+
+
+ cout << " T0 offset selected for this sample : " << fT0Offset << 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];
+
+ for (Int_t jtrk=0; jtrk< ngoodtrk; jtrk++) {
+ AliESDtrack *t=tracks[jtrk];
+ Double_t time=t->GetTOFsignal();
+
+ 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<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>nmaxtracksinsetCurrent) {nset= (Int_t)(ngoodtrkt0/ntracksinset)+1;}
+
+ cout << " number of sets = " << nset << endl;
+
+ if(strstr(option,"tim") || strstr(option,"all"))gBenchmark->Start("TOFT0v1");
+
+ // 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 time=t->GetTOFsignal();
+
+ 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+fT0Offset;
+ tracktoflen[j]=toflen+fLOffset;
+ 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++) {
+
+ cout << "starting t0 calculation for current set" << endl;
+ 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.;
+
+ for (Int_t itz=0; itz<ntracksinsetmy;itz++) {
+ 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]
+
+ // 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.;
+
+ 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]
+
+ // 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){
+
+ cout << "T0 all " << t0bestallSel << " +/- " << Et0bestallSel << "Number of tracks used: "<<ngoodtrktrulyused<<endl;
+ }
+
+ t0def=t0bestallSel;
+ deltat0def=Et0bestallSel;
+ if ((t0bestallSel==0)&&(Et0bestallSel==0)){
+ t0def=-999; deltat0def=0.600;
+ }
+
+ fT0SigmaT0def[0]=t0def;
+ fT0SigmaT0def[1]=TMath::Sqrt(deltat0def*deltat0def*(ngoodtrktrulyused/(ngoodtrktrulyused-1)));
+ fT0SigmaT0def[2]=ngoodtrkt0;//ngoodtrktrulyused;
+ fT0SigmaT0def[3]=ngoodtrktrulyused;
+ }
+ }
+
+ if(strstr(option,"tim") || strstr(option,"all")){
+ gBenchmark->Stop("TOFT0v1");
+ cout << "AliTOFT0v1:" << endl ;
+ cout << " took " << gBenchmark->GetCpuTime("TOFT0v1") << " seconds in order to calculate T0 " << endl ;
+ }
+ printf("T0 from TOF = %f ns\n",fT0SigmaT0def[0]);
+
+ return fT0SigmaT0def;
+ }
+//__________________________________________________________________
+Double_t * AliTOFT0v1::DefineT0RawCorrection(Option_t *option)
+{
+ 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 = fCalib->GetCalibMap(AliTOFcalibHisto::kDDL,chan);
+ if(crate == 63 || crate == 62){
+ time += 9200;
+ }
+
+ Int_t strip = fCalib->GetCalibMap(AliTOFcalibHisto::kSectorStrip,chan);
+
+ time*=1.E-3; // tof given in nanoseconds
+ if (!(mom<=fUpperMomBound && mom>=fLowerMomBound))continue;
+
+ AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts();
+ Bool_t tpcRefit = kTRUE;
+ Double_t nSigma = 4;
+ Bool_t sigmaToVertex = kTRUE;
+ esdTrackCuts->SetRequireSigmaToVertex(sigmaToVertex);
+ if (sigmaToVertex) {
+ esdTrackCuts->SetMaxNsigmaToVertex(nSigma);
+ }
+ else{
+ esdTrackCuts->SetMaxDCAToVertexZ(3.0);
+ esdTrackCuts->SetMaxDCAToVertexXY(3.0);
+ }
+ esdTrackCuts->SetRequireTPCRefit(tpcRefit);
+ esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
+ esdTrackCuts->SetMinNClustersTPC(50);
+ esdTrackCuts->SetMaxChi2PerClusterTPC(3.5);
+ Bool_t accepted;
+ accepted=esdTrackCuts->AcceptTrack(t);
+ if(!accepted) 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 << " T0 offset selected for this sample : " << fT0Offset << 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];
+
+ 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 create = fCalib->GetCalibMap(AliTOFcalibHisto::kDDL,chan);
+ if(create == 63 || create == 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;}
+
+ cout << " number of sets = " << nset << endl;
+
+ if(strstr(option,"tim") || strstr(option,"all"))gBenchmark->Start("TOFT0v1");
+
+ // 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 create = fCalib->GetCalibMap(AliTOFcalibHisto::kDDL,chan);
+ if(create == 63 || create == 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+fT0Offset;
+ tracktoflen[j]=toflen+fLOffset;
+ 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++) {
+
+ cout << "starting t0 calculation for current set" << endl;
+ 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){
+
+ cout << "T0 all " << t0bestallSel << " +/- " << Et0bestallSel << "Number of tracks used: "<<ngoodtrktrulyused<<endl;
+ }
+
+ t0def=t0bestallSel;
+ deltat0def=Et0bestallSel;
+ if ((t0bestallSel==0)&&(Et0bestallSel==0)){
+ t0def=-999; deltat0def=0.600;
+ }
+
+ fT0SigmaT0def[0]=t0def;
+ fT0SigmaT0def[1]=TMath::Sqrt(deltat0def*deltat0def*(ngoodtrktrulyused/(ngoodtrktrulyused-1)));
+ fT0SigmaT0def[2]=ngoodtrkt0;//ngoodtrktrulyused;
+ fT0SigmaT0def[3]=ngoodtrktrulyused;
+ }
+ }
+
+ if(strstr(option,"tim") || strstr(option,"all")){
+ gBenchmark->Stop("TOFT0v1");
+ cout << "AliTOFT0v1:" << endl ;
+ cout << " took " << gBenchmark->GetCpuTime("TOFT0v1") << " seconds in order to calculate T0 " << endl ;
+ }
+ printf("T0 from TOF = %f ns\n",fT0SigmaT0def[0]);
+
+ return fT0SigmaT0def;
+ }
+//__________________________________________________________________
+void AliTOFT0v1::SetTZeroFile(char * file ){
+ cout << "Destination file : " << file << endl ;
+ fT0File=file;
+}
+//__________________________________________________________________
+void AliTOFT0v1::Print(Option_t* /*option*/)const
+{
+ cout << "------------------- "<< GetName() << " -------------" << endl ;
+ if(!fT0File.IsNull())
+ cout << " Writing T0 Distribution to file " << (char*) fT0File.Data() << endl ;
+}
+
+//__________________________________________________________________
+Bool_t AliTOFT0v1::operator==( AliTOFT0v1 const &tzero )const
+{
+ // Equal operator.
+ //
+ if( (fTimeResolution==tzero.fTimeResolution)&&(fLowerMomBound==tzero.fLowerMomBound)&&(fUpperMomBound==tzero.fUpperMomBound))
+ return kTRUE ;
+ else
+ return kFALSE ;
+}
+
+
+//__________________________________________________________________
+Float_t AliTOFT0v1::GetMomError(Int_t index, Float_t mom, Float_t texp)
+{
+ static const Double_t kMasses[]={
+ 0.000511, 0.105658, 0.139570, 0.493677, 0.938272, 1.875613
+ };
+
+ Double_t mass=kMasses[index+2];
+ Double_t dpp=0.01; //mean relative pt resolution;
+ Double_t sigma=dpp*texp*1E3/(1.+ mom*mom/(mass*mass));
+
+ sigma =TMath::Sqrt(sigma*sigma);
+
+ return sigma;
+}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff