//-----------------------------------------------------------------
+//
// Implementation of the TPC seed class
// This class is used by the AliTPCtrackerMI class
// Origin: Marian Ivanov, CERN, Marian.Ivanov@cern.ch
//-----------------------------------------------------------------
#include "TClonesArray.h"
+#include "TGraphErrors.h"
#include "AliTPCseed.h"
#include "AliTPCReconstructor.h"
#include "AliTPCClusterParam.h"
return mean;
}
-Float_t AliTPCseed::CookdEdxAnalytical(Double_t low, Double_t up, Int_t type, Int_t i1, Int_t i2, Int_t returnVal){
+Float_t AliTPCseed::CookdEdxAnalytical(Double_t low, Double_t up, Int_t type, Int_t i1, Int_t i2, Int_t returnVal, Int_t rowThres, Int_t mode){
//
// calculates dedx using the cluster
// i1-i2 - the pad-row range used for calculation
//
// posNorm - usage of pos normalization
- // returnVal - 0 return mean
- // - 1 return RMS
- // - 2 return number of clusters
+ // returnVal - 0 return mean
+ // - 1 return RMS
+ // - 2 return number of clusters
+ // - 3 ratio
+ // - 4 mean upper half
+ // - 5 mean - lower half
+ // - 6 third moment
+ // mode - 0 - linear
+ // - 1 - logatithmic
+ // rowThres - number of rows before and after given pad row to check for clusters below threshold
//
// normalization parametrization taken from AliTPCClusterParam
//
Float_t amp[160];
Int_t indexes[160];
Int_t ncl=0;
+ Int_t nclBelowThr = 0; // counts number of clusters below threshold
//
//
Float_t gainGG = 1; // gas gain factor -always enabled
Float_t gainPad = 1; // gain map - used always
Float_t corrPos = 1; // local position correction - if posNorm enabled
-
//
//
//
if (AliTPCcalibDB::Instance()->GetParameters()){
gainGG= AliTPCcalibDB::Instance()->GetParameters()->GetGasGain()/20000; //relative gas gain
}
-
+ //
+ // extract time-dependent correction for pressure and temperature variations
+ //
+ UInt_t runNumber = 1;
+ Float_t corrTimeGain = 1;
+ TObjArray * timeGainSplines = 0x0;
+ TGraphErrors * grPadEqual = 0x0;
+ //
+ AliTPCTransform * trans = AliTPCcalibDB::Instance()->GetTransform();
+ const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
+ //
+ if (recoParam->GetNeighborRowsDedx() == 0) rowThres = 0;
+ //
+ if (trans) {
+ runNumber = trans->GetCurrentRunNumber();
+ //AliTPCcalibDB::Instance()->SetRun(runNumber);
+ timeGainSplines = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(runNumber);
+ if (timeGainSplines && recoParam->GetUseGainCorrectionTime()>0) {
+ UInt_t time = trans->GetCurrentTimeStamp();
+ AliSplineFit * fitMIP = (AliSplineFit *) timeGainSplines->At(0);
+ AliSplineFit * fitFPcosmic = (AliSplineFit *) timeGainSplines->At(1);
+ if (fitMIP) {
+ corrTimeGain = AliTPCcalibDButil::EvalGraphConst(fitMIP, time); /*fitMIP->Eval(time);*/
+ } else {
+ if (fitFPcosmic) corrTimeGain = AliTPCcalibDButil::EvalGraphConst(fitFPcosmic, time); /*fitFPcosmic->Eval(time); */
+ }
+ //
+ if (type==1) grPadEqual = (TGraphErrors * ) timeGainSplines->FindObject("TGRAPHERRORS_MEANQMAX_PADREGIONGAIN_BEAM_ALL");
+ if (type==0) grPadEqual = (TGraphErrors * ) timeGainSplines->FindObject("TGRAPHERRORS_MEANQTOT_PADREGIONGAIN_BEAM_ALL");
+ }
+ }
+
+ const Float_t kClusterShapeCut = 1.5; // IMPPRTANT TO DO: move value to AliTPCRecoParam
const Float_t ktany = TMath::Tan(TMath::DegToRad()*10);
const Float_t kedgey =3.;
//
//
for (Int_t irow=i1; irow<i2; irow++){
AliTPCclusterMI* cluster = GetClusterPointer(irow);
+ if (!cluster && irow > 1 && irow < 157) {
+ Bool_t isClBefore = kFALSE;
+ Bool_t isClAfter = kFALSE;
+ for(Int_t ithres = 1; ithres <= rowThres; ithres++) {
+ AliTPCclusterMI * clusterBefore = GetClusterPointer(irow - ithres);
+ if (clusterBefore) isClBefore = kTRUE;
+ AliTPCclusterMI * clusterAfter = GetClusterPointer(irow + ithres);
+ if (clusterAfter) isClAfter = kTRUE;
+ }
+ if (isClBefore && isClAfter) nclBelowThr++;
+ }
if (!cluster) continue;
+ //
+ //
if (TMath::Abs(cluster->GetY())>cluster->GetX()*ktany-kedgey) continue; // edge cluster
+ //
+ AliTPCTrackerPoint * point = GetTrackPoint(irow);
+ if (point==0) continue;
+ Float_t rsigmay = TMath::Sqrt(point->GetSigmaY());
+ if (rsigmay > kClusterShapeCut) continue;
+ //
+ if (cluster->IsUsed(11)) continue; // remove shared clusters for PbPb
+ //
Float_t charge= (type%2)? cluster->GetMax():cluster->GetQ();
Int_t ipad= 0;
if (irow>=row0) ipad=1;
} else { // OROC
factor = roc->GetValue(irow - row0, TMath::Nint(cluster->GetPad()));
}
- if (factor>0.5) gainPad=factor;
+ if (factor>0.3) gainPad=factor;
}
-
//
// Do position normalization - relative distance to
// center of pad- time bin
- AliTPCTrackerPoint * point = GetTrackPoint(irow);
Float_t ty = TMath::Abs(point->GetAngleY());
Float_t tz = TMath::Abs(point->GetAngleZ()*TMath::Sqrt(1+ty*ty));
Float_t yres0 = parcl->GetRMS0(0,ipad,0,0)/param->GetPadPitchWidth(cluster->GetDetector());
corrPos*=(1+parcl->GetQnormCorr(ipad, type+2,2)*tz*tz);
//
}
-
+ //
+ // pad region equalization outside of cluster param
+ //
+ Float_t gainEqualPadRegion = 1;
+ if (grPadEqual) gainEqualPadRegion = grPadEqual->Eval(ipad);
//
amp[ncl]=charge;
amp[ncl]/=gainGG;
amp[ncl]/=gainPad;
amp[ncl]/=corrPos;
+ amp[ncl]/=gainEqualPadRegion;
//
ncl++;
}
- if (type>3) return ncl;
+ if (type==2) return ncl;
TMath::Sort(ncl,amp, indexes, kFALSE);
-
+ //
if (ncl<10) return 0;
-
+ //
+ Double_t * ampWithBelow = new Double_t[ncl + nclBelowThr];
+ for(Int_t iCl = 0; iCl < ncl + nclBelowThr; iCl++) {
+ if (iCl < nclBelowThr) {
+ ampWithBelow[iCl] = amp[indexes[0]];
+ } else {
+ ampWithBelow[iCl] = amp[indexes[iCl - nclBelowThr]];
+ }
+ }
+ //printf("DEBUG: %i shit %f", nclBelowThr, amp[indexes[0]]);
+ //
Float_t suma=0;
Float_t suma2=0;
+ Float_t suma3=0;
+ Float_t sumaS=0;
Float_t sumn=0;
- Int_t icl0=TMath::Nint(ncl*low);
- Int_t icl1=TMath::Nint(ncl*up);
- for (Int_t icl=icl0; icl<icl1;icl++){
- suma+=amp[indexes[icl]];
- suma2+=amp[indexes[icl]]*amp[indexes[icl]];
- sumn++;
- }
- Float_t mean =suma/sumn;
- Float_t rms =TMath::Sqrt(TMath::Abs(suma2/sumn-mean*mean));
+ // upper,and lower part statistic
+ Float_t sumL=0, sumL2=0, sumLN=0;
+ Float_t sumD=0, sumD2=0, sumDN=0;
+
+ Int_t icl0=TMath::Nint((ncl + nclBelowThr)*low);
+ Int_t icl1=TMath::Nint((ncl + nclBelowThr)*up);
+ Int_t iclm=TMath::Nint((ncl + nclBelowThr)*(low +(up+low)*0.5));
//
- // do time-dependent correction for pressure and temperature variations
- UInt_t runNumber = 1;
- Float_t corrTimeGain = 1;
- AliTPCTransform * trans = AliTPCcalibDB::Instance()->GetTransform();
- const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
- if (trans && recoParam->GetUseGainCorrectionTime()>0) {
- runNumber = trans->GetCurrentRunNumber();
- //AliTPCcalibDB::Instance()->SetRun(runNumber);
- TObjArray * timeGainSplines = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(runNumber);
- if (timeGainSplines) {
- UInt_t time = trans->GetCurrentTimeStamp();
- AliSplineFit * fitMIP = (AliSplineFit *) timeGainSplines->At(0);
- AliSplineFit * fitFPcosmic = (AliSplineFit *) timeGainSplines->At(1);
- if (fitMIP) {
- corrTimeGain = fitMIP->Eval(time);
- } else {
- if (fitFPcosmic) corrTimeGain = fitFPcosmic->Eval(time); // This value describes the ratio FP-to-MIP, hardwired for the moment
+ for (Int_t icl=icl0; icl<icl1;icl++){
+ if (ampWithBelow[icl]<0.1) continue;
+ Double_t camp=ampWithBelow[icl]/corrTimeGain;
+ if (mode==1) camp= TMath::Log(camp);
+ if (icl<icl1){
+ suma+=camp;
+ suma2+=camp*camp;
+ suma3+=camp*camp*camp;
+ sumaS+=TMath::Power(TMath::Abs(camp),1./3.);
+ sumn++;
+ }
+ if (icl>iclm){
+ sumL+=camp;
+ sumL2+=camp*camp;
+ sumLN++;
}
+ if (icl<=iclm){
+ sumD+=camp;
+ sumD2+=camp*camp;
+ sumDN++;
}
}
- mean /= corrTimeGain;
- rms /= corrTimeGain;
+ //
+ Float_t mean = 0;
+ Float_t meanL = 0;
+ Float_t meanD = 0; // lower half mean
+ if (sumn > 1e-30) mean =suma/sumn;
+ if (sumLN > 1e-30) meanL =sumL/sumLN;
+ if (sumDN > 1e-30) meanD =(sumD/sumDN);
+ /*
+ Float_t mean =suma/sumn;
+ Float_t meanL = sumL/sumLN;
+ Float_t meanD =(sumD/sumDN); // lower half mean
+ */
+
+ Float_t rms = 0;
+ Float_t mean2=0;
+ Float_t mean3=0;
+ Float_t meanS=0;
+
+ if(sumn>0){
+ rms = TMath::Sqrt(TMath::Abs(suma2/sumn-mean*mean));
+ mean2=suma2/sumn;
+ mean3=suma3/sumn;
+ meanS=sumaS/sumn;
+ }
+
+ if (mean2>0) mean2=TMath::Power(TMath::Abs(mean2),1./2.);
+ if (mean3>0) mean3=TMath::Power(TMath::Abs(mean3),1./3.);
+ if (meanS>0) meanS=TMath::Power(TMath::Abs(meanS),3.);
+ //
+ if (mode==1) mean=TMath::Exp(mean);
+ if (mode==1) meanL=TMath::Exp(meanL); // upper truncation
+ if (mode==1) meanD=TMath::Exp(meanD); // lower truncation
+ //
+ delete [] ampWithBelow;
+
+
//
if (returnVal==1) return rms;
if (returnVal==2) return ncl;
+ if (returnVal==3) return Double_t(nclBelowThr)/Double_t(nclBelowThr+ncl);
+ if (returnVal==4) return meanL;
+ if (returnVal==5) return meanD;
+ if (returnVal==6) return mean2;
+ if (returnVal==7) return mean3;
+ if (returnVal==8) return meanS;
return mean;
}