- 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){