From 46e897933a84394e980cc941064170a9d61dbfac Mon Sep 17 00:00:00 2001 From: marian Date: Tue, 15 Feb 2011 14:45:32 +0000 Subject: [PATCH] AliTPCcalibTimeGain.cxx - Adding the Gamma conversion selected electorns 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 --- TPC/AliTPCCorrection.cxx | 687 ++++++++-- TPC/AliTPCCorrection.h | 6 +- TPC/AliTPCExBBShape.cxx | 25 + TPC/AliTPCExBBShape.h | 1 + TPC/AliTPCExBEffectiveSector.cxx | 27 +- TPC/AliTPCExBEffectiveSector.h | 2 +- TPC/AliTPCcalibTimeGain.cxx | 90 +- TPC/AliTPCcalibTracks.cxx | 2002 ++++++------------------------ TPC/AliTPCcalibTracks.h | 42 +- TPC/AliTPCcalibV0.cxx | 932 +++++++++++--- TPC/AliTPCcalibV0.h | 52 +- TPC/AliTPCkalmanAlign.cxx | 108 ++ TPC/AliTPCkalmanAlign.h | 6 + 13 files changed, 2024 insertions(+), 1956 deletions(-) diff --git a/TPC/AliTPCCorrection.cxx b/TPC/AliTPCCorrection.cxx index cd2cdf5e8cc..a4317b80e37 100644 --- a/TPC/AliTPCCorrection.cxx +++ b/TPC/AliTPCCorrection.cxx @@ -207,6 +207,41 @@ void AliTPCCorrection::DistortPoint(Float_t x[],const Short_t roc) { 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 @@ -1345,6 +1380,7 @@ AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackPara 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(); @@ -1368,7 +1404,8 @@ AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackPara 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())kRTPC1) continue; AliTrackPoint pIn0; // space point AliTrackPoint pIn1; Int_t sector= (xyz[2]>0)? 0:18; @@ -1408,7 +1445,11 @@ AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackPara 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); @@ -1426,6 +1467,8 @@ AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackPara 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())kRTPC1) continue; // if (!AliTrackerBase::PropagateTrackTo(track0,prot0.GetX(),kMass,5,kFALSE,kMaxSnp)) break; if (!AliTrackerBase::PropagateTrackTo(track1,prot0.GetX(),kMass,5,kFALSE,kMaxSnp)) break; @@ -1545,7 +1588,7 @@ TTree* AliTPCCorrection::CreateDistortionTree(Double_t step){ -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) @@ -1571,9 +1614,11 @@ void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t // 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); @@ -1583,7 +1628,7 @@ void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t 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(); @@ -1598,7 +1643,7 @@ void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t if (dtype==3) { dir=1;} if (dtype==4) { dir=-1;} // - for (Int_t ientry=0; ientryGetEntry(ientry); if (TMath::Abs(snp)>kMaxSnp) continue; tPar[0]=0; @@ -1631,9 +1676,12 @@ void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t 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="<_.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; ientryGetEntry(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 (entries1) isOK=kFALSE; + // + Double_t dRrec=0; // dummy value - needed for points - e.g for laser + // + (*pcstream)<<"fit"<< + "run="<GetMatrixArray()[irow]>cutErrY||veceZ->GetMatrixArray()[irow]>cutErrZ) nentries=0; - if (veceY->GetMatrixArray()[irow]GetMatrixArray()[irow]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="<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]GetMatrixArray()[irow]GetVecSec())[irow]>=0 && err) { + for (Int_t jrow=first3; jrow<=last3; jrow++){ + if ((*ltr->GetVecSec())[irow]!= (*ltr->GetVecSec())[jrow]) continue; + if ((*err)[jrow]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="<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; ientryGetEntry(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])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])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])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="<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]; +} diff --git a/TPC/AliTPCExBBShape.h b/TPC/AliTPCExBBShape.h index defcce7fbb8..f2763f8a037 100644 --- a/TPC/AliTPCExBBShape.h +++ b/TPC/AliTPCExBBShape.h @@ -65,6 +65,7 @@ public: 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) diff --git a/TPC/AliTPCExBEffectiveSector.cxx b/TPC/AliTPCExBEffectiveSector.cxx index c3954a492a3..fcd1b9cc59f 100644 --- a/TPC/AliTPCExBEffectiveSector.cxx +++ b/TPC/AliTPCExBEffectiveSector.cxx @@ -190,7 +190,7 @@ void AliTPCExBEffectiveSector::Print(const Option_t* option) const { } } -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: @@ -201,7 +201,7 @@ void AliTPCExBEffectiveSector::MakeResidualMap(THnSparse * hisInput, const char 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(), @@ -314,13 +314,36 @@ void AliTPCExBEffectiveSector::MakeResidualMap(THnSparse * hisInput, const char 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="<GetTrack(i); if (!track) continue; @@ -392,8 +393,16 @@ void AliTPCcalibTimeGain::ProcessBeamEvent(AliESDEvent *event) { // 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; @@ -402,19 +411,70 @@ void AliTPCcalibTimeGain::ProcessBeamEvent(AliESDEvent *event) { if ((seed=dynamic_cast(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(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); + } + } + } } diff --git a/TPC/AliTPCcalibTracks.cxx b/TPC/AliTPCcalibTracks.cxx index 6687d07f1c5..700bfaef615 100644 --- a/TPC/AliTPCcalibTracks.cxx +++ b/TPC/AliTPCcalibTracks.cxx @@ -34,25 +34,6 @@ 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(); -*/ - // /////////////////////////////////////////////////////////////////////////////// @@ -108,7 +89,6 @@ using namespace std; #include "AliTPCcalibTracks.h" #include "AliTPCClusterParam.h" #include "AliTPCcalibTracksCuts.h" -#include "AliTPCCalPadRegion.h" #include "AliTPCCalPad.h" #include "AliTPCCalROC.h" #include "TText.h" @@ -117,6 +97,7 @@ using namespace std; #include "TStatToolkit.h" #include "TCut.h" #include "THnSparse.h" +#include "AliRieman.h" @@ -133,25 +114,16 @@ AliTPCcalibTracks::AliTPCcalibTracks(): 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) @@ -174,25 +146,16 @@ AliTPCcalibTracks::AliTPCcalibTracks(const AliTPCcalibTracks& calibTracks): 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) @@ -206,14 +169,6 @@ AliTPCcalibTracks::AliTPCcalibTracks(const AliTPCcalibTracks& calibTracks): 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); @@ -239,20 +194,9 @@ AliTPCcalibTracks::AliTPCcalibTracks(const AliTPCcalibTracks& calibTracks): 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(); @@ -286,25 +230,16 @@ AliTPCcalibTracks::AliTPCcalibTracks(const Text_t *name, const Text_t *title, Al 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) @@ -339,61 +274,14 @@ AliTPCcalibTracks::AliTPCcalibTracks(const Text_t *name, const Text_t *title, Al 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); @@ -456,19 +344,6 @@ AliTPCcalibTracks::AliTPCcalibTracks(const Text_t *name, const Text_t *title, Al } } - 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; @@ -483,16 +358,7 @@ AliTPCcalibTracks::~AliTPCcalibTracks() { 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(); @@ -518,31 +384,17 @@ AliTPCcalibTracks::~AliTPCcalibTracks() { } } - 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 @@ -664,666 +516,272 @@ void AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){ // 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())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="< 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="<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="< 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="<GetPad()); + if (cstream){ + Float_t zdrift = 250 - TMath::Abs(cluster0->GetZ()); + (*cstream)<<"Resol0"<< + "run="<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 { @@ -1340,67 +798,6 @@ 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){ // @@ -1411,471 +808,10 @@ 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"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; - hisResZ02->Fit(fres, "q"); - fresol<<"Z\t"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; - // - fresol<<"Pad 1.00 cm"<<1<<"\n"; - hisResY12->Fit(fres, "q"); // valgrind - fresol<<"Y\t"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; - hisResZ12->Fit(fres, "q"); - fresol<<"Z\t"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; - // - fresol<<"Pad 1.50 cm"<<0<<"\n"; - hisResY22->Fit(fres, "q"); - fresol<<"Y\t"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; - hisResZ22->Fit(fres, "q"); - fresol<<"Z\t"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<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"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; - hisResZ02->Fit(frms, "qn0"); // valgrind 23 blocks indirectly lost - filerms<<"Z\t"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<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"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; - hisResZ12->Fit(frms, "qn0"); // valgrind 66,036 bytes in 3 blocks are still reachable - filerms<<"Z\t"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<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"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; - hisResZ22->Fit(frms, "qn0"); - filerms<<"Z\t"<GetParameter(0)<<"\t"<GetParameter(1)<<"\t"<GetParameter(2)<<"\n"; - - TH1::AddDirectory(kFALSE); - gSystem->ChangeDirectory(".."); - ps->Close(); - delete ps; - delete frms; -} void AliTPCcalibTracks::MakeResPlotsQTree(Int_t minEntries, const char* pathName){ @@ -2209,17 +1145,12 @@ Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) { 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; @@ -2236,10 +1167,6 @@ Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) { while ( (calibTracks = dynamic_cast (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()); @@ -2248,13 +1175,8 @@ Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) { 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()); @@ -2267,49 +1189,15 @@ Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) { // 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 @@ -2373,52 +1261,7 @@ Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) { 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; - } - } @@ -2504,179 +1347,6 @@ Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) { } -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(){ @@ -2703,46 +1373,61 @@ 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); @@ -2761,3 +1446,128 @@ void AliTPCcalibTracks::AddHistos(AliTPCcalibTracks* calib){ 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)<GetMagneticField()); + if (TMath::Abs(AliTracker::GetBz())<1) return; + DumpToTree(esd); + DumpToTreeHPT(esd); + return; + // MakeV0s(); if (stack) { fStack = stack; @@ -96,6 +114,304 @@ void AliTPCcalibV0::ProcessESD(AliESDEvent *esd, AliStack *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(esd->FindListObject("AliESDfriend")); + if (!esdFriend) { + Printf("ERROR: esdFriend not available"); + return; + } + // + Int_t ntracks=esd->GetNumberOfTracks(); + for (Int_t i=0;iGetTrack(i); + if (track->GetTPCncls()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(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(esd->FindListObject("AliESDfriend")); + if (!esdFriend) { + Printf("ERROR: esdFriend not available"); + return; + } + // + + for (Int_t ivertex=0; ivertexGetV0(ivertex); + AliESDtrack * track0 = fESD->GetTrack(v0->GetIndex(0)); // negative track + AliESDtrack * track1 = fESD->GetTrack(v0->GetIndex(1)); // positive track + if (track0->GetTPCNcls()GetKinkIndex(0)>0) continue; + if (track1->GetTPCNcls()GetKinkIndex(0)>0) continue; + if (v0->GetOnFlyStatus()==kFALSE) continue; + // + if (TMath::Min(track0->Pt(),track1->Pt())Pt(),track1->Pt())>kMinPt) isOK=kTRUE; + if (gRandom->Rndm()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(calibObject))) break; + } + for (Int_t l=0;(calibObject=ftrack1->GetCalibObject(l));++l) { + if ((seed1=dynamic_cast(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; iSetBranchAddress("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())Pt(),track1->Pt())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; iSetBranchAddress("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 @@ -257,11 +573,430 @@ void AliTPCcalibV0::MakeMC(){ } } +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; itrackGetEntry(itrack); + if (!track) continue; + if (seed->Pt()Pt()GetTPCncls()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=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.="<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; iv0GetEntry(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()Pt()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=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="<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= (xstartGetClusterPointer(irow); //cluster in local system + if (!cluster) continue; + if (cluster->GetX()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; @@ -333,166 +1068,12 @@ void AliTPCcalibV0::MakeV0s(){ -// 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; ivertexAt(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 (chi2K0CGetP()/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="<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(); @@ -613,10 +1194,6 @@ void AliTPCcalibV0::ProcessV0(Int_t ftype){ 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()); - } } // @@ -750,13 +1327,6 @@ AliKFParticle * AliTPCcalibV0::Fit(AliKFVertex & /*primVtx*/, AliESDv0 *v0, Int_ return V0; } -TH2F * AliTPCcalibV0::GetHistograms() { - // - // - // - return fTPCdEdx; -} - diff --git a/TPC/AliTPCcalibV0.h b/TPC/AliTPCcalibV0.h index 0bbba452dfa..c4628c179a9 100644 --- a/TPC/AliTPCcalibV0.h +++ b/TPC/AliTPCcalibV0.h @@ -21,6 +21,8 @@ class AliESDv0; class TArrayI; class TTree; class AliStack; +class AliExternalTrackParam; +class AliTPCCorrection; class AliTPCcalibV0 : public AliTPCcalibBase { public : @@ -28,47 +30,53 @@ 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); }; diff --git a/TPC/AliTPCkalmanAlign.cxx b/TPC/AliTPCkalmanAlign.cxx index 513afa63e7a..e6049f05b6b 100644 --- a/TPC/AliTPCkalmanAlign.cxx +++ b/TPC/AliTPCkalmanAlign.cxx @@ -969,3 +969,111 @@ void AliTPCkalmanAlign::MakeAliasCE(TTree * chain){ 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;ielGetEntries(); i++){paramM(i,0)=param(i);} + + for (Int_t i=0; iGetEntries(); i++){ + Bool_t isOK=kTRUE; + TString str(array0->At(i)->GetName()); + for (Int_t j=0; jGetEntries(); 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; iGetEntries(); 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; iGetEntries(); i++){ + Bool_t isOK=kTRUE; + TString str(array0->At(i)->GetName()); + for (Int_t j=0; jGetEntries(); 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; +} + + diff --git a/TPC/AliTPCkalmanAlign.h b/TPC/AliTPCkalmanAlign.h index 09e54fcba32..50a1d1564c6 100644 --- a/TPC/AliTPCkalmanAlign.h +++ b/TPC/AliTPCkalmanAlign.h @@ -6,6 +6,7 @@ #include "TNamed.h" #include "TMatrixD.h" +#include "TString.h" class TTreeSRedirector; class TObjArray; class AliTPCcalibAlign; @@ -24,6 +25,11 @@ public: 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); -- 2.43.0