AliTPCcalibTracks.cxx AliTPCcalibTracks.h - Use THnSparse for cluster residuals - removed obsolete histograms
AliTPCcalibV0.cxx AliTPCcalibV0.h - high pt V0 filtering
AliTPCCorrection.cxx AliTPCCorrection.h - Additional arguemnt in function to create track distortion fit tree
AliTPCExBBShape.cxx AliTPCExBBShape.h - Adding visualization function for magnetic field
AliTPCExBEffectiveSector.cxx AliTPCExBEffectiveSector.h - Modified MakeResidualMap function
AliTPCkalmanAlign.cxx AliTPCkalmanAlign.h - Fit constrain functions
AliTPCReconstructor.cxx AliTPCReconstructor.h - Possibility to use the Altro emulator in reconstruction
Marian
for (Int_t j=0;j<3;++j) x[j]+=dx[j];
}
+void AliTPCCorrection::DistortPointLocal(Float_t x[],const Short_t roc) {
+ //
+ // Distorts the initial coordinates x (cartesian coordinates)
+ // according to the given effect (inherited classes)
+ // roc represents the TPC read out chamber (offline numbering convention)
+ //
+ Float_t gxyz[3]={0,0,0};
+ Double_t alpha = TMath::Pi()*(roc%18+0.5)/18;
+ Double_t ca=TMath::Cos(alpha), sa= TMath::Sin(alpha);
+ gxyz[0]= ca*x[0]+sa*x[1];
+ gxyz[1]= -sa*x[0]+ca*x[1];
+ gxyz[2]= x[2];
+ DistortPoint(gxyz,roc);
+ x[0]= ca*gxyz[0]-sa*gxyz[1];
+ x[1]= +sa*gxyz[0]+ca*gxyz[1];
+ x[2]= gxyz[2];
+}
+void AliTPCCorrection::CorrectPointLocal(Float_t x[],const Short_t roc) {
+ //
+ // Distorts the initial coordinates x (cartesian coordinates)
+ // according to the given effect (inherited classes)
+ // roc represents the TPC read out chamber (offline numbering convention)
+ //
+ Float_t gxyz[3]={0,0,0};
+ Double_t alpha = TMath::Pi()*(roc%18+0.5)/18;
+ Double_t ca=TMath::Cos(alpha), sa= TMath::Sin(alpha);
+ gxyz[0]= ca*x[0]+sa*x[1];
+ gxyz[1]= -sa*x[0]+ca*x[1];
+ gxyz[2]= x[2];
+ CorrectPoint(gxyz,roc);
+ x[0]= ca*gxyz[0]-sa*gxyz[1];
+ x[1]= sa*gxyz[0]+ca*gxyz[1];
+ x[2]= gxyz[2];
+}
+
void AliTPCCorrection::DistortPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
//
// Distorts the initial coordinates x (cartesian coordinates) and stores the new
const Double_t kSigmaZ=0.1;
const Double_t kMaxR=500;
const Double_t kMaxZ=500;
+
const Double_t kMaxZ0=220;
const Double_t kZcut=3;
const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
xyz[1]+=gRandom->Gaus(0,0.000005);
xyz[2]+=gRandom->Gaus(0,0.000005);
if (TMath::Abs(track.GetZ())>kMaxZ0) continue;
- if (TMath::Abs(track.GetX())>kMaxR) break;
+ if (TMath::Abs(track.GetX())<kRTPC0) continue;
+ if (TMath::Abs(track.GetX())>kRTPC1) continue;
AliTrackPoint pIn0; // space point
AliTrackPoint pIn1;
Int_t sector= (xyz[2]>0)? 0:18;
pointArray0.GetPoint(point1,npoints-21);
pointArray0.GetPoint(point2,npoints-11);
pointArray0.GetPoint(point3,npoints-1);
- }
+ }
+ if ((TMath::Abs(point1.GetX()-point3.GetX())+TMath::Abs(point1.GetY()-point3.GetY()))<10){
+ printf("fit points not properly initialized\n");
+ return 0;
+ }
track0 = AliTrackerBase::MakeSeed(point1, point2, point3);
track1 = AliTrackerBase::MakeSeed(point1, point2, point3);
track0->ResetCovariance(10);
pointArray1.GetPoint(pIn1,ipoint);
AliTrackPoint prot0 = pIn0.Rotate(track0->GetAlpha()); // rotate to the local frame - non distoted point
AliTrackPoint prot1 = pIn1.Rotate(track1->GetAlpha()); // rotate to the local frame - distorted point
+ if (TMath::Abs(prot0.GetX())<kRTPC0) continue;
+ if (TMath::Abs(prot0.GetX())>kRTPC1) continue;
//
if (!AliTrackerBase::PropagateTrackTo(track0,prot0.GetX(),kMass,5,kFALSE,kMaxSnp)) break;
if (!AliTrackerBase::PropagateTrackTo(track1,prot0.GetX(),kMass,5,kFALSE,kMaxSnp)) break;
-void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Bool_t debug ){
+void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Int_t offset, Bool_t debug ){
//
// Make a fit tree:
// For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
//
const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
// const Double_t kB2C=-0.299792458e-3;
- const Int_t kMinEntries=10;
+ const Int_t kMinEntries=20;
Double_t phi,theta, snp, mean,rms, entries,sector,dsec;
- Float_t refX;
+ Float_t refX;
+ Int_t run;
+ tinput->SetBranchAddress("run",&run);
tinput->SetBranchAddress("theta",&theta);
tinput->SetBranchAddress("phi", &phi);
tinput->SetBranchAddress("snp",&snp);
tinput->SetBranchAddress("sector",§or);
tinput->SetBranchAddress("dsec",&dsec);
tinput->SetBranchAddress("refX",&refX);
- TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortion%d_%d.root",dtype,ptype));
+ TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortion%d_%d_%d.root",dtype,ptype,offset));
//
Int_t nentries=tinput->GetEntries();
Int_t ncorr=corrArray->GetEntries();
if (dtype==3) { dir=1;}
if (dtype==4) { dir=-1;}
//
- for (Int_t ientry=0; ientry<nentries; ientry+=step){
+ for (Int_t ientry=offset; ientry<nentries; ientry+=step){
tinput->GetEntry(ientry);
if (TMath::Abs(snp)>kMaxSnp) continue;
tPar[0]=0;
Double_t xyz[3];
track.GetXYZ(xyz);
Int_t id=0;
+ Double_t pt=1./tPar[4];
Double_t dRrec=0; // dummy value - needed for points - e.g for laser
//if (ptype==4 &&bz<0) mean*=-1; // interpret as curvature -- COMMENTED out - in lookup signed 1/pt used
+ Double_t refXD=refX;
(*pcstream)<<"fit"<<
+ "run="<<run<< // run number
"bz="<<bz<< // magnetic filed used
"dtype="<<dtype<< // detector match type
"ptype="<<ptype<< // parameter type
"rms="<<rms<< // rms
"sector="<<sector<<
"dsec="<<dsec<<
- "refX="<<refX<< // referece X
+ "refX="<<refXD<< // referece X as double
"gx="<<xyz[0]<< // global position at reference
"gy="<<xyz[1]<< // global position at reference
"gz="<<xyz[2]<< // global position at reference
"dRrec="<<dRrec<< // delta Radius in reconstruction
+ "pt="<<pt<< // pt
"id="<<id<< // track id
"entries="<<entries;// number of entries in bin
//
Bool_t isOK=kTRUE;
+ if (entries<kMinEntries) isOK=kFALSE;
+ //
if (dtype!=4) for (Int_t icorr=0; icorr<ncorr; icorr++) {
AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
corrections[icorr]=0;
}
//if (ptype==4 &&bz<0) corrections[icorr]*=-1; // interpret as curvature - commented out
}
- Double_t dRdummy=0;
(*pcstream)<<"fit"<<
- Form("%s=",corr->GetName())<<corrections[icorr]<< // dump correction value
- Form("dR%s=",corr->GetName())<<dRdummy; // dump dummy correction value not needed for tracks
- // for points it is neccessary
+ Form("%s=",corr->GetName())<<corrections[icorr]; // dump correction value
}
if (dtype==4) for (Int_t icorr=0; icorr<ncorr; icorr++) {
AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
corrections[icorr]=0;
if (entries>kMinEntries){
- AliExternalTrackParam trackIn0(refX,phi,tPar,cov);
- AliExternalTrackParam trackIn1(refX,phi,tPar,cov);
+ AliExternalTrackParam trackIn0(refX,phi,tPar,cov); //Outer - direction to vertex
+ AliExternalTrackParam trackIn1(refX,phi,tPar,cov); //Inner - direction magnet
AliExternalTrackParam *trackOut0 = 0;
AliExternalTrackParam *trackOut1 = 0;
//
//
corrections[icorr] = (trackOut0->GetParameter()[ptype]-trackIn0.GetParameter()[ptype]);
corrections[icorr]-= (trackOut1->GetParameter()[ptype]-trackIn1.GetParameter()[ptype]);
+ if (isOK)
+ if ((TMath::Abs(trackOut0->GetX()-trackOut1->GetX())>0.1)||
+ (TMath::Abs(trackOut0->GetX()-trackIn1.GetX())>0.1)||
+ (TMath::Abs(trackOut0->GetAlpha()-trackOut1->GetAlpha())>0.00001)||
+ (TMath::Abs(trackOut0->GetAlpha()-trackIn1.GetAlpha())>0.00001)||
+ (TMath::Abs(trackIn0.GetTgl()-trackIn1.GetTgl())>0.0001)||
+ (TMath::Abs(trackIn0.GetSnp()-trackIn1.GetSnp())>0.0001)
+ ){
+ isOK=kFALSE;
+ }
delete trackOut0;
- delete trackOut1;
+ delete trackOut1;
}else{
corrections[icorr]=0;
isOK=kFALSE;
//
//if (ptype==4 &&bz<0) corrections[icorr]*=-1; // interpret as curvature - commented out no in lookup
}
- Double_t dRdummy=0;
(*pcstream)<<"fit"<<
- Form("%s=",corr->GetName())<<corrections[icorr]<< // dump correction value
- Form("dR%s=",corr->GetName())<<dRdummy; // dump dummy correction value not needed for tracks
- // for points it is neccessary
+ Form("%s=",corr->GetName())<<corrections[icorr]; // dump correction value
}
//
(*pcstream)<<"fit"<<"isOK="<<isOK<<"\n";
-void AliTPCCorrection::MakeLaserDistortionTree(TTree* tree, TObjArray *corrArray, Int_t itype){
+void AliTPCCorrection::MakeSectorDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Int_t offset, Bool_t debug ){
+ //
+ // Make a fit tree:
+ // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
+ // calculates partial distortions
+ // Partial distortion is stored in the resulting tree
+ // Output is storred in the file distortion_<dettype>_<partype>.root
+ // Partial distortion is stored with the name given by correction name
+ //
+ //
+ // Parameters of function:
+ // input - input tree
+ // dtype - distortion type 10 - IROC-OROC
+ // ppype - parameter type
+ // corrArray - array with partial corrections
+ // step - skipe entries - if 1 all entries processed - it is slow
+ // debug 0 if debug on also space points dumped - it is slow
+
+ const Double_t kMaxSnp = 0.8;
+ const Int_t kMinEntries=200;
+ // AliTPCROC *tpcRoc =AliTPCROC::Instance();
+ //
+ const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
+ // const Double_t kB2C=-0.299792458e-3;
+ Double_t phi,theta, snp, mean,rms, entries,sector,dsec,globalZ;
+ Int_t isec1, isec0;
+ Double_t refXD;
+ Float_t refX;
+ Int_t run;
+ tinput->SetBranchAddress("run",&run);
+ tinput->SetBranchAddress("theta",&theta);
+ tinput->SetBranchAddress("phi", &phi);
+ tinput->SetBranchAddress("snp",&snp);
+ tinput->SetBranchAddress("mean",&mean);
+ tinput->SetBranchAddress("rms",&rms);
+ tinput->SetBranchAddress("entries",&entries);
+ tinput->SetBranchAddress("sector",§or);
+ tinput->SetBranchAddress("dsec",&dsec);
+ tinput->SetBranchAddress("refX",&refXD);
+ tinput->SetBranchAddress("z",&globalZ);
+ tinput->SetBranchAddress("isec0",&isec0);
+ tinput->SetBranchAddress("isec1",&isec1);
+ TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortionSector%d_%d_%d.root",dtype,ptype,offset));
+ //
+ Int_t nentries=tinput->GetEntries();
+ Int_t ncorr=corrArray->GetEntries();
+ Double_t corrections[100]={0}; //
+ Double_t tPar[5];
+ Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
+ Int_t dir=0;
+ //
+ for (Int_t ientry=offset; ientry<nentries; ientry+=step){
+ tinput->GetEntry(ientry);
+ refX=refXD;
+ Int_t id=-1;
+ if (TMath::Abs(TMath::Abs(isec0%18)-TMath::Abs(isec1%18))==0) id=1; // IROC-OROC - opposite side
+ if (TMath::Abs(TMath::Abs(isec0%36)-TMath::Abs(isec1%36))==0) id=2; // IROC-OROC - same side
+ if (dtype==10 && id==-1) continue;
+ //
+ dir=-1;
+ tPar[0]=0;
+ tPar[1]=globalZ;
+ tPar[2]=snp;
+ tPar[3]=theta;
+ tPar[4]=(gRandom->Rndm()-0.1)*0.2; //
+ Double_t pt=1./tPar[4];
+ //
+ printf("%f\t%f\t%f\t%f\t%f\t%f\n",entries, sector,theta,snp, mean,rms);
+ Double_t bz=AliTrackerBase::GetBz();
+ AliExternalTrackParam track(refX,phi,tPar,cov);
+ Double_t xyz[3],xyzIn[3],xyzOut[3];
+ track.GetXYZ(xyz);
+ track.GetXYZAt(85,bz,xyzIn);
+ track.GetXYZAt(245,bz,xyzOut);
+ Double_t phiIn = TMath::ATan2(xyzIn[1],xyzIn[0]);
+ Double_t phiOut = TMath::ATan2(xyzOut[1],xyzOut[0]);
+ Double_t phiRef = TMath::ATan2(xyz[1],xyz[0]);
+ Int_t sectorRef = TMath::Nint(9.*phiRef/TMath::Pi()-0.5);
+ Int_t sectorIn = TMath::Nint(9.*phiIn/TMath::Pi()-0.5);
+ Int_t sectorOut = TMath::Nint(9.*phiOut/TMath::Pi()-0.5);
+ //
+ Bool_t isOK=kTRUE;
+ if (sectorIn!=sectorOut) isOK=kFALSE; // requironment - cluster in the same sector
+ if (sectorIn!=sectorRef) isOK=kFALSE; // requironment - cluster in the same sector
+ if (entries<kMinEntries/(1+TMath::Abs(globalZ/100.))) isOK=kFALSE; // requironment - minimal amount of tracks in bin
+ // Do downscale
+ if (TMath::Abs(theta)>1) isOK=kFALSE;
+ //
+ Double_t dRrec=0; // dummy value - needed for points - e.g for laser
+ //
+ (*pcstream)<<"fit"<<
+ "run="<<run<< //run
+ "bz="<<bz<< // magnetic filed used
+ "dtype="<<dtype<< // detector match type
+ "ptype="<<ptype<< // parameter type
+ "theta="<<theta<< // theta
+ "phi="<<phi<< // phi
+ "snp="<<snp<< // snp
+ "mean="<<mean<< // mean dist value
+ "rms="<<rms<< // rms
+ "sector="<<sector<<
+ "dsec="<<dsec<<
+ "refX="<<refXD<< // referece X
+ "gx="<<xyz[0]<< // global position at reference
+ "gy="<<xyz[1]<< // global position at reference
+ "gz="<<xyz[2]<< // global position at reference
+ "dRrec="<<dRrec<< // delta Radius in reconstruction
+ "pt="<<pt<< //pt
+ "id="<<id<< // track id
+ "entries="<<entries;// number of entries in bin
+ //
+ AliExternalTrackParam *trackOut0 = 0;
+ AliExternalTrackParam *trackOut1 = 0;
+ AliExternalTrackParam *ptrackIn0 = 0;
+ AliExternalTrackParam *ptrackIn1 = 0;
+
+ for (Int_t icorr=0; icorr<ncorr; icorr++) {
+ //
+ // special case of the TPC tracks crossing the CE
+ //
+ AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
+ corrections[icorr]=0;
+ if (entries>kMinEntries &&isOK){
+ AliExternalTrackParam trackIn0(refX,phi,tPar,cov);
+ AliExternalTrackParam trackIn1(refX,phi,tPar,cov);
+ ptrackIn1=&trackIn0;
+ ptrackIn0=&trackIn1;
+ //
+ if (debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,pcstream);
+ if (!debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,0);
+ if (debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,pcstream);
+ if (!debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,0);
+ //
+ if (trackOut0 && trackOut1){
+ //
+ if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kTRUE,kMaxSnp)) isOK=kFALSE;
+ if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
+ // rotate all tracks to the same frame
+ if (!trackOut0->Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
+ if (!trackIn1.Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
+ if (!trackOut1->Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
+ //
+ if (!AliTrackerBase::PropagateTrackTo(trackOut0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
+ if (!AliTrackerBase::PropagateTrackTo(&trackIn1,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
+ if (!AliTrackerBase::PropagateTrackTo(trackOut1,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
+ //
+ corrections[icorr] = (trackOut0->GetParameter()[ptype]-trackIn0.GetParameter()[ptype]);
+ corrections[icorr]-= (trackOut1->GetParameter()[ptype]-trackIn1.GetParameter()[ptype]);
+ (*pcstream)<<"fitDebug"<< // just to debug the correction
+ "mean="<<mean<<
+ "pIn0.="<<ptrackIn0<<
+ "pIn1.="<<ptrackIn1<<
+ "pOut0.="<<trackOut0<<
+ "pOut1.="<<trackOut1<<
+ "refX="<<refXD<<
+ "\n";
+ delete trackOut0;
+ delete trackOut1;
+ }else{
+ corrections[icorr]=0;
+ isOK=kFALSE;
+ }
+ }
+ (*pcstream)<<"fit"<<
+ Form("%s=",corr->GetName())<<corrections[icorr]; // dump correction value
+ }
+ //
+ (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n";
+ }
+ delete pcstream;
+}
+
+
+
+void AliTPCCorrection::MakeLaserDistortionTreeOld(TTree* tree, TObjArray *corrArray, Int_t itype){
//
// Make a laser fit tree for global minimization
//
const Double_t cutErrY=0.1;
const Double_t cutErrZ=0.1;
const Double_t kEpsilon=0.00000001;
+ const Double_t kMaxDist=1.; // max distance - space correction
+ const Double_t kMaxRMS=0.05; // max distance -between point and local mean
TVectorD *vecdY=0;
TVectorD *vecdZ=0;
TVectorD *veceY=0;
}
TVectorD * delta= (itype==0)? vecdY:vecdZ;
TVectorD * err= (itype==0)? veceY:veceZ;
-
- for (Int_t irow=0; irow<159; irow++){
- Int_t nentries = 1000;
- if (veceY->GetMatrixArray()[irow]>cutErrY||veceZ->GetMatrixArray()[irow]>cutErrZ) nentries=0;
- if (veceY->GetMatrixArray()[irow]<kEpsilon||veceZ->GetMatrixArray()[irow]<kEpsilon) nentries=0;
- Int_t dtype=5;
- Double_t phi =(*ltr->GetVecPhi())[irow];
- Double_t theta =ltr->GetTgl();
- Double_t mean=delta->GetMatrixArray()[irow];
- Double_t gx=0,gy=0,gz=0;
- Double_t snp = (*ltr->GetVecP2())[irow];
- Double_t rms = 0.1+err->GetMatrixArray()[irow];
- gx = (*ltr->GetVecGX())[irow];
- gy = (*ltr->GetVecGY())[irow];
- gz = (*ltr->GetVecGZ())[irow];
- Int_t bundle= ltr->GetBundle();
- Double_t dRrec=0;
- //
- // get delta R used in reconstruction
- AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
- AliTPCCorrection * correction = calib->GetTPCComposedCorrection(AliTrackerBase::GetBz());
- // const AliTPCRecoParam * recoParam = calib->GetTransform()->GetCurrentRecoParam();
- //Double_t xyz0[3]={gx,gy,gz};
- Double_t oldR=TMath::Sqrt(gx*gx+gy*gy);
- Double_t fphi = TMath::ATan2(gy,gx);
- Double_t fsector = 9.*fphi/TMath::Pi();
- if (fsector<0) fsector+=18;
- Double_t dsec = fsector-Int_t(fsector)-0.5;
- Double_t refX=0;
- //
- if (1 && oldR>1) {
- Float_t xyz1[3]={gx,gy,gz};
- Int_t sector=(gz>0)?0:18;
- correction->CorrectPoint(xyz1, sector);
- refX=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
- dRrec=oldR-refX;
- }
-
- (*pcstream)<<"fit"<<
- "bz="<<bz<< // magnetic filed used
- "dtype="<<dtype<< // detector match type
- "ptype="<<itype<< // parameter type
- "theta="<<theta<< // theta
- "phi="<<phi<< // phi
- "snp="<<snp<< // snp
- "mean="<<mean<< // mean dist value
- "rms="<<rms<< // rms
- "sector="<<fsector<<
- "dsec="<<dsec<<
+ TLinearFitter fitter(2,"pol1");
+ for (Int_t iter=0; iter<2; iter++){
+ Double_t kfit0=0, kfit1=0;
+ Int_t npoints=fitter.GetNpoints();
+ if (npoints>80){
+ fitter.Eval();
+ kfit0=fitter.GetParameter(0);
+ kfit1=fitter.GetParameter(1);
+ }
+ for (Int_t irow=0; irow<159; irow++){
+ Bool_t isOK=kTRUE;
+ Int_t isOKF=0;
+ Int_t nentries = 1000;
+ if (veceY->GetMatrixArray()[irow]>cutErrY||veceZ->GetMatrixArray()[irow]>cutErrZ) nentries=0;
+ if (veceY->GetMatrixArray()[irow]<kEpsilon||veceZ->GetMatrixArray()[irow]<kEpsilon) nentries=0;
+ Int_t dtype=5;
+ Double_t array[10];
+ Int_t first3=TMath::Max(irow-3,0);
+ Int_t last3 =TMath::Min(irow+3,159);
+ Int_t counter=0;
+ if ((*ltr->GetVecSec())[irow]>=0 && err) {
+ for (Int_t jrow=first3; jrow<=last3; jrow++){
+ if ((*ltr->GetVecSec())[irow]!= (*ltr->GetVecSec())[jrow]) continue;
+ if ((*err)[jrow]<kEpsilon) continue;
+ array[counter]=(*delta)[jrow];
+ counter++;
+ }
+ }
+ Double_t rms3 = 0;
+ Double_t mean3 = 0;
+ if (counter>2){
+ rms3 = TMath::RMS(counter,array);
+ mean3 = TMath::Mean(counter,array);
+ }else{
+ isOK=kFALSE;
+ }
+ Double_t phi =(*ltr->GetVecPhi())[irow];
+ Double_t theta =ltr->GetTgl();
+ Double_t mean=delta->GetMatrixArray()[irow];
+ Double_t gx=0,gy=0,gz=0;
+ Double_t snp = (*ltr->GetVecP2())[irow];
+ Double_t dRrec=0;
+ // Double_t rms = err->GetMatrixArray()[irow];
//
- "refX="<<refX<< // reference radius
- "gx="<<gx<< // global position
- "gy="<<gy<< // global position
- "gz="<<gz<< // global position
- "dRrec="<<dRrec<< // delta Radius in reconstruction
- "id="<<bundle<< //bundle
- "entries="<<nentries;// number of entries in bin
- //
- //
- Double_t ky = TMath::Tan(TMath::ASin(snp));
- Int_t ncorr = corrArray->GetEntries();
- Double_t r0 = TMath::Sqrt(gx*gx+gy*gy);
- Double_t phi0 = TMath::ATan2(gy,gx);
- Double_t distortions[1000]={0};
- Double_t distortionsR[1000]={0};
- for (Int_t icorr=0; icorr<ncorr; icorr++) {
- AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
- Float_t distPoint[3]={gx,gy,gz};
- Int_t sector= (gz>0)? 0:18;
- if (r0>80){
- corr->DistortPoint(distPoint, sector);
+ gx = (*ltr->GetVecGX())[irow];
+ gy = (*ltr->GetVecGY())[irow];
+ gz = (*ltr->GetVecGZ())[irow];
+ //
+ // get delta R used in reconstruction
+ AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
+ AliTPCCorrection * correction = calib->GetTPCComposedCorrection(AliTrackerBase::GetBz());
+ // const AliTPCRecoParam * recoParam = calib->GetTransform()->GetCurrentRecoParam();
+ //Double_t xyz0[3]={gx,gy,gz};
+ Double_t oldR=TMath::Sqrt(gx*gx+gy*gy);
+ Double_t fphi = TMath::ATan2(gy,gx);
+ Double_t fsector = 9.*fphi/TMath::Pi();
+ if (fsector<0) fsector+=18;
+ Double_t dsec = fsector-Int_t(fsector)-0.5;
+ Double_t refX=0;
+ Int_t id= ltr->GetId();
+ Double_t pt=0;
+ //
+ if (1 && oldR>1) {
+ Float_t xyz1[3]={gx,gy,gz};
+ Int_t sector=(gz>0)?0:18;
+ correction->CorrectPoint(xyz1, sector);
+ refX=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
+ dRrec=oldR-refX;
+ }
+ if (TMath::Abs(rms3)>kMaxRMS) isOK=kFALSE;
+ if (TMath::Abs(mean-mean3)>kMaxRMS) isOK=kFALSE;
+ if (counter<4) isOK=kFALSE;
+ if (npoints<90) isOK=kFALSE;
+ if (isOK){
+ fitter.AddPoint(&refX,mean);
}
- // Double_t value=distPoint[2]-gz;
- if (itype==0 && r0>1){
- Double_t r1 = TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
- Double_t phi1 = TMath::ATan2(distPoint[1],distPoint[0]);
- Double_t drphi= r0*(phi1-phi0);
- Double_t dr = r1-r0;
- distortions[icorr] = drphi-ky*dr;
- distortionsR[icorr] = dr;
+ Double_t deltaF=kfit0+kfit1*refX;
+ if (iter==1){
+ (*pcstream)<<"fitFull"<< // dumpe also intermediate results
+ "bz="<<bz<< // magnetic filed used
+ "dtype="<<dtype<< // detector match type
+ "ptype="<<itype<< // parameter type
+ "theta="<<theta<< // theta
+ "phi="<<phi<< // phi
+ "snp="<<snp<< // snp
+ "mean="<<mean3<< // mean dist value
+ "rms="<<rms3<< // rms
+ "deltaF="<<deltaF<<
+ "npoints="<<npoints<< //number of points
+ "mean3="<<mean3<< // mean dist value
+ "rms3="<<rms3<< // rms
+ "counter="<<counter<<
+ "sector="<<fsector<<
+ "dsec="<<dsec<<
+ //
+ "refX="<<refX<< // reference radius
+ "gx="<<gx<< // global position
+ "gy="<<gy<< // global position
+ "gz="<<gz<< // global position
+ "dRrec="<<dRrec<< // delta Radius in reconstruction
+ "id="<<id<< //bundle
+ "entries="<<nentries<<// number of entries in bin
+ "\n";
+ }
+ if (iter==1) (*pcstream)<<"fit"<< // dump valus for fit
+ "bz="<<bz<< // magnetic filed used
+ "dtype="<<dtype<< // detector match type
+ "ptype="<<itype<< // parameter type
+ "theta="<<theta<< // theta
+ "phi="<<phi<< // phi
+ "snp="<<snp<< // snp
+ "mean="<<mean3<< // mean dist value
+ "rms="<<rms3<< // rms
+ "sector="<<fsector<<
+ "dsec="<<dsec<<
+ //
+ "refX="<<refX<< // reference radius
+ "gx="<<gx<< // global position
+ "gy="<<gy<< // global position
+ "gz="<<gz<< // global position
+ "dRrec="<<dRrec<< // delta Radius in reconstruction
+ "pt="<<pt<< //pt
+ "id="<<id<< //bundle
+ "entries="<<nentries;// number of entries in bin
+ //
+ //
+ Double_t ky = TMath::Tan(TMath::ASin(snp));
+ Int_t ncorr = corrArray->GetEntries();
+ Double_t r0 = TMath::Sqrt(gx*gx+gy*gy);
+ Double_t phi0 = TMath::ATan2(gy,gx);
+ Double_t distortions[1000]={0};
+ Double_t distortionsR[1000]={0};
+ if (iter==1){
+ for (Int_t icorr=0; icorr<ncorr; icorr++) {
+ AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
+ Float_t distPoint[3]={gx,gy,gz};
+ Int_t sector= (gz>0)? 0:18;
+ if (r0>80){
+ corr->DistortPoint(distPoint, sector);
+ }
+ // Double_t value=distPoint[2]-gz;
+ if (itype==0 && r0>1){
+ Double_t r1 = TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
+ Double_t phi1 = TMath::ATan2(distPoint[1],distPoint[0]);
+ Double_t drphi= r0*(phi1-phi0);
+ Double_t dr = r1-r0;
+ distortions[icorr] = drphi-ky*dr;
+ distortionsR[icorr] = dr;
+ }
+ if (TMath::Abs(distortions[icorr])>kMaxDist) {isOKF=icorr+1; isOK=kFALSE; }
+ if (TMath::Abs(distortionsR[icorr])>kMaxDist) {isOKF=icorr+1; isOK=kFALSE;}
+ (*pcstream)<<"fit"<<
+ Form("%s=",corr->GetName())<<distortions[icorr]; // dump correction value
+ }
+ (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n";
}
- (*pcstream)<<"fit"<<
- Form("%s=",corr->GetName())<<distortions[icorr]<< // dump correction value
- Form("dR%s=",corr->GetName())<<distortionsR[icorr]; // dump correction R value
}
- (*pcstream)<<"fit"<<"\n";
}
}
delete pcstream;
if (axisType==2) return distPoint[2]-gz;
return phi1-phi0;
}
+
+
+
+
+
+void AliTPCCorrection::MakeLaserDistortionTree(TTree* tree, TObjArray *corrArray, Int_t itype){
+ //
+ // Make a laser fit tree for global minimization
+ //
+ AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
+ AliTPCCorrection * correction = calib->GetTPCComposedCorrection();
+ if (!correction) correction = calib->GetTPCComposedCorrection(AliTrackerBase::GetBz());
+ correction->AddVisualCorrection(correction,0); //register correction
+
+ AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
+ AliTPCParam *param = AliTPCcalibDB::Instance()->GetParameters();
+ //
+ const Double_t cutErrY=0.05;
+ const Double_t kSigmaCut=4;
+ // const Double_t cutErrZ=0.03;
+ const Double_t kEpsilon=0.00000001;
+ const Double_t kMaxDist=1.; // max distance - space correction
+ TVectorD *vecdY=0;
+ TVectorD *vecdZ=0;
+ TVectorD *veceY=0;
+ TVectorD *veceZ=0;
+ AliTPCLaserTrack *ltr=0;
+ AliTPCLaserTrack::LoadTracks();
+ tree->SetBranchAddress("dY.",&vecdY);
+ tree->SetBranchAddress("dZ.",&vecdZ);
+ tree->SetBranchAddress("eY.",&veceY);
+ tree->SetBranchAddress("eZ.",&veceZ);
+ tree->SetBranchAddress("LTr.",<r);
+ Int_t entries= tree->GetEntries();
+ TTreeSRedirector *pcstream= new TTreeSRedirector("distortionLaser_0.root");
+ Double_t bz=AliTrackerBase::GetBz();
+ //
+ Double_t globalXYZ[3];
+ Double_t globalXYZCorr[3];
+ for (Int_t ientry=0; ientry<entries; ientry++){
+ tree->GetEntry(ientry);
+ if (!ltr->GetVecGX()){
+ ltr->UpdatePoints();
+ }
+ //
+ TVectorD fit10(5);
+ TVectorD fit5(5);
+ printf("Entry\t%d\n",ientry);
+ for (Int_t irow0=0; irow0<158; irow0+=1){
+ //
+ TLinearFitter fitter10(4,"hyp3");
+ TLinearFitter fitter5(2,"hyp1");
+ Int_t sector= (Int_t)(*ltr->GetVecSec())[irow0];
+ if (sector<0) continue;
+ //if (TMath::Abs(vecdY->GetMatrixArray()[irow0])<kEpsilon) continue;
+
+ Double_t refX= (*ltr->GetVecLX())[irow0];
+ Int_t firstRow1 = TMath::Max(irow0-10,0);
+ Int_t lastRow1 = TMath::Min(irow0+10,158);
+ Double_t padWidth=(irow0<64)?0.4:0.6;
+ // make long range fit
+ for (Int_t irow1=firstRow1; irow1<=lastRow1; irow1++){
+ if (TMath::Abs((*ltr->GetVecSec())[irow1]-sector)>kEpsilon) continue;
+ if (veceY->GetMatrixArray()[irow1]>cutErrY) continue;
+ if (TMath::Abs(vecdY->GetMatrixArray()[irow1])<kEpsilon) continue;
+ Double_t idealX= (*ltr->GetVecLX())[irow1];
+ Double_t idealY= (*ltr->GetVecLY())[irow1];
+ Double_t idealZ= (*ltr->GetVecLZ())[irow1];
+ Double_t gx= (*ltr->GetVecGX())[irow1];
+ Double_t gy= (*ltr->GetVecGY())[irow1];
+ Double_t gz= (*ltr->GetVecGZ())[irow1];
+ Double_t measY=(*vecdY)[irow1]+idealY;
+ Double_t deltaR = GetCorrXYZ(gx, gy, gz, 0,0);
+ // deltaR = R distorted -R ideal
+ Double_t xxx[4]={idealX+deltaR-refX,TMath::Cos(idealY/padWidth), TMath::Sin(idealY/padWidth)};
+ fitter10.AddPoint(xxx,measY,1);
+ }
+ Bool_t isOK=kTRUE;
+ Double_t rms10=0;//TMath::Sqrt(fitter10.GetChisquare()/(fitter10.GetNpoints()-4));
+ Double_t mean10 =0;// fitter10.GetParameter(0);
+ Double_t slope10 =0;// fitter10.GetParameter(0);
+ Double_t cosPart10 = 0;// fitter10.GetParameter(2);
+ Double_t sinPart10 =0;// fitter10.GetParameter(3);
+
+ if (fitter10.GetNpoints()>10){
+ fitter10.Eval();
+ rms10=TMath::Sqrt(fitter10.GetChisquare()/(fitter10.GetNpoints()-4));
+ mean10 = fitter10.GetParameter(0);
+ slope10 = fitter10.GetParameter(1);
+ cosPart10 = fitter10.GetParameter(2);
+ sinPart10 = fitter10.GetParameter(3);
+ //
+ // make short range fit
+ //
+ for (Int_t irow1=firstRow1+5; irow1<=lastRow1-5; irow1++){
+ if (TMath::Abs((*ltr->GetVecSec())[irow1]-sector)>kEpsilon) continue;
+ if (veceY->GetMatrixArray()[irow1]>cutErrY) continue;
+ if (TMath::Abs(vecdY->GetMatrixArray()[irow1])<kEpsilon) continue;
+ Double_t idealX= (*ltr->GetVecLX())[irow1];
+ Double_t idealY= (*ltr->GetVecLY())[irow1];
+ Double_t idealZ= (*ltr->GetVecLZ())[irow1];
+ Double_t gx= (*ltr->GetVecGX())[irow1];
+ Double_t gy= (*ltr->GetVecGY())[irow1];
+ Double_t gz= (*ltr->GetVecGZ())[irow1];
+ Double_t measY=(*vecdY)[irow1]+idealY;
+ Double_t deltaR = GetCorrXYZ(gx, gy, gz, 0,0);
+ // deltaR = R distorted -R ideal
+ Double_t expY= mean10+slope10*(idealX+deltaR-refX);
+ if (TMath::Abs(measY-expY)>kSigmaCut*rms10) continue;
+ //
+ Double_t corr=cosPart10*TMath::Cos(idealY/padWidth)+sinPart10*TMath::Sin(idealY/padWidth);
+ Double_t xxx[4]={idealX+deltaR-refX,TMath::Cos(idealY/padWidth), TMath::Sin(idealY/padWidth)};
+ fitter5.AddPoint(xxx,measY-corr,1);
+ }
+ }else{
+ isOK=kFALSE;
+ }
+ if (fitter5.GetNpoints()<8) isOK=kFALSE;
+
+ Double_t rms5=0;//TMath::Sqrt(fitter5.GetChisquare()/(fitter5.GetNpoints()-4));
+ Double_t offset5 =0;// fitter5.GetParameter(0);
+ Double_t slope5 =0;// fitter5.GetParameter(0);
+ if (isOK){
+ fitter5.Eval();
+ rms5=TMath::Sqrt(fitter5.GetChisquare()/(fitter5.GetNpoints()-4));
+ offset5 = fitter5.GetParameter(0);
+ slope5 = fitter5.GetParameter(0);
+ }
+ //
+ Double_t dtype=5;
+ Double_t ptype=0;
+ Double_t phi =(*ltr->GetVecPhi())[irow0];
+ Double_t theta =ltr->GetTgl();
+ Double_t mean=(vecdY)->GetMatrixArray()[irow0];
+ Double_t gx=0,gy=0,gz=0;
+ Double_t snp = (*ltr->GetVecP2())[irow0];
+ Int_t bundle= ltr->GetBundle();
+ Int_t id= ltr->GetId();
+ // Double_t rms = err->GetMatrixArray()[irow];
+ //
+ gx = (*ltr->GetVecGX())[irow0];
+ gy = (*ltr->GetVecGY())[irow0];
+ gz = (*ltr->GetVecGZ())[irow0];
+ Double_t dRrec = GetCorrXYZ(gx, gy, gz, 0,0);
+ fitter10.GetParameters(fit10);
+ fitter5.GetParameters(fit5);
+ Double_t idealY= (*ltr->GetVecLY())[irow0];
+ Double_t measY=(*vecdY)[irow0]+idealY;
+ Double_t corr=cosPart10*TMath::Cos(idealY/padWidth)+sinPart10*TMath::Sin(idealY/padWidth);
+ if (TMath::Max(rms5,rms10)>0.06) isOK=kFALSE;
+ //
+ (*pcstream)<<"fitFull"<< // dumpe also intermediate results
+ "bz="<<bz<< // magnetic filed used
+ "dtype="<<dtype<< // detector match type
+ "ptype="<<ptype<< // parameter type
+ "theta="<<theta<< // theta
+ "phi="<<phi<< // phi
+ "snp="<<snp<< // snp
+ "sector="<<sector<<
+ "bundle="<<bundle<<
+// // "dsec="<<dsec<<
+ "refX="<<refX<< // reference radius
+ "gx="<<gx<< // global position
+ "gy="<<gy<< // global position
+ "gz="<<gz<< // global position
+ "dRrec="<<dRrec<< // delta Radius in reconstruction
+ "id="<<id<< //bundle
+ "rms10="<<rms10<<
+ "rms5="<<rms5<<
+ "fit10.="<<&fit10<<
+ "fit5.="<<&fit5<<
+ "measY="<<measY<<
+ "mean="<<mean<<
+ "idealY="<<idealY<<
+ "corr="<<corr<<
+ "isOK="<<isOK<<
+ "\n";
+ }
+ }
+ delete pcstream;
+}
// functions to correct a space point
void CorrectPoint ( Float_t x[],const Short_t roc);
+ void CorrectPointLocal(Float_t x[],const Short_t roc);
void CorrectPoint (const Float_t x[],const Short_t roc,Float_t xp[]);
virtual void GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]);
// functions to distort a space point
void DistortPoint ( Float_t x[],const Short_t roc);
+ void DistortPointLocal(Float_t x[],const Short_t roc);
void DistortPoint (const Float_t x[],const Short_t roc,Float_t xp[]);
virtual void GetDistortion(const Float_t x[],const Short_t roc,Float_t dx[]);
virtual void SetOmegaTauT1T2(Float_t omegaTau,Float_t t1,Float_t t2);
AliExternalTrackParam * FitDistortedTrack(AliExternalTrackParam & trackIn, Double_t refX, Int_t dir,TTreeSRedirector *pcstream);
void StoreInOCDB(Int_t startRun, Int_t endRun, const char *comment=0);
- static void MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step=1, Bool_t debug=0);
+ static void MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step=1, Int_t offset=0, Bool_t debug=0);
+ static void MakeSectorDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step=1, Int_t offset=0, Bool_t debug=0);
+ static void MakeLaserDistortionTreeOld(TTree* tree, TObjArray *corrArray, Int_t itype);
static void MakeLaserDistortionTree(TTree* tree, TObjArray *corrArray, Int_t itype);
void FastSimDistortedVertex(Double_t orgVertex[3], Int_t nTracks, AliESDVertex &aV, AliESDVertex &avOrg, AliESDVertex &cV, AliESDVertex &cvOrg, TTreeSRedirector * const pcstream, Double_t etaCuts);
printf(" - C1: %1.4f, C2: %1.4f \n",fC1,fC2);
}
}
+
+Double_t AliTPCExBBShape::GetBFieldXYZ(Double_t gx, Double_t gy, Double_t gz, Int_t axisType){
+ //
+ // return B field at given x,y,z
+ //
+ AliMagF* field = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
+ if (!field) return 0;
+ Double_t xyz[3]={gx,gy,gz};
+ Double_t bxyz[3]={0};
+ field->Field(xyz,bxyz);
+ //
+ Double_t r=TMath::Sqrt(gx*gx+gy*gy);
+ Double_t b=TMath::Sqrt(bxyz[0]*bxyz[0]+bxyz[1]*bxyz[1]);
+ if (axisType==0) {
+ return (xyz[0]*bxyz[1]-xyz[1]*bxyz[0])/(bxyz[2]*r);
+ }
+ if (axisType==1){
+ return (xyz[0]*bxyz[0]+xyz[1]*bxyz[1])/(bxyz[2]*r);
+ }
+ if (axisType==2) return bxyz[2];
+ if (axisType==3) return bxyz[0];
+ if (axisType==4) return bxyz[1];
+ if (axisType==5) return bxyz[2];
+ return bxyz[2];
+}
void GetBxAndByOverBz(const Float_t x[],const Short_t roc,Float_t BxByOverBz[]);
virtual void Print(Option_t* option="") const;
+ static Double_t GetBFieldXYZ(Double_t gx, Double_t gy, Double_t gz, Int_t axisType);
private:
Float_t fC1; // coefficient C1 (compare Jim Thomas's notes for definitions)
}
}
-void AliTPCExBEffectiveSector::MakeResidualMap(THnSparse * hisInput, const char *sname){
+void AliTPCExBEffectiveSector::MakeResidualMap(THnSparse * hisInput, const char *sname, Int_t ptype, Int_t dtype){
//
// Make cluster residual map from the n-dimensional histogram
// hisInput supposed to have given format:
Int_t nbins1=hisInput->GetAxis(1)->GetNbins();
Int_t nbins2=hisInput->GetAxis(2)->GetNbins();
Int_t nbins3=hisInput->GetAxis(3)->GetNbins();
-
+ TF1 *fgaus=0;
TH3F * hisResMap3D =
new TH3F("his3D","his3D",
nbins1,hisInput->GetAxis(1)->GetXmin(), hisInput->GetAxis(1)->GetXmax(),
hisResMap3D->SetBinContent(ibin1,ibin2,ibin3, mean);
Double_t phi=TMath::Pi()*sector/9;
if (phi>TMath::Pi()) phi+=TMath::Pi();
+ Double_t meanG=0;
+ Double_t rmsG=0;
+ if (entries>50){
+ if (!fgaus) {
+ his->Fit("gaus","Q","goff");
+ fgaus= (TF1*)((his->GetListOfFunctions()->FindObject("gaus"))->Clone());
+ fgaus->SetName("gausk");
+ }
+ if (fgaus) his->Fit(fgaus,"Q","goff");
+ if (fgaus){
+ meanG=fgaus->GetParameter(1);
+ rmsG=fgaus->GetParameter(2);
+ }
+ }
+ Double_t dsec=sector-Int_t(sector)-0.5;
+ Double_t snp=dsec*TMath::Pi()/9.;
(*pcstream)<<"delta"<<
+ "ptype="<<ptype<<
+ "dtype="<<dtype<<
"sector="<<sector<<
+ "dsec="<<dsec<<
+ "snp="<<snp<<
"phi="<<phi<<
"localX="<<localX<<
"kZ="<<kZ<<
+ "theta="<<kZ<<
"mean="<<mean<<
"rms="<<rms<<
+ "meanG="<<meanG<<
+ "rmsG="<<rmsG<<
"entries="<<entries<<
"meanA="<<meanA<<
"rmsA="<<rmsA<<
Float_t GetC1() const {return fC1;}
Float_t GetC0() const {return fC0;}
void Print(const Option_t* option) const;
- static void MakeResidualMap(THnSparse * hisInput, const char *sname);
+ static void MakeResidualMap(THnSparse * hisInput, const char *sname, Int_t ptype, Int_t dtype=0);
public:
virtual void GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]);
public:
#include "AliTracker.h"
#include "AliMagF.h"
#include "AliTPCCalROC.h"
+#include "AliESDv0.h"
#include "AliLog.h"
Double_t deltaTime = EndTime - StartTime;
- // main histogram for time dependence: dE/dx, time, type (1-muon cosmic,2-pion beam data), meanDriftlength, momenta (only filled if enough space is available), run number
+ // main histogram for time dependence: dE/dx, time, type (1-muon cosmic,2-pion beam data, 3&4 - proton points at higher dE/dx), meanDriftlength, momenta (only filled if enough space is available), run number, eta
Int_t timeBins = TMath::Nint(deltaTime/deltaIntegrationTimeGain);
- Int_t binsGainTime[6] = {150, timeBins, 2, 25, 200, 10000000};
- Double_t xminGainTime[6] = {0.5, StartTime, 0.5, 0, 0.1, -0.5};
- Double_t xmaxGainTime[6] = { 8, EndTime, 2.5, 250, 50, 9999999.5};
- fHistGainTime = new THnSparseF("HistGainTime","dEdx time dep.;dEdx,time,type,driftlength,momenta,run number;dEdx",6,binsGainTime,xminGainTime,xmaxGainTime);
+ Int_t binsGainTime[7] = {150, timeBins, 4, 25, 200, 10000000, 20};
+ Double_t xminGainTime[7] = {0.5, StartTime, 0.5, 0, 0.1, -0.5, -1};
+ Double_t xmaxGainTime[7] = { 8, EndTime, 4.5, 250, 50, 9999999.5, 1};
+ fHistGainTime = new THnSparseF("HistGainTime","dEdx time dep.;dEdx,time,type,driftlength,momenta,run number, eta;dEdx",7,binsGainTime,xminGainTime,xmaxGainTime);
BinLogX(fHistGainTime, 4);
//
fHistDeDxTotal = new TH2F("DeDx","dEdx; momentum p (GeV); TPC signal (a.u.)",250,0.01,100.,1000,0.,8);
//
// track loop
//
- for (Int_t i=0;i<ntracks;++i) {
+ for (Int_t i=0;i<ntracks;++i) { // begin track loop
AliESDtrack *track = event->GetTrack(i);
if (!track) continue;
// exclude tracks which do not look like primaries or are simply too short or on wrong sectors
if (nclsDeDx < 60) continue;
- if (TMath::Abs(trackIn->GetTgl()) > 1) continue;
- if (TMath::Abs(trackIn->GetSnp()) > 0.6) continue;
+ //if (TMath::Abs(trackIn->GetTgl()) > 1) continue;
+ //if (TMath::Abs(trackIn->GetSnp()) > 0.6) continue;
+ if (TMath::Abs(trackIn->Eta()) > 1) continue;
+ UInt_t status = track->GetStatus();
+ if ((status&AliESDtrack::kTPCrefit)==0) continue;
+ if (track->GetNcls(0) < 3) continue; // ITS clusters
+ Float_t dca[2], cov[3];
+ track->GetImpactParameters(dca,cov);
+ if (dca[0] > 7 || dca[0] < 0.000001) continue; // cut in xy
+ Double_t eta = trackIn->Eta();
// Get seeds
TObject *calibObject;
if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
}
- if (seed) {
- if (fLowMemoryConsumption) {
- if (meanP > 0.5 || meanP < 0.4) continue;
- meanP = 0.45; // set all momenta to one in order to save memory
- }
+ if (seed) {
+ Int_t particleCase = 0;
+ if (meanP > 0.5 || meanP < 0.4) particleCase = 2; // MIP pions
+ if (meanP > 0.56 || meanP < 0.57) particleCase = 3; // protons 1
+ if (meanP > 0.64 || meanP < 0.66) particleCase = 4; // protons 2
+ //
+ if (fLowMemoryConsumption && particleCase == 0) continue;
+ //
Double_t tpcSignal = GetTPCdEdx(seed);
fHistDeDxTotal->Fill(meanP, tpcSignal);
//
- //dE/dx, time, type (1-muon cosmic,2-pion beam data), momenta
- Double_t vec[6] = {tpcSignal,time,2,meanDrift,meanP,runNumber};
+ //dE/dx, time, type (1-muon cosmic,2-pion beam data, 3&4 protons), momenta, runNumner, eta
+ Double_t vec[7] = {tpcSignal,time,particleCase,meanDrift,meanP,runNumber, eta};
fHistGainTime->Fill(vec);
}
+ } // end track loop
+ //
+ // V0 loop -- in beam events the cosmic part of the histogram is filled with GammaConversions
+ //
+ for(Int_t iv0 = 0; iv0 < event->GetNumberOfV0s(); iv0++) {
+ AliESDv0 * v0 = event->GetV0(iv0);
+ if (!v0->GetOnFlyStatus()) continue;
+ if (v0->GetEffMass(0,0) > 0.02) continue; // select low inv. mass
+ Double_t xyz[3];
+ v0->GetXYZ(xyz[0], xyz[1], xyz[2]);
+ if (TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]) < 3 || TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]) > 30) continue;
+ //
+ // "loop over daughters"
+ //
+ for(Int_t idaughter = 0; idaughter < 2; idaughter++) { // daughter loop
+ Int_t index = idaughter == 0 ? v0->GetPindex() : v0->GetNindex();
+ AliESDtrack * trackP = event->GetTrack(index);
+ AliESDfriendTrack *friendTrackP = esdFriend->GetTrack(index);
+ if (!friendTrackP) continue;
+ const AliExternalTrackParam * trackPIn = trackP->GetInnerParam();
+ const AliExternalTrackParam * trackPOut = friendTrackP->GetTPCOut();
+ if (!trackPIn) continue;
+ if (!trackPOut) continue;
+ // calculate necessary track parameters
+ Double_t meanP = trackPIn->GetP();
+ Double_t meanDrift = 250 - 0.5*TMath::Abs(trackPIn->GetZ() + trackPOut->GetZ());
+ Int_t nclsDeDx = trackP->GetTPCNcls();
+ // exclude tracks which do not look like primaries or are simply too short or on wrong sectors
+ if (nclsDeDx < 60) continue;
+ if (TMath::Abs(trackPIn->GetTgl()) > 1) continue;
+ //
+ TObject *calibObject;
+ AliTPCseed *seed = 0;
+ for (Int_t l=0;(calibObject=friendTrackP->GetCalibObject(l));++l) {
+ if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
+ }
+ if (seed) {
+ if (fLowMemoryConsumption) {
+ if (meanP > 0.5 || meanP < 0.4) continue;
+ meanP = 0.45; // set all momenta to one in order to save memory
+ }
+ Double_t tpcSignal = GetTPCdEdx(seed);
+ //dE/dx, time, type (1-muon cosmic,2-pion beam data), momenta
+ Double_t vec[6] = {tpcSignal,time,1,meanDrift,meanP,runNumber};
+ fHistGainTime->Fill(vec);
+ }
+ }
+
}
}
Offline/HLT Offline/HLT OCDB entries (AliTPCClusterParam)
*/
-/*
-
-How to retrive it from file (created using calibration task):
-
-gSystem->Load("libANALYSIS");
-gSystem->Load("libTPCcalib");
-TFile fcalib("CalibObjects.root");
-TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
-AliTPCcalibTracks * calibTracks = ( AliTPCcalibTracks *)array->FindObject("calibTracks");
-
-
-//USAGE of debug stream example
- gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
- gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
- AliXRDPROOFtoolkit tool;
- TChain * chainres = tool.MakeChain("tracks.txt","ResolCl",0,10200);
- chainres->Lookup();
-*/
-
//
///////////////////////////////////////////////////////////////////////////////
#include "AliTPCcalibTracks.h"
#include "AliTPCClusterParam.h"
#include "AliTPCcalibTracksCuts.h"
-#include "AliTPCCalPadRegion.h"
#include "AliTPCCalPad.h"
#include "AliTPCCalROC.h"
#include "TText.h"
#include "TStatToolkit.h"
#include "TCut.h"
#include "THnSparse.h"
+#include "AliRieman.h"
fHisRMSZ(0), // THnSparse - rms Z
fHisQmax(0), // THnSparse - qmax
fHisQtot(0), // THnSparse - qtot
- fArrayAmpRow(0),
- fArrayAmp(0),
fArrayQDY(0),
fArrayQDZ(0),
fArrayQRMSY(0),
fArrayQRMSZ(0),
- fArrayChargeVsDriftlength(0),
- fcalPadRegionChargeVsDriftlength(0),
- fDeltaY(0),
- fDeltaZ(0),
fResolY(0),
fResolZ(0),
fRMSY(0),
fRMSZ(0),
fCuts(0),
- fHclus(0),
fRejectedTracksHisto(0),
- fHclusterPerPadrow(0),
- fHclusterPerPadrowRaw(0),
fClusterCutHisto(0),
fCalPadClusterPerPad(0),
fCalPadClusterPerPadRaw(0)
fHisRMSZ(0), // THnSparse - rms Z
fHisQmax(0), // THnSparse - qmax
fHisQtot(0), // THnSparse - qtot
- fArrayAmpRow(0),
- fArrayAmp(0),
fArrayQDY(0),
fArrayQDZ(0),
fArrayQRMSY(0),
fArrayQRMSZ(0),
- fArrayChargeVsDriftlength(0),
- fcalPadRegionChargeVsDriftlength(0),
- fDeltaY(0),
- fDeltaZ(0),
fResolY(0),
fResolZ(0),
fRMSY(0),
fRMSZ(0),
fCuts(0),
- fHclus(0),
fRejectedTracksHisto(0),
- fHclusterPerPadrow(0),
- fHclusterPerPadrowRaw(0),
fClusterCutHisto(0),
fCalPadClusterPerPad(0),
fCalPadClusterPerPadRaw(0)
TH1::AddDirectory(kFALSE);
Int_t length = -1;
- // backward compatibility: if the data member doesn't yet exist, it will not be merged
- (calibTracks.fArrayAmpRow) ? length = calibTracks.fArrayAmpRow->GetEntriesFast() : length = -1;
- fArrayAmpRow = new TObjArray(length);
- fArrayAmp = new TObjArray(length);
- for (Int_t i = 0; i < length; i++) {
- fArrayAmpRow->AddAt( (TProfile*)calibTracks.fArrayAmpRow->At(i)->Clone(), i);
- fArrayAmp->AddAt( ((TProfile*)calibTracks.fArrayAmp->At(i)->Clone()), i);
- }
(calibTracks.fArrayQDY) ? length = calibTracks.fArrayQDY->GetEntriesFast() : length = -1;
fArrayQDY= new TObjArray(length);
fRMSZ->AddAt( ((TH1F*)calibTracks.fRMSZ->At(i)->Clone()), i);
}
- (calibTracks.fArrayChargeVsDriftlength) ? length = calibTracks.fArrayChargeVsDriftlength->GetEntriesFast() : length = -1;
- (calibTracks.fArrayChargeVsDriftlength) ? fArrayChargeVsDriftlength = new TObjArray(length) : fArrayChargeVsDriftlength = 0;
- for (Int_t i = 0; i < length; i++) {
- fArrayChargeVsDriftlength->AddAt( ((TProfile*)calibTracks.fArrayChargeVsDriftlength->At(i)->Clone()), i);
- }
- fDeltaY = (TH1F*)calibTracks.fDeltaY->Clone();
- fDeltaZ = (TH1F*)calibTracks.fDeltaZ->Clone();
- fHclus = (TH1I*)calibTracks.fHclus->Clone();
fClusterCutHisto = (TH2I*)calibTracks.fClusterCutHisto->Clone();
fRejectedTracksHisto = (TH1I*)calibTracks.fRejectedTracksHisto->Clone();
- fHclusterPerPadrow = (TH1I*)calibTracks.fHclusterPerPadrow->Clone();
- fHclusterPerPadrowRaw = (TH1I*)calibTracks.fHclusterPerPadrowRaw->Clone();
- fcalPadRegionChargeVsDriftlength = (AliTPCCalPadRegion*)calibTracks.fcalPadRegionChargeVsDriftlength->Clone();
fCalPadClusterPerPad = (AliTPCCalPad*)calibTracks.fCalPadClusterPerPad->Clone();
fCalPadClusterPerPadRaw = (AliTPCCalPad*)calibTracks.fCalPadClusterPerPadRaw->Clone();
fHisRMSZ(0), // THnSparse - rms Z
fHisQmax(0), // THnSparse - qmax
fHisQtot(0), // THnSparse - qtot
- fArrayAmpRow(0),
- fArrayAmp(0),
fArrayQDY(0),
fArrayQDZ(0),
fArrayQRMSY(0),
fArrayQRMSZ(0),
- fArrayChargeVsDriftlength(0),
- fcalPadRegionChargeVsDriftlength(0),
- fDeltaY(0),
- fDeltaZ(0),
fResolY(0),
fResolZ(0),
fRMSY(0),
fRMSZ(0),
fCuts(0),
- fHclus(0),
fRejectedTracksHisto(0),
- fHclusterPerPadrow(0),
- fHclusterPerPadrowRaw(0),
fClusterCutHisto(0),
fCalPadClusterPerPad(0),
fCalPadClusterPerPadRaw(0)
TH1::AddDirectory(kFALSE);
- char chname[1000];
- TProfile * prof1=0;
- TH1F * his1 =0;
- fHclus = new TH1I("hclus","Number of clusters per track",160, 0, 160); // valgrind 3
fRejectedTracksHisto = new TH1I("RejectedTracksHisto", "Rejected tracks, sorted by failed cut", 100, -1, 10);
- fHclusterPerPadrow = new TH1I("fHclusterPerPadrow", " clusters per padRow, used for the resolution tree", 160, 0, 160);
- fHclusterPerPadrowRaw = new TH1I("fHclusterPerPadrowRaw", " clusters per padRow, before cutting clusters", 160, 0, 160);
fCalPadClusterPerPad = new AliTPCCalPad("fCalPadClusterPerPad", "clusters per pad");
fCalPadClusterPerPadRaw = new AliTPCCalPad("fCalPadClusterPerPadRaw", "clusters per pad, before cutting clusters");
fClusterCutHisto = new TH2I("fClusterCutHisto", "Cutted cluster over padRow; Cut Criterium; PadRow", 5,1,5, 160,0,159);
- // Amplitude - sector - row histograms
- fArrayAmpRow = new TObjArray(72);
- fArrayAmp = new TObjArray(72);
- fArrayChargeVsDriftlength = new TObjArray(72);
-
- for (Int_t i = 0; i < 36; i++){
- snprintf(chname,100,"Amp_row_Sector%d",i);
- prof1 = new TProfile(chname,chname,63,0,64); // valgrind 3 193,536 bytes in 354 blocks are still reachable
- prof1->SetXTitle("Pad row");
- prof1->SetYTitle("Mean Max amplitude");
- fArrayAmpRow->AddAt(prof1,i);
- snprintf(chname,100,"Amp_row_Sector%d",i+36);
- prof1 = new TProfile(chname,chname,96,0,97); // valgrind 3 3,912 bytes in 6 blocks are possibly lost
- prof1->SetXTitle("Pad row");
- prof1->SetYTitle("Mean Max amplitude");
- fArrayAmpRow->AddAt(prof1,i+36);
-
- // amplitude
- sprintf(chname,"Amp_Sector%d",i);
- his1 = new TH1F(chname,chname,100,0,32); // valgrind
- his1->SetXTitle("Max Amplitude (ADC)");
- fArrayAmp->AddAt(his1,i);
- sprintf(chname,"Amp_Sector%d",i+36);
- his1 = new TH1F(chname,chname,100,0,32); // valgrind 3 13,408,208 bytes in 229 blocks are still reachable
- his1->SetXTitle("Max Amplitude (ADC)");
- fArrayAmp->AddAt(his1,i+36);
-
- // driftlength
- snprintf(chname,100, "driftlengt vs. charge, ROC %i", i);
- prof1 = new TProfile(chname, chname, 25, 0, 250);
- prof1->SetYTitle("Charge");
- prof1->SetXTitle("Driftlength");
- fArrayChargeVsDriftlength->AddAt(prof1,i);
- snprintf(chname,100, "driftlengt vs. charge, ROC %i", i+36);
- prof1 = new TProfile(chname, chname, 25, 0, 250);
- prof1->SetYTitle("Charge");
- prof1->SetXTitle("Driftlength");
- fArrayChargeVsDriftlength->AddAt(prof1,i+36);
- }
TH1::AddDirectory(kFALSE);
- fDeltaY = new TH1F("DeltaY","DeltaY",100,-1,1);
- fDeltaZ = new TH1F("DeltaZ","DeltaZ",100,-1,1);
fResolY = new TObjArray(3);
fResolZ = new TObjArray(3);
}
}
- fcalPadRegionChargeVsDriftlength = new AliTPCCalPadRegion("fcalPadRegionChargeVsDriftlength", "TProfiles with charge vs driftlength for each pad region");
- TProfile *tempProf;
- for (UInt_t padSize = 0; padSize < 3; padSize++) {
- for (UInt_t isector = 0; isector < 36; isector++) {
- if (padSize == 0) sprintf(chname, "driftlengt vs. charge, sector %i, short pads", isector);
- if (padSize == 1) sprintf(chname, "driftlengt vs. charge, sector %i, medium pads", isector);
- if (padSize == 2) sprintf(chname, "driftlengt vs. charge, sector %i, long pads", isector);
- tempProf = new TProfile(chname, chname, 500, 0, 250);
- tempProf->SetYTitle("Charge");
- tempProf->SetXTitle("Driftlength");
- fcalPadRegionChargeVsDriftlength->SetObject(tempProf, isector, padSize);
- }
- }
if (GetDebugLevel() > 1) cout << "AliTPCcalibTracks object sucessfully constructed: " << GetName() << endl;
if (GetDebugLevel() > 0) cout << "AliTPCcalibTracks' destuctor called." << endl;
Int_t length = 0;
- if (fArrayAmpRow) length = fArrayAmpRow->GetEntriesFast();
- for (Int_t i = 0; i < length; i++){
- delete fArrayAmpRow->At(i);
- delete fArrayAmp->At(i);
- }
- delete fArrayAmpRow;
- delete fArrayAmp;
- delete fDeltaY;
- delete fDeltaZ;
if (fResolY) {
length = fResolY->GetEntriesFast();
}
}
- if (fArrayChargeVsDriftlength) {
- length = fArrayChargeVsDriftlength->GetEntriesFast();
- for (Int_t i = 0; i < length; i++){
- delete fArrayChargeVsDriftlength->At(i);
- }
- }
delete fArrayQDY;
delete fArrayQDZ;
delete fArrayQRMSY;
delete fArrayQRMSZ;
- delete fArrayChargeVsDriftlength;
- delete fHclus;
delete fRejectedTracksHisto;
delete fClusterCutHisto;
- delete fHclusterPerPadrow;
- delete fHclusterPerPadrowRaw;
if (fCalPadClusterPerPad) delete fCalPadClusterPerPad;
if (fCalPadClusterPerPadRaw) delete fCalPadClusterPerPadRaw;
- if(fcalPadRegionChargeVsDriftlength) {
- fcalPadRegionChargeVsDriftlength->Delete();
- delete fcalPadRegionChargeVsDriftlength;
- }
delete fHisDeltaY; // THnSparse - delta Y
delete fHisDeltaZ; // THnSparse - delta Z
delete fHisRMSY; // THnSparse - rms Y
// if this chi2 is bigger than a given threshold, assume that the current cluster is
// a kink an goto next padrow
// if not:
- // fill fArrayAmpRow, array with amplitudes vs. row for given sector
- // fill fArrayAmp, array with amplitude histograms for give sector
- // fill fRMSY, fRMSZ, fArrayQRMSY and fArrayQRMSZ, fDeltaY, fDeltaZ, fResolY, fResolZ, fArrayQDY, fArrayQDY
+ // fill fRMSY, fRMSZ, fArrayQRMSY and fArrayQRMSZ, fResolY, fResolZ, fArrayQDY, fArrayQDY
//
// write debug information to 'TPCSelectorDebug.root'
// only for every kDeltaWriteDebugStream'th padrow to reduce data volume
// and to avoid redundant data
//
- static TLinearFitter fFitterLinY1(2,"pol1"); //
- static TLinearFitter fFitterLinZ1(2,"pol1"); //
- static TLinearFitter fFitterLinY2(2,"pol1"); //
- static TLinearFitter fFitterLinZ2(2,"pol1"); //
static TLinearFitter fFitterParY(3,"pol2"); //
static TLinearFitter fFitterParZ(3,"pol2"); //
-
- fFitterLinY1.StoreData(kFALSE);
- fFitterLinZ1.StoreData(kFALSE);
- fFitterLinY2.StoreData(kFALSE);
- fFitterLinZ2.StoreData(kFALSE);
- fFitterParY.StoreData(kFALSE);
- fFitterParZ.StoreData(kFALSE);
-
-
+ static TLinearFitter fFitterParYWeight(3,"pol2"); //
+ static TLinearFitter fFitterParZWeight(3,"pol2"); //
+ fFitterParY.StoreData(kTRUE);
+ fFitterParZ.StoreData(kTRUE);
+ fFitterParYWeight.StoreData(kTRUE);
+ fFitterParZWeight.StoreData(kTRUE);
if (GetDebugLevel() > 5) Info("FillResolutionHistoLocal"," ***** Start of FillResolutionHistoLocal *****");
- const Int_t kDelta = 10; // delta rows to fit
- const Float_t kMinRatio = 0.75; // minimal ratio
- // const Float_t kCutChi2 = 6.; // cut chi2 - left right - kink removal
- const Float_t kErrorFraction = 0.5; // use only clusters with small interpolation error - for error param
- const Int_t kFirstLargePad = 127; // medium pads -> long pads
- const Float_t kLargePadSize = 1.5; // factor between medium and long pads' area
- const Int_t kDeltaWriteDebugStream = 5; // only for every kDeltaWriteDebugStream'th padrow debug information is calulated and written to debugstream
- TVectorD paramY0(2);
- TVectorD paramZ0(2);
- TVectorD paramY1(2);
- TVectorD paramZ1(2);
- TVectorD paramY2(3);
- TVectorD paramZ2(3);
- TMatrixD matrixY0(2,2);
- TMatrixD matrixZ0(2,2);
- TMatrixD matrixY1(2,2);
- TMatrixD matrixZ1(2,2);
-
- // estimate mean error
- Int_t nTrackletsAll = 0; // number of tracklets for given track
- Float_t csigmaY = 0; // mean sigma for tracklet refit in Y direction
- Float_t csigmaZ = 0; // mean sigma for tracklet refit in Z direction
- Int_t nClusters = 0; // working variable, number of clusters per tracklet
- Int_t sectorG = -1; // working variable, sector of tracklet, has to stay constant for one tracklet
-
- fHclus->Fill(track->GetNumberOfClusters()); // for statistics overview
- // ---------------------------------------------------------------------
- for (Int_t irow = 0; irow < 159; irow++){
- // loop over all rows along the track
- // fit tracklets (length: 13 rows) with pol2 in Y and Z direction
- // calculate mean chi^2 for this track-fit in Y and Z direction
- AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow);
- if (!cluster0) continue; // no cluster found
- Int_t sector = cluster0->GetDetector();
- fHclusterPerPadrowRaw->Fill(irow);
-
- Int_t ipad = TMath::Nint(cluster0->GetPad());
- Float_t value = fCalPadClusterPerPadRaw->GetCalROC(sector)->GetValue((sector<36)?irow:irow-64, TMath::Nint(cluster0->GetPad()));
- fCalPadClusterPerPadRaw->GetCalROC(sector)->SetValue((sector<36)?irow:irow-64, ipad, value + 1 );
-
- if (sector != sectorG){
- // track leaves sector before it crossed enough rows to fit / initialization
- nClusters = 0;
- fFitterParY.ClearPoints();
- fFitterParZ.ClearPoints();
- sectorG = sector;
- }
- else {
- nClusters++;
- Double_t x = cluster0->GetX();
- fFitterParY.AddPoint(&x, cluster0->GetY(), 1);
- fFitterParZ.AddPoint(&x, cluster0->GetZ(), 1);
- //
- if ( nClusters >= kDelta + 3 ){
- // if more than 13 (kDelta+3) clusters were added to the fitters
- // fit the tracklet, increase trackletCounter
- fFitterParY.Eval();
- fFitterParZ.Eval();
- nTrackletsAll++;
- csigmaY += fFitterParY.GetChisquare() / (nClusters - 3.);
- csigmaZ += fFitterParZ.GetChisquare() / (nClusters - 3.);
- nClusters = -1;
- fFitterParY.ClearPoints();
- fFitterParZ.ClearPoints();
- }
- }
- } // for (Int_t irow = 0; irow < 159; irow++)
- // mean chi^2 for all tracklet fits in Y and in Z direction:
- csigmaY = TMath::Sqrt(TMath::Abs(csigmaY) / (nTrackletsAll+0.1));
- csigmaZ = TMath::Sqrt(TMath::Abs(csigmaZ) / (nTrackletsAll+0.1));
- // ---------------------------------------------------------------------
-
- for (Int_t irow = 0; irow < 159; irow++){
- // loop again over all rows along the track
- // do analysis
- //
- Int_t nclFound = 0; // number of clusters in the neighborhood
- Int_t ncl0 = 0; // number of clusters in rows < rowOfCenterCluster
- Int_t ncl1 = 0; // number of clusters in rows > rowOfCenterCluster
- AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow);
- if (!cluster0) continue;
- Int_t sector = cluster0->GetDetector();
- Float_t xref = cluster0->GetX();
-
- // Make Fit
+ const Int_t kDelta = 10; // delta rows to fit
+ const Float_t kMinRatio = 0.75; // minimal ratio
+ const Float_t kChi2Cut = 10.; // cut chi2 - left right
+ const Float_t kSigmaCut = 3.; //sigma cut
+ const Float_t kdEdxCut=300;
+ const Float_t kPtCut=0.300;
+
+ if (track->GetTPCsignal()>kdEdxCut) return;
+ if (TMath::Abs(AliTracker::GetBz())>0.1 &&TMath::Abs(track->Pt())<kPtCut) return;
+
+ // estimate mean error
+ Int_t nTrackletsAll = 0; // number of tracklets for given track
+ Float_t csigmaY = 0; // mean sigma for tracklet refit in Y direction
+ Float_t csigmaZ = 0; // mean sigma for tracklet refit in Z direction
+ Int_t nClusters = 0; // working variable, number of clusters per tracklet
+ Int_t sectorG = -1; // working variable, sector of tracklet, has to stay constant for one tracklet
+ Double_t refX=0;
+ // ---------------------------------------------------------------------
+ for (Int_t irow = 0; irow < 159; irow++){
+ // loop over all rows along the track
+ // fit tracklets (length: 13 rows) with pol2 in Y and Z direction
+ // calculate mean chi^2 for this track-fit in Y and Z direction
+ AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow);
+ if (!cluster0) continue; // no cluster found
+ Int_t sector = cluster0->GetDetector();
+
+ Int_t ipad = TMath::Nint(cluster0->GetPad());
+ Float_t value = fCalPadClusterPerPadRaw->GetCalROC(sector)->GetValue((sector<36)?irow:irow-64, TMath::Nint(cluster0->GetPad()));
+ fCalPadClusterPerPadRaw->GetCalROC(sector)->SetValue((sector<36)?irow:irow-64, ipad, value + 1 );
+
+ if (sector != sectorG){
+ // track leaves sector before it crossed enough rows to fit / initialization
+ nClusters = 0;
fFitterParY.ClearPoints();
fFitterParZ.ClearPoints();
- fFitterLinY1.ClearPoints();
- fFitterLinZ1.ClearPoints();
- fFitterLinY2.ClearPoints();
- fFitterLinZ2.ClearPoints();
-
- // fit tracklet (clusters in given padrow +- kDelta padrows)
- // with polynom of 2nd order and two polynoms of 1st order
- // take both polynoms of 1st order, calculate difference of their parameters
- // add covariance matrixes and calculate chi2 of this difference
- // if this chi2 is bigger than a given threshold, assume that the current cluster is
- // a kink an goto next padrow
-
- for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){
- // loop over irow +- kDelta rows (neighboured rows)
- //
- //
- if (idelta == 0) continue; // don't use center cluster
- if (idelta + irow < 0 || idelta + irow > 159) continue; // don't go out of ROC
- AliTPCclusterMI * currentCluster = track->GetClusterPointer(irow + idelta);
- if (!currentCluster) continue;
- if (currentCluster->GetType() < 0) continue;
- if (currentCluster->GetDetector() != sector) continue;
- Double_t x = currentCluster->GetX() - xref; // x = differece: current cluster - cluster @ irow
- nclFound++;
- if (idelta < 0){
- ncl0++;
- fFitterLinY1.AddPoint(&x, currentCluster->GetY(), csigmaY);
- fFitterLinZ1.AddPoint(&x, currentCluster->GetZ(), csigmaZ);
- }
- if (idelta > 0){
- ncl1++;
- fFitterLinY2.AddPoint(&x, currentCluster->GetY(), csigmaY);
- fFitterLinZ2.AddPoint(&x, currentCluster->GetZ(), csigmaZ);
- }
- fFitterParY.AddPoint(&x, currentCluster->GetY(), csigmaY);
- fFitterParZ.AddPoint(&x, currentCluster->GetZ(), csigmaZ);
- } // loop over neighbourhood for fitter filling
-
-
-
- if (nclFound < kDelta * kMinRatio) fRejectedTracksHisto->Fill(10);
- if (nclFound < kDelta * kMinRatio) fClusterCutHisto->Fill(1, irow);
- if (nclFound < kDelta * kMinRatio) continue; // if not enough clusters (7.5) found in neighbourhood goto next padrow
- fFitterParY.Eval();
- fFitterParZ.Eval();
- Double_t chi2 = (fFitterParY.GetChisquare() + fFitterParZ.GetChisquare()) / (2. * nclFound - 6.);
- TTreeSRedirector *cstream = GetDebugStreamer();
- if (cstream){
- (*cstream)<<"Cut9"<<
- "chi2="<<chi2<<
- "\n";
- }
- // REMOVE KINK
- // only when there are enough clusters (4) in each direction
- if (ncl0 > 4){
- fFitterLinY1.Eval();
- fFitterLinZ1.Eval();
- }
- if (ncl1 > 4){
- fFitterLinY2.Eval();
- fFitterLinZ2.Eval();
- }
-
- if (ncl0 > 4 && ncl1 > 4){
- fFitterLinY1.GetCovarianceMatrix(matrixY0);
- fFitterLinY2.GetCovarianceMatrix(matrixY1);
- fFitterLinZ1.GetCovarianceMatrix(matrixZ0);
- fFitterLinZ2.GetCovarianceMatrix(matrixZ1);
- fFitterLinY2.GetParameters(paramY1);
- fFitterLinZ2.GetParameters(paramZ1);
- fFitterLinY1.GetParameters(paramY0);
- fFitterLinZ1.GetParameters(paramZ0);
- paramY0 -= paramY1;
- paramZ0 -= paramZ1;
- matrixY0 += matrixY1;
- matrixZ0 += matrixZ1;
- Double_t cchi2 = 0;
-
- TMatrixD difY(2, 1, paramY0.GetMatrixArray());
- TMatrixD difYT(1, 2, paramY0.GetMatrixArray());
- matrixY0.Invert();
- TMatrixD mulY(matrixY0, TMatrixD::kMult, difY);
- TMatrixD chi2Y(difYT, TMatrixD::kMult, mulY);
- cchi2 += chi2Y(0, 0);
-
- TMatrixD difZ(2, 1, paramZ0.GetMatrixArray());
- TMatrixD difZT(1, 2, paramZ0.GetMatrixArray());
- matrixZ0.Invert();
- TMatrixD mulZ(matrixZ0, TMatrixD::kMult, difZ);
- TMatrixD chi2Z(difZT, TMatrixD::kMult, mulZ);
- cchi2 += chi2Z(0, 0);
-
- if (cstream){
- (*cstream)<<"Cut8"<<
- "chi2="<<cchi2<<
- "\n";
- }
- }
-
- // current padrow has no kink
-
- // get fit parameters from pol2 fit:
- Double_t paramY[4], paramZ[4];
- paramY[0] = fFitterParY.GetParameter(0);
- paramY[1] = fFitterParY.GetParameter(1);
- paramY[2] = fFitterParY.GetParameter(2);
- paramZ[0] = fFitterParZ.GetParameter(0);
- paramZ[1] = fFitterParZ.GetParameter(1);
- paramZ[2] = fFitterParZ.GetParameter(2);
-
- Double_t tracky = paramY[0];
- Double_t trackz = paramZ[0];
- Float_t deltay = tracky - cluster0->GetY();
- Float_t deltaz = trackz - cluster0->GetZ();
- Float_t angley = paramY[1] - paramY[0] / xref;
- Float_t anglez = paramZ[1];
-
- Float_t max = cluster0->GetMax();
- UInt_t isegment = cluster0->GetDetector() % 36;
- Int_t padSize = 0; // short pads
- if (cluster0->GetDetector() >= 36) {
- padSize = 1; // medium pads
- if (cluster0->GetRow() > 63) padSize = 2; // long pads
+ sectorG = sector;
+ }
+ else {
+ nClusters++;
+ if (refX<1) refX=cluster0->GetX()+kDelta*0.5;
+ Double_t x = cluster0->GetX()-refX;
+ fFitterParY.AddPoint(&x, cluster0->GetY(), 1);
+ fFitterParZ.AddPoint(&x, cluster0->GetZ(), 1);
+ //
+ if ( nClusters >= kDelta + 3 ){
+ // if more than 13 (kDelta+3) clusters were added to the fitters
+ // fit the tracklet, increase trackletCounter
+ fFitterParY.Eval();
+ fFitterParZ.Eval();
+ nTrackletsAll++;
+ csigmaY += fFitterParY.GetChisquare() / (nClusters - 3.);
+ csigmaZ += fFitterParZ.GetChisquare() / (nClusters - 3.);
+ nClusters = -1;
+ fFitterParY.ClearPoints();
+ fFitterParZ.ClearPoints();
+ refX=0;
}
+ }
+ } // for (Int_t irow = 0; irow < 159; irow++)
+ // mean chi^2 for all tracklet fits in Y and in Z direction:
+ csigmaY = TMath::Sqrt(TMath::Abs(csigmaY) / (nTrackletsAll+0.1));
+ csigmaZ = TMath::Sqrt(TMath::Abs(csigmaZ) / (nTrackletsAll+0.1));
+ // ---------------------------------------------------------------------
+ //
+ //
- // =========================================
- // wirte collected information to histograms
- // =========================================
-
- TProfile *profAmpRow = (TProfile*)fArrayAmpRow->At(sector);
- if ( irow >= kFirstLargePad) max /= kLargePadSize;
- Double_t smax = TMath::Sqrt(max);
- profAmpRow->Fill( (Double_t)cluster0->GetRow(), smax );
- TH1F *hisAmp = (TH1F*)fArrayAmp->At(sector);
- hisAmp->Fill(smax);
-
- // remove the following two lines one day:
- TProfile *profDriftLength = (TProfile*)fArrayChargeVsDriftlength->At(sector);
- profDriftLength->Fill( 250.-(Double_t)TMath::Abs(cluster0->GetZ()), smax );
-
- TProfile *profDriftLengthTmp = (TProfile*)(fcalPadRegionChargeVsDriftlength->GetObject(isegment, padSize));
- profDriftLengthTmp->Fill( 250.-(Double_t)TMath::Abs(cluster0->GetZ()), smax );
-
- fHclusterPerPadrow->Fill(irow); // fill histogram showing clusters per padrow
- Int_t ipad = TMath::Nint(cluster0->GetPad());
- Float_t value = fCalPadClusterPerPad->GetCalROC(sector)->GetValue((sector<36)?irow:irow-64, TMath::Nint(cluster0->GetPad()));
- fCalPadClusterPerPad->GetCalROC(sector)->SetValue((sector<36)?irow:irow-64, ipad, value + 1 );
-
-
- TH3F * his3 = 0;
- his3 = (TH3F*)fRMSY->At(padSize);
- if (his3) his3->Fill(250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(cluster0->GetSigmaY2()) );
- his3 = (TH3F*)fRMSZ->At(padSize);
- if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(cluster0->GetSigmaZ2()) );
-
- his3 = (TH3F*)fArrayQRMSY->At(GetBin(cluster0->GetMax(), padSize));
- if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(cluster0->GetSigmaY2()) );
- his3 = (TH3F*)fArrayQRMSZ->At(GetBin(cluster0->GetMax(), padSize));
- if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(cluster0->GetSigmaZ2()) );
-
-
- // Fill resolution histograms
- Bool_t useForResol = kTRUE;
- if (fFitterParY.GetParError(0) > kErrorFraction * csigmaY) useForResol = kFALSE;
-
- if (cstream){
- Float_t zdrift = 250 - TMath::Abs(cluster0->GetZ());
- Float_t sy = cluster0->GetSigmaY2();
- Float_t sz = cluster0->GetSigmaZ2();
- (*cstream)<<"Resol0"<<
- "run="<<fRun<< // run number
- "event="<<fEvent<< // event number
- "time="<<fTime<< // time stamp of event
- "trigger="<<fTrigger<< // trigger
- "mag="<<fMagF<< // magnetic field
- "padSize="<<padSize<<
- "angley="<<angley<<
- "anglez="<<anglez<<
- "zdr="<<zdrift<<
- "dy="<<deltay<<
- "dz="<<deltaz<<
- "sy="<<sy<<
- "sz="<<sz<<
- "\n";
+ for (Int_t irow = kDelta; irow < 159-kDelta; irow++){
+ // loop again over all rows along the track
+ // do analysis
+ //
+ Int_t nclFound = 0; // number of clusters in the neighborhood
+ Int_t ncl0 = 0; // number of clusters in rows < rowOfCenterCluster
+ Int_t ncl1 = 0; // number of clusters in rows > rowOfCenterCluster
+ AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow);
+ if (!cluster0) continue;
+ Int_t sector = cluster0->GetDetector();
+ Float_t xref = cluster0->GetX();
+
+ // Make Fit
+ fFitterParY.ClearPoints();
+ fFitterParZ.ClearPoints();
+ fFitterParYWeight.ClearPoints();
+ fFitterParZWeight.ClearPoints();
+ // fit tracklet (clusters in given padrow +- kDelta padrows)
+ // with polynom of 2nd order and two polynoms of 1st order
+ // take both polynoms of 1st order, calculate difference of their parameters
+ // add covariance matrixes and calculate chi2 of this difference
+ // if this chi2 is bigger than a given threshold, assume that the current cluster is
+ // a kink an goto next padrow
+ AliRieman riemanFit(2*kDelta);
+ AliRieman riemanFitW(2*kDelta);
+ for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){
+ // loop over irow +- kDelta rows (neighboured rows)
+ //
+ //
+ if (idelta == 0) continue; // don't use center cluster
+ if (idelta + irow < 0 || idelta + irow > 159) continue; // don't go out of ROC
+ AliTPCclusterMI * currentCluster = track->GetClusterPointer(irow + idelta);
+ if (!currentCluster) continue;
+ if (currentCluster->GetType() < 0) continue;
+ if (currentCluster->GetDetector() != sector) continue;
+ nclFound++;
+ if (idelta < 0){
+ ncl0++;
}
-
- if (useForResol){
- fDeltaY->Fill(deltay);
- fDeltaZ->Fill(deltaz);
- his3 = (TH3F*)fResolY->At(padSize);
- if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), deltay );
- his3 = (TH3F*)fResolZ->At(padSize);
- if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), deltaz );
- his3 = (TH3F*)fArrayQDY->At(GetBin(cluster0->GetMax(), padSize));
- if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(angley), deltay );
- his3 = (TH3F*)fArrayQDZ->At(GetBin(cluster0->GetMax(), padSize));
- if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(anglez), deltaz );
+ if (idelta > 0){
+ ncl1++;
}
-
- //=============================================================================================
-
- if (useForResol && nclFound > 2 * kMinRatio * kDelta
- && irow % kDeltaWriteDebugStream == 0 && GetDebugLevel() > 4){
- if (GetDebugLevel() > 20) Info("FillResolutionHistoLocal","Filling 'TPCSelectorDebug.root', irow = %i", irow);
- FillResolutionHistoLocalDebugPart(track, cluster0, irow, angley, anglez, nclFound, kDelta);
- } // if (useForResol && nclFound > 2 * kMinRatio * kDelta)
- //
- // Fill THN histograms
- //
- Double_t xvar[9];
- xvar[1]=padSize;
- xvar[2]=1.-TMath::Abs(cluster0->GetZ())/250.;
- xvar[3]=cluster0->GetMax();
- xvar[5]=angley;
- xvar[6]=anglez;
-
- xvar[4]=cluster0->GetPad()-TMath::Nint(cluster0->GetPad());
- xvar[0]=deltay;
- fHisDeltaY->Fill(xvar);
- xvar[0]=TMath::Sqrt(cluster0->GetSigmaY2());
- fHisRMSY->Fill(xvar);
-
- xvar[4]=cluster0->GetTimeBin()-TMath::Nint(cluster0->GetTimeBin());
- xvar[0]=deltaz;
- fHisDeltaZ->Fill(xvar);
- xvar[0]=TMath::Sqrt(cluster0->GetSigmaZ2());
- fHisRMSZ->Fill(xvar);
-
- } // loop over all padrows along the track: for (Int_t irow = 0; irow < 159; irow++)
-} // FillResolutionHistoLocal(...)
-
-
-
-void AliTPCcalibTracks::FillResolutionHistoLocalDebugPart(AliTPCseed *track, AliTPCclusterMI *cluster0, Int_t irow, Float_t angley, Float_t anglez, Int_t nclFound, Int_t kDelta) {
- //
- // - debug part of FillResolutionHistoLocal -
- // called only for every kDeltaWriteDebugStream'th padrow, to avoid to much redundant data
- // called only for GetStreamLevel() > 4
- // fill resolution trees
- //
-
- Int_t sector = cluster0->GetDetector();
- Float_t xref = cluster0->GetX();
- Int_t padSize = 0; // short pads
- if (cluster0->GetDetector() >= 36) {
+ riemanFit.AddPoint(currentCluster->GetX(), currentCluster->GetY(),currentCluster->GetZ(), csigmaY,csigmaZ);
+ riemanFitW.AddPoint(currentCluster->GetX(), currentCluster->GetY(),currentCluster->GetZ(), csigmaY*TMath::Sqrt(1+TMath::Abs(idelta)),csigmaZ*TMath::Sqrt(1+TMath::Abs(idelta)));
+ } // loop over neighbourhood for fitter filling
+ if (nclFound < kDelta * kMinRatio) fRejectedTracksHisto->Fill(10);
+ if (nclFound < kDelta * kMinRatio) fClusterCutHisto->Fill(1, irow);
+ if (nclFound < kDelta * kMinRatio) continue; // if not enough clusters (7.5) found in neighbourhood goto next padrow
+ riemanFit.Update();
+ riemanFitW.Update();
+ Double_t chi2R=TMath::Sqrt(riemanFit.CalcChi2()/(2*nclFound-5));
+ Double_t chi2RW=TMath::Sqrt(riemanFitW.CalcChi2()/(2*nclFound-5));
+ Double_t paramYR[3], paramZR[3];
+ Double_t paramYRW[3], paramZRW[3];
+ //
+ paramYR[0] = riemanFit.GetYat(xref);
+ paramZR[0] = riemanFit.GetZat(xref);
+ paramYRW[0] = riemanFitW.GetYat(xref);
+ paramZRW[0] = riemanFitW.GetZat(xref);
+ //
+ paramYR[1] = riemanFit.GetDYat(xref);
+ paramZR[1] = riemanFit.GetDZat(xref);
+ paramYRW[1] = riemanFitW.GetDYat(xref);
+ paramZRW[1] = riemanFitW.GetDZat(xref);
+ //
+ Int_t reject=0;
+ if (chi2R>kChi2Cut) reject+=1;
+ if (chi2RW>kChi2Cut) reject+=2;
+ if (TMath::Abs(paramYR[0]-paramYRW[0])>kSigmaCut*csigmaY) reject+=4;
+ if (TMath::Abs(paramZR[0]-paramZRW[0])>kSigmaCut*csigmaZ) reject+=8;
+ if (TMath::Abs(paramYR[1]-paramYRW[1])>csigmaY) reject+=16;
+ if (TMath::Abs(paramZR[1]-paramZRW[1])>csigmaZ) reject+=32;
+ //
+ TTreeSRedirector *cstream = GetDebugStreamer();
+ // get fit parameters from pol2 fit:
+ Double_t tracky = paramYR[0];
+ Double_t trackz = paramZR[0];
+ Float_t deltay = cluster0->GetY()-tracky;
+ Float_t deltaz = cluster0->GetZ()-trackz;
+ Float_t angley = paramYR[1];
+ Float_t anglez = paramZR[1];
+
+ Int_t padSize = 0; // short pads
+ if (cluster0->GetDetector() >= 36) {
padSize = 1; // medium pads
if (cluster0->GetRow() > 63) padSize = 2; // long pads
- }
-
- static TLinearFitter fitY0(3, "pol2");
- static TLinearFitter fitZ0(3, "pol2");
- static TLinearFitter fitY2(5, "hyp4");
- static TLinearFitter fitZ2(5, "hyp4");
- static TLinearFitter fitY2Q(5, "hyp4");
- static TLinearFitter fitZ2Q(5, "hyp4");
- static TLinearFitter fitY2S(5, "hyp4");
- static TLinearFitter fitZ2S(5, "hyp4");
- fitY0.ClearPoints();
- fitZ0.ClearPoints();
- fitY2.ClearPoints();
- fitZ2.ClearPoints();
- fitY2Q.ClearPoints();
- fitZ2Q.ClearPoints();
- fitY2S.ClearPoints();
- fitZ2S.ClearPoints();
-
- for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){
- // loop over irow +- kDelta rows (neighboured rows)
- //
- //
- if (idelta == 0) continue;
- if (idelta + irow < 0 || idelta + irow > 159) continue; // don't go out of ROC
- AliTPCclusterMI * cluster = track->GetClusterPointer(irow + idelta);
- if (!cluster) continue;
- if (cluster->GetType() < 0) continue;
- if (cluster->GetDetector() != sector) continue;
- Double_t x = cluster->GetX() - xref;
- Double_t sigmaY0 = fClusterParam->GetError0Par( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(angley) );
- Double_t sigmaZ0 = fClusterParam->GetError0Par( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(anglez) );
- //
- Double_t sigmaYQ = fClusterParam->GetErrorQPar( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(angley), TMath::Abs(cluster->GetMax()) );
- Double_t sigmaZQ = fClusterParam->GetErrorQPar( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(anglez), TMath::Abs(cluster->GetMax()) );
- Double_t sigmaYS = fClusterParam->GetErrorQParScaled( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())),
- TMath::Abs(angley), TMath::Abs(cluster->GetMax()) );
- Double_t sigmaZS = fClusterParam->GetErrorQParScaled( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())),
- TMath::Abs(anglez), TMath::Abs(cluster->GetMax()) );
- Float_t rmsYFactor = fClusterParam->GetShapeFactor( 0, padSize,(250.0 - TMath::Abs(cluster->GetZ())),
- TMath::Abs(anglez), TMath::Abs(cluster->GetMax()),
- TMath::Sqrt(cluster0->GetSigmaY2()), 0 );
- Float_t rmsZFactor = fClusterParam->GetShapeFactor(0, padSize,(250.0 - TMath::Abs(cluster->GetZ())),
- TMath::Abs(anglez), TMath::Abs(cluster->GetMax()),
- TMath::Sqrt(cluster0->GetSigmaZ2()),0 );
- sigmaYS = TMath::Sqrt(sigmaYS * sigmaYS + rmsYFactor * rmsYFactor / 12.);
- sigmaZS = TMath::Sqrt(sigmaZS * sigmaZS + rmsZFactor * rmsZFactor / 12. + rmsYFactor * rmsYFactor / 24.);
- //
- if (kDelta != 0){
- fitY0.AddPoint(&x, cluster->GetY(), sigmaY0);
- fitZ0.AddPoint(&x, cluster->GetZ(), sigmaZ0);
- }
- Double_t xxx[4];
- xxx[0] = ( (idelta+irow) % 2 == 0 ) ? 1 : 0;
- xxx[1] = x;
- xxx[2] = ( (idelta+irow) % 2 == 0 ) ? x : 0;
- xxx[3] = x * x;
- fitY2.AddPoint(xxx, cluster->GetY(), sigmaY0);
- fitY2Q.AddPoint(xxx, cluster->GetY(), sigmaYQ);
- fitY2S.AddPoint(xxx, cluster->GetY(), sigmaYS);
- fitZ2.AddPoint(xxx, cluster->GetZ(), sigmaZ0);
- fitZ2Q.AddPoint(xxx, cluster->GetZ(), sigmaZQ);
- fitZ2S.AddPoint(xxx, cluster->GetZ(), sigmaZS);
- //
- } // neigbouhood-loop
- //
- fitY0.Eval();
- fitZ0.Eval();
- fitY2.Eval();
- fitZ2.Eval();
- fitY2Q.Eval();
- fitZ2Q.Eval();
- fitY2S.Eval();
- fitZ2S.Eval();
- Float_t chi2Y0 = fitY0.GetChisquare() / (nclFound-3.);
- Float_t chi2Z0 = fitZ0.GetChisquare() / (nclFound-3.);
- Float_t chi2Y2 = fitY2.GetChisquare() / (nclFound-5.);
- Float_t chi2Z2 = fitZ2.GetChisquare() / (nclFound-5.);
- Float_t chi2Y2Q = fitY2Q.GetChisquare() / (nclFound-5.);
- Float_t chi2Z2Q = fitZ2Q.GetChisquare() / (nclFound-5.);
- Float_t chi2Y2S = fitY2S.GetChisquare() / (nclFound-5.);
- Float_t chi2Z2S = fitZ2S.GetChisquare() / (nclFound-5.);
- //
- static TVectorD parY0(3);
- static TMatrixD matY0(3, 3);
- static TVectorD parZ0(3);
- static TMatrixD matZ0(3, 3);
- fitY0.GetParameters(parY0);
- fitY0.GetCovarianceMatrix(matY0);
- fitZ0.GetParameters(parZ0);
- fitZ0.GetCovarianceMatrix(matZ0);
- //
- static TVectorD parY2(5);
- static TMatrixD matY2(5,5);
- static TVectorD parZ2(5);
- static TMatrixD matZ2(5,5);
- fitY2.GetParameters(parY2);
- fitY2.GetCovarianceMatrix(matY2);
- fitZ2.GetParameters(parZ2);
- fitZ2.GetCovarianceMatrix(matZ2);
- //
- static TVectorD parY2Q(5);
- static TMatrixD matY2Q(5,5);
- static TVectorD parZ2Q(5);
- static TMatrixD matZ2Q(5,5);
- fitY2Q.GetParameters(parY2Q);
- fitY2Q.GetCovarianceMatrix(matY2Q);
- fitZ2Q.GetParameters(parZ2Q);
- fitZ2Q.GetCovarianceMatrix(matZ2Q);
- static TVectorD parY2S(5);
- static TMatrixD matY2S(5,5);
- static TVectorD parZ2S(5);
- static TMatrixD matZ2S(5,5);
- fitY2S.GetParameters(parY2S);
- fitY2S.GetCovarianceMatrix(matY2S);
- fitZ2S.GetParameters(parZ2S);
- fitZ2S.GetCovarianceMatrix(matZ2S);
- Float_t sigmaY0 = TMath::Sqrt(matY0(0,0));
- Float_t sigmaZ0 = TMath::Sqrt(matZ0(0,0));
- Float_t sigmaDY0 = TMath::Sqrt(matY0(1,1));
- Float_t sigmaDZ0 = TMath::Sqrt(matZ0(1,1));
- Float_t sigmaY2 = TMath::Sqrt(matY2(1,1));
- Float_t sigmaZ2 = TMath::Sqrt(matZ2(1,1));
- Float_t sigmaDY2 = TMath::Sqrt(matY2(3,3));
- Float_t sigmaDZ2 = TMath::Sqrt(matZ2(3,3));
- Float_t sigmaY2Q = TMath::Sqrt(matY2Q(1,1));
- Float_t sigmaZ2Q = TMath::Sqrt(matZ2Q(1,1));
- Float_t sigmaDY2Q = TMath::Sqrt(matY2Q(3,3));
- Float_t sigmaDZ2Q = TMath::Sqrt(matZ2Q(3,3));
- Float_t sigmaY2S = TMath::Sqrt(matY2S(1,1));
- Float_t sigmaZ2S = TMath::Sqrt(matZ2S(1,1));
- Float_t sigmaDY2S = TMath::Sqrt(matY2S(3,3));
- Float_t sigmaDZ2S = TMath::Sqrt(matZ2S(3,3));
-
- // Error parameters
- Float_t csigmaY0 = fClusterParam->GetError0Par(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),TMath::Abs(angley));
- Float_t csigmaZ0 = fClusterParam->GetError0Par(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),TMath::Abs(anglez));
- //
- Float_t csigmaYQ = fClusterParam->GetErrorQPar(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
- TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
- Float_t csigmaZQ = fClusterParam->GetErrorQPar(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
- TMath::Abs(anglez),TMath::Abs(cluster0->GetMax()));
- Float_t csigmaYS = fClusterParam->GetErrorQParScaled(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
- TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
- Float_t csigmaZS = fClusterParam->GetErrorQParScaled(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
- TMath::Abs(anglez),TMath::Abs(cluster0->GetMax()));
-
-
- // RMS parameters
- Float_t meanRMSY = 0;
- Float_t meanRMSZ = 0;
- Int_t nclRMS = 0;
- for (Int_t idelta = -2; idelta <= 2; idelta++){
- // loop over neighbourhood
- if (idelta+irow < 0 || idelta+irow > 159) continue;
- // if (idelta+irow>159) continue;
- AliTPCclusterMI * cluster = track->GetClusterPointer(irow+idelta);
- if (!cluster) continue;
- meanRMSY += TMath::Sqrt(cluster->GetSigmaY2());
- meanRMSZ += TMath::Sqrt(cluster->GetSigmaZ2());
- nclRMS++;
- }
- meanRMSY /= nclRMS;
- meanRMSZ /= nclRMS;
-
- Float_t rmsY = TMath::Sqrt(cluster0->GetSigmaY2());
- Float_t rmsZ = TMath::Sqrt(cluster0->GetSigmaZ2());
- Float_t rmsYT = fClusterParam->GetRMSQ(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
- TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
- Float_t rmsZT = fClusterParam->GetRMSQ(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
- TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
- Float_t rmsYT0 = fClusterParam->GetRMS0(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
- TMath::Abs(angley));
- Float_t rmsZT0 = fClusterParam->GetRMS0(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
- TMath::Abs(anglez));
- Float_t rmsYSigma = fClusterParam->GetRMSSigma(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
- TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
- Float_t rmsZSigma = fClusterParam->GetRMSSigma(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
- TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
- Float_t rmsYFactor = fClusterParam->GetShapeFactor(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
- TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()),
- rmsY,meanRMSY);
- Float_t rmsZFactor = fClusterParam->GetShapeFactor(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
- TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()),
- rmsZ,meanRMSZ);
- //
- // cluster debug
- TTreeSRedirector *cstream = GetDebugStreamer();
- if (cstream){
- (*cstream)<<"ResolCl"<< // valgrind 3 40,000 bytes in 1 blocks are possibly
- "run="<<fRun<< // run number
- "event="<<fEvent<< // event number
- "time="<<fTime<< // time stamp of event
- "trigger="<<fTrigger<< // trigger
- "mag="<<fMagF<< // magnetic field
- "Sector="<<sector<<
- "Cl.="<<cluster0<<
- "CSigmaY0="<<csigmaY0<< // cluster errorY
- "CSigmaYQ="<<csigmaYQ<< // cluster errorY - q scaled
- "CSigmaYS="<<csigmaYS<< // cluster errorY - q scaled
- "CSigmaZ0="<<csigmaZ0<< //
- "CSigmaZQ="<<csigmaZQ<<
- "CSigmaZS="<<csigmaZS<<
- "shapeYF="<<rmsYFactor<<
- "shapeZF="<<rmsZFactor<<
- "rmsY="<<rmsY<<
- "rmsZ="<<rmsZ<<
- "rmsYM="<<meanRMSY<<
- "rmsZM="<<meanRMSZ<<
- "rmsYT="<<rmsYT<<
- "rmsZT="<<rmsZT<<
- "rmsYT0="<<rmsYT0<<
- "rmsZT0="<<rmsZT0<<
- "rmsYS="<<rmsYSigma<<
- "rmsZS="<<rmsZSigma<<
- "padSize="<<padSize<<
- "Ncl="<<nclFound<<
- "PY0.="<<&parY0<<
- "PZ0.="<<&parZ0<<
- "SigmaY0="<<sigmaY0<<
- "SigmaZ0="<<sigmaZ0<<
- "angley="<<angley<<
- "anglez="<<anglez<<
- "\n";
-
- // tracklet dubug
- (*cstream)<<"ResolTr"<<
- "run="<<fRun<< // run number
- "event="<<fEvent<< // event number
- "time="<<fTime<< // time stamp of event
- "trigger="<<fTrigger<< // trigger
- "mag="<<fMagF<< // magnetic field
- "padSize="<<padSize<<
- "IPad="<<padSize<<
- "Sector="<<sector<<
- "Ncl="<<nclFound<<
- "chi2Y0="<<chi2Y0<<
- "chi2Z0="<<chi2Z0<<
- "chi2Y2="<<chi2Y2<<
- "chi2Z2="<<chi2Z2<<
- "chi2Y2Q="<<chi2Y2Q<<
- "chi2Z2Q="<<chi2Z2Q<<
- "chi2Y2S="<<chi2Y2S<<
- "chi2Z2S="<<chi2Z2S<<
- "PY0.="<<&parY0<<
- "PZ0.="<<&parZ0<<
- "PY2.="<<&parY2<<
- "PZ2.="<<&parZ2<<
- "PY2Q.="<<&parY2Q<<
- "PZ2Q.="<<&parZ2Q<<
- "PY2S.="<<&parY2S<<
- "PZ2S.="<<&parZ2S<<
- "SigmaY0="<<sigmaY0<<
- "SigmaZ0="<<sigmaZ0<<
- "SigmaDY0="<<sigmaDY0<<
- "SigmaDZ0="<<sigmaDZ0<<
- "SigmaY2="<<sigmaY2<<
- "SigmaZ2="<<sigmaZ2<<
- "SigmaDY2="<<sigmaDY2<<
- "SigmaDZ2="<<sigmaDZ2<<
- "SigmaY2Q="<<sigmaY2Q<<
- "SigmaZ2Q="<<sigmaZ2Q<<
- "SigmaDY2Q="<<sigmaDY2Q<<
- "SigmaDZ2Q="<<sigmaDZ2Q<<
- "SigmaY2S="<<sigmaY2S<<
- "SigmaZ2S="<<sigmaZ2S<<
- "SigmaDY2S="<<sigmaDY2S<<
- "SigmaDZ2S="<<sigmaDZ2S<<
- "angley="<<angley<<
- "anglez="<<anglez<<
- "\n";
- }
-}
+ }
+ Int_t ipad = TMath::Nint(cluster0->GetPad());
+ if (cstream){
+ Float_t zdrift = 250 - TMath::Abs(cluster0->GetZ());
+ (*cstream)<<"Resol0"<<
+ "run="<<fRun<< // run number
+ "event="<<fEvent<< // event number
+ "time="<<fTime<< // time stamp of event
+ "trigger="<<fTrigger<< // trigger
+ "mag="<<fMagF<< // magnetic field
+ "padSize="<<padSize<<
+ //
+ "reject="<<reject<<
+ "cl.="<<cluster0<< // cluster info
+ "xref="<<xref<< // cluster reference X
+ //rieman fit
+ "yr="<<paramYR[0]<< // track position y - rieman fit
+ "zr="<<paramZR[0]<< // track position z - rieman fit
+ "yrW="<<paramYRW[0]<< // track position y - rieman fit - weight
+ "zrW="<<paramZRW[0]<< // track position z - rieman fit - weight
+ "dyr="<<paramYR[1]<< // track position y - rieman fit
+ "dzr="<<paramZR[1]<< // track position z - rieman fit
+ "dyrW="<<paramYRW[1]<< // track position y - rieman fit - weight
+ "dzrW="<<paramZRW[1]<< // track position z - rieman fit - weight
+ //
+ "angley="<<angley<< // angle par fit
+ "anglez="<<anglez<< // angle par fit
+ "zdr="<<zdrift<< //
+ "dy="<<deltay<<
+ "dz="<<deltaz<<
+ "sy="<<csigmaY<<
+ "sz="<<csigmaZ<<
+ "chi2R="<<chi2R<<
+ "chi2RW="<<chi2RW<<
+ "\n";
+ }
+ if (reject) continue;
+
+
+ // =========================================
+ // wirte collected information to histograms
+ // =========================================
+
+ Float_t value = fCalPadClusterPerPad->GetCalROC(sector)->GetValue((sector<36)?irow:irow-64, TMath::Nint(cluster0->GetPad()));
+ fCalPadClusterPerPad->GetCalROC(sector)->SetValue((sector<36)?irow:irow-64, ipad, value + 1 );
+
+
+ TH3F * his3 = 0;
+ his3 = (TH3F*)fRMSY->At(padSize);
+ if (his3) his3->Fill(250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(cluster0->GetSigmaY2()) );
+ his3 = (TH3F*)fRMSZ->At(padSize);
+ if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(cluster0->GetSigmaZ2()) );
+
+ his3 = (TH3F*)fArrayQRMSY->At(GetBin(cluster0->GetMax(), padSize));
+ if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(cluster0->GetSigmaY2()) );
+ his3 = (TH3F*)fArrayQRMSZ->At(GetBin(cluster0->GetMax(), padSize));
+ if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(cluster0->GetSigmaZ2()) );
+
+
+ his3 = (TH3F*)fResolY->At(padSize);
+ if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), deltay );
+ his3 = (TH3F*)fResolZ->At(padSize);
+ if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), deltaz );
+ his3 = (TH3F*)fArrayQDY->At(GetBin(cluster0->GetMax(), padSize));
+ if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(angley), deltay );
+ his3 = (TH3F*)fArrayQDZ->At(GetBin(cluster0->GetMax(), padSize));
+ if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(anglez), deltaz );
+ //=============================================================================================
+ //
+ // Fill THN histograms
+ //
+ Double_t xvar[9];
+ xvar[1]=padSize;
+ xvar[2]=cluster0->GetZ();
+ xvar[3]=cluster0->GetMax();
+
+ xvar[0]=deltay;
+ xvar[4]=cluster0->GetPad()-Int_t(cluster0->GetPad());
+ xvar[5]=angley;
+ fHisDeltaY->Fill(xvar);
+ xvar[0]=TMath::Sqrt(cluster0->GetSigmaY2());
+ fHisRMSY->Fill(xvar);
+
+ xvar[0]=deltaz;
+ xvar[4]=cluster0->GetTimeBin()-Int_t(cluster0->GetTimeBin());
+ xvar[5]=anglez;
+ fHisDeltaZ->Fill(xvar);
+ xvar[0]=TMath::Sqrt(cluster0->GetSigmaZ2());
+ fHisRMSZ->Fill(xvar);
+
+ } // loop over all padrows along the track: for (Int_t irow = 0; irow < 159; irow++)
+} // FillResolutionHistoLocal(...)
+
-TH2D * AliTPCcalibTracks::MakeDiff(TH2D * hfit, TF2 * func){
- //
- // creates a new histogram which contains the difference between
- // the histogram hfit and the function func
- //
- TH2D * result = (TH2D*)hfit->Clone(); // valgrind 3 40,139 bytes in 11 blocks are still reachable
- result->SetTitle(Form("%s fit residuals",result->GetTitle()));
- result->SetName(Form("%s fit residuals",result->GetName()));
- TAxis *xaxis = hfit->GetXaxis();
- TAxis *yaxis = hfit->GetYaxis();
- Double_t x[2];
- for (Int_t biny = 0; biny <= yaxis->GetNbins(); biny++) {
- x[1] = yaxis->GetBinCenter(biny);
- for (Int_t binx = 0; binx <= xaxis->GetNbins(); binx++) {
- x[0] = xaxis->GetBinCenter(binx);
- Int_t bin = hfit->GetBin(binx, biny);
- Double_t val = hfit->GetBinContent(bin);
-// result->SetBinContent( bin, (val - func->Eval(x[0], x[1])) / func->Eval(x[0], x[1]) );
- result->SetBinContent( bin, (val / func->Eval(x[0], x[1])) - 1 );
- }
- }
- return result;
-}
void AliTPCcalibTracks::SetStyle() const {
}
-void AliTPCcalibTracks::Draw(Option_t* opt){
- //
- // draw-function of AliTPCcalibTracks
- // will draws some exemplaric pictures
- //
-
- if (GetDebugLevel() > 6) Info("Draw", "Drawing an exemplaric picture.");
- SetStyle();
- Double_t min = 0;
- Double_t max = 0;
- TCanvas *c1 = new TCanvas();
- c1->Divide(0, 3);
- TVirtualPad *upperThird = c1->GetPad(1);
- TVirtualPad *middleThird = c1->GetPad(2);
- TVirtualPad *lowerThird = c1->GetPad(3);
- upperThird->Divide(2,0);
- TVirtualPad *upleft = upperThird->GetPad(1);
- TVirtualPad *upright = upperThird->GetPad(2);
- middleThird->Divide(2,0);
- TVirtualPad *middleLeft = middleThird->GetPad(1);
- TVirtualPad *middleRight = middleThird->GetPad(2);
- lowerThird->Divide(2,0);
- TVirtualPad *downLeft = lowerThird->GetPad(1);
- TVirtualPad *downRight = lowerThird->GetPad(2);
-
-
- upleft->cd(0);
- min = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())-20;
- max = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())+20;
- fDeltaY->SetAxisRange(min, max);
- fDeltaY->Fit("gaus","q","",min, max); // valgrind 3 7 block possibly lost 2,400 bytes in 1 blocks are still reachable
- c1->Update();
-
- upright->cd(0);
- max = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())+20;
- min = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())-20;
- fDeltaZ->SetAxisRange(min, max);
- fDeltaZ->Fit("gaus","q","",min, max);
- c1->Update();
-
- middleLeft->cd();
- fHclus->Draw(opt);
-
- middleRight->cd();
- fRejectedTracksHisto->Draw(opt);
- TPaveText *pt = new TPaveText(0.6,0.6, 0.8,0.8, "NDC");
- TText *t1 = pt->AddText("1: kEdgeThetaCutNoise");
- TText *t2 = pt->AddText("2: kMinClusters");
- TText *t3 = pt->AddText("3: kMinRatio");
- TText *t4 = pt->AddText("4: kMax1pt");
- t1 = t1; t2 = t2; t3 = t3; t4 = t4; // avoid compiler warnings
- pt->SetToolTipText("Legend for failed cuts");
- pt->Draw();
-
- downLeft->cd();
- fHclusterPerPadrowRaw->Draw(opt);
-
- downRight->cd();
- fHclusterPerPadrow->Draw(opt);
-}
-
void AliTPCcalibTracks::MakeReport(Int_t stat, const char* pathName){
//
//
if (GetDebugLevel() > 0) Info("MakeReport","Writing plots and trees to '%s'.", pathName);
- MakeAmpPlots(stat, pathName);
- MakeDeltaPlots(pathName);
- FitResolutionNew(pathName);
- FitRMSNew(pathName);
- MakeChargeVsDriftLengthPlots(pathName);
-// MakeResPlotsQ(1, 1);
- MakeResPlotsQTree(stat, pathName);
+ MakeResPlotsQTree(stat, pathName);
}
-void AliTPCcalibTracks::MakeAmpPlots(Int_t stat, const char* pathName){
- //
- // creates several plots:
- // fArrayAmp.ps, fArrayAmpRow.ps and DeltaYZ.ps
- // fArrayAmp.ps: one histogram per sector, the histogram shows the charge per cluster
- // fArrayAmpRow.ps: one histogram per sector, mean max. amplitude vs. pad row with landau fit
- // DeltaYZ.ps: DeltaY and DeltaZ histogram with gaus fit
- // Empty histograms (sectors without data) are not written to file
- // the ps-files are written to the directory 'pathName', that is created if it does not exist
- // 'stat': only histograms with more than 'stat' entries are written to file.
- //
-
- SetStyle();
- gSystem->MakeDirectory(pathName);
- gSystem->ChangeDirectory(pathName);
-
- TCanvas* c1 = new TCanvas(); // valgrind 3 ??? 634 bytes in 28 blocks are still reachable
- TPostScript *ps;
- // histograms with accumulated amplitude for all IROCs and OROCs
- TH1F *allAmpHisIROC = ((TH1F*)(fArrayAmp->At(0))->Clone());
- allAmpHisIROC->SetName("Amp all IROCs");
- allAmpHisIROC->SetTitle("Amp all IROCs");
- TH1F *allAmpHisOROC = ((TH1F*)(fArrayAmp->At(36))->Clone());
- allAmpHisOROC->SetName("Amp all OROCs");
- allAmpHisOROC->SetTitle("Amp all OROCs");
-
- ps = new TPostScript("fArrayAmp.ps", 112);
- if (GetDebugLevel() > -1) cout << "creating fArrayAmp.ps..." << endl;
- for (Int_t i = 0; i < fArrayAmp->GetEntriesFast(); i++){
- if ( ((TH1F*)fArrayAmp->At(i))->GetEntries() < stat ) continue;
- ps->NewPage();
- ((TH1F*)fArrayAmp->At(i))->Draw();
- c1->Update(); // valgrind 3
- if (i > 0 && i < 36) {
- allAmpHisIROC->Add(((TH1F*)fArrayAmp->At(i)));
- allAmpHisOROC->Add(((TH1F*)fArrayAmp->At(i+36)));
- }
- }
- ps->NewPage();
- allAmpHisIROC->Draw();
- c1->Update(); // valgrind
- ps->NewPage();
- allAmpHisOROC->Draw();
- c1->Update();
- ps->Close();
- delete ps;
-
- TH1F *his = 0;
- Double_t min = 0;
- Double_t max = 0;
- ps = new TPostScript("fArrayAmpRow.ps", 112);
- if (GetDebugLevel() > -1) cout << "creating fArrayAmpRow.ps..." << endl;
- for (Int_t i = 0; i < fArrayAmpRow->GetEntriesFast(); i++){
- his = (TH1F*)fArrayAmpRow->At(i);
- if (his->GetEntries() < stat) continue;
- ps->NewPage();
- min = TMath::Max( his->GetBinCenter(his->GetMaximumBin() )-100., 0.);
- max = his->GetBinCenter(5*his->GetMaximumBin()) + 100;
- his->SetAxisRange(min, max);
- his->Fit("pol3", "q", "", min, max);
- // his->Draw("error"); // don't use this line when you don't want to have empty pages in the ps-file
- c1->Update();
- }
- ps->Close();
- delete ps;
- delete c1;
- gSystem->ChangeDirectory("..");
-}
-
-
-void AliTPCcalibTracks::MakeDeltaPlots(const char* pathName){
- //
- // creates several plots:
- // DeltaYZ.ps: DeltaY and DeltaZ histogram with gaus fit
- // the ps-files are written to the directory 'pathName', that is created if it does not exist
- //
-
- SetStyle();
- gSystem->MakeDirectory(pathName);
- gSystem->ChangeDirectory(pathName);
-
- TCanvas* c1 = new TCanvas(); // valgrind 3 ??? 634 bytes in 28 blocks are still reachable
- TPostScript *ps;
- Double_t min = 0;
- Double_t max = 0;
-
- ps = new TPostScript("DeltaYZ.ps", 112);
- if (GetDebugLevel() > -1) cout << "creating DeltaYZ.ps..." << endl;
- min = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())-20;
- max = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())+20;
- fDeltaY->SetAxisRange(min, max);
- ps->NewPage();
- fDeltaY->Fit("gaus","q","",min, max); // valgrind 3 7 block possibly lost 2,400 bytes in 1 blocks are still reachable
- c1->Update();
- ps->NewPage();
- max = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())+20;
- min = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())-20;
- fDeltaZ->SetAxisRange(min, max);
- fDeltaZ->Fit("gaus","q","",min, max);
- c1->Update();
- ps->Close();
- delete ps;
- delete c1;
- gSystem->ChangeDirectory("..");
-}
-
-
-void AliTPCcalibTracks::MakeChargeVsDriftLengthPlotsOld(const char* pathName){
- //
- // creates charge vs. driftlength plots, one TProfile for each ROC
- // is not correct like this, should be one TProfile for each sector and padsize
- //
-
- SetStyle();
- gSystem->MakeDirectory(pathName);
- gSystem->ChangeDirectory(pathName);
-
- TCanvas* c1 = new TCanvas(); // valgrind 3 ??? 634 bytes in 28 blocks are still reachable
- TPostScript *ps;
- ps = new TPostScript("chargeVsDriftlengthOld.ps", 112);
- if (GetDebugLevel() > -1) cout << "creating chargeVsDriftlength.ps..." << endl;
- TProfile *chargeVsDriftlengthAllIROCs = ((TProfile*)fArrayChargeVsDriftlength->At(0)->Clone());
- TProfile *chargeVsDriftlengthAllOROCs = ((TProfile*)fArrayChargeVsDriftlength->At(36)->Clone());
- chargeVsDriftlengthAllIROCs->SetName("allAmpHisIROC");
- chargeVsDriftlengthAllIROCs->SetTitle("charge vs. driftlength, all IROCs");
- chargeVsDriftlengthAllOROCs->SetName("allAmpHisOROC");
- chargeVsDriftlengthAllOROCs->SetTitle("charge vs. driftlength, all OROCs");
-
- for (Int_t i = 0; i < fArrayChargeVsDriftlength->GetEntriesFast(); i++) {
- ((TProfile*)fArrayChargeVsDriftlength->At(i))->Draw();
- c1->Update();
- if (i > 0 && i < 36) {
- chargeVsDriftlengthAllIROCs->Add(((TProfile*)fArrayChargeVsDriftlength->At(i)));
- chargeVsDriftlengthAllOROCs->Add(((TProfile*)fArrayChargeVsDriftlength->At(i+36)));
- }
- ps->NewPage();
- }
- chargeVsDriftlengthAllIROCs->Draw();
- c1->Update(); // valgrind
- ps->NewPage();
- chargeVsDriftlengthAllOROCs->Draw();
- c1->Update();
- ps->Close();
- delete ps;
- delete c1;
- gSystem->ChangeDirectory("..");
-}
-
-
-void AliTPCcalibTracks::MakeChargeVsDriftLengthPlots(const char* pathName){
- //
- // creates charge vs. driftlength plots, one TProfile for each ROC
- // under development....
- //
-
- SetStyle();
- gSystem->MakeDirectory(pathName);
- gSystem->ChangeDirectory(pathName);
-
- TCanvas* c1 = new TCanvas("c1", "c1", 700,(Int_t)(TMath::Sqrt(2)*700)); // valgrind 3 ??? 634 bytes in 28 blocks are still reachable
-// TCanvas c1("c1", "c1", 500,(sqrt(2)*500))
- c1->Divide(0,3);
- TPostScript *ps;
- ps = new TPostScript("chargeVsDriftlength.ps", 111);
- if (GetDebugLevel() > -1) cout << "creating chargeVsDriftlengthNew.ps..." << endl;
-
- TProfile *chargeVsDriftlengthAllShortPads = ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,0)->Clone());
- TProfile *chargeVsDriftlengthAllMediumPads = ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,1)->Clone());
- TProfile *chargeVsDriftlengthAllLongPads = ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,2)->Clone());
- chargeVsDriftlengthAllShortPads->SetName("allAmpHisShortPads");
- chargeVsDriftlengthAllShortPads->SetTitle("charge vs. driftlength, all sectors, short pads");
- chargeVsDriftlengthAllMediumPads->SetName("allAmpHisMediumPads");
- chargeVsDriftlengthAllMediumPads->SetTitle("charge vs. driftlength, all sectors, medium pads");
- chargeVsDriftlengthAllLongPads->SetName("allAmpHisLongPads");
- chargeVsDriftlengthAllLongPads->SetTitle("charge vs. driftlength, all sectors, long pads");
-
- for (Int_t i = 0; i < 36; i++) {
- c1->cd(1)->SetGridx();
- c1->cd(1)->SetGridy();
- ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(i,0))->Draw();
- c1->cd(2)->SetGridx();
- c1->cd(2)->SetGridy();
- ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(i,1))->Draw();
- c1->cd(3)->SetGridx();
- c1->cd(3)->SetGridy();
- ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(i,2))->Draw();
- c1->Update();
- chargeVsDriftlengthAllShortPads->Add( (TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,0));
- chargeVsDriftlengthAllMediumPads->Add((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,1));
- chargeVsDriftlengthAllLongPads->Add( (TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,2));
- ps->NewPage();
- }
- c1->cd(1)->SetGridx();
- c1->cd(1)->SetGridy();
- chargeVsDriftlengthAllShortPads->Draw();
- c1->cd(2)->SetGridx();
- c1->cd(2)->SetGridy();
- chargeVsDriftlengthAllMediumPads->Draw();
- c1->cd(3)->SetGridx();
- c1->cd(3)->SetGridy();
- chargeVsDriftlengthAllLongPads->Draw();
- c1->Update(); // valgrind
-// ps->NewPage();
- ps->Close();
- delete ps;
- delete c1;
- gSystem->ChangeDirectory("..");
-}
-
-
-
-void AliTPCcalibTracks::FitResolutionNew(const char* pathName){
- //
- // calculates different resulution fits in Y and Z direction
- // the histograms are written to 'ResolutionYZ.ps'
- // writes calculated resolution to 'resol.txt'
- // all files are stored in the directory pathName
- //
-
- SetStyle();
- gSystem->MakeDirectory(pathName);
- gSystem->ChangeDirectory(pathName);
-
- TCanvas c;
- c.Divide(2,1);
- if (GetDebugLevel() > -1) cout << "creating ResolutionYZ.ps..." << endl;
- TPostScript *ps = new TPostScript("ResolutionYZ.ps", 112);
- TF2 *fres = new TF2("fres","TMath::Sqrt([0]*[0]+[1]*[1]*x+[2]*[2]*y*y)",0,250,0,1);
- fres->SetParameter(0,0.02);
- fres->SetParameter(1,0.0054);
- fres->SetParameter(2,0.13);
-
- TH1::AddDirectory(kTRUE); // TH3F::FitSlicesZ() writes histograms into the current directory
-
- // create histogramw for Y-resolution
- TH3F * hisResY0 = (TH3F*)fResolY->At(0);
- hisResY0->FitSlicesZ();
- TH2D * hisResY02 = (TH2D*)gDirectory->Get("Resol Y0_2");
- TH3F * hisResY1 = (TH3F*)fResolY->At(1);
- hisResY1->FitSlicesZ();
- TH2D * hisResY12 = (TH2D*)gDirectory->Get("Resol Y1_2");
- TH3F * hisResY2 = (TH3F*)fResolY->At(2);
- hisResY2->FitSlicesZ();
- TH2D * hisResY22 = (TH2D*)gDirectory->Get("Resol Y2_2");
- //
- ps->NewPage();
- c.cd(1);
- hisResY02->Fit(fres, "q"); // valgrind 132,072 bytes in 6 blocks are indirectly lost
- hisResY02->Draw("surf1");
- c.cd(2);
- MakeDiff(hisResY02,fres)->Draw("surf1");
- c.Update();
- // c.SaveAs("ResolutionYPad0.eps");
- ps->NewPage();
- c.cd(1);
- hisResY12->Fit(fres, "q");
- hisResY12->Draw("surf1");
- c.cd(2);
- MakeDiff(hisResY12,fres)->Draw("surf1");
- c.Update();
- // c.SaveAs("ResolutionYPad1.eps");
- ps->NewPage();
- c.cd(1);
- hisResY22->Fit(fres, "q");
- hisResY22->Draw("surf1");
- c.cd(2);
- MakeDiff(hisResY22,fres)->Draw("surf1");
- c.Update();
-// c.SaveAs("ResolutionYPad2.eps");
-
- // create histogramw for Z-resolution
- TH3F * hisResZ0 = (TH3F*)fResolZ->At(0);
- hisResZ0->FitSlicesZ();
- TH2D * hisResZ02 = (TH2D*)gDirectory->Get("Resol Z0_2");
- TH3F * hisResZ1 = (TH3F*)fResolZ->At(1);
- hisResZ1->FitSlicesZ();
- TH2D * hisResZ12 = (TH2D*)gDirectory->Get("Resol Z1_2");
- TH3F * hisResZ2 = (TH3F*)fResolZ->At(2);
- hisResZ2->FitSlicesZ();
- TH2D * hisResZ22 = (TH2D*)gDirectory->Get("Resol Z2_2");
-
- ps->NewPage();
- c.cd(1);
- hisResZ02->Fit(fres, "q");
- hisResZ02->Draw("surf1");
- c.cd(2);
- MakeDiff(hisResZ02,fres)->Draw("surf1");
- c.Update();
-// c.SaveAs("ResolutionZPad0.eps");
- ps->NewPage();
- c.cd(1);
- hisResZ12->Fit(fres, "q");
- hisResZ12->Draw("surf1");
- c.cd(2);
- MakeDiff(hisResZ12,fres)->Draw("surf1");
- c.Update();
-// c.SaveAs("ResolutionZPad1.eps");
- ps->NewPage();
- c.cd(1);
- hisResZ22->Fit(fres, "q");
- hisResZ22->Draw("surf1");
- c.cd(2);
- MakeDiff(hisResZ22,fres)->Draw("surf1");
- c.Update();
-// c.SaveAs("ResolutionZPad2.eps");
- ps->Close();
- delete ps;
-
- // write calculated resoltuions to 'resol.txt'
- ofstream fresol("resol.txt");
- fresol<<"Pad 0.75 cm"<<"\n";
- hisResY02->Fit(fres, "q"); // valgrind
- fresol<<"Y\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
- hisResZ02->Fit(fres, "q");
- fresol<<"Z\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
- //
- fresol<<"Pad 1.00 cm"<<1<<"\n";
- hisResY12->Fit(fres, "q"); // valgrind
- fresol<<"Y\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
- hisResZ12->Fit(fres, "q");
- fresol<<"Z\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
- //
- fresol<<"Pad 1.50 cm"<<0<<"\n";
- hisResY22->Fit(fres, "q");
- fresol<<"Y\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
- hisResZ22->Fit(fres, "q");
- fresol<<"Z\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
-
- TH1::AddDirectory(kFALSE);
- gSystem->ChangeDirectory("..");
- delete fres;
-}
-
-
-void AliTPCcalibTracks::FitRMSNew(const char* pathName){
- //
- // calculates different resulution-rms fits in Y and Z direction
- // the histograms are written to 'RMS_YZ.ps'
- // writes calculated resolution rms to 'rms.txt'
- // all files are stored in the directory pathName
- //
-
- SetStyle();
- gSystem->MakeDirectory(pathName);
- gSystem->ChangeDirectory(pathName);
-
- TCanvas c; // valgrind 3 42,120 bytes in 405 blocks are still reachable 23,816 bytes in 229 blocks are still reachable
- c.Divide(2,1);
- if (GetDebugLevel() > -1) cout << "creating RMS_YZ.ps..." << endl;
- TPostScript *ps = new TPostScript("RMS_YZ.ps", 112);
- TF2 *frms = new TF2("fres","TMath::Sqrt([0]*[0]+[1]*[1]*x+[2]*[2]*y*y)",0,250,0,1);
- frms->SetParameter(0,0.02);
- frms->SetParameter(1,0.0054);
- frms->SetParameter(2,0.13);
-
- TH1::AddDirectory(kTRUE); // TH3F::FitSlicesZ() writes histograms into the current directory
-
- // create histogramw for Y-RMS
- TH3F * hisResY0 = (TH3F*)fRMSY->At(0);
- hisResY0->FitSlicesZ();
- TH2D * hisResY02 = (TH2D*)gDirectory->Get("RMS Y0_1");
- TH3F * hisResY1 = (TH3F*)fRMSY->At(1);
- hisResY1->FitSlicesZ();
- TH2D * hisResY12 = (TH2D*)gDirectory->Get("RMS Y1_1");
- TH3F * hisResY2 = (TH3F*)fRMSY->At(2);
- hisResY2->FitSlicesZ();
- TH2D * hisResY22 = (TH2D*)gDirectory->Get("RMS Y2_1");
- //
- ps->NewPage();
- c.cd(1);
- hisResY02->Fit(frms, "qn0");
- hisResY02->Draw("surf1");
- c.cd(2);
- MakeDiff(hisResY02,frms)->Draw("surf1");
- c.Update();
-// c.SaveAs("RMSYPad0.eps");
- ps->NewPage();
- c.cd(1);
- hisResY12->Fit(frms, "qn0"); // valgrind several blocks possibly lost
- hisResY12->Draw("surf1");
- c.cd(2);
- MakeDiff(hisResY12,frms)->Draw("surf1");
- c.Update();
-// c.SaveAs("RMSYPad1.eps");
- ps->NewPage();
- c.cd(1);
- hisResY22->Fit(frms, "qn0");
- hisResY22->Draw("surf1");
- c.cd(2);
- MakeDiff(hisResY22,frms)->Draw("surf1");
- c.Update();
-// c.SaveAs("RMSYPad2.eps");
-
- // create histogramw for Z-RMS
- TH3F * hisResZ0 = (TH3F*)fRMSZ->At(0);
- hisResZ0->FitSlicesZ();
- TH2D * hisResZ02 = (TH2D*)gDirectory->Get("RMS Z0_1");
- TH3F * hisResZ1 = (TH3F*)fRMSZ->At(1);
- hisResZ1->FitSlicesZ();
- TH2D * hisResZ12 = (TH2D*)gDirectory->Get("RMS Z1_1");
- TH3F * hisResZ2 = (TH3F*)fRMSZ->At(2);
- hisResZ2->FitSlicesZ();
- TH2D * hisResZ22 = (TH2D*)gDirectory->Get("RMS Z2_1");
- //
- ps->NewPage();
- c.cd(1);
- hisResZ02->Fit(frms, "qn0"); // valgrind
- hisResZ02->Draw("surf1");
- c.cd(2);
- MakeDiff(hisResZ02,frms)->Draw("surf1");
- c.Update();
-// c.SaveAs("RMSZPad0.eps");
- ps->NewPage();
- c.cd(1);
- hisResZ12->Fit(frms, "qn0");
- hisResZ12->Draw("surf1");
- c.cd(2);
- MakeDiff(hisResZ12,frms)->Draw("surf1");
- c.Update();
-// c.SaveAs("RMSZPad1.eps");
- ps->NewPage();
- c.cd(1);
- hisResZ22->Fit(frms, "qn0"); // valgrind 1 block possibly lost
- hisResZ22->Draw("surf1");
- c.cd(2);
- MakeDiff(hisResZ22,frms)->Draw("surf1");
- c.Update();
-// c.SaveAs("RMSZPad2.eps");
-
- // write calculated resoltuion rms to 'rms.txt'
- ofstream filerms("rms.txt");
- filerms<<"Pad 0.75 cm"<<"\n";
- hisResY02->Fit(frms, "qn0"); // valgrind 23 blocks indirectly lost
- filerms<<"Y\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
- hisResZ02->Fit(frms, "qn0"); // valgrind 23 blocks indirectly lost
- filerms<<"Z\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
- //
- filerms<<"Pad 1.00 cm"<<1<<"\n";
- hisResY12->Fit(frms, "qn0"); // valgrind 3,256 bytes in 22 blocks are indirectly lost
- filerms<<"Y\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
- hisResZ12->Fit(frms, "qn0"); // valgrind 66,036 bytes in 3 blocks are still reachable
- filerms<<"Z\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
- //
- filerms<<"Pad 1.50 cm"<<0<<"\n";
- hisResY22->Fit(frms, "qn0"); // valgrind 40,139 bytes in 11 blocks are still reachable 330,180 bytes in 15 blocks are possibly lost
- filerms<<"Y\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
- hisResZ22->Fit(frms, "qn0");
- filerms<<"Z\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
-
- TH1::AddDirectory(kFALSE);
- gSystem->ChangeDirectory("..");
- ps->Close();
- delete ps;
- delete frms;
-}
void AliTPCcalibTracks::MakeResPlotsQTree(Int_t minEntries, const char* pathName){
TList* deltaZList = new TList;
TList* arrayAmpRowList = new TList;
TList* rejectedTracksList = new TList;
- TList* hclusList = new TList;
TList* clusterCutHistoList = new TList;
TList* arrayAmpList = new TList;
TList* arrayQDYList = new TList;
TList* arrayQDZList = new TList;
TList* arrayQRMSYList = new TList;
TList* arrayQRMSZList = new TList;
- TList* arrayChargeVsDriftlengthList = new TList;
- TList* calPadRegionChargeVsDriftlengthList = new TList;
- TList* hclusterPerPadrowList = new TList;
- TList* hclusterPerPadrowRawList = new TList;
TList* resolYList = new TList;
TList* resolZList = new TList;
TList* rMSYList = new TList;
while ( (calibTracks = dynamic_cast<AliTPCcalibTracks*> (listIterator->Next())) ){
// loop over all entries in the collectionList and get dataMembers into lists
- deltaYList->Add( calibTracks->GetfDeltaY() );
- deltaZList->Add( calibTracks->GetfDeltaZ() );
- arrayAmpRowList->Add(calibTracks->GetfArrayAmpRow());
- arrayAmpList->Add(calibTracks->GetfArrayAmp());
arrayQDYList->Add(calibTracks->GetfArrayQDY());
arrayQDZList->Add(calibTracks->GetfArrayQDZ());
arrayQRMSYList->Add(calibTracks->GetfArrayQRMSY());
resolZList->Add(calibTracks->GetfResolZ());
rMSYList->Add(calibTracks->GetfRMSY());
rMSZList->Add(calibTracks->GetfRMSZ());
- arrayChargeVsDriftlengthList->Add(calibTracks->GetfArrayChargeVsDriftlength());
- calPadRegionChargeVsDriftlengthList->Add(calibTracks->GetCalPadRegionchargeVsDriftlength());
- hclusList->Add(calibTracks->GetfHclus());
rejectedTracksList->Add(calibTracks->GetfRejectedTracksHisto());
clusterCutHistoList->Add(calibTracks->GetfClusterCutHisto());
- hclusterPerPadrowList->Add(calibTracks->GetfHclusterPerPadrow());
- hclusterPerPadrowRawList->Add(calibTracks->GetfHclusterPerPadrowRaw());
//
if (fCalPadClusterPerPad && calibTracks->GetfCalPadClusterPerPad())
fCalPadClusterPerPad->Add(calibTracks->GetfCalPadClusterPerPad());
// merge data members
if (GetDebugLevel() > 0) cout << "histogram's merge-functins are called... " << endl;
- fDeltaY->Merge(deltaYList);
- fDeltaZ->Merge(deltaZList);
- fHclus->Merge(hclusList);
fClusterCutHisto->Merge(clusterCutHistoList);
fRejectedTracksHisto->Merge(rejectedTracksList);
- fHclusterPerPadrow->Merge(hclusterPerPadrowList);
- fHclusterPerPadrowRaw->Merge(hclusterPerPadrowRawList);
TObjArray* objarray = 0;
TH1* hist = 0;
TList* histList = 0;
TIterator *objListIterator = 0;
- if (GetDebugLevel() > 0) cout << "merging fArrayAmpRows..." << endl;
- // merge fArrayAmpRows
- for (Int_t i = 0; i < fArrayAmpRow->GetEntriesFast(); i++ ) { // loop over data member, i<72
- objListIterator = arrayAmpRowList->MakeIterator();
- histList = new TList;
- while (( objarray = (TObjArray*)objListIterator->Next() )) {
- // loop over arrayAmpRowList, get TObjArray, get object at position i, cast it into TProfile
- hist = (TProfile*)(objarray->At(i));
- histList->Add(hist);
- }
- ((TProfile*)(fArrayAmpRow->At(i)))->Merge(histList);
- delete histList;
- delete objListIterator;
- }
-
- if (GetDebugLevel() > 0) cout << "merging fArrayAmps..." << endl;
- // merge fArrayAmps
- for (Int_t i = 0; i < fArrayAmp->GetEntriesFast(); i++ ) { // loop over data member, i<72
- TIterator *cobjListIterator = arrayAmpList->MakeIterator();
- histList = new TList;
- while (( objarray = (TObjArray*)cobjListIterator->Next() )) {
- // loop over arrayAmpList, get TObjArray, get object at position i, cast it into TH1F
- hist = (TH1F*)(objarray->At(i));
- histList->Add(hist);
- }
- ((TH1F*)(fArrayAmp->At(i)))->Merge(histList);
- delete histList;
- delete cobjListIterator;
- }
-
+
if (GetDebugLevel() > 0) cout << "merging fArrayQDY..." << endl;
// merge fArrayQDY
for (Int_t i = 0; i < fArrayQDY->GetEntriesFast(); i++) { // loop over data member, i < 300
delete objListIterator;
}
- if (GetDebugLevel() > 0) cout << "merging fArrayChargeVsDriftlength..." << endl;
- // merge fArrayChargeVsDriftlength
- for (Int_t i = 0; i < fArrayChargeVsDriftlength->GetEntriesFast(); i++) { // loop over data member, i < 300
- objListIterator = arrayChargeVsDriftlengthList->MakeIterator();
- histList = new TList;
- while (( objarray = (TObjArray*)objListIterator->Next() )) {
- // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TProfile
- if (!objarray) continue;
- hist = (TProfile*)(objarray->At(i));
- histList->Add(hist);
- }
- ((TProfile*)(fArrayChargeVsDriftlength->At(i)))->Merge(histList);
- delete histList;
- delete objListIterator;
- }
-
- if (GetDebugLevel() > 0) cout << "merging fcalPadRegionChargeVsDriftlength..." << endl;
- // merge fcalPadRegionChargeVsDriftlength
- AliTPCCalPadRegion *cpr = 0x0;
-
- /*
- TIterator *regionIterator = fcalPadRegionChargeVsDriftlength->MakeIterator();
- while (hist = (TProfile*)regionIterator->Next()) {
- // loop over all calPadRegion's in destination calibTracks object
- objListIterator = calPadRegionChargeVsDriftlengthList->MakeIterator();
- while (( cpr = (AliTPCCalPadRegion*)objListIterator->Next() )) {
-
-
- hist->Merge(...);
- }
- */
- for (UInt_t isec = 0; isec < 36; isec++) {
- for (UInt_t padSize = 0; padSize < 3; padSize++){
- objListIterator = calPadRegionChargeVsDriftlengthList->MakeIterator();
- histList = new TList;
- while (( cpr = (AliTPCCalPadRegion*)objListIterator->Next() )) {
- // loop over calPadRegionChargeVsDriftlengthList, get AliTPCCalPadRegion, get object
- hist = (TProfile*)cpr->GetObject(isec, padSize);
- histList->Add(hist);
- }
- ((TProfile*)(fcalPadRegionChargeVsDriftlength->GetObject(isec, padSize)))->Merge(histList);
- delete histList;
- delete objListIterator;
- }
- }
}
-AliTPCcalibTracks* AliTPCcalibTracks::TestMerge(AliTPCcalibTracks *ct, AliTPCClusterParam *clusterParam, Int_t nCalTracks){
- //
- // function to test AliTPCcalibTrack::Merge:
- // in the file 'f' is a AliTPCcalibTrack object with name "calibTracks"
- // this object is appended 'nCalTracks' times to a TList
- // A new AliTPCcalibTrack object is created which merges the list
- // this object is returned
- //
- /*
- // .L AliTPCcalibTracks.cxx+g
- .L libTPCcalib.so
- TFile f("Output.root");
- AliTPCcalibTracks* calTracks = (AliTPCcalibTracks*)f.Get("calibTracks");
- //f.Close();
- TFile clusterParamFile("/u/lbozyk/calibration/workdir/calibTracks/TPCClusterParam.root");
- AliTPCClusterParam *clusterParam = (AliTPCClusterParam *) clusterParamFile.Get("Param");
- clusterParamFile.Close();
-
- AliTPCcalibTracks::TestMerge(calTracks, clusterParam);
- */
- TList *list = new TList();
- if (ct == 0 || clusterParam == 0) return 0;
- cout << "making list with " << nCalTracks << " AliTPCcalibTrack objects" << endl;
- for (Int_t i = 0; i < nCalTracks; i++) {
- if (i%10==0) cout << "Adding element " << i << " of " << nCalTracks << endl;
- list->Add(new AliTPCcalibTracks(*ct));
- }
-
- // only for check at the end
- AliTPCcalibTracks* cal1 = new AliTPCcalibTracks(*ct);
- Double_t cal1Entries = ((TH1F*)cal1->GetfArrayAmpRow()->At(5))->GetEntries();
-// Double_t cal1Entries = 5; //((TH1F*)ct->GetfArrayAmpRow()->At(5))->GetEntries();
-
- cout << "The list contains " << list->GetEntries() << " entries. " << endl;
-
-
- AliTPCcalibTracksCuts *cuts = new AliTPCcalibTracksCuts(20, 0.4, 0.5, 0.13, 0.018);
- AliTPCcalibTracks* cal = new AliTPCcalibTracks("calTracksMerged", "calTracksMerged", clusterParam, cuts, 5);
- cal->Merge(list);
-
- cout << "cal->GetfArrayAmpRow()->At(5)->Print():" << endl;
- cal->GetfArrayAmpRow()->At(5)->Print();
- Double_t calEntries = ((TH1F*)cal->GetfArrayAmpRow()->At(5))->GetEntries();
-
- cout << "cal1->GetfArrayAmpRow()->At(5))->GetEntries() = " << cal1Entries << endl;
- cout << " cal->GetfArrayAmpRow()->At(5))->GetEntries() = " << calEntries << endl;
- printf("That means there were %f / %f = %f AliTPCcalibTracks-Objects merged. \n",
- calEntries, cal1Entries, ((Double_t)calEntries/cal1Entries));
-
- return cal;
-
-}
-
-
-void AliTPCcalibTracks::MakeQPosNormAll(TTree * chainres, AliTPCClusterParam * param, Int_t maxPoints){
- //
- // Make position corrections
- // for the moment Only using debug streamer
- // chainres - debug tree
- // param - parameters to be updated
- // maxPoints - maximal number of points using for fit
- // verbose - print info flag
- //
- // Current implementation - only using debug streamers
- //
-
- /*
- //Defaults
- Int_t maxPoints=100000;
- */
- /*
- Usage:
- //0. Load libraries
- gSystem->Load("libANALYSIS");
- gSystem->Load("libSTAT");
- gSystem->Load("libTPCcalib");
-
-
- //1. Load Parameters to be modified:
- //e.g:
-
- AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
- AliCDBManager::Instance()->SetRun(0)
- AliTPCClusterParam * param = AliTPCcalibDB::Instance()->GetClusterParam();
-
- //2. Load chain from debug streamers
- //
- //e.g
- gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
- gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
- AliXRDPROOFtoolkit tool;
- TChain * chainres = tool.MakeChain("tracks.txt","ResolCl",0,10200);
- chainres->Lookup();
- //3. Do fits and store results
- //
- AliTPCcalibTracks::MakeQPosNormAll(chainres,param,200000,0) ;
- TFile f("paramout.root","recreate");
- param->Write("clusterParam");
- f.Close();
- //4. Verification
- TFile f2("paramout.root");
- AliTPCClusterParam *param2 = (AliTPCClusterParam*)f2.Get("clusterParam");
- param2->SetInstance(param2);
- chainres->Draw("fitZ0:AliTPCClusterParam::SPosCorrection(1,0,Cl.fPad,Cl.fTimeBin,Cl.fZ,Cl.fSigmaY2,Cl.fSigmaZ2,Cl.fMax)","Cl.fDetector<36","",10000); // should be line
-
- */
-
-
- TStatToolkit toolkit;
- Double_t chi2;
- TVectorD fitParamY0;
- TVectorD fitParamY1;
- TVectorD fitParamZ0;
- TVectorD fitParamZ1;
- TMatrixD covMatrix;
- Int_t npoints;
-
- chainres->SetAlias("dp","(-1+(Cl.fZ>0)*2)*((Cl.fPad-int(Cl.fPad))-0.5)");
- chainres->SetAlias("dt","(-1+(Cl.fZ>0)*2)*((Cl.fTimeBin-0.66-int(Cl.fTimeBin-0.66))-0.5)");
- chainres->SetAlias("sp","(sin(dp*pi)-dp*pi)");
- chainres->SetAlias("st","(sin(dt)-dt)");
- //
- chainres->SetAlias("di","sqrt(1.-abs(Cl.fZ/250.))");
- //
- //
- //
- TCut cutA("1");
- TString fstringY="";
- //
- fstringY+="(dp)++"; //1
- fstringY+="(dp)*di++"; //2
- fstringY+="(sp)++"; //3
- fstringY+="(sp)*di++"; //4
- TString fstringZ="";
- fstringZ+="(dt)++"; //1
- fstringZ+="(dt)*di++"; //2
- fstringZ+="(st)++"; //3
- fstringZ+="(st)*di++"; //4
- //
- // Z corrections
- //
- TString *strZ0 = toolkit.FitPlane(chainres,"(Cl.fZ-PZ0.fElements[0]):CSigmaZ0",fstringZ.Data(), "Cl.fDetector<36"+cutA, chi2,npoints,fitParamZ0,covMatrix,-1,0,maxPoints);
- printf("Z0 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
- param->PosZcor(0) = (TVectorD*) fitParamZ0.Clone();
- //
- TString *strZ1 = toolkit.FitPlane(chainres,"(Cl.fZ-PZ0.fElements[0]):CSigmaZ0",fstringZ.Data(), "Cl.fDetector>36"+cutA, chi2,npoints,fitParamZ1,covMatrix,-1,0,maxPoints);
- //
- printf("Z1 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
- param->PosZcor(1) = (TVectorD*) fitParamZ1.Clone();
- param->PosZcor(2) = (TVectorD*) fitParamZ1.Clone();
- //
- // Y corrections
- //
- TString *strY0 = toolkit.FitPlane(chainres,"(Cl.fY-PY0.fElements[0]):CSigmaY0",fstringY.Data(), "Cl.fDetector<36"+cutA, chi2,npoints,fitParamY0,covMatrix,-1,0,maxPoints);
- printf("Y0 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
- param->PosYcor(0) = (TVectorD*) fitParamY0.Clone();
-
-
- TString *strY1 = toolkit.FitPlane(chainres,"(Cl.fY-PY0.fElements[0]):CSigmaY0",fstringY.Data(), "Cl.fDetector>36"+cutA, chi2,npoints,fitParamY1,covMatrix,-1,0,maxPoints);
- //
- printf("Y1 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
- param->PosYcor(1) = (TVectorD*) fitParamY1.Clone();
- param->PosYcor(2) = (TVectorD*) fitParamY1.Clone();
- //
- //
- //
- chainres->SetAlias("fitZ0",strZ0->Data());
- chainres->SetAlias("fitZ1",strZ1->Data());
- chainres->SetAlias("fitY0",strY0->Data());
- chainres->SetAlias("fitY1",strY1->Data());
- // chainres->Draw("Cl.fZ-PZ0.fElements[0]","CSigmaY0<0.7&&CSigmaZ0<0.7"+cutA,"",10000);
-}
-
void AliTPCcalibTracks::MakeHistos(){
TString axisName[8];
//
- binsTrack[0]=100;
+ binsTrack[0]=200;
axisName[0] ="var";
binsTrack[1] =3;
xminTrack[1] =0; xmaxTrack[1]=3;
axisName[1] ="pad type";
//
- binsTrack[2] =10;
- xminTrack[2] =0; xmaxTrack[2]=1;
- axisName[2] ="drift length";
+ binsTrack[2] =20;
+ xminTrack[2] =-250; xmaxTrack[2]=250;
+ axisName[2] ="z";
//
binsTrack[3] =10;
xminTrack[3] =1; xmaxTrack[3]=400;
axisName[3] ="Qmax";
//
- binsTrack[4] =10;
+ binsTrack[4] =20;
xminTrack[4] =0; xmaxTrack[4]=1;
axisName[4] ="cog";
//
- binsTrack[5] =10;
- xminTrack[5] =0; xmaxTrack[5]=2;
- axisName[5] ="tan(phi)";
+ binsTrack[5] =15;
+ xminTrack[5] =-1.5; xmaxTrack[5]=1.5;
+ axisName[5] ="tan(angle)";
//
- binsTrack[6] =10;
- xminTrack[6] =0; xmaxTrack[6]=2;
- axisName[6] ="tan(theta)";
//
- xminTrack[0] =-0.5; xmaxTrack[0]=0.5;
- fHisDeltaY=new THnSparseS("#Delta_{y} (cm)","#Delta_{y} (cm)", 7, binsTrack,xminTrack, xmaxTrack);
- xminTrack[0] =-0.5; xmaxTrack[0]=0.5;
- fHisDeltaZ=new THnSparseS("#Delta_{z} (cm)","#Delta_{z} (cm)", 7, binsTrack,xminTrack, xmaxTrack);
+ xminTrack[0] =-1.5; xmaxTrack[0]=1.5;
+ fHisDeltaY=new THnSparseF("#Delta_{y} (cm)","#Delta_{y} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
+ xminTrack[0] =-1.5; xmaxTrack[0]=1.5;
+ fHisDeltaZ=new THnSparseF("#Delta_{z} (cm)","#Delta_{z} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
xminTrack[0] =0.; xmaxTrack[0]=0.5;
- fHisRMSY=new THnSparseS("#RMS_{y} (cm)","#RMS_{y} (cm)", 7, binsTrack,xminTrack, xmaxTrack);
+ fHisRMSY=new THnSparseF("#RMS_{y} (cm)","#RMS_{y} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
xminTrack[0] =0.; xmaxTrack[0]=0.5;
- fHisRMSZ=new THnSparseS("#RMS_{z} (cm)","#RMS_{z} (cm)", 7, binsTrack,xminTrack, xmaxTrack);
+ fHisRMSZ=new THnSparseF("#RMS_{z} (cm)","#RMS_{z} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
xminTrack[0] =0.; xmaxTrack[0]=100;
- fHisQmax=new THnSparseS("Qmax (ADC)","Qmax (ADC)", 7, binsTrack,xminTrack, xmaxTrack);
+ fHisQmax=new THnSparseF("Qmax (ADC)","Qmax (ADC)", 6, binsTrack,xminTrack, xmaxTrack);
xminTrack[0] =0.; xmaxTrack[0]=250;
- fHisQtot=new THnSparseS("Qtot (ADC)","Qtot (ADC)", 7, binsTrack,xminTrack, xmaxTrack);
+ fHisQtot=new THnSparseF("Qtot (ADC)","Qtot (ADC)", 6, binsTrack,xminTrack, xmaxTrack);
+
+
+ for (Int_t ivar=0;ivar<6;ivar++){
+ fHisDeltaY->GetAxis(ivar)->SetName(axisName[ivar].Data());
+ fHisDeltaZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
+ fHisRMSY->GetAxis(ivar)->SetName(axisName[ivar].Data());
+ fHisRMSZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
+ fHisQmax->GetAxis(ivar)->SetName(axisName[ivar].Data());
+ fHisQtot->GetAxis(ivar)->SetName(axisName[ivar].Data());
+ fHisDeltaY->GetAxis(ivar)->SetTitle(axisName[ivar].Data());
+ fHisDeltaZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
+ fHisRMSY->GetAxis(ivar)->SetName(axisName[ivar].Data());
+ fHisRMSZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
+ fHisQmax->GetAxis(ivar)->SetName(axisName[ivar].Data());
+ fHisQtot->GetAxis(ivar)->SetName(axisName[ivar].Data());
+ }
+
+
BinLogX(fHisDeltaY,3);
BinLogX(fHisDeltaZ,3);
BinLogX(fHisRMSY,3);
if (calib->fHisRMSY) fHisRMSY->Add(calib->fHisRMSY);
if (calib->fHisRMSZ) fHisRMSZ->Add(calib->fHisRMSZ);
}
+
+
+
+void AliTPCcalibTracks::MakeSummaryTree(THnSparse *hisInput, TTreeSRedirector *pcstream, Int_t ptype){
+ //
+ // Dump summary info
+ //
+ // 0.OBJ: TAxis var
+ // 1.OBJ: TAxis pad type
+ // 2.OBJ: TAxis z
+ // 3.OBJ: TAxis Qmax
+ // 4.OBJ: TAxis cog
+ // 5.OBJ: TAxis tan(angle)
+ //
+ if (ptype>3) return;
+ Int_t idim[6]={0,1,2,3,4,5};
+ TString hname[4]={"dy","dz","rmsy","rmsz"};
+ //
+ Int_t nbins5=hisInput->GetAxis(5)->GetNbins();
+ Int_t first5=hisInput->GetAxis(5)->GetFirst();
+ Int_t last5 =hisInput->GetAxis(5)->GetLast();
+ //
+ for (Int_t ibin5=first5-1; ibin5<=last5; ibin5+=1){ // axis 5 - local angle
+ Bool_t bin5All=kTRUE;
+ if (ibin5>=first5){
+ hisInput->GetAxis(5)->SetRange(TMath::Max(ibin5-1,first5),TMath::Min(ibin5+1,nbins5));
+ if (ibin5==first5) hisInput->GetAxis(5)->SetRange(TMath::Max(ibin5,first5),TMath::Min(ibin5,nbins5));
+ bin5All=kFALSE;
+ }
+ Double_t x5= hisInput->GetAxis(5)->GetBinCenter(ibin5);
+ THnSparse * his5= hisInput->Projection(5,idim); //projected histogram according selection 5
+ Int_t nbins4=his5->GetAxis(4)->GetNbins();
+ Int_t first4=his5->GetAxis(4)->GetFirst();
+ Int_t last4 =his5->GetAxis(4)->GetLast();
+
+ for (Int_t ibin4=first4-1; ibin4<=last4; ibin4+=1){ // axis 4 - cog
+ Bool_t bin4All=kTRUE;
+ if (ibin4>=first4){
+ his5->GetAxis(4)->SetRange(TMath::Max(ibin4+1,first4),TMath::Min(ibin4-1,nbins4));
+ if (ibin4==first4||ibin4==last4) his5->GetAxis(4)->SetRange(TMath::Max(ibin4,first4),TMath::Min(ibin4,nbins4));
+ bin4All=kFALSE;
+ }
+ Double_t x4= his5->GetAxis(4)->GetBinCenter(ibin4);
+ THnSparse * his4= his5->Projection(4,idim); //projected histogram according selection 4
+ //
+ Int_t nbins3=his4->GetAxis(3)->GetNbins();
+ Int_t first3=his4->GetAxis(3)->GetFirst();
+ Int_t last3 =his4->GetAxis(3)->GetLast();
+ //
+ for (Int_t ibin3=first3-1; ibin3<=last3; ibin3+=1){ // axis 3 - Qmax
+ Bool_t bin3All=kTRUE;
+ if (ibin3>=first3){
+ his4->GetAxis(3)->SetRange(TMath::Max(ibin3,first3),TMath::Min(ibin3,nbins3));
+ bin3All=kFALSE;
+ }
+ Double_t x3= his4->GetAxis(3)->GetBinCenter(ibin3);
+ THnSparse * his3= his4->Projection(3,idim); //projected histogram according selection 3
+ //
+ Int_t nbins2 = his3->GetAxis(2)->GetNbins();
+ Int_t first2 = his3->GetAxis(2)->GetFirst();
+ Int_t last2 = his3->GetAxis(2)->GetLast();
+ //
+ for (Int_t ibin2=first2-1; ibin2<=last2; ibin2+=1){ // axis 2 - z
+ Bool_t bin2All=kTRUE;
+ Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
+ if (ibin2>=first2){
+ his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,first2),TMath::Min(ibin2+1,nbins2));
+ if (ibin2==first2||ibin2==last2||TMath::Abs(x2)<20) his3->GetAxis(2)->SetRange(TMath::Max(ibin2,first2),TMath::Min(ibin2,nbins2));
+ bin2All=kFALSE;
+ }
+ THnSparse * his2= his3->Projection(2,idim); //projected histogram according selection 2
+ //
+ Int_t nbins1 = his2->GetAxis(1)->GetNbins();
+ Int_t first1 = his2->GetAxis(1)->GetFirst();
+ Int_t last1 = his2->GetAxis(1)->GetLast();
+ for (Int_t ibin1=first1-1; ibin1<=last1; ibin1++){ //axis 1 - pad type
+ Bool_t bin1All=kTRUE;
+ if (ibin1>=first1){
+ his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1));
+ bin1All=kFALSE;
+ }
+ Double_t x1= TMath::Nint(his2->GetAxis(1)->GetBinCenter(ibin1)-0.5);
+ TH1 * hisDelta = his2->Projection(0);
+ Double_t entries = hisDelta->GetEntries();
+ Double_t mean=0, rms=0;
+ if (entries>10){
+ mean = hisDelta->GetMean();
+ rms = hisDelta->GetRMS();
+ hisDelta->GetXaxis()->SetRangeUser(mean-4*rms,mean+4*rms);
+ mean = hisDelta->GetMean();
+ rms = hisDelta->GetRMS();
+ }
+ //
+ (*pcstream)<<hname[ptype].Data()<<
+ // flag - intgrated values for given bin
+ "angleA="<<bin5All<<
+ "cogA="<<bin4All<<
+ "qmaxA="<<bin3All<<
+ "zA="<<bin2All<<
+ "ipadA="<<bin1All<<
+ // center of bin value
+ "angle="<<x5<<
+ "cog="<<x4<<
+ "qmax="<<x3<<
+ "z="<<x2<<
+ "ipad="<<x1<<
+ // mean values
+ //
+ "entries="<<entries<<
+ "mean="<<mean<<
+ "rms="<<rms<<
+ "ptype="<<ptype<< //
+ "\n";
+ delete hisDelta;
+ printf("%f\t%f\t%f\t%f\t%f\t%f\t%f\n",x5,x4,x3,x2,x1, entries,mean);
+ }//loop z
+ delete his2;
+ } // loop Q max
+ delete his3;
+ } // loop COG
+ delete his4;
+ }//loop local angle
+ delete his5;
+ }
+}
///////////////////////////////////////////////////////////////////////////////
// //
// Class to analyse tracks for calibration //
-// to be used as a component in AliTPCSelectorTracks //
+// to be used as a component in AliTPCSelectorTracks //
// In the constructor you have to specify name and title //
// to get the Object out of a file. //
// The parameter 'clusterParam', a AliTPCClusterParam object //
class AliESDtrack;
class AliTPCclusterMI;
class AliTPCcalibTracksCuts;
-class AliTPCCalPadRegion;
class AliTPCCalPad;
class TChain;
class TTree;
virtual Long64_t Merge(TCollection *li);
void AddHistos(AliTPCcalibTracks* calib);
void MakeResPlotsQTree(Int_t minEntries = 100, const char* pathName = "plots");
- static void MakeQPosNormAll(TTree * chain, AliTPCClusterParam * param, Int_t maxPoints=1000000);
void MakeReport(Int_t stat, const char* pathName = "plots"); // calls all functions that procude pictures, results are written to pathName, stat is the minimal statistic threshold
- //
-
-
+ //
Int_t AcceptTrack(AliTPCseed * track);
void FillResolutionHistoLocal(AliTPCseed * track); // the MAIN-FUNCTION, called for each track to fill the histograms, called by Process(...)
- static TH2D* MakeDiff(TH2D * hfit, TF2 * func);
- static AliTPCcalibTracks* TestMerge(AliTPCcalibTracks *ct, AliTPCClusterParam *clusterParam, Int_t nCalTracks = 50);
void SetStyle() const;
- void Draw(Option_t* opt); // draw some exemplaric histograms for fast result-check
- void MakeAmpPlots(Int_t stat, const char* pathName = "plots");
- void MakeDeltaPlots(const char* pathName = "plots");
- void MakeChargeVsDriftLengthPlotsOld(const char* pathName = "plots");
- void MakeChargeVsDriftLengthPlots(const char* pathName = "plots");
- void FitResolutionNew(const char* pathName = "plots");
- void FitRMSNew(const char* pathName = "plots");
- TObjArray* GetfArrayAmpRow() const {return fArrayAmpRow;}
- TObjArray* GetfArrayAmp() const {return fArrayAmp;}
TObjArray* GetfArrayQDY() const {return fArrayQDY;}
TObjArray* GetfArrayQDZ() const {return fArrayQDZ;}
TObjArray* GetfArrayQRMSY() const {return fArrayQRMSY;}
TObjArray* GetfArrayQRMSZ() const {return fArrayQRMSZ;}
- TObjArray* GetfArrayChargeVsDriftlength() const {return fArrayChargeVsDriftlength;}
- TH1F* GetfDeltaY() const {return fDeltaY;}
- TH1F* GetfDeltaZ() const {return fDeltaZ;}
TObjArray* GetfResolY() const {return fResolY;}
TObjArray* GetfResolZ() const {return fResolZ;}
TObjArray* GetfRMSY() const {return fRMSY;}
TObjArray* GetfRMSZ() const {return fRMSZ;}
- TH1I* GetfHclus() const {return fHclus;}
TH1I* GetfRejectedTracksHisto() const {return fRejectedTracksHisto;}
- TH1I* GetfHclusterPerPadrow() const {return fHclusterPerPadrow;}
- TH1I* GetfHclusterPerPadrowRaw() const {return fHclusterPerPadrowRaw;}
TH2I* GetfClusterCutHisto() const {return fClusterCutHisto;}
AliTPCCalPad* GetfCalPadClusterPerPad() const {return fCalPadClusterPerPad; }
AliTPCCalPad* GetfCalPadClusterPerPadRaw() const {return fCalPadClusterPerPadRaw;}
- AliTPCCalPadRegion* GetCalPadRegionchargeVsDriftlength() const {return fcalPadRegionChargeVsDriftlength;}
AliTPCcalibTracksCuts* GetCuts() {return fCuts;}
void MakeHistos(); //make THnSparse
+ static void MakeSummaryTree(THnSparse *hisInput, TTreeSRedirector *pcstream, Int_t ptype);
protected:
private:
static Int_t GetBin(Int_t iq, Int_t pad);
static Float_t GetQ(Int_t bin);
static Float_t GetPad(Int_t bin);
- void FillResolutionHistoLocalDebugPart(AliTPCseed *track, AliTPCclusterMI *cluster0, Int_t irow, Float_t angley, Float_t anglez, Int_t nclFound, Int_t kDelta);
AliTPCClusterParam *fClusterParam; // pointer to cluster parameterization
AliTPCROC *fROC; //!
+public:
THnSparse *fHisDeltaY; // THnSparse - delta Y
THnSparse *fHisDeltaZ; // THnSparse - delta Z
THnSparse *fHisRMSY; // THnSparse - rms Y
THnSparse *fHisQmax; // THnSparse - qmax
THnSparse *fHisQtot; // THnSparse - qtot
-
- TObjArray *fArrayAmpRow; // array with amplitudes versus row for given sector
- TObjArray *fArrayAmp; // array with amplitude for sectors
+private:
TObjArray *fArrayQDY; // q binned delta Y histograms
TObjArray *fArrayQDZ; // q binned delta Z histograms
TObjArray *fArrayQRMSY; // q binned delta Y histograms
TObjArray *fArrayQRMSZ; // q binned delta Z histograms
- TObjArray *fArrayChargeVsDriftlength; // array of arrays of TProfiles with charge vs. driftlength for each padsize and sector
- AliTPCCalPadRegion *fcalPadRegionChargeVsDriftlength; // CalPadRegion, one TProfile for charge vs. driftlength for each padsize and sector
- TH1F *fDeltaY; // integrated delta y histo
- TH1F *fDeltaZ; // integrated delta z histo
TObjArray *fResolY; // array of resolution histograms Y
TObjArray *fResolZ; // array of resolution histograms Z
TObjArray *fRMSY; // array of RMS histograms Y
TObjArray *fRMSZ; // array of RMS histograms Z
AliTPCcalibTracksCuts *fCuts; // object with cuts, that is passed to the constructor
- TH1I *fHclus; // number of clusters per track
TH1I *fRejectedTracksHisto; // histogram of rejecteced tracks, the number coresponds to the failed cut
- TH1I *fHclusterPerPadrow; // histogram showing the number of clusters per padRow
- TH1I *fHclusterPerPadrowRaw;// histogram showing the number of clusters per padRow before cuts on clusters are applied
TH2I *fClusterCutHisto; // histogram showing in which padRow the clusters were cutted by which criterium
AliTPCCalPad *fCalPadClusterPerPad; // AliTPCCalPad showing the number of clusters per Pad
AliTPCCalPad *fCalPadClusterPerPadRaw; // AliTPCCalPad showing the number of clusters per Pad before cuts on clusters are applied
- ClassDef(AliTPCcalibTracks,1)
+ ClassDef(AliTPCcalibTracks,2)
};
#include "TCint.h"
#include "AliTPCcalibV0.h"
#include "AliV0.h"
+#include "TRandom.h"
+#include "TTreeStream.h"
+#include "AliTPCcalibDB.h"
+#include "AliTPCCorrection.h"
+#include "AliGRPObject.h"
+#include "AliTPCTransform.h"
AliTPCcalibV0::AliTPCcalibV0() :
AliTPCcalibBase(),
+ fV0Tree(0),
+ fHPTTree(0),
fStack(0),
fESD(0),
fPdg(0),
fParticles(0),
fV0s(0),
- fGammas(0),
- fV0type(0),
- fTPCdEdx(0),
- fTPCdEdxPi(0),
- fTPCdEdxEl(0),
- fTPCdEdxP(0)
+ fGammas(0)
+{
+
+}
+AliTPCcalibV0::AliTPCcalibV0(const Text_t *name, const Text_t *title):
+ AliTPCcalibBase(),
+ fV0Tree(0),
+ fHPTTree(0),
+ fStack(0),
+ fESD(0),
+ fPdg(0),
+ fParticles(0),
+ fV0s(0),
+ fGammas(0)
{
fPdg = new TDatabasePDG;
// create output histograms
- fTPCdEdx = new TH2F("TPCdEdX", "dE/dX; BetaGamma; TPC signal (a.u.)", 1000, 0.1, 10000, 300, 0, 300);
- BinLogX(fTPCdEdx);
- fTPCdEdxPi = new TH2F("TPCdEdXPi","dE/dX; BetaGamma; TPC signal (a.u.)", 1000, 0.1, 10000, 300, 0, 300);
- fTPCdEdxEl = new TH2F("TPCdEdXEl","dE/dX; BetaGamma; TPC signal (a.u.)", 1000, 0.1, 10000, 300, 0, 300);
- fTPCdEdxP = new TH2F("TPCdEdXP", "dE/dX; BetaGamma; TPC signal (a.u.)", 1000, 0.1, 10000, 300, 0, 300);
-
-
+ SetName(name);
+ SetTitle(title);
}
AliTPCcalibV0::~AliTPCcalibV0(){
//
//
//
+ delete fV0Tree;
+ delete fHPTTree;
}
//
fESD = esd;
AliKFParticle::SetField(esd->GetMagneticField());
+ if (TMath::Abs(AliTracker::GetBz())<1) return;
+ DumpToTree(esd);
+ DumpToTreeHPT(esd);
+ return;
+ //
MakeV0s();
if (stack) {
fStack = stack;
}
}
+void AliTPCcalibV0::DumpToTreeHPT(AliESDEvent *esd){
+ //
+ // Dump V0s fith full firend information to the
+ //
+ if (TMath::Abs(AliTracker::GetBz())<1) return;
+ const Int_t kMinCluster=110;
+ const Float_t kMinPt =3.;
+ AliESDfriend *esdFriend=static_cast<AliESDfriend*>(esd->FindListObject("AliESDfriend"));
+ if (!esdFriend) {
+ Printf("ERROR: esdFriend not available");
+ return;
+ }
+ //
+ Int_t ntracks=esd->GetNumberOfTracks();
+ for (Int_t i=0;i<ntracks;++i) {
+ Bool_t isOK=kFALSE;
+ AliESDtrack *track = esd->GetTrack(i);
+ if (track->GetTPCncls()<kMinCluster) continue;
+ AliESDfriendTrack *friendTrack = esdFriend->GetTrack(i);
+ if (!friendTrack) continue;
+ if (TMath::Abs(AliTracker::GetBz())>1){ // cut on momenta if measured
+ if (track->Pt()>kMinPt) isOK=kTRUE;
+ }
+ if (TMath::Abs(AliTracker::GetBz())<1){ // require primary track for the B field OFF data
+ Bool_t isAccepted=kTRUE;
+ if (!track->IsOn(AliESDtrack::kITSrefit)) isAccepted=kFALSE;
+ if (!track->IsOn(AliESDtrack::kTPCrefit)) isAccepted=kFALSE;
+ if (!track->IsOn(AliESDtrack::kTOFout)) isAccepted=kFALSE;
+ Float_t dvertex[2],cvertex[3];
+ track->GetImpactParametersTPC(dvertex,cvertex);
+ if (TMath::Abs(dvertex[0]/TMath::Sqrt(cvertex[0]+0.01))>20) isAccepted=kFALSE;
+ if (TMath::Abs(dvertex[1]/TMath::Sqrt(TMath::Abs(cvertex[2]+0.01)))>20) isAccepted=kFALSE;
+ track->GetImpactParameters(dvertex,cvertex);
+ if (TMath::Abs(dvertex[0]/TMath::Sqrt(cvertex[0]+0.01))>10) isAccepted=kFALSE;
+ if (TMath::Abs(dvertex[1]/TMath::Sqrt(TMath::Abs(cvertex[2]+0.01)))>10) isAccepted=kFALSE;
+ if (!isAccepted) isOK=kFALSE;
+ }
+
+ if (!isOK) continue;
+ //
+
+ TObject *calibObject;
+ AliTPCseed *seed = 0;
+ for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
+ if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
+ }
+ if (!seed) continue;
+ if (!fHPTTree) {
+ fHPTTree = new TTree("HPT","HPT");
+ fHPTTree->SetDirectory(0);
+ }
+ if (fHPTTree->GetEntries()==0){
+ //
+ fHPTTree->SetDirectory(0);
+ fHPTTree->Branch("t.",&track);
+ fHPTTree->Branch("ft.",&friendTrack);
+ fHPTTree->Branch("s.",&seed);
+ }else{
+ fHPTTree->SetBranchAddress("t.",&track);
+ fHPTTree->SetBranchAddress("ft.",&friendTrack);
+ fHPTTree->SetBranchAddress("s.",&seed);
+ }
+ fHPTTree->Fill();
+ //
+ }
+}
+
+
+
+void AliTPCcalibV0::DumpToTree(AliESDEvent *esd){
+ //
+ // Dump V0s fith full firend information to the
+ //
+ Int_t nV0s = fESD->GetNumberOfV0s();
+ const Int_t kMinCluster=110;
+ const Double_t kDownscale=0.01;
+ const Float_t kMinR =0;
+ const Float_t kMinPt =1.0;
+ const Float_t kMinMinPt =0.7;
+ AliESDfriend *esdFriend=static_cast<AliESDfriend*>(esd->FindListObject("AliESDfriend"));
+ if (!esdFriend) {
+ Printf("ERROR: esdFriend not available");
+ return;
+ }
+ //
+
+ for (Int_t ivertex=0; ivertex<nV0s; ivertex++){
+ Bool_t isOK=kFALSE;
+ AliESDv0 * v0 = (AliESDv0*) esd->GetV0(ivertex);
+ AliESDtrack * track0 = fESD->GetTrack(v0->GetIndex(0)); // negative track
+ AliESDtrack * track1 = fESD->GetTrack(v0->GetIndex(1)); // positive track
+ if (track0->GetTPCNcls()<kMinCluster) continue;
+ if (track0->GetKinkIndex(0)>0) continue;
+ if (track1->GetTPCNcls()<kMinCluster) continue;
+ if (track1->GetKinkIndex(0)>0) continue;
+ if (v0->GetOnFlyStatus()==kFALSE) continue;
+ //
+ if (TMath::Min(track0->Pt(),track1->Pt())<kMinMinPt) continue;
+ //
+ //
+ if (TMath::Max(track0->Pt(),track1->Pt())>kMinPt) isOK=kTRUE;
+ if (gRandom->Rndm()<kDownscale) isOK=kTRUE;
+ if (!isOK) continue;
+ //
+ AliESDfriendTrack *ftrack0 = esdFriend->GetTrack(v0->GetIndex(0));
+ if (!ftrack0) continue;
+ AliESDfriendTrack *ftrack1 = esdFriend->GetTrack(v0->GetIndex(1));
+ if (!ftrack1) continue;
+ //
+ TObject *calibObject;
+ AliTPCseed *seed0 = 0;
+ AliTPCseed *seed1 = 0;
+ for (Int_t l=0;(calibObject=ftrack0->GetCalibObject(l));++l) {
+ if ((seed0=dynamic_cast<AliTPCseed*>(calibObject))) break;
+ }
+ for (Int_t l=0;(calibObject=ftrack1->GetCalibObject(l));++l) {
+ if ((seed1=dynamic_cast<AliTPCseed*>(calibObject))) break;
+ }
+ if (!seed0) continue;
+ if (!seed1) continue;
+ AliExternalTrackParam * paramIn0 = (AliExternalTrackParam *)track0->GetInnerParam();
+ AliExternalTrackParam * paramIn1 = (AliExternalTrackParam *)track1->GetInnerParam();
+ if (!paramIn0) continue;
+ if (!paramIn1) continue;
+ //
+ //
+ if (!fV0Tree) {
+ fV0Tree = new TTree("V0s","V0s");
+ fV0Tree->SetDirectory(0);
+ }
+ if (fV0Tree->GetEntries()==0){
+ //
+ fV0Tree->SetDirectory(0);
+ fV0Tree->Branch("v0.",&v0);
+ fV0Tree->Branch("t0.",&track0);
+ fV0Tree->Branch("t1.",&track1);
+ fV0Tree->Branch("ft0.",&ftrack0);
+ fV0Tree->Branch("ft1.",&ftrack1);
+ fV0Tree->Branch("s0.",&seed0);
+ fV0Tree->Branch("s1.",&seed1);
+ }else{
+ fV0Tree->SetBranchAddress("v0.",&v0);
+ fV0Tree->SetBranchAddress("t0.",&track0);
+ fV0Tree->SetBranchAddress("t1.",&track1);
+ fV0Tree->SetBranchAddress("ft0.",&ftrack0);
+ fV0Tree->SetBranchAddress("ft1.",&ftrack1);
+ fV0Tree->SetBranchAddress("s0.",&seed0);
+ fV0Tree->SetBranchAddress("s1.",&seed1);
+ }
+ fV0Tree->Fill();
+ }
+}
+
+
+Long64_t AliTPCcalibV0::Merge(TCollection *const li) {
+
+ TIterator* iter = li->MakeIterator();
+ AliTPCcalibV0* cal = 0;
+
+ while ((cal = (AliTPCcalibV0*)iter->Next())) {
+ if (cal->fV0Tree){
+ if (!fV0Tree) {
+ fV0Tree = new TTree("V0s","V0s");
+ fV0Tree->SetDirectory(0);
+ }
+ if (cal->fV0Tree->GetEntries()>0) AliTPCcalibV0::AddTree(cal->fV0Tree);
+ if (cal->fHPTTree->GetEntries()>0) AliTPCcalibV0::AddTreeHPT(cal->fHPTTree);
+ }
+ }
+ return 0;
+}
+
+
+void AliTPCcalibV0::AddTree(TTree * treeInput){
+ //
+ // Add the content of tree:
+ // Notice automatic copy of tree in ROOT does not work for such complicated tree
+ //
+ AliESDv0 * v0 = new AliESDv0;
+ Double_t kMinPt=0.8;
+ AliESDtrack * track0 = 0; // negative track
+ AliESDtrack * track1 = 0; // positive track
+ AliESDfriendTrack *ftrack0 = 0;
+ AliESDfriendTrack *ftrack1 = 0;
+ AliTPCseed *seed0 = 0;
+ AliTPCseed *seed1 = 0;
+ treeInput->SetBranchStatus("ft0.",kFALSE);
+ treeInput->SetBranchStatus("ft1.",kFALSE);
+ TDatabasePDG pdg;
+ Double_t massK0= pdg.GetParticle("K0")->Mass();
+ Double_t massLambda= pdg.GetParticle("Lambda0")->Mass();
+
+ Int_t entries= treeInput->GetEntries();
+ for (Int_t i=0; i<entries; i++){
+ treeInput->SetBranchAddress("v0.",&v0);
+ treeInput->SetBranchAddress("t0.",&track0);
+ treeInput->SetBranchAddress("t1.",&track1);
+ treeInput->SetBranchAddress("ft0.",&ftrack0);
+ treeInput->SetBranchAddress("ft1.",&ftrack1);
+ treeInput->SetBranchAddress("s0.",&seed0);
+ treeInput->SetBranchAddress("s1.",&seed1);
+ if (fV0Tree->GetEntries()==0){
+ fV0Tree->SetDirectory(0);
+ fV0Tree->Branch("v0.",&v0);
+ fV0Tree->Branch("t0.",&track0);
+ fV0Tree->Branch("t1.",&track1);
+ fV0Tree->Branch("ft0.",&ftrack0);
+ fV0Tree->Branch("ft1.",&ftrack1);
+ fV0Tree->Branch("s0.",&seed0);
+ fV0Tree->Branch("s1.",&seed1);
+ }else{
+ fV0Tree->SetBranchAddress("v0.",&v0);
+ fV0Tree->SetBranchAddress("t0.",&track0);
+ fV0Tree->SetBranchAddress("t1.",&track1);
+ fV0Tree->SetBranchAddress("ft0.",&ftrack0);
+ fV0Tree->SetBranchAddress("ft1.",&ftrack1);
+ fV0Tree->SetBranchAddress("s0.",&seed0);
+ fV0Tree->SetBranchAddress("s1.",&seed1);
+ }
+ //
+ treeInput->GetEntry(i);
+ ftrack0->GetCalibContainer()->SetOwner(kTRUE);
+ ftrack1->GetCalibContainer()->SetOwner(kTRUE);
+ Bool_t isOK=kTRUE;
+ if (v0->GetOnFlyStatus()==kFALSE) isOK=kFALSE;
+ if (track0->GetTPCncls()<100) isOK=kFALSE;
+ if (track1->GetTPCncls()<100) isOK=kFALSE;
+ if (TMath::Min(seed0->Pt(),seed1->Pt())<kMinPt) isOK=kFALSE;
+ if (TMath::Min(track0->Pt(),track1->Pt())<kMinPt) isOK=kFALSE;
+ Bool_t isV0=kFALSE;
+ if (TMath::Abs(v0->GetEffMass(2,2)-massK0)<0.05) isV0=kTRUE;
+ if (TMath::Abs(v0->GetEffMass(4,2)-massLambda)<0.05) isV0=kTRUE;
+ if (TMath::Abs(v0->GetEffMass(2,4)-massLambda)<0.05) isV0=kTRUE;
+ if (TMath::Abs(v0->GetEffMass(0,0))<0.02) isV0=kFALSE; //reject electrons
+ if (!isV0) isOK=kFALSE;
+ if (isOK) fV0Tree->Fill();
+ delete v0;
+ delete track0;
+ delete track1;
+ delete ftrack0;
+ delete ftrack1;
+ delete seed0;
+ delete seed1;
+ v0=0;
+ track0=0;
+ track1=0;
+ ftrack0=0;
+ ftrack1=0;
+ seed0=0;
+ seed1=0;
+ }
+}
+
+void AliTPCcalibV0::AddTreeHPT(TTree * treeInput){
+ //
+ // Add the content of tree:
+ // Notice automatic copy of tree in ROOT does not work for such complicated tree
+ //
+ AliESDtrack *track = 0;
+ AliESDfriendTrack *friendTrack = 0;
+ AliTPCseed *seed = 0;
+ if (!treeInput) return;
+ if (treeInput->GetEntries()==0) return;
+ //
+ Int_t entries= treeInput->GetEntries();
+ //
+ for (Int_t i=0; i<entries; i++){
+ track=0;
+ friendTrack=0;
+ seed=0;
+ //
+ treeInput->SetBranchAddress("t.",&track);
+ treeInput->SetBranchAddress("ft.",&friendTrack);
+ treeInput->SetBranchAddress("s.",&seed);
+ treeInput->GetEntry(i);
+ //
+ if (fHPTTree->GetEntries()==0){
+ fHPTTree->SetDirectory(0);
+ fHPTTree->Branch("t.",&track);
+ fHPTTree->Branch("ft.",&friendTrack);
+ fHPTTree->Branch("s.",&seed);
+ }else{
+ fHPTTree->SetBranchAddress("t.",&track);
+ fHPTTree->SetBranchAddress("ft.",&friendTrack);
+ fHPTTree->SetBranchAddress("s.",&seed);
+ }
+ Bool_t isOK=kTRUE;
+ if (!track->IsOn(AliESDtrack::kITSrefit)) isOK=kFALSE;
+ if (!track->IsOn(AliESDtrack::kTOFout)) isOK=kFALSE;
+ if (isOK) fHPTTree->Fill();
+ //
+ delete track;
+ delete friendTrack;
+ delete seed;
+ }
+}
+
+
void AliTPCcalibV0::MakeMC(){
//
// MC comparison
}
}
+void AliTPCcalibV0::MakeFitTreeTrack(const TObjArray * corrArray, Double_t ptCut, Int_t run){
+ //
+ // Make a fit tree
+ //
+ // 0. Loop over selected tracks
+ // 1. Loop over all transformation - refit the track with and without the
+ // transformtation
+ // 2. Dump the matching paramaeters to the debugStremer
+ //
+
+ //Connect input
+ const Int_t kMinNcl=120;
+ TFile f("TPCV0Objects.root");
+ AliTPCcalibV0 *v0TPC = (AliTPCcalibV0*) f.Get("v0TPC");
+ TTree * treeInput = v0TPC->GetHPTTree();
+ TTreeSRedirector *pcstream = new TTreeSRedirector("fitHPT.root");
+ AliESDtrack *track = 0;
+ AliESDfriendTrack *friendTrack = 0;
+ AliTPCseed *seed = 0;
+ if (!treeInput) return;
+ if (treeInput->GetEntries()==0) return;
+ //
+ treeInput->SetBranchAddress("t.",&track);
+ treeInput->SetBranchAddress("ft.",&friendTrack);
+ treeInput->SetBranchAddress("s.",&seed);
+ //
+ Int_t ncorr=0;
+ if (corrArray) ncorr = corrArray->GetEntries();
+ AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
+ AliTPCParam *param = AliTPCcalibDB::Instance()->GetParameters();
+ AliGRPObject* grp = AliTPCcalibDB::Instance()->GetGRP(run);
+ Double_t time=0.5*(grp->GetTimeStart() +grp->GetTimeEnd());
+ //
+ //
+ //
+ Int_t ntracks= treeInput->GetEntries();
+ for (Int_t itrack=0; itrack<ntracks; itrack++){
+ treeInput->GetEntry(itrack);
+ if (!track) continue;
+ if (seed->Pt()<ptCut) continue;
+ if (track->Pt()<ptCut) continue;
+ if (track->GetTPCncls()<kMinNcl) continue;
+ //
+ // Reapply transformation
+ //
+ for (Int_t irow=0; irow<159; irow++){
+ AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
+ if (cluster &&cluster->GetX()>10){
+ Double_t x0[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
+ Int_t index0[1]={cluster->GetDetector()};
+ transform->Transform(x0,index0,0,1);
+ cluster->SetX(x0[0]);
+ cluster->SetY(x0[1]);
+ cluster->SetZ(x0[2]);
+ //
+ }
+ }
+ //
+ Double_t xref=134;
+ AliExternalTrackParam* paramInner=0;
+ AliExternalTrackParam* paramOuter=0;
+ AliExternalTrackParam* paramIO=0;
+ Bool_t isOK=kTRUE;
+ for (Int_t icorr=-1; icorr<ncorr; icorr++){
+ //
+ AliTPCCorrection *corr = 0;
+ if (icorr>=0) corr = (AliTPCCorrection*)corrArray->At(icorr);
+ AliExternalTrackParam * trackInner = RefitTrack(seed, corr,85,134,0.1);
+ AliExternalTrackParam * trackIO = RefitTrack(seed, corr,245,85,0.1);
+ AliExternalTrackParam * trackOuter = RefitTrack(seed, corr,245,134,0.1 );
+ if (trackInner&&trackOuter&&trackIO){
+ trackOuter->Rotate(trackInner->GetAlpha());
+ trackOuter->PropagateTo(trackInner->GetX(),AliTracker::GetBz());
+ if (icorr<0) {
+ paramInner=trackInner;
+ paramOuter=trackOuter;
+ paramIO=trackIO;
+ paramIO->Rotate(seed->GetAlpha());
+ paramIO->PropagateTo(seed->GetX(),AliTracker::GetBz());
+ }
+ }else{
+ isOK=kFALSE;
+ }
+
+ }
+ if (paramOuter&& paramInner) {
+ // Bool_t isOK=kTRUE;
+ if (paramInner->GetSigmaY2()>0.01) isOK&=kFALSE;
+ if (paramOuter->GetSigmaY2()>0.01) isOK&=kFALSE;
+ if (paramInner->GetSigmaZ2()>0.01) isOK&=kFALSE;
+ if (paramOuter->GetSigmaZ2()>0.01) isOK&=kFALSE;
+ (*pcstream)<<"fit"<<
+ "s.="<<seed<<
+ "io.="<<paramIO<<
+ "pIn.="<<paramInner<<
+ "pOut.="<<paramOuter;
+ }
+ //
+ }
+ delete pcstream;
+ /*
+ .x ~/rootlogon.C
+ Int_t run=117112;
+ .x ../ConfigCalibTrain.C(run)
+ .L ../AddTaskTPCCalib.C
+ ConfigOCDB(run)
+ TFile fexb("../../RegisterCorrectionExB.root");
+ AliTPCComposedCorrection *cc= (AliTPCComposedCorrection*) fexb.Get("ComposedExB");
+ cc->Init();
+ cc->Print("DA"); // Print used correction classes
+ TObjArray *array = cc->GetCorrections()
+ AliTPCcalibV0::MakeFitTreeTrack(array,2,run);
+
+ */
+}
+void AliTPCcalibV0::MakeFitTreeV0(const TObjArray * corrArray, Double_t ptCut, Int_t run){
+ //
+ // Make a fit tree
+ //
+ // 0. Loop over selected tracks
+ // 1. Loop over all transformation - refit the track with and without the
+ // transformtation
+ // 2. Dump the matching paramaeters to the debugStremer
+ //
+
+ //Connect input
+ const Int_t kMinNcl=120;
+ TFile f("TPCV0Objects.root");
+ AliTPCcalibV0 *v0TPC = (AliTPCcalibV0*) f.Get("v0TPC");
+ TTree * treeInput = v0TPC->GetV0Tree();
+ TTreeSRedirector *pcstream = new TTreeSRedirector("fitV0.root");
+ AliESDv0 *v0 = 0;
+ AliESDtrack *track0 = 0;
+ AliESDfriendTrack *friendTrack0 = 0;
+ AliTPCseed *seed0 = 0;
+ AliTPCseed *s0 = 0;
+ AliESDtrack *track1 = 0;
+ AliESDfriendTrack *friendTrack1 = 0;
+ AliTPCseed *s1 = 0;
+ AliTPCseed *seed1 = 0;
+ if (!treeInput) return;
+ if (treeInput->GetEntries()==0) return;
+ //
+ treeInput->SetBranchAddress("v0.",&v0);
+ treeInput->SetBranchAddress("t0.",&track0);
+ treeInput->SetBranchAddress("ft0.",&friendTrack0);
+ treeInput->SetBranchAddress("s0.",&s0);
+ treeInput->SetBranchAddress("t1.",&track1);
+ treeInput->SetBranchAddress("ft1.",&friendTrack1);
+ treeInput->SetBranchAddress("s1.",&s1);
+ //
+ TDatabasePDG pdg;
+ Int_t ncorr=0;
+ if (corrArray) ncorr = corrArray->GetEntries();
+ AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
+ AliTPCParam *param = AliTPCcalibDB::Instance()->GetParameters();
+ AliGRPObject* grp = AliTPCcalibDB::Instance()->GetGRP(run);
+ Double_t time=0.5*(grp->GetTimeStart() +grp->GetTimeEnd());
+ Double_t massK0= pdg.GetParticle("K0")->Mass();
+ Double_t massLambda= pdg.GetParticle("Lambda0")->Mass();
+ Double_t massPion=pdg.GetParticle("pi+")->Mass();
+ Double_t massProton=pdg.GetParticle("proton")->Mass();
+ Double_t pdgPion=pdg.GetParticle("pi+")->PdgCode();
+ Double_t pdgProton=pdg.GetParticle("proton")->PdgCode();
+ Double_t mass0=0;
+ Double_t mass1=0;
+ Double_t massV0=0;
+ Int_t pdg0=0;
+ Int_t pdg1=0;
+ //
+ //
+ //
+ Int_t nv0s= treeInput->GetEntries();
+ for (Int_t iv0=0; iv0<nv0s; iv0++){
+ Int_t v0Type=0;
+ Int_t isK0=0;
+ Int_t isLambda=0;
+ Int_t isAntiLambda=0;
+ treeInput->GetEntry(iv0);
+ if (TMath::Abs(v0->GetEffMass(2,2)-massK0)<0.03) {isK0=1; v0Type=1;} //select K0s
+ if (TMath::Abs(v0->GetEffMass(4,2)-massLambda)<0.01) {isLambda=1; v0Type=2;} //select Lambda
+ if (TMath::Abs(v0->GetEffMass(2,4)-massLambda)<0.01) {isAntiLambda=1;v0Type=3;} //select Anti Lambda
+ if (isK0+isLambda+isAntiLambda!=1) continue;
+ mass0=massPion;
+ mass1=massPion;
+ pdg0=pdgPion;
+ pdg1=pdgPion;
+ if (isLambda) {mass0=massProton; pdg0=pdgProton;}
+ if (isAntiLambda) {mass1=massProton; pdg1=pdgProton;}
+ massV0=massK0;
+ if (isK0==0) massV0=massLambda;
+ //
+ if (!s0) continue;
+ seed0=(s0->GetSign()>0)?s0:s1;
+ seed1=(s0->GetSign()>0)?s1:s0;
+ if (seed0->GetZ()*seed1->GetZ()<0) continue; //remove membrane crossed tracks
+ if (seed0->Pt()<ptCut) continue;
+ if (seed1->Pt()<ptCut) continue;
+ //
+ // Reapply transformation
+ //
+ for (Int_t itype=0; itype<2; itype++){
+ AliTPCseed * seed = (itype==0) ? seed0: seed1;
+ for (Int_t irow=0; irow<159; irow++){
+ AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
+ if (cluster &&cluster->GetX()>10){
+ Double_t x0[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
+ Int_t index0[1]={cluster->GetDetector()};
+ transform->Transform(x0,index0,0,1);
+ cluster->SetX(x0[0]);
+ cluster->SetY(x0[1]);
+ cluster->SetZ(x0[2]);
+ //
+ }
+ }
+ }
+ Bool_t isOK=kTRUE;
+ Double_t radius = v0->GetRr();
+ Double_t xyz[3];
+ v0->GetXYZ(xyz[0],xyz[1],xyz[2]);
+ Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
+ TObjArray arrayV0in(ncorr+1);
+ TObjArray arrayV0io(ncorr+1);
+ TObjArray arrayT0(ncorr+1);
+ TObjArray arrayT1(ncorr+1);
+ arrayV0in.SetOwner(kTRUE);
+ arrayV0io.SetOwner(kTRUE);
+ //
+ for (Int_t icorr=-1; icorr<ncorr; icorr++){
+ AliTPCCorrection *corr =0;
+ if (icorr>=0) corr = (AliTPCCorrection*)corrArray->At(icorr);
+ //
+ AliExternalTrackParam * trackInner0 = RefitTrack(seed0, corr,160,85,mass0);
+ AliExternalTrackParam * trackIO0 = RefitTrack(seed0, corr,245,85,mass0);
+ AliExternalTrackParam * trackInner1 = RefitTrack(seed1, corr,160,85,mass1);
+ AliExternalTrackParam * trackIO1 = RefitTrack(seed1, corr,245,85,mass1);
+ if (!trackInner0) isOK=kFALSE;
+ if (!trackInner1) isOK=kFALSE;
+ if (!trackIO0) isOK=kFALSE;
+ if (!trackIO1) isOK=kFALSE;
+ if (isOK){
+ if (!trackInner0->Rotate(alpha)) isOK=kFALSE;
+ if (!trackInner1->Rotate(alpha)) isOK=kFALSE;
+ if (!trackIO0->Rotate(alpha)) isOK=kFALSE;
+ if (!trackIO1->Rotate(alpha)) isOK=kFALSE;
+ //
+ if (!AliTracker::PropagateTrackToBxByBz(trackInner0, radius, mass0, 1, kFALSE)) isOK=kFALSE;
+ if (!AliTracker::PropagateTrackToBxByBz(trackInner1, radius, mass1, 1, kFALSE)) isOK=kFALSE;
+ if (!AliTracker::PropagateTrackToBxByBz(trackIO0, radius, mass0, 1, kFALSE)) isOK=kFALSE;
+ if (!AliTracker::PropagateTrackToBxByBz(trackIO1, radius, mass1, 1, kFALSE)) isOK=kFALSE;
+ if (!isOK) continue;
+ arrayT0.AddAt(trackIO0->Clone(),icorr+1);
+ arrayT1.AddAt(trackIO1->Clone(),icorr+1);
+ Int_t charge=(trackIO0->GetSign());
+ AliKFParticle pin0( *trackInner0, pdg0*charge);
+ AliKFParticle pin1( *trackInner1, -pdg1*charge);
+ AliKFParticle pio0( *trackIO0, pdg0*charge);
+ AliKFParticle pio1( *trackIO1, -pdg1*charge);
+ AliKFParticle v0in;
+ AliKFParticle v0io;
+ v0in+=pin0;
+ v0in+=pin1;
+ v0io+=pio0;
+ v0io+=pio1;
+ arrayV0in.AddAt(v0in.Clone(),icorr+1);
+ arrayV0io.AddAt(v0io.Clone(),icorr+1);
+ }
+ }
+ if (!isOK) continue;
+ //
+ //AliKFParticle* pin0= (AliKFParticle*)arrayV0in.At(0);
+ AliKFParticle* pio0= (AliKFParticle*)arrayV0io.At(0);
+ AliExternalTrackParam *param0=(AliExternalTrackParam *)arrayT0.At(0);
+ AliExternalTrackParam *param1=(AliExternalTrackParam *)arrayT1.At(0);
+ Double_t mass0=0, mass0E=0;
+ pio0->GetMass( mass0,mass0E);
+ //
+ Double_t mean=mass0-massV0;
+ if (TMath::Abs(mean)>0.05) continue;
+ Double_t mass[10000];
+ //
+ Int_t dtype=30; // id for V0
+ Int_t ptype=5; // id for invariant mass
+ // Int_t id=TMath::Nint(100.*(param0->Pt()-param1->Pt())/(param0->Pt()+param1->Pt())); // K0s V0 asymetry
+ Int_t id=1000.*(param0->Pt()-param1->Pt()); // K0s V0 asymetry
+ Double_t gx,gy,gz, px,py,pz;
+ Double_t pt = v0->Pt();
+ v0->GetXYZ(gx,gy,gz);
+ v0->GetPxPyPz(px,py,pz);
+ Double_t theta=pz/TMath::Sqrt(px*px+py*py);
+ Double_t phi=TMath::ATan2(py,px);
+ Double_t snp=0.5*(seed0->GetSnp()+seed1->GetSnp());
+ Double_t sector=9*phi/TMath::Pi();
+ Double_t dsec=sector-TMath::Nint(sector);
+ Double_t refX=TMath::Sqrt(gx*gx+gy*gy);
+ //Int_t nentries=v0Type;
+ Double_t bz=AliTracker::GetBz();
+ Double_t dRrec=0;
+ (*pcstream)<<"fitDebug"<<
+ "id="<<id<<
+ "v0.="<<v0<<
+ "mean="<<mean<<
+ "rms="<<mass0E<<
+ "pio0.="<<pio0<<
+ "p0.="<<param0<<
+ "p1.="<<param1;
+ (*pcstream)<<"fit"<< // dump valus for fit
+ "run="<<run<< //run number
+ "bz="<<bz<< // magnetic filed used
+ "dtype="<<dtype<< // detector match type 30
+ "ptype="<<ptype<< // parameter type
+ "theta="<<theta<< // theta
+ "phi="<<phi<< // phi
+ "snp="<<snp<< // snp
+ "mean="<<mean<< // mean dist value
+ "rms="<<mass0E<< // rms
+ "sector="<<sector<<
+ "dsec="<<dsec<<
+ //
+ "refX="<<refX<< // reference radius
+ "gx="<<gx<< // global position
+ "gy="<<gy<< // global position
+ "gz="<<gz<< // global position
+ "dRrec="<<dRrec<< // delta Radius in reconstruction
+ "pt="<<pt<< // pt of the particle
+ "id="<<id<< //delta of the momenta
+ "entries="<<v0Type;// type of the V0
+ for (Int_t icorr=0; icorr<ncorr; icorr++){
+ AliTPCCorrection *corr =0;
+ if (icorr>=0) corr = (AliTPCCorrection*)corrArray->At(icorr);
+ AliKFParticle* pin= (AliKFParticle*)arrayV0in.At(icorr+1);
+ AliKFParticle* pio= (AliKFParticle*)arrayV0io.At(icorr+1);
+ AliExternalTrackParam *par0=(AliExternalTrackParam *)arrayT0.At(icorr+1);
+ AliExternalTrackParam *par1=(AliExternalTrackParam *)arrayT1.At(icorr+1);
+ Double_t massE=0;
+ pio->GetMass( mass[icorr],massE);
+ mass[icorr]-=mass0;
+ (*pcstream)<<"fit"<<
+ Form("%s=",corr->GetName())<<mass[icorr];
+ (*pcstream)<<"fitDebug"<<
+ Form("%s=",corr->GetName())<<mass[icorr]<<
+ Form("%sp0.=",corr->GetName())<<par0<<
+ Form("%sp1=",corr->GetName())<<par1;
+ }
+ (*pcstream)<<"fit"<< "isOK="<<isOK<<"\n";
+ (*pcstream)<<"fitDebug"<< "isOK="<<isOK<<"\n";
+ }
+ delete pcstream;
+}
-void AliTPCcalibV0::MakeV0s(){
+AliExternalTrackParam * AliTPCcalibV0::RefitTrack(AliTPCseed *seed, AliTPCCorrection * corr, Double_t xstart, Double_t xstop, Double_t mass){
//
+ // Refit the track:
+ // seed - tpc track with cluster
+ // corr - distrotion/correction class - apllied to the points
+ // xstart - radius to start propagate/update
+ // xstop - radius to stop propagate/update
+ //
+ const Double_t kResetCov=20.;
+ const Double_t kSigma=5.;
+ Double_t covar[15];
+ for (Int_t i=0;i<15;i++) covar[i]=0;
+ covar[0]=kSigma*kSigma;
+ covar[2]=kSigma*kSigma;
+ covar[5]=kSigma*kSigma/Float_t(150*150);
+ covar[9]=kSigma*kSigma/Float_t(150*150);
+ covar[14]=1*1;
+ //
+ AliExternalTrackParam *refit = new AliExternalTrackParam(*seed);
+ refit->PropagateTo(xstart, AliTracker::GetBz());
+ refit->AddCovariance(covar);
+ refit->ResetCovariance(kResetCov);
+ Double_t xmin = TMath::Min(xstart,xstop);
+ Double_t xmax = TMath::Max(xstart,xstop);
+ Int_t ncl=0;
//
+ Bool_t isOK=kTRUE;
+ for (Int_t index0=0; index0<159; index0++){
+ Int_t irow= (xstart<xstop)? index0:159-index0;
+ AliTPCclusterMI *cluster=seed->GetClusterPointer(irow); //cluster in local system
+ if (!cluster) continue;
+ if (cluster->GetX()<xmin) continue;
+ if (cluster->GetX()>xmax) continue;
+ Double_t alpha = TMath::Pi()*(cluster->GetDetector()+0.5)/9.;
+ if (!refit->Rotate(alpha)) isOK=kFALSE;
+ Double_t x = cluster->GetX();
+ Double_t y = cluster->GetY();
+ Double_t z = cluster->GetZ();
+ if (corr){
+ Float_t xyz[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()}; // original position
+ corr->DistortPointLocal(xyz,cluster->GetDetector());
+ x=xyz[0];
+ y=xyz[1];
+ z=xyz[2];
+ }
+ if (!AliTracker::PropagateTrackToBxByBz(refit, x,mass,1.,kFALSE)) isOK=kFALSE;
+ if (!isOK) continue;
+ Double_t cov[3]={0.01,0.,0.01};
+ Double_t yz[2]={y,z};
+ if (!refit->Update(yz,cov)) isOK=kFALSE;
+ ncl++;
+ }
+ if (!AliTracker::PropagateTrackToBxByBz(refit, xstop, mass,1.,kTRUE)) isOK=kFALSE;
+ //
+ if (!isOK) {
+ delete refit;
+ return 0;
+ }
+ return refit;
+}
+
+
+
+
+
+//
+// Obsolete part
+//
+
+
+void AliTPCcalibV0::MakeV0s(){
+ //
+ //
//
const Int_t kMinCluster=40;
const Float_t kMinR =0;
-// void AliTPCcalibV0::ProcessV0(Int_t ftype){
-// //
-// //
-// const Double_t ktimeK0 = 2.684;
-// const Double_t ktimeLambda = 7.89;
-
-
-// if (! fGammas) fGammas = new TObjArray(10);
-// fGammas->Clear();
-// Int_t nV0s = fV0s->GetEntries();
-// if (nV0s==0) return;
-// AliKFVertex primVtx(*(fESD->GetPrimaryVertex()));
-// //
-// for (Int_t ivertex=0; ivertex<nV0s; ivertex++){
-// AliESDv0 * v0 = (AliESDv0*)fV0s->At(ivertex);
-// AliESDtrack * trackN = fESD->GetTrack(v0->GetIndex(0));
-// AliESDtrack * trackP = fESD->GetTrack(v0->GetIndex(1));
-// //
-// //
-// //
-// AliKFParticle *v0K0 = Fit(primVtx,v0,211,211);
-// AliKFParticle *v0Gamma = Fit(primVtx,v0,11,-11);
-// AliKFParticle *v0Lambda42 = Fit(primVtx,v0,2212,211);
-// AliKFParticle *v0Lambda24 = Fit(primVtx,v0,211,2212);
-// //Set production vertex
-// v0K0->SetProductionVertex( primVtx );
-// v0Gamma->SetProductionVertex( primVtx );
-// v0Lambda42->SetProductionVertex( primVtx );
-// v0Lambda24->SetProductionVertex( primVtx );
-// Double_t massK0, massGamma, massLambda42,massLambda24, massSigma;
-// v0K0->GetMass( massK0,massSigma);
-// v0Gamma->GetMass( massGamma,massSigma);
-// v0Lambda42->GetMass( massLambda42,massSigma);
-// v0Lambda24->GetMass( massLambda24,massSigma);
-// Float_t chi2K0 = v0K0->GetChi2()/v0K0->GetNDF();
-// Float_t chi2Gamma = v0Gamma->GetChi2()/v0Gamma->GetNDF();
-// Float_t chi2Lambda42 = v0Lambda42->GetChi2()/v0Lambda42->GetNDF();
-// Float_t chi2Lambda24 = v0Lambda24->GetChi2()/v0Lambda24->GetNDF();
-// //
-// // Mass Contrained params
-// //
-// AliKFParticle *v0K0C = Fit(primVtx,v0,211,211);
-// AliKFParticle *v0GammaC = Fit(primVtx,v0,11,-11);
-// AliKFParticle *v0Lambda42C = Fit(primVtx,v0,2212,211);
-// AliKFParticle *v0Lambda24C = Fit(primVtx,v0,211,2212);
-// //
-// v0K0C->SetProductionVertex( primVtx );
-// v0GammaC->SetProductionVertex( primVtx );
-// v0Lambda42C->SetProductionVertex( primVtx );
-// v0Lambda24C->SetProductionVertex( primVtx );
-
-// v0K0C->SetMassConstraint(fPdg->GetParticle(310)->Mass());
-// v0GammaC->SetMassConstraint(0);
-// v0Lambda42C->SetMassConstraint(fPdg->GetParticle(3122)->Mass());
-// v0Lambda24C->SetMassConstraint(fPdg->GetParticle(-3122)->Mass());
-// //
-// Double_t timeK0, sigmaTimeK0;
-// Double_t timeLambda42, sigmaTimeLambda42;
-// Double_t timeLambda24, sigmaTimeLambda24;
-// v0K0C->GetLifeTime(timeK0, sigmaTimeK0);
-// //v0K0Gamma->GetLifeTime(timeK0, sigmaTimeK0);
-// v0Lambda42C->GetLifeTime(timeLambda42, sigmaTimeLambda42);
-// v0Lambda24C->GetLifeTime(timeLambda24, sigmaTimeLambda24);
-
-
-// //
-// Float_t chi2K0C = v0K0C->GetChi2()/v0K0C->GetNDF();
-// if (chi2K0C<0) chi2K0C=100;
-// Float_t chi2GammaC = v0GammaC->GetChi2()/v0GammaC->GetNDF();
-// if (chi2GammaC<0) chi2GammaC=100;
-// Float_t chi2Lambda42C = v0Lambda42C->GetChi2()/v0Lambda42C->GetNDF();
-// if (chi2Lambda42C<0) chi2Lambda42C=100;
-// Float_t chi2Lambda24C = v0Lambda24C->GetChi2()/v0Lambda24C->GetNDF();
-// if (chi2Lambda24C<0) chi2Lambda24C=100;
-// //
-// Float_t minChi2C=99;
-// Int_t type =-1;
-// if (chi2K0C<minChi2C) { minChi2C= chi2K0C; type=0;}
-// if (chi2GammaC<minChi2C) { minChi2C= chi2GammaC; type=1;}
-// if (chi2Lambda42C<minChi2C) { minChi2C= chi2Lambda42C; type=2;}
-// if (chi2Lambda24C<minChi2C) { minChi2C= chi2Lambda24C; type=3;}
-// Float_t minChi2=99;
-// Int_t type0 =-1;
-// if (chi2K0<minChi2) { minChi2= chi2K0; type0=0;}
-// if (chi2Gamma<minChi2) { minChi2= chi2Gamma; type0=1;}
-// if (chi2Lambda42<minChi2) { minChi2= chi2Lambda42; type0=2;}
-// if (chi2Lambda24<minChi2) { minChi2= chi2Lambda24; type0=3;}
-// Float_t betaGammaP = trackN->GetP()/fPdg->GetParticle(-2212)->Mass();
-// Float_t betaGammaPi = trackN->GetP()/fPdg->GetParticle(-211)->Mass();
-// Float_t betaGammaEl = trackN->GetP()/fPdg->GetParticle(11)->Mass();
-// Float_t dedxTeorP = BetheBlochAleph(betaGammaP);
-// Float_t dedxTeorPi = BetheBlochAleph(betaGammaPi);;
-// Float_t dedxTeorEl = BetheBlochAleph(betaGammaEl);;
-// //
-// //
-// if (minChi2>50) continue;
-// (*fDebugStream)<<"V0"<<
-// "ftype="<<ftype<<
-// "v0.="<<v0<<
-// "trackN.="<<trackN<<
-// "trackP.="<<trackP<<
-// //
-// "dedxTeorP="<<dedxTeorP<<
-// "dedxTeorPi="<<dedxTeorPi<<
-// "dedxTeorEl="<<dedxTeorEl<<
-// //
-// "type="<<type<<
-// "chi2C="<<minChi2C<<
-// "v0K0.="<<v0K0<<
-// "v0Gamma.="<<v0Gamma<<
-// "v0Lambda42.="<<v0Lambda42<<
-// "v0Lambda24.="<<v0Lambda24<<
-// //
-// "chi20K0.="<<chi2K0<<
-// "chi2Gamma.="<<chi2Gamma<<
-// "chi2Lambda42.="<<chi2Lambda42<<
-// "chi2Lambda24.="<<chi2Lambda24<<
-// //
-// "chi20K0c.="<<chi2K0C<<
-// "chi2Gammac.="<<chi2GammaC<<
-// "chi2Lambda42c.="<<chi2Lambda42C<<
-// "chi2Lambda24c.="<<chi2Lambda24C<<
-// //
-// "v0K0C.="<<v0K0C<<
-// "v0GammaC.="<<v0GammaC<<
-// "v0Lambda42C.="<<v0Lambda42C<<
-// "v0Lambda24C.="<<v0Lambda24C<<
-// //
-// "massK0="<<massK0<<
-// "massGamma="<<massGamma<<
-// "massLambda42="<<massLambda42<<
-// "massLambda24="<<massLambda24<<
-// //
-// "timeK0="<<timeK0<<
-// "timeLambda42="<<timeLambda42<<
-// "timeLambda24="<<timeLambda24<<
-// "\n";
-// if (type==1) fGammas->AddLast(v0);
-// //
-// //
-// //
-// delete v0K0;
-// delete v0Gamma;
-// delete v0Lambda42;
-// delete v0Lambda24;
-// delete v0K0C;
-// delete v0GammaC;
-// delete v0Lambda42C;
-// delete v0Lambda24C;
-// }
-// ProcessPI0();
-// }
-
-
void AliTPCcalibV0::ProcessV0(Int_t ftype){
//
- //
-
+ // Obsolete
+ //
if (! fGammas) fGammas = new TObjArray(10);
fGammas->Clear();
Int_t nV0s = fV0s->GetEntries();
if (chi2Lambda42C < 10*minChi2C) numCand++;
if (chi2Lambda24C < 10*minChi2C) numCand++;
//
- if (numCand < 2) {
- if (paramInNeg->GetP() > 0.4) fTPCdEdx->Fill(betaGamma0, trackN->GetTPCsignal());
- if (paramInPos->GetP() > 0.4) fTPCdEdx->Fill(betaGamma1, trackP->GetTPCsignal());
- }
}
//
return V0;
}
-TH2F * AliTPCcalibV0::GetHistograms() {
- //
- //
- //
- return fTPCdEdx;
-}
-
class TArrayI;
class TTree;
class AliStack;
+class AliExternalTrackParam;
+class AliTPCCorrection;
class AliTPCcalibV0 : public AliTPCcalibBase {
public :
// List of branches
AliTPCcalibV0();
+ AliTPCcalibV0(const Text_t *name, const Text_t *title);
virtual ~AliTPCcalibV0();
virtual void Process(AliESDEvent *event) {return ProcessESD(event,0);}
-
+ Long64_t Merge(TCollection *const li);
+ void AddTree(TTree * treeInput);
+ void AddTreeHPT(TTree * treeInput);
+ static AliExternalTrackParam * RefitTrack(AliTPCseed *seed, AliTPCCorrection * corr, Double_t xstart, Double_t xstop, Double_t mass);
+ static void MakeFitTreeTrack(const TObjArray * corrArray, Double_t ptCut, Int_t run);
+ static void MakeFitTreeV0(const TObjArray * corrArray, Double_t ptCut, Int_t run);
//
//
//
void ProcessESD(AliESDEvent *esd, AliStack *stack=0);
+ void DumpToTree(AliESDEvent *esd);
+ void DumpToTreeHPT(AliESDEvent *esd);
+ TTree * GetV0Tree(){return fV0Tree;}
+ TTree * GetHPTTree(){return fHPTTree;}
+ //
void MakeMC();
void MakeV0s();
void ProcessV0(Int_t ftype);
void ProcessPI0();
- TH2F * GetHistograms();
void BinLogX(TH2F * h);
//
//
//
static AliKFParticle * Fit(AliKFVertex & primVtx, AliESDv0 *v0, Int_t PDG1, Int_t PDG2);
- void Process(AliESDtrack *track, Int_t runNo=-1){AliTPCcalibBase::Process(track,runNo);};
- void Process(AliTPCseed *track){return AliTPCcalibBase::Process(track);}
protected:
private:
- AliTPCcalibV0(const AliTPCcalibV0&); // Not implemented
- AliTPCcalibV0& operator=(const AliTPCcalibV0&); // Not implemented
-
-
- AliStack *fStack; // pointer to kinematic tree
- AliESDEvent *fESD; //! current ED to proccess - NOT OWNER
- TDatabasePDG *fPdg; // particle database
- TObjArray *fParticles; // array of selected MC particles
- TObjArray *fV0s; // array of V0s
- TObjArray *fGammas; // gamma conversion candidates
- //
- TArrayI *fV0type; // array of types for V0s
- TH2F *fTPCdEdx; // dEdx spectra
- TH2F *fTPCdEdxPi; // dEdx spectra - pion anti-pion
- TH2F *fTPCdEdxEl; // dEdx spectra - electroen -positrons from gamma
- TH2F *fTPCdEdxP; // dEdx spectra - proton antiproton - lambda - antilambda
- //
- ClassDef(AliTPCcalibV0,1);
+ AliTPCcalibV0(const AliTPCcalibV0&); // Not implemented
+ AliTPCcalibV0& operator=(const AliTPCcalibV0&); // Not implemented
+ TTree *fV0Tree; // tree with full V0 information
+ TTree *fHPTTree; // tree with high mometa tracks - full calib info
+ //
+ AliStack *fStack; // pointer to kinematic tree
+ AliESDEvent *fESD; //! current ED to proccess - NOT OWNER
+ TDatabasePDG *fPdg; //! particle database
+ TObjArray *fParticles; // array of selected MC particles
+ TObjArray *fV0s; // array of V0s
+ TObjArray *fGammas; // gamma conversion candidates
+ //
+ void Process(AliESDtrack *track, Int_t runNo=-1){AliTPCcalibBase::Process(track,runNo);};
+ void Process(AliTPCseed *track){return AliTPCcalibBase::Process(track);}
+ //
+ ClassDef(AliTPCcalibV0,3);
};
chain->SetAlias("La","(ly.fElements/lx.fElements/0.155)");
chain->SetAlias("deltaT","(CETmean.fElements-PulserTmean.fElements)");
}
+
+
+
+
+void AliTPCkalmanAlign::Update1D(Double_t delta, Double_t sigma, Int_t s1, TMatrixD &vecXk, TMatrixD &covXk){
+ //
+ // Update parameters and covariance - with one measurement
+ //
+ const Int_t knMeas=1;
+ Int_t knElem=vecXk.GetNrows();
+
+ TMatrixD mat1(knElem,knElem); // update covariance matrix
+ TMatrixD matHk(1,knElem); // vector to mesurement
+ TMatrixD vecYk(knMeas,1); // Innovation or measurement residual
+ TMatrixD matHkT(knElem,knMeas); // helper matrix Hk transpose
+ TMatrixD matSk(knMeas,knMeas); // Innovation (or residual) covariance
+ TMatrixD matKk(knElem,knMeas); // Optimal Kalman gain
+ TMatrixD covXk2(knElem,knElem); // helper matrix
+ TMatrixD covXk3(knElem,knElem); // helper matrix
+ TMatrixD vecZk(1,1);
+ TMatrixD measR(1,1);
+ vecZk(0,0)=delta;
+ measR(0,0)=sigma*sigma;
+ //
+ // reset matHk
+ for (Int_t iel=0;iel<knElem;iel++)
+ for (Int_t ip=0;ip<knMeas;ip++) matHk(ip,iel)=0;
+ //mat1
+ for (Int_t iel=0;iel<knElem;iel++) {
+ for (Int_t jel=0;jel<knElem;jel++) mat1(iel,jel)=0;
+ mat1(iel,iel)=1;
+ }
+ //
+ matHk(0, s1)=1;
+ vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual
+ matHkT=matHk.T(); matHk.T();
+ matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance
+ matSk.Invert();
+ matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain
+ vecXk += matKk*vecYk; // updated vector
+ covXk2= (mat1-(matKk*matHk));
+ covXk3 = covXk2*covXk;
+ covXk = covXk3;
+ Int_t nrows=covXk3.GetNrows();
+
+ for (Int_t irow=0; irow<nrows; irow++)
+ for (Int_t icol=0; icol<nrows; icol++){
+ // rounding problems - make matrix again symteric
+ covXk(irow,icol)=(covXk3(irow,icol)+covXk3(icol,irow))*0.5;
+ }
+}
+
+
+void AliTPCkalmanAlign::Update1D(TString &input, TString filter, TVectorD ¶m, TMatrixD & covar, Double_t mean, Double_t sigma){
+ //
+ // Update Parameters
+ // input - variable name description
+ // filter - filter string
+ // param - parameter vector
+ // covar - covariance
+ // mean - value to set
+ // sigma - sigma of measurement
+ TObjArray *array0= input.Tokenize("++");
+ TObjArray *array1= filter.Tokenize("++");
+ TMatrixD paramM(param.GetNrows(),1);
+ for (Int_t i=0; i<array0->GetEntries(); i++){paramM(i,0)=param(i);}
+
+ for (Int_t i=0; i<array0->GetEntries(); i++){
+ Bool_t isOK=kTRUE;
+ TString str(array0->At(i)->GetName());
+ for (Int_t j=0; j<array1->GetEntries(); j++){
+ if (str.Contains(array1->At(j)->GetName())==0) isOK=kFALSE;
+ }
+ if (isOK) {
+ AliTPCkalmanAlign::Update1D(mean, sigma, i, paramM, covar);
+ }
+ }
+ for (Int_t i=0; i<array0->GetEntries(); i++){
+ param(i)=paramM(i,0);
+ }
+}
+
+
+TString AliTPCkalmanAlign::FilterFit(TString &input, TString filter, TVectorD ¶m, TMatrixD & covar){
+ //
+ //
+ //
+ TObjArray *array0= input.Tokenize("++");
+ TObjArray *array1= filter.Tokenize("++");
+ //TString *presult=new TString("(0");
+ TString result="(0.0";
+ for (Int_t i=0; i<array0->GetEntries(); i++){
+ Bool_t isOK=kTRUE;
+ TString str(array0->At(i)->GetName());
+ for (Int_t j=0; j<array1->GetEntries(); j++){
+ if (str.Contains(array1->At(j)->GetName())==0) isOK=kFALSE;
+ }
+ if (isOK) {
+ result+="+"+str;
+ result+=Form("*(%f)",param[i+1]);
+ printf("%f\t%f\t%s\n",param[i+1], TMath::Sqrt(covar(i+1,i+1)),str.Data());
+ }
+ }
+ result+="-0.)";
+ return result;
+}
+
+
#include "TNamed.h"
#include "TMatrixD.h"
+#include "TString.h"
class TTreeSRedirector;
class TObjArray;
class AliTPCcalibAlign;
void UpdateOCDBTime0( AliTPCCalPad *pad, Int_t ustartRun, Int_t uendRun, const char* storagePath );
static void UpdateAlign1D(Double_t delta, Double_t sigma, Int_t s1, Int_t s2, TMatrixD ¶m, TMatrixD &covar);
static void UpdateAlign1D(Double_t delta, Double_t sigma, Int_t s1, TMatrixD ¶m, TMatrixD &covar);
+ //
+ static void Update1D(Double_t delta, Double_t sigma, Int_t s1, TMatrixD ¶m, TMatrixD &covar);
+ static void Update1D(TString &input, TString filter, TVectorD ¶m, TMatrixD & covar, Double_t mean, Double_t sigma);
+ static TString FilterFit(TString &input, TString filter, TVectorD ¶m, TMatrixD & covar);
+ //
static void BookAlign1D(TMatrixD ¶m, TMatrixD &covar, Double_t sigma, Double_t mean);
void DumpOldAlignment(TTreeSRedirector *pcstream);
void MakeNewAlignment(Bool_t add,TTreeSRedirector *pcstream=0);
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff