AliTPCcalibCalib.cxx - bigger error close to the TPC boundaried in calibration refit
AliTPCcalibCosmic.cxx AliTPCcalibCosmic.h - high pt cosmic filtering to tree added
AliTPCCalibGlobalMisalignment.cxx AliTPCCalibGlobalMisalignment.h - adding global delta A side C side misaalignment
AliTPCcalibLaser.cxx AliTPCcalibLaser.h - histograming the cluster-survey trak residuals in the THnSparse
fQuadrantRQ1(0), //OROC long pads -delta ly+ - ly - rotation
fQuadrantRQ2(0), //OROC long pads -rotation
fMatrixGlobal(0), // global Alignment common
+ fMatrixGlobalDelta(0), // global Alignment Delta A side-c side
fArraySector(0) // fArraySector
{
//
delete fQuadrantRQ1; //OROC long pads -delta ly+ - ly - rotation
delete fQuadrantRQ2; //OROC long pads -rotation
delete fMatrixGlobal; // global matrix
+ delete fMatrixGlobal; // global matrix
delete fArraySector; // sector matrices
}
if (matrixGlobal) fMatrixGlobal = new TGeoHMatrix(*matrixGlobal);
}
+void AliTPCCalibGlobalMisalignment::SetAlignGlobalDelta(const TGeoMatrix * matrixGlobalDelta){
+ //
+ // Set global misalignment
+ // Object is OWNER
+ //
+ if (fMatrixGlobalDelta) delete fMatrixGlobalDelta;
+ fMatrixGlobalDelta=0;
+ if (matrixGlobalDelta) fMatrixGlobalDelta = new TGeoHMatrix(*matrixGlobalDelta);
+}
+
void AliTPCCalibGlobalMisalignment::SetAlignSectors(const TObjArray *arraySector){
//
// Set misalignment TObjArray of TGeoMatrices - for each sector
//
static AliTPCROC *tpcRoc =AliTPCROC::Instance();
Double_t xref = ( tpcRoc->GetPadRowRadii(0,0)+tpcRoc->GetPadRowRadii(36,tpcRoc->GetNRows(36)-1))*0.5;
-
+ Double_t xquadrant = tpcRoc->GetPadRowRadii(36,53); // row 53 from uli
+ Double_t xIO = ( tpcRoc->GetPadRowRadii(0,tpcRoc->GetNRows(0)-1)+tpcRoc->GetPadRowRadii(36,0))*0.5;
Double_t r=0, phi=0;
r = TMath::Sqrt( x[0]*x[0] + x[1]*x[1] );
phi = TMath::ATan2(x[1],x[0]);
Double_t posQG[3]={x[0],x[1],x[2]};
{
Double_t dly=0;
- Bool_t isQ0 = TMath::Abs(pos[0]-161)<28;
- Bool_t isQ1 = (pos[0]>189);
+ Bool_t isQ0 = (pos[0]<xquadrant)&&(pos[0]>xIO);
+ Bool_t isQ1 = (pos[0]>xquadrant);
Double_t sign = (pos[1]>0)? 1.: -1.;
if (isQ0){
if (fQuadrantQ0) dly+=sign*(*fQuadrantQ0)[isec%36]; // shift in cm
dx[1]+=pposC[1]-ppos[1];
dx[2]+=pposC[2]-ppos[2];
}
+ if (fMatrixGlobalDelta){
+ // apply global alignment matrix A-C Side side
+ Double_t ppos[3]={x[0],x[1],x[2]};
+ Double_t pposC[3]={x[0],x[1],x[2]};
+ fMatrixGlobalDelta->LocalToMaster(ppos,pposC);
+ Double_t ssign=(roc%36<18) ? 1.:-1.;
+ dx[0]+=ssign*(pposC[0]-ppos[0]);
+ dx[1]+=ssign*(pposC[1]-ppos[1]);
+ dx[2]+=ssign*(pposC[2]-ppos[2]);
+ }
if (fArraySector){
// apply global alignment matrix
// Alignment manipulation using TGeoMatrix
void SetAlignGlobal(const TGeoMatrix * matrixGlobal);
+ void SetAlignGlobalDelta(const TGeoMatrix * matrixGlobalDelta);
void SetAlignSectors(const TObjArray *arraySector);
TGeoMatrix* GetAlignGlobal() const {return fMatrixGlobal;}
+ TGeoMatrix* GetAlignGlobalDelta() const {return fMatrixGlobalDelta;}
TObjArray * GetAlignSectors() const {return fArraySector;}
//
static AliTPCCalibGlobalMisalignment* CreateOCDBAlign();
// Global alignment - use native ROOT representation
//
TGeoMatrix * fMatrixGlobal; // global Alignment common
+ TGeoMatrix * fMatrixGlobalDelta; // global Alignment common A side-C side
TObjArray * fArraySector; // local Alignmnet Sector
//
AliTPCCalibGlobalMisalignment& operator=(const AliTPCCalibGlobalMisalignment&);
AliTPCCalibGlobalMisalignment(const AliTPCCalibGlobalMisalignment&);
- ClassDef(AliTPCCalibGlobalMisalignment,2);
+ ClassDef(AliTPCCalibGlobalMisalignment,3);
};
#endif
if (TMath::Abs(parLine1[2])<0.8 &&TMath::Abs(parLine1[2])<0.8 ){ //angular cut
FillHisto(parLine1,parLine2,s1,s2);
ProcessAlign(parLine1,parLine2,s1,s2);
+ FillHisto((AliExternalTrackParam*)&tp1,(AliExternalTrackParam*)&tp2,s1,s2);
+ FillHisto((AliExternalTrackParam*)&tp2,(AliExternalTrackParam*)&tp1,s2,s1);
//UpdateKalman(s1,s2,par1, cov1, par2, cov2); - OBSOLETE to be removed - 50 % of time here
}
}
// 4 - local kz
//
axisName[0]="delta"; axisTitle[0]="#Delta (cm)";
- binsTrack[0]=60; xminTrack[0]=-0.6; xmaxTrack[0]=0.6;
+ if (TMath::Abs(AliTracker::GetBz())<0.01){
+ binsTrack[0]=60; xminTrack[0]=-1.2; xmaxTrack[0]=1.2;
+ }else{
+ binsTrack[0]=60; xminTrack[0]=-0.6; xmaxTrack[0]=0.6;
+ }
//
axisName[1]="sector"; axisTitle[1]="Sector Number";
binsTrack[1]=180; xminTrack[1]=0; xmaxTrack[1]=18;
// 3 - local ky
// 4 - local kz
// 5 - sector 1
- // 5 - sector 0
+ // 6 - sector 0
+ // 7 - z position 0
axisName[0]="delta"; axisTitle[0]="#Delta (cm)";
binsTrack[0]=60; xminTrack[0]=-0.6; xmaxTrack[0]=0.6;
//
axisName[6]="is0"; axisTitle[6]="is0";
binsTrack[6]=72; xminTrack[6]=0; xmaxTrack[6]=72;
+ //
+ axisName[7]="z"; axisTitle[7]="z(cm)";
+ binsTrack[7]=12; xminTrack[7]=-240; xmaxTrack[7]=240;
+ //
+ axisName[8]="IsPrimary"; axisTitle[8]="Is Primary";
+ binsTrack[8]=2; xminTrack[8]=-0.1; xmaxTrack[8]=1.1;
//
- xminTrack[0]=-0.3; xmaxTrack[0]=0.3;
- fTrackletDelta[0] = new THnSparseF("#Delta_{Y} (cm)","#Delta_{Y} (cm)", 7, binsTrack,xminTrack, xmaxTrack);
+ xminTrack[0]=-0.25; xmaxTrack[0]=0.25;
+ fTrackletDelta[0] = new THnSparseF("#Delta_{Y} (cm)","#Delta_{Y} (cm)", 9, binsTrack,xminTrack, xmaxTrack);
xminTrack[0]=-0.5; xmaxTrack[0]=0.5;
- fTrackletDelta[1] = new THnSparseF("#Delta_{Z} (cm)","#Delta_{Z} (cm)", 7, binsTrack,xminTrack, xmaxTrack);
+ fTrackletDelta[1] = new THnSparseF("#Delta_{Z} (cm)","#Delta_{Z} (cm)", 9, binsTrack,xminTrack, xmaxTrack);
xminTrack[0]=-0.005; xmaxTrack[0]=0.005;
- fTrackletDelta[2] = new THnSparseF("#Delta_{kY}","#Delta_{kY}", 7, binsTrack,xminTrack, xmaxTrack);
- xminTrack[0]=-0.005; xmaxTrack[0]=0.005;
- fTrackletDelta[3] = new THnSparseF("#Delta_{kZ}","#Delta_{kZ}", 7, binsTrack,xminTrack, xmaxTrack);
+ fTrackletDelta[2] = new THnSparseF("#Delta_{kY}","#Delta_{kY}", 9, binsTrack,xminTrack, xmaxTrack);
+ xminTrack[0]=-0.008; xmaxTrack[0]=0.008;
+ fTrackletDelta[3] = new THnSparseF("#Delta_{kZ}","#Delta_{kZ}", 9, binsTrack,xminTrack, xmaxTrack);
//
//
//
for (Int_t ivar=0;ivar<4;ivar++){
- for (Int_t ivar2=0;ivar2<7;ivar2++){
+ for (Int_t ivar2=0;ivar2<9;ivar2++){
fTrackletDelta[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
fTrackletDelta[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
}
//
// Fill residual histograms
// Track2-Track1
+ if (s2<s1) return;//
const Double_t kEpsilon=0.001;
Double_t x[8]={0,0,0,0,0,0,0,0};
AliExternalTrackParam p1(*tp1);
x[4]=0.5*(p1.GetTgl()+p2.GetTgl()); // mean tgl
x[5]=s2;
x[6]=s1;
-
+ x[7]=0.5*(p1.GetZ()+p2.GetZ());
+ // is primary ?
+ Int_t isPrimary = (TMath::Abs(p1.GetTgl()-p1.GetZ()/p1.GetX())<0.1) ? 1:0;
+ x[8]= isPrimary;
+ //
x[0]=p2.GetY()-p1.GetY();
fTrackletDelta[0]->Fill(x);
x[0]=p2.GetZ()-p1.GetZ();
if (!fClusterDelta[0]) MakeResidualHistos();
for (Int_t i=0; i<2; i++){
- if (align->fClusterDelta[i]) fClusterDelta[i]->Add(align->fClusterDelta[i]);
+ if (align->fClusterDelta[i]){
+ fClusterDelta[i]->Add(align->fClusterDelta[i]);
+ // align->fClusterDelta[i]->GetAxis(0)->SetRangeUser(-0.87,0.87);
+// align->fClusterDelta[i]->GetAxis(3)->SetRangeUser(-0.87,0.87);
+// fClusterDelta[i]->GetAxis(0)->SetRangeUser(-0.87,0.87);
+// fClusterDelta[i]->GetAxis(3)->SetRangeUser(-0.87,0.87);
+// Int_t idim[4]={0,1,2,3};
+// THnSparse *htemp=align->fClusterDelta[i]->Projection(4,idim);
+// THnSparse *htemp1=fClusterDelta[i]->Projection(4,idim);
+// htemp1->Add(htemp);
+// delete fClusterDelta[i];
+// fClusterDelta[i]=htemp1;
+// delete htemp;
+ }
}
-
+
for (Int_t i=0; i<4; i++){
if (!fTrackletDelta[i] && align->fTrackletDelta[i]) {
fTrackletDelta[i]= (THnSparse*)(align->fTrackletDelta[i]->Clone());
continue;
}
- if (align->fTrackletDelta[i]) fTrackletDelta[i]->Add(align->fTrackletDelta[i]);
+ if (align->fTrackletDelta[i]) {
+ fTrackletDelta[i]->Add(align->fTrackletDelta[i]);
+ //
+ // align->fTrackletDelta[i]->GetAxis(3)->SetRangeUser(-0.36,0.36);
+// align->fTrackletDelta[i]->GetAxis(4)->SetRangeUser(-0.87,0.87);
+// fTrackletDelta[i]->GetAxis(3)->SetRangeUser(-0.36,0.36);
+// fTrackletDelta[i]->GetAxis(4)->SetRangeUser(-0.87,0.87);
+// //
+// Int_t idim[9]={0,1,2,3,4,5,6,7,8};
+// THnSparse *htemp=align->fTrackletDelta[i]->Projection(9,idim);
+// THnSparse *htemp1=fTrackletDelta[i]->Projection(9,idim);
+// htemp1->Add(htemp);
+// delete fTrackletDelta[i];
+// fTrackletDelta[i]=htemp1;
+// delete htemp;
+ }
}
}
//
// 1. Apply selection
// 2. Refit the track - in-out
- // - update the cluster delta in upper part
// 3. Refit the track - out-in
- // - update the cluster delta histo lower part
+ // 4. Combine In and Out track - - fil cluster residuals
//
const Double_t kPtCut=1.0; // pt
const Double_t kSnpCut=0.2; // snp cut
const Double_t kVertexCut=1;
const Double_t kMaxDist=0.5; // max distance between tracks and cluster
const Double_t kEdgeCut = 2.5;
+ const Double_t kDelta2=0.2*0.2; // initial increase in covar matrix
+ const Double_t kSigma=0.3; // error increase towards edges of TPC
+ const Double_t kSkipBoundary=7.5; // skip track updates in the boundary IFC,OFC, IO
+ //
if (!fCurrentTrack) return;
if (!fCurrentFriendTrack) return;
Float_t vertexXY=0,vertexZ=0;
if (seed->GetNumberOfClusters()<kNclCut) return;
if (TMath::Abs(seed->GetSnp())>kSnpCut) return;
if (!fClusterDelta[0]) MakeResidualHistos();
-
+ //
+ AliExternalTrackParam fitIn[160];
+ AliExternalTrackParam fitOut[160];
+ AliTPCROC * roc = AliTPCROC::Instance();
+ Double_t xmiddle = ( roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
+ Double_t xDiff = ( -roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
+ Double_t xIFC = ( roc->GetPadRowRadii(0,0));
+ Double_t xOFC = ( roc->GetPadRowRadii(36,roc->GetNRows(36)-1));
+ //
Int_t detector=-1;
//
//
AliExternalTrackParam trackIn = *(fCurrentTrack->GetInnerParam());
AliExternalTrackParam trackOut = *(fCurrentFriendTrack->GetTPCOut());
+ trackIn.ResetCovariance(10);
+ trackOut.ResetCovariance(10);
+ Double_t *covarIn = (Double_t*)trackIn.GetCovariance();
+ Double_t *covarOut = (Double_t*)trackOut.GetCovariance();
+ covarIn[0]+=kDelta2; covarIn[2]+=kDelta2;
+ covarIn[5]+=kDelta2/(100.*100.); covarIn[9]=kDelta2/(100.*100.);
+ covarIn[14]+=kDelta2/(5.*5.);
+ covarOut[0]+=kDelta2; covarOut[2]+=kDelta2;
+ covarOut[5]+=kDelta2/(100.*100.); covarOut[9]=kDelta2/(100.*100.);
+ covarOut[14]+=kDelta2/(5.*5.);
+ //
static Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
//
Int_t ncl=0;
ncl++;
}
if (ncl<kNclCut) return;
-
Int_t nclIn=0,nclOut=0;
Double_t xyz[3];
//
if (detector<0) detector=cl->GetDetector()%36;
Int_t sector = cl->GetDetector();
Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackOut.GetAlpha();
+ if (cl->GetDetector()%36!=detector) continue;
if (TMath::Abs(dalpha)>0.01){
if (!trackOut.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break;
}
Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
- Double_t cov[3]={0.01,0.,0.01};
- AliTPCseed::GetError(cl, &trackOut,cov[0],cov[2]);
+ Double_t cov[3]={0.1,0.,0.1};
Double_t dedge = cl->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(trackOut.GetY());
- cov[0]+=1./(irow+1.); // bigger error at boundary
- cov[0]+=1./(160.-irow); // bigger error at boundary
- cov[2]+=1./(irow+1.); // bigger error at boundary
- cov[2]+=1./(160.-irow); // bigger error at boundary
- cov[0]+=0.5/dedge; // bigger error close to the boundary
- cov[2]+=0.5/dedge; // bigger error close to the boundary
+ Double_t dmiddle = TMath::Abs(cl->GetX()-xmiddle)/xDiff;
+ dmiddle*=dmiddle;
+ //
+ cov[0]+=kSigma*dmiddle; // bigger error at boundary
+ cov[0]+=kSigma*dmiddle; // bigger error at boundary
+ cov[2]+=kSigma*dmiddle; // bigger error at boundary
+ cov[2]+=kSigma*dmiddle; // bigger error at boundary
+ cov[0]+=kSigma/dedge; // bigger error close to the boundary
+ cov[2]+=kSigma/dedge; // bigger error close to the boundary
cov[0]*=cov[0];
cov[2]*=cov[2];
- if (!AliTracker::PropagateTrackToBxByBz(&trackOut, r[0],mass,1.,kFALSE)) continue;
-
+ if (!AliTracker::PropagateTrackToBxByBz(&trackOut, r[0],mass,1.,kFALSE)) continue;
if (TMath::Abs(dedge)<kEdgeCut) continue;
-
+ //
+ Bool_t doUpdate=kTRUE;
+ if (TMath::Abs(cl->GetX()-xIFC)<kSkipBoundary) doUpdate=kFALSE;
+ if (TMath::Abs(cl->GetX()-xOFC)<kSkipBoundary) doUpdate=kFALSE;
+ if (TMath::Abs(cl->GetX()-fXIO)<kSkipBoundary) doUpdate=kFALSE;
+ //
if (TMath::Abs(cl->GetY()-trackOut.GetY())<kMaxDist){
nclOut++;
- trackOut.Update(&r[1],cov);
+ if (doUpdate) trackOut.Update(&r[1],cov);
}
- if (nclOut<kNclCut/2) continue;
- if (cl->GetDetector()%36!=detector) continue;
- //
- // fill residual histogram
- //
- Double_t resVector[5];
- trackOut.GetXYZ(xyz);
- resVector[1]= 9.*TMath::ATan2(xyz[1],xyz[0])/TMath::Pi();
- if (resVector[1]<0) resVector[1]+=18;
- resVector[2]= TMath::Sqrt(cl->GetX()*cl->GetX()+cl->GetY()*cl->GetY());
- resVector[3]= cl->GetZ()/resVector[2];
- //
- resVector[0]= cl->GetY()-trackOut.GetY();
- fClusterDelta[0]->Fill(resVector);
- resVector[0]= cl->GetZ()-trackOut.GetZ();
- fClusterDelta[1]->Fill(resVector);
+ fitOut[irow]=trackOut;
}
+
//
- // Refit in - store residual maps
+ // Refit In - store residual maps
//
for (Int_t irow=159; irow>=0; irow--){
AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
if (detector<0) detector=cl->GetDetector()%36;
Int_t sector = cl->GetDetector();
Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackIn.GetAlpha();
+ if (cl->GetDetector()%36!=detector) continue;
if (TMath::Abs(dalpha)>0.01){
if (!trackIn.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break;
}
Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
- Double_t cov[3]={0.01,0.,0.01};
- AliTPCseed::GetError(cl, &trackIn,cov[0],cov[2]);
+ Double_t cov[3]={0.1,0.,0.1};
Double_t dedge = cl->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(trackIn.GetY());
- cov[0]+=1./(irow+1.); // bigger error at boundary
- cov[0]+=1./(160.-irow); // bigger error at boundary
- cov[2]+=1./(irow+1.); // bigger error at boundary
- cov[2]+=1./(160.-irow); // bigger error at boundary
- cov[0]+=0.5/dedge; // bigger error close to the boundary +-
- cov[2]+=0.5/dedge; // bigger error close to the boundary +-
+ Double_t dmiddle = TMath::Abs(cl->GetX()-xmiddle)/xDiff;
+ dmiddle*=dmiddle;
+ //
+ cov[0]+=kSigma*dmiddle; // bigger error at boundary
+ cov[0]+=kSigma*dmiddle; // bigger error at boundary
+ cov[2]+=kSigma*dmiddle; // bigger error at boundary
+ cov[2]+=kSigma*dmiddle; // bigger error at boundary
+ cov[0]+=kSigma/dedge; // bigger error close to the boundary
+ cov[2]+=kSigma/dedge; // bigger error close to the boundary
cov[0]*=cov[0];
cov[2]*=cov[2];
- if (!AliTracker::PropagateTrackToBxByBz(&trackIn, r[0],mass,1.,kFALSE)) continue;
+ if (!AliTracker::PropagateTrackToBxByBz(&trackIn, r[0],mass,1.,kFALSE)) continue;
if (TMath::Abs(dedge)<kEdgeCut) continue;
-
-
+ Bool_t doUpdate=kTRUE;
+ if (TMath::Abs(cl->GetX()-xIFC)<kSkipBoundary) doUpdate=kFALSE;
+ if (TMath::Abs(cl->GetX()-xOFC)<kSkipBoundary) doUpdate=kFALSE;
+ if (TMath::Abs(cl->GetX()-fXIO)<kSkipBoundary) doUpdate=kFALSE;
if (TMath::Abs(cl->GetY()-trackIn.GetY())<kMaxDist){
nclIn++;
- trackIn.Update(&r[1],cov);
+ if (doUpdate) trackIn.Update(&r[1],cov);
}
- if (nclIn<kNclCut/2) continue;
- if (cl->GetDetector()%36!=detector) continue;
+ fitIn[irow]=trackIn;
+ }
+ //
+ //
+ for (Int_t irow=159; irow>=0; irow--){
//
- // fill residual histogram
+ // Update kalman - +- direction
+ // Store cluster residuals
+ AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
+ if (!cl) continue;
+ if (cl->GetX()<80) continue;
+ if (detector<0) detector=cl->GetDetector()%36;
+ if (cl->GetDetector()%36!=detector) continue;
+ AliExternalTrackParam trackSmooth = fitIn[irow];
+ AliTrackerBase::UpdateTrack(trackSmooth, fitOut[irow]);
//
Double_t resVector[5];
- trackIn.GetXYZ(xyz);
+ trackSmooth.GetXYZ(xyz);
resVector[1]= 9.*TMath::ATan2(xyz[1],xyz[0])/TMath::Pi();
if (resVector[1]<0) resVector[1]+=18;
resVector[2]= TMath::Sqrt(cl->GetX()*cl->GetX()+cl->GetY()*cl->GetY());
resVector[3]= cl->GetZ()/resVector[2];
//
- resVector[0]= cl->GetY()-trackIn.GetY();
+ resVector[0]= cl->GetY()-trackSmooth.GetY();
fClusterDelta[0]->Fill(resVector);
- resVector[0]= cl->GetZ()-trackIn.GetZ();
+ resVector[0]= cl->GetZ()-trackSmooth.GetZ();
fClusterDelta[1]->Fill(resVector);
}
-
}
//
if (TMath::Abs(AliTracker::GetBz())>0.5) return;
if (!fClusterDelta[0]) MakeResidualHistos();
- const Int_t kMinClusterF=40;
+ // const Int_t kMinClusterF=40;
const Int_t kMinClusterFit=10;
const Int_t kMinClusterQ=10;
//
}
if (TMath::Abs(x[0])<10){
fyf.AddPoint(x,c->GetY(),0.1); //use only middle rows+-10cm
+ fzf.AddPoint(x,c->GetZ(),0.1);
}
- fzf.AddPoint(x,c->GetZ(),0.1);
}
nf = fyf.GetNpoints();
if (fyf.GetNpoints()<kMinClusterFit) return; // not enough points - skip
- if (fzf.GetNpoints()<kMinClusterF) return; // not enough points - skip
+ if (fzf.GetNpoints()<kMinClusterFit) return; // not enough points - skip
fyf.Eval();
fyf.GetParameters(pyf);
fyf.GetErrors(peyf);
//
// 0 - Setup transform object
//
+ const Double_t kxIFC = 83.; // position of IFC
+ const Double_t kxOFC = 250.; // position of OFC
+ const Double_t kaFC = 1.; // amplitude
+ const Double_t ktFC = 5.0; // slope of error
+ //cov[0]+= kaFC*(TMath::Exp(-TMath::Abs(cl->GetX()-kxIFC)/ktFC)+TMath::Exp(-TMath::Abs(cl->GetX()-kxOFC)/ktFC));
+ //cov[2]+= kaFC*(TMath::Exp(-TMath::Abs(cl->GetX()-kxIFC)/ktFC)+TMath::Exp(-TMath::Abs(cl->GetX()-kxOFC)/ktFC));
+
static Int_t streamCounter=0;
streamCounter++;
AliESDfriendTrack *friendTrack = fCurrentFriendTrack;
AliTPCseed::GetError(cl, &trackIn,cov[0],cov[2]);
cov[0]*=cov[0];
cov[2]*=cov[2];
+ cov[0]+= kaFC*(TMath::Exp(-TMath::Abs(cl->GetX()-kxIFC)/ktFC)+TMath::Exp(-TMath::Abs(cl->GetX()-kxOFC)/ktFC));
+ cov[2]+= kaFC*(TMath::Exp(-TMath::Abs(cl->GetX()-kxIFC)/ktFC)+TMath::Exp(-TMath::Abs(cl->GetX()-kxOFC)/ktFC));
trackIn.GetXYZ(xyz);
// Double_t bz = AliTracker::GetBz(xyz);
AliTPCseed::GetError(cl, &trackOut,cov[0],cov[2]);
cov[0]*=cov[0];
cov[2]*=cov[2];
+ cov[0]+= kaFC*(TMath::Exp(-TMath::Abs(cl->GetX()-kxIFC)/ktFC)+TMath::Exp(-TMath::Abs(cl->GetX()-kxOFC)/ktFC));
+ cov[2]+= kaFC*(TMath::Exp(-TMath::Abs(cl->GetX()-kxIFC)/ktFC)+TMath::Exp(-TMath::Abs(cl->GetX()-kxOFC)/ktFC));
trackOut.GetXYZ(xyz);
//Double_t bz = AliTracker::GetBz(xyz);
// if (!trackOut.PropagateTo(r[0],bz)) continue;
AliTPCseed::GetError(cl, &trackIn,cov[0],cov[2]);
cov[0]*=cov[0];
cov[2]*=cov[2];
+ cov[0]+= kaFC*(TMath::Exp(-TMath::Abs(cl->GetX()-kxIFC)/ktFC)+TMath::Exp(-TMath::Abs(cl->GetX()-kxOFC)/ktFC));
+ cov[2]+= kaFC*(TMath::Exp(-TMath::Abs(cl->GetX()-kxIFC)/ktFC)+TMath::Exp(-TMath::Abs(cl->GetX()-kxOFC)/ktFC));
+
trackIn.GetXYZ(xyz);
//Double_t bz = AliTracker::GetBz(xyz);
+
+
/**************************************************************************
* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
* *
#include "AliTracker.h"
#include "AliMagF.h"
#include "AliTPCCalROC.h"
-
+#include "AliTPCParam.h"
#include "AliLog.h"
#include "AliTPCcalibCosmic.h"
#include "TTreeStream.h"
#include "AliTPCTracklet.h"
//#include "AliESDcosmic.h"
-
-
+#include "AliRecoParam.h"
+#include "AliMultiplicity.h"
+#include "AliTPCTransform.h"
+#include "AliTPCcalibDB.h"
+#include "AliTPCseed.h"
+#include "AliGRPObject.h"
+#include "AliTPCCorrection.h"
ClassImp(AliTPCcalibCosmic)
fCutMaxD(5), // maximal distance in rfi ditection
fCutMaxDz(40), // maximal distance in z ditection
fCutTheta(0.03), // maximal distan theta
- fCutMinDir(-0.99) // direction vector products
+ fCutMinDir(-0.99), // direction vector products
+ fCosmicTree(0) // tree with cosmic data
{
AliInfo("Default Constructor");
for (Int_t ihis=0; ihis<6;ihis++){
fCutMaxD(5), // maximal distance in rfi ditection
fCutMaxDz(40), // maximal distance in z ditection
fCutTheta(0.03), // maximal distan theta
- fCutMinDir(-0.99) // direction vector products
+ fCutMinDir(-0.99), // direction vector products
+ fCosmicTree(0) // tree with cosmic data
{
SetName(name);
SetTitle(title);
// 6 - P4 - 1/pt mean
// 7 - pt - pt mean
// 8 - alpha
-
- Double_t xminTrack[9], xmaxTrack[9];
- Int_t binsTrack[9];
- TString axisName[9];
+ // 9 - is corss indicator
+ Int_t ndim=10;
+ Double_t xminTrack[10], xmaxTrack[10];
+ Int_t binsTrack[10];
+ TString axisName[10];
//
binsTrack[0] =100;
axisName[0] ="#Delta";
xminTrack[6] =-2; xmaxTrack[6]=2; //
axisName[6] ="1/pt (1/GeV)";
//
- binsTrack[7] =40;
- xminTrack[7] =0.2; xmaxTrack[7]=50; //
+ binsTrack[7] =50;
+ xminTrack[7] =1; xmaxTrack[7]=1000; //
axisName[7] ="pt (GeV)";
//
binsTrack[8] =18;
xminTrack[8] =0; xmaxTrack[8]=TMath::Pi(); //
axisName[8] ="alpha";
//
+ binsTrack[9] =3;
+ xminTrack[9] =-0.1; xmaxTrack[9]=2.1; //
+ axisName[9] ="cross";
+ //
// delta y
xminTrack[0] =-1; xmaxTrack[0]=1; //
- fHistoDelta[0] = new THnSparseS("#Delta_{Y} (cm)","#Delta_{Y} (cm)", 9, binsTrack,xminTrack, xmaxTrack);
+ fHistoDelta[0] = new THnSparseS("#Delta_{Y} (cm)","#Delta_{Y} (cm)", ndim, binsTrack,xminTrack, xmaxTrack);
xminTrack[0] =-5; xmaxTrack[0]=5; //
- fHistoPull[0] = new THnSparseS("#Delta_{Y} (unit)","#Delta_{Y} (unit)", 9, binsTrack,xminTrack, xmaxTrack);
+ fHistoPull[0] = new THnSparseS("#Delta_{Y} (unit)","#Delta_{Y} (unit)", ndim, binsTrack,xminTrack, xmaxTrack);
//
// delta z
xminTrack[0] =-1; xmaxTrack[0]=1; //
- fHistoDelta[1] = new THnSparseS("#Delta_{Z} (cm)","#Delta_{Z} (cm)", 9, binsTrack,xminTrack, xmaxTrack);
+ fHistoDelta[1] = new THnSparseS("#Delta_{Z} (cm)","#Delta_{Z} (cm)", ndim, binsTrack,xminTrack, xmaxTrack);
xminTrack[0] =-5; xmaxTrack[0]=5; //
- fHistoPull[1] = new THnSparseS("#Delta_{Z} (unit)","#Delta_{Z} (unit)", 9, binsTrack,xminTrack, xmaxTrack);
+ fHistoPull[1] = new THnSparseS("#Delta_{Z} (unit)","#Delta_{Z} (unit)", ndim, binsTrack,xminTrack, xmaxTrack);
//
// delta P2
xminTrack[0] =-10; xmaxTrack[0]=10; //
- fHistoDelta[2] = new THnSparseS("#Delta_{#phi} (mrad)","#Delta_{#phi} (mrad)", 9, binsTrack,xminTrack, xmaxTrack);
+ fHistoDelta[2] = new THnSparseS("#Delta_{#phi} (mrad)","#Delta_{#phi} (mrad)", ndim, binsTrack,xminTrack, xmaxTrack);
xminTrack[0] =-5; xmaxTrack[0]=5; //
- fHistoPull[2] = new THnSparseS("#Delta_{#phi} (unit)","#Delta_{#phi} (unit)", 9, binsTrack,xminTrack, xmaxTrack);
+ fHistoPull[2] = new THnSparseS("#Delta_{#phi} (unit)","#Delta_{#phi} (unit)", ndim, binsTrack,xminTrack, xmaxTrack);
//
// delta P3
xminTrack[0] =-10; xmaxTrack[0]=10; //
- fHistoDelta[3] = new THnSparseS("#Delta_{#theta} (mrad)","#Delta_{#theta} (mrad)", 9, binsTrack,xminTrack, xmaxTrack);
+ fHistoDelta[3] = new THnSparseS("#Delta_{#theta} (mrad)","#Delta_{#theta} (mrad)", ndim, binsTrack,xminTrack, xmaxTrack);
xminTrack[0] =-5; xmaxTrack[0]=5; //
- fHistoPull[3] = new THnSparseS("#Delta_{#theta} (unit)","#Delta_{#theta} (unit)", 9, binsTrack,xminTrack, xmaxTrack);
+ fHistoPull[3] = new THnSparseS("#Delta_{#theta} (unit)","#Delta_{#theta} (unit)", ndim, binsTrack,xminTrack, xmaxTrack);
//
// delta P4
xminTrack[0] =-0.2; xmaxTrack[0]=0.2; //
- fHistoDelta[4] = new THnSparseS("#Delta_{1/pt} (1/GeV)","#Delta_{1/pt} (1/GeV)", 9, binsTrack,xminTrack, xmaxTrack);
+ fHistoDelta[4] = new THnSparseS("#Delta_{1/pt} (1/GeV)","#Delta_{1/pt} (1/GeV)", ndim, binsTrack,xminTrack, xmaxTrack);
xminTrack[0] =-5; xmaxTrack[0]=5; //
- fHistoPull[4] = new THnSparseS("#Delta_{1/pt} (unit)","#Delta_{1/pt} (unit)", 9, binsTrack,xminTrack, xmaxTrack);
+ fHistoPull[4] = new THnSparseS("#Delta_{1/pt} (unit)","#Delta_{1/pt} (unit)", ndim, binsTrack,xminTrack, xmaxTrack);
//
// delta Pt
xminTrack[0] =-0.5; xmaxTrack[0]=0.5; //
- fHistoDelta[5] = new THnSparseS("#Delta_{pt}/p_{t}","#Delta_{pt}/p_{t}", 9, binsTrack,xminTrack, xmaxTrack);
+ fHistoDelta[5] = new THnSparseS("#Delta_{pt}/p_{t}","#Delta_{pt}/p_{t}", ndim, binsTrack,xminTrack, xmaxTrack);
xminTrack[0] =-5; xmaxTrack[0]=5; //
- fHistoPull[5] = new THnSparseS("#Delta_{pt}/p_{t} (unit)","#Delta_{pt}/p_{t} (unit)", 9, binsTrack,xminTrack, xmaxTrack);
+ fHistoPull[5] = new THnSparseS("#Delta_{pt}/p_{t} (unit)","#Delta_{pt}/p_{t} (unit)", ndim, binsTrack,xminTrack, xmaxTrack);
//
for (Int_t idedx=0;idedx<4;idedx++){
fHistodEdxMax[idedx] = new THnSparseS(Form("dEdx_{MaxUp}/dEdx_{MaxDown}_Pad%d",idedx),
Form("dEdx_{MaxUp}/dEdx_{MaxDown}_Pad%d",idedx),
- 9, binsTrack,xminTrack, xmaxTrack);
+ ndim, binsTrack,xminTrack, xmaxTrack);
fHistodEdxTot[idedx] = new THnSparseS(Form("dEdx_{TotUp}/dEdx_{TotDown}_Pad%d",idedx),
Form("dEdx_{TotUp}/dEdx_{TotDown}_Pad%d",idedx),
- 9, binsTrack,xminTrack, xmaxTrack);
+ ndim, binsTrack,xminTrack, xmaxTrack);
}
for (Int_t ivar=0;ivar<6;ivar++){
- for (Int_t ivar2=0;ivar2<9;ivar2++){
+ for (Int_t ivar2=0;ivar2<ndim;ivar2++){
fHistoDelta[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
fHistoDelta[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
fHistoPull[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
fHistodEdxTot[ivar]->Add(cosmic->fHistodEdxTot[ivar]);
}
}
+ if (cosmic->fCosmicTree){
+ if (!fCosmicTree) {
+ fCosmicTree = new TTree("pairs","pairs");
+ fCosmicTree->SetDirectory(0);
+ }
+ AliTPCcalibCosmic::AddTree(fCosmicTree,cosmic->fCosmicTree);
+ }
}
return;
}
+ //
+ //Int_t isOK=kTRUE;
+ // COSMIC not signed properly
+ // UInt_t specie = event->GetEventSpecie(); // select only cosmic events
+ //if (specie==AliRecoParam::kCosmic || specie==AliRecoParam::kCalib) {
+ // isOK = kTRUE;
+ //}
+ //if (!isOK) return;
+ // Work around
+ FindCosmicPairs(event);
+ return;
+ const AliMultiplicity *multiplicity = event->GetMultiplicity();
+ Int_t ntracklets = multiplicity->GetNumberOfTracklets();
+ if (ntracklets>6) return; // filter out "normal" event with high multiplicity
+ const TString &trigger = event->GetFiredTriggerClasses();
+ if (trigger.Contains("C0OB0")==0) return;
+
FindPairs(event); // nearly everything takes place in find pairs...
if (GetDebugLevel()>20) printf("Hallo world: Im here and processing an event\n");
Int_t ntracks=event->GetNumberOfTracks();
fHistNTracks->Fill(ntracks);
- if (ntracks==0) return;
- // AliESDcosmic cosmicESD;
-// TTreeSRedirector * cstream = GetDebugStreamer();
-// cosmicESD.SetDebugStreamer(cstream);
-// cosmicESD.ProcessEvent(event);
-// if (cstream) cosmicESD.DumpToTree();
-
}
-void AliTPCcalibCosmic::FillHistoPerformance(const AliExternalTrackParam *par0, const AliExternalTrackParam *par1, const AliExternalTrackParam *inner0, const AliExternalTrackParam */*inner1*/, AliTPCseed *seed0, AliTPCseed *seed1, const AliExternalTrackParam *param0Combined ){
+void AliTPCcalibCosmic::FillHistoPerformance(const AliExternalTrackParam *par0, const AliExternalTrackParam *par1, const AliExternalTrackParam *inner0, const AliExternalTrackParam */*inner1*/, AliTPCseed *seed0, AliTPCseed *seed1, const AliExternalTrackParam *param0Combined , Int_t cross){
//
// par0,par1 - parameter of tracks at DCA 0
// inner0,inner1 - parameter of tracks at the TPC entrance
Int_t kMinCldEdx =20;
Int_t ncl0 = seed0->GetNumberOfClusters();
Int_t ncl1 = seed1->GetNumberOfClusters();
-
const Double_t kpullCut = 10;
- Double_t x[9];
+ Double_t x[10];
Double_t xyz0[3];
Double_t xyz1[3];
par0->GetXYZ(xyz0);
x[6] = param0Combined->GetSigned1Pt();
x[7] = param0Combined->Pt();
x[8] = alpha;
+ x[9] = cross;
// deltas
Double_t delta[6];
Double_t sigma[6];
//
if (isOK) for (Int_t ivar=0;ivar<6;ivar++){
- x[0]= delta[ivar]/TMath::Sqrt(2);
- if (ivar==2 || ivar ==3) x[0]*=1000;
+ x[0]= delta[ivar]; // Modifiation 10.10 use not normalized deltas
+ if (ivar==2 || ivar ==3) x[0]*=1000; // angles in radian
fHistoDelta[ivar]->Fill(x);
if (sigma[ivar]>0){
x[0]= delta[ivar]/sigma[ivar];
Bool_t isPair = IsPair(¶m0,¶m1);
//
if (isPair) FillAcordeHist(track0);
+ if (isPair &¶m0.Pt()>1) {
+ const TString &trigger = event->GetFiredTriggerClasses();
+ UInt_t specie = event->GetEventSpecie();
+ printf("COSMIC ?\t%s\t%d\t%f\t%f\n", trigger.Data(),specie, param0.GetZ(), param1.GetZ());
+ }
//
// combined track params
//
AliExternalTrackParam *par0U=MakeCombinedTrack(¶m0,¶m1);
AliExternalTrackParam *par1U=MakeCombinedTrack(¶m1,¶m0);
-
-
+
//
if (fStreamLevel>0){
TTreeSRedirector * cstream = GetDebugStreamer();
//
//
//
- FillHistoPerformance(¶m0, ¶m1, ip0, ip1, seed0, seed1,par0U);
+ Int_t cross =0; // 0 no cross, 2 cross on both sides
+ if (isCrossI) cross+=1;
+ if (isCrossO) cross+=1;
+ FillHistoPerformance(¶m0, ¶m1, ip0, ip1, seed0, seed1,par0U, cross);
MaterialBudgetDump(¶m0, ¶m1, ip0, ip1, seed0, seed1,par0U,par1U);
if (cstream) {
(*cstream) << "Track0" <<
+void AliTPCcalibCosmic::FindCosmicPairs(AliESDEvent * event){
+ //
+ // find cosmic pairs trigger by random trigger
+ //
+ //
+ AliESDVertex *vertexSPD = (AliESDVertex *)event->GetPrimaryVertexSPD();
+ AliESDVertex *vertexTPC = (AliESDVertex *)event->GetPrimaryVertexTPC();
+ AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
+ const Double_t kMinPt=1;
+ const Double_t kMinPtMax=0.8;
+ const Double_t kMinNcl=50;
+ const Double_t kMaxDelta[5]={2,600,0.02,0.02,0.1};
+ Int_t ntracks=event->GetNumberOfTracks();
+ Float_t dcaTPC[2]={0,0};
+ Float_t covTPC[3]={0,0,0};
+
+ UInt_t specie = event->GetEventSpecie(); // skip laser events
+ if (specie==AliRecoParam::kCalib) return;
+
+
+
+ for (Int_t itrack0=0;itrack0<ntracks;itrack0++) {
+ AliESDtrack *track0 = event->GetTrack(itrack0);
+ if (!track0) continue;
+ if (!track0->IsOn(AliESDtrack::kTPCrefit)) continue;
+
+ if (TMath::Abs(AliTracker::GetBz())>1&&track0->Pt()<kMinPt) continue;
+ if (track0->GetTPCncls()<kMinNcl) continue;
+ if (TMath::Abs(track0->GetY())<kMaxDelta[0]) continue;
+ if (track0->GetKinkIndex(0)>0) continue;
+ const Double_t * par0=track0->GetParameter(); //track param at rhe DCA
+ //rm primaries
+ //
+ //track0->GetImpactParametersTPC(dcaTPC,covTPC);
+ //if (TMath::Abs(dcaTPC[0])<kMaxDelta[0]) continue;
+ //if (TMath::Abs(dcaTPC[1])<kMaxDelta[0]*2) continue;
+ const AliExternalTrackParam * trackIn0 = track0->GetInnerParam();
+ for (Int_t itrack1=itrack0+1;itrack1<ntracks;itrack1++) {
+ AliESDtrack *track1 = event->GetTrack(itrack1);
+ if (!track1) continue;
+ if (!track1->IsOn(AliESDtrack::kTPCrefit)) continue;
+ if (track1->GetKinkIndex(0)>0) continue;
+ if (TMath::Abs(AliTracker::GetBz())>1&&track1->Pt()<kMinPt) continue;
+ if (track1->GetTPCncls()<kMinNcl) continue;
+ if (TMath::Abs(AliTracker::GetBz())>1&&TMath::Max(track1->Pt(), track0->Pt())<kMinPtMax) continue;
+ if (TMath::Abs(track1->GetY())<kMaxDelta[0]) continue;
+ //track1->GetImpactParametersTPC(dcaTPC,covTPC);
+ // if (TMath::Abs(dcaTPC[0])<kMaxDelta[0]) continue;
+ //if (TMath::Abs(dcaTPC[1])<kMaxDelta[0]*2) continue;
+ //
+ const Double_t* par1=track1->GetParameter(); //track param at rhe DCA
+ //
+ Bool_t isPair=kTRUE;
+ for (Int_t ipar=0; ipar<5; ipar++){
+ if (ipar==4&&TMath::Abs(AliTracker::GetBz())<1) continue; // 1/pt not defined for B field off
+ if (TMath::Abs(TMath::Abs(par0[ipar])-TMath::Abs(par1[ipar]))>kMaxDelta[ipar]) isPair=kFALSE;
+ }
+ if (!isPair) continue;
+ if (TMath::Abs(TMath::Abs(track0->GetAlpha()-track1->GetAlpha())-TMath::Pi())>kMaxDelta[2]) isPair=kFALSE;
+ //delta with correct sign
+ /*
+ TCut cut0="abs(t1.fP[0]+t0.fP[0])<2"
+ TCut cut3="abs(t1.fP[3]+t0.fP[3])<0.02"
+ TCut cut4="abs(t1.fP[4]+t0.fP[4])<0.2"
+ */
+ if (TMath::Abs(par0[0]+par1[0])>kMaxDelta[0]) isPair=kFALSE; //delta y opposite sign
+ if (TMath::Abs(par0[3]+par1[3])>kMaxDelta[3]) isPair=kFALSE; //delta tgl opposite sign
+ if (TMath::Abs(AliTracker::GetBz())>1 && TMath::Abs(par0[4]+par1[4])>kMaxDelta[4]) isPair=kFALSE; //delta 1/pt opposite sign
+ if (!isPair) continue;
+ const AliExternalTrackParam * trackIn1 = track1->GetInnerParam();
+ //
+ //
+ TTreeSRedirector * pcstream = GetDebugStreamer();
+ Int_t ntracksSPD = vertexSPD->GetNContributors();
+ Int_t ntracksTPC = vertexTPC->GetNContributors();
+ //
+ AliESDfriendTrack *friendTrack0 = esdFriend->GetTrack(itrack0);
+ if (!friendTrack0) continue;
+ AliESDfriendTrack *friendTrack1 = esdFriend->GetTrack(itrack1);
+ if (!friendTrack1) continue;
+ TObject *calibObject;
+ AliTPCseed *seed0 = 0;
+ AliTPCseed *seed1 = 0;
+ //
+ for (Int_t l=0;(calibObject=friendTrack0->GetCalibObject(l));++l) {
+ if ((seed0=dynamic_cast<AliTPCseed*>(calibObject))) break;
+ }
+ for (Int_t l=0;(calibObject=friendTrack1->GetCalibObject(l));++l) {
+ if ((seed1=dynamic_cast<AliTPCseed*>(calibObject))) break;
+ }
+ //
+ if (pcstream){
+ (*pcstream)<<"pairs"<<
+ "run="<<fRun<< // run number
+ "event="<<fEvent<< // event number
+ "time="<<fTime<< // time stamp of event
+ "trigger="<<fTrigger<< // trigger
+ "triggerClass="<<&fTriggerClass<< // trigger
+ "mag="<<fMagF<< // magnetic field
+ //
+ "nSPD="<<ntracksSPD<<
+ "nTPC="<<ntracksTPC<<
+ "vSPD.="<<vertexSPD<< //primary vertex -SPD
+ "vTPC.="<<vertexTPC<< //primary vertex -TPC
+ "t0.="<<track0<< //track0
+ "t1.="<<track1<< //track1
+ "ft0.="<<friendTrack0<< //track0
+ "ft1.="<<friendTrack1<< //track1
+ "s0.="<<seed0<< //track0
+ "s1.="<<seed1<< //track1
+ "\n";
+ }
+ if (!fCosmicTree) {
+ fCosmicTree = new TTree("pairs","pairs");
+ fCosmicTree->SetDirectory(0);
+ }
+ if (fCosmicTree->GetEntries()==0){
+ //
+ fCosmicTree->SetDirectory(0);
+ fCosmicTree->Branch("t0.",&track0);
+ fCosmicTree->Branch("t1.",&track1);
+ fCosmicTree->Branch("ft0.",&friendTrack0);
+ fCosmicTree->Branch("ft1.",&friendTrack1);
+ }else{
+ fCosmicTree->SetBranchAddress("t0.",&track0);
+ fCosmicTree->SetBranchAddress("t1.",&track1);
+ fCosmicTree->SetBranchAddress("ft0.",&friendTrack0);
+ fCosmicTree->SetBranchAddress("ft1.",&friendTrack1);
+ }
+ fCosmicTree->Fill();
+ }
+ }
+}
+
+
+void AliTPCcalibCosmic::Terminate(){
+ //
+ // copy the cosmic tree to memory resident tree
+ //
+ static Int_t counter=0;
+ printf("AliTPCcalibCosmic::Terminate\t%d\n",counter);
+ counter++;
+ AliTPCcalibBase::Terminate();
+}
+
+void AliTPCcalibCosmic::AddTree(TTree* treeOutput, TTree * treeInput){
+ //
+ // Add the content of tree:
+ // Notice automatic copy of tree in ROOT does not work for such complicated tree
+ //
+ AliESDtrack *track0=new AliESDtrack;
+ AliESDtrack *track1=new AliESDtrack;
+ AliESDfriendTrack *ftrack0=new AliESDfriendTrack;
+ AliESDfriendTrack *ftrack1=new AliESDfriendTrack;
+ treeInput->SetBranchAddress("t0.",&track0);
+ treeInput->SetBranchAddress("t1.",&track1);
+ treeInput->SetBranchAddress("ft0.",&ftrack0);
+ treeInput->SetBranchAddress("ft1.",&ftrack1);
+ if (treeOutput->GetEntries()==0){
+ //
+ treeOutput->SetDirectory(0);
+ treeOutput->Branch("t0.",&track0);
+ treeOutput->Branch("t1.",&track1);
+ treeOutput->Branch("ft0.",&ftrack0);
+ treeOutput->Branch("ft1.",&ftrack1);
+ }else{
+ treeOutput->SetBranchAddress("t0.",&track0);
+ treeOutput->SetBranchAddress("t1.",&track1);
+ treeOutput->SetBranchAddress("ft0.",&ftrack0);
+ treeOutput->SetBranchAddress("ft1.",&ftrack1);
+ }
+ Int_t entries= treeInput->GetEntries();
+ for (Int_t i=0; i<entries; i++){
+ treeInput->GetEntry(i);
+ treeOutput->Fill();
+ }
+}
+
+
+
+void AliTPCcalibCosmic::MakeFitTree(TTree * treeInput, TTreeSRedirector *pcstream, const TObjArray * corrArray, Int_t step, Int_t run){
+ //
+ // Make fit tree
+ // refit the tracks with original points + corrected points for each correction
+ // Input:
+ // treeInput - tree with cosmic tracks
+ // pcstream - debug output
+
+ // Algorithm:
+ // Loop over pair of cosmic tracks:
+ // 1. Find trigger offset between cosmic event and "physic" trigger
+ // 2. Refit tracks with current transformation
+ // 3. Refit tracks using additional "primitive" distortion on top
+ // Best correction estimated as linear combination of corrections
+ // minimizing the observed distortions
+ // Observed distortions - matching in the y,z, snp, theta and 1/pt
+ //
+ const Double_t kResetCov=20.;
+ const Double_t kMaxDelta[5]={1,40,0.03,0.01,0.2};
+ const Double_t kSigma=2.;
+ const Double_t kMaxTime=1050;
+ const Double_t kMaxSnp=0.8;
+ Int_t 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());
+ transform->SetCurrentRun(run);
+ transform->SetCurrentTimeStamp(time);
+ Double_t covar[15];
+ for (Int_t i=0;i<15;i++) covar[i]=0;
+ covar[0]=kSigma*kSigma;
+ covar[2]=kSigma*kSigma;
+ covar[5]=kSigma*kSigma/Float_t(150*150);
+ covar[9]=kSigma*kSigma/Float_t(150*150);
+ covar[14]=0.2*0.2;
+ Double_t *distortions = new Double_t[ncorr+1];
+
+ AliESDtrack *track0=new AliESDtrack;
+ AliESDtrack *track1=new AliESDtrack;
+ AliESDfriendTrack *ftrack0=new AliESDfriendTrack;
+ AliESDfriendTrack *ftrack1=new AliESDfriendTrack;
+ treeInput->SetBranchAddress("t0.",&track0);
+ treeInput->SetBranchAddress("t1.",&track1);
+ treeInput->SetBranchAddress("ft0.",&ftrack0);
+ treeInput->SetBranchAddress("ft1.",&ftrack1);
+ Int_t entries= treeInput->GetEntries();
+ for (Int_t i=0; i<entries; i+=step){
+ treeInput->GetEntry(i);
+ if (i%20==0) printf("%d\n",i);
+ if (!ftrack0->GetTPCOut()) continue;
+ if (!ftrack1->GetTPCOut()) continue;
+ AliTPCseed *seed0=0;
+ AliTPCseed *seed1=0;
+ TObject *calibObject;
+ for (Int_t l=0;(calibObject=ftrack0->GetCalibObject(l));++l) {
+ if ((seed0=dynamic_cast<AliTPCseed*>(calibObject))) break;
+ }
+ for (Int_t l=0;(calibObject=ftrack1->GetCalibObject(l));++l) {
+ if ((seed1=dynamic_cast<AliTPCseed*>(calibObject))) break;
+ }
+ if (!seed0) continue;
+ if (!seed1) continue;
+ if (TMath::Abs(seed0->GetSnp())>kMaxSnp) continue;
+ if (TMath::Abs(seed1->GetSnp())>kMaxSnp) continue;
+ //
+ //
+ Int_t nclA0=0, nclC0=0; // number of clusters
+ Int_t nclA1=0, nclC1=0; // number of clusters
+ Int_t ncl0=0,ncl1=0;
+ Double_t rmin0=300, rmax0=-300; // variables to estimate the time0 - trigger offset
+ Double_t rmin1=300, rmax1=-300;
+ Double_t tmin0=2000, tmax0=-2000;
+ Double_t tmin1=2000, tmax1=-2000;
+ //
+ //
+ // calculate trigger offset usig "missing clusters"
+ for (Int_t irow=0; irow<159; irow++){
+ AliTPCclusterMI *cluster0=seed0->GetClusterPointer(irow);
+ if (cluster0 &&cluster0->GetX()>10){
+ if (cluster0->GetX()<rmin0) { rmin0=cluster0->GetX(); tmin0=cluster0->GetTimeBin();}
+ if (cluster0->GetX()>rmax0) { rmax0=cluster0->GetX(); tmax0=cluster0->GetTimeBin();}
+ ncl0++;
+ if (cluster0->GetDetector()%36<18) nclA0++;
+ if (cluster0->GetDetector()%36>=18) nclC0++;
+ }
+ AliTPCclusterMI *cluster1=seed1->GetClusterPointer(irow);
+ if (cluster1&&cluster1->GetX()>10){
+ if (cluster1->GetX()<rmin1) { rmin1=cluster1->GetX(); tmin1=cluster1->GetTimeBin();}
+ if (cluster1->GetX()>rmax1) { rmax1=cluster1->GetX(); tmax1=cluster1->GetTimeBin();}
+ ncl1++;
+ if (cluster1->GetDetector()%36<18) nclA1++;
+ if (cluster1->GetDetector()%36>=18) nclC1++;
+ }
+ }
+ Int_t cosmicType=0; // types of cosmic topology
+ if ((nclA0>nclC0) && (nclA1>nclC1)) cosmicType=0; // AA side
+ if ((nclA0<nclC0) && (nclA1<nclC1)) cosmicType=1; // CC side
+ if ((nclA0>nclC0) && (nclA1<nclC1)) cosmicType=2; // AC side
+ if ((nclA0<nclC0) && (nclA1>nclC1)) cosmicType=3; // CA side
+ //if ((nclA0>nclC0) && (nclA1<nclC1)) cosmicType=6; // AC side out of time
+ //if ((nclA0>nclC0) && (nclA1<nclC1)) cosmicType=7; // CA side out of time
+ //
+ Double_t deltaTime = 0; // correction for trigger not in time - equalizing the time arival
+ if ((cosmicType>1)&&TMath::Abs(track1->GetZ()-track0->GetZ())>4){
+ cosmicType+=4;
+ deltaTime=0.5*(track1->GetZ()-track0->GetZ())/param->GetZWidth();
+ if (nclA0>nclC0) deltaTime*=-1; // if A side track
+ }
+ //
+ TVectorD vectorDT(8);
+ Int_t crossCounter=0;
+ Double_t deltaTimeCross = AliTPCcalibCosmic::GetDeltaTime(rmin0, rmax0, rmin1, rmax1, tmin0, tmax0, tmin1, tmax1, TMath::Abs(track0->GetY()),vectorDT);
+ Bool_t isOKTrigger=kTRUE;
+ for (Int_t ic=0; ic<6;ic++) {
+ if (TMath::Abs(vectorDT[ic])>0) {
+ if (vectorDT[ic]+vectorDT[6]<0) isOKTrigger=kFALSE;
+ if (vectorDT[ic]+vectorDT[7]>kMaxTime) isOKTrigger=kFALSE;
+ if (isOKTrigger){
+ crossCounter++;
+ }
+ }
+ }
+ Double_t deltaTimeCluster=deltaTime;
+ if ((cosmicType==0 || cosmicType==1) && crossCounter>0){
+ deltaTimeCluster=deltaTimeCross;
+ cosmicType+=8;
+ }
+ if (nclA0*nclC0>0 || nclA1*nclC1>0) cosmicType+=16; // mixed A side C side - bad for visualization
+ //
+ // Apply current transformation
+ //
+ //
+ for (Int_t irow=0; irow<159; irow++){
+ AliTPCclusterMI *cluster0=seed0->GetClusterPointer(irow);
+ if (cluster0 &&cluster0->GetX()>10){
+ Double_t x0[3]={cluster0->GetRow(),cluster0->GetPad(),cluster0->GetTimeBin()+deltaTimeCluster};
+ Int_t index0[1]={cluster0->GetDetector()};
+ transform->Transform(x0,index0,0,1);
+ cluster0->SetX(x0[0]);
+ cluster0->SetY(x0[1]);
+ cluster0->SetZ(x0[2]);
+ //
+ }
+ AliTPCclusterMI *cluster1=seed1->GetClusterPointer(irow);
+ if (cluster1&&cluster1->GetX()>10){
+ Double_t x1[3]={cluster1->GetRow(),cluster1->GetPad(),cluster1->GetTimeBin()+deltaTimeCluster};
+ Int_t index1[1]={cluster1->GetDetector()};
+ transform->Transform(x1,index1,0,1);
+ cluster1->SetX(x1[0]);
+ cluster1->SetY(x1[1]);
+ cluster1->SetZ(x1[2]);
+ }
+ }
+ //
+ //
+ Double_t alpha=track0->GetAlpha(); // rotation frame
+ Double_t cos = TMath::Cos(alpha);
+ Double_t sin = TMath::Sin(alpha);
+ Double_t mass = TDatabasePDG::Instance()->GetParticle("mu+")->Mass();
+ AliExternalTrackParam btrack0=*(ftrack0->GetTPCOut());
+ AliExternalTrackParam btrack1=*(ftrack1->GetTPCOut());
+ btrack0.Rotate(alpha);
+ btrack1.Rotate(alpha);
+ // change the sign for track 1
+ Double_t* par1=(Double_t*)btrack0.GetParameter();
+ par1[3]*=-1;
+ par1[4]*=-1;
+ btrack0.AddCovariance(covar);
+ btrack1.AddCovariance(covar);
+ btrack0.ResetCovariance(kResetCov);
+ btrack1.ResetCovariance(kResetCov);
+ Bool_t isOK=kTRUE;
+ Bool_t isOKT=kTRUE;
+ TObjArray tracks0(ncorr+1);
+ TObjArray tracks1(ncorr+1);
+ //
+ Double_t dEdx0Tot=seed0->CookdEdxAnalytical(0.02,0.6,kTRUE);
+ Double_t dEdx1Tot=seed1->CookdEdxAnalytical(0.02,0.6,kTRUE);
+ Double_t dEdx0Max=seed0->CookdEdxAnalytical(0.02,0.6,kFALSE);
+ Double_t dEdx1Max=seed1->CookdEdxAnalytical(0.02,0.6,kFALSE);
+ //if (TMath::Abs((dEdx0Max+1)/(dEdx0Tot+1)-1.)>0.1) isOK=kFALSE;
+ //if (TMath::Abs((dEdx1Max+1)/(dEdx1Tot+1)-1.)>0.1) isOK=kFALSE;
+ ncl0=0; ncl1=0;
+ for (Int_t icorr=-1; icorr<ncorr; icorr++){
+ AliExternalTrackParam rtrack0=btrack0;
+ AliExternalTrackParam rtrack1=btrack1;
+ AliTPCCorrection *corr = 0;
+ if (icorr>=0) corr = (AliTPCCorrection*)corrArray->At(icorr);
+ //
+ for (Int_t irow=159; irow>0; irow--){
+ AliTPCclusterMI *cluster=seed0->GetClusterPointer(irow);
+ if (!cluster) continue;
+ if (!isOKT) break;
+ Double_t rD[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
+ transform->RotatedGlobal2Global(cluster->GetDetector()%36,rD); // transform to global
+ Float_t r[3]={rD[0],rD[1],rD[2]};
+ if (corr){
+ corr->DistortPoint(r, cluster->GetDetector());
+ }
+ //
+ Double_t cov[3]={0.01,0.,0.01};
+ Double_t lx =cos*r[0]+sin*r[1];
+ Double_t ly =-sin*r[0]+cos*r[1];
+ rD[1]=ly; rD[0]=lx; rD[2]=r[2]; //transform to track local
+ if (!AliTracker::PropagateTrackToBxByBz(&rtrack0, lx,mass,1.,kFALSE)) isOKT=kFALSE;;
+ if (!rtrack0.Update(&rD[1],cov)) isOKT =kFALSE;
+ if (icorr<0) ncl0++;
+ }
+ //
+ for (Int_t irow=159; irow>0; irow--){
+ AliTPCclusterMI *cluster=seed1->GetClusterPointer(irow);
+ if (!cluster) continue;
+ if (!isOKT) break;
+ Double_t rD[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
+ transform->RotatedGlobal2Global(cluster->GetDetector()%36,rD);
+ Float_t r[3]={rD[0],rD[1],rD[2]};
+ if (corr){
+ corr->DistortPoint(r, cluster->GetDetector());
+ }
+ //
+ Double_t cov[3]={0.01,0.,0.01};
+ Double_t lx =cos*r[0]+sin*r[1];
+ Double_t ly =-sin*r[0]+cos*r[1];
+ rD[1]=ly; rD[0]=lx; rD[2]=r[2];
+ if (!AliTracker::PropagateTrackToBxByBz(&rtrack1, lx,mass,1.,kFALSE)) isOKT=kFALSE;
+ if (!rtrack1.Update(&rD[1],cov)) isOKT=kFALSE;
+ if (icorr<0) ncl1++;
+ }
+ if (!AliTracker::PropagateTrackToBxByBz(&rtrack0, 0,mass,10.,kFALSE)) isOKT=kFALSE;
+ if (!AliTracker::PropagateTrackToBxByBz(&rtrack1, 0,mass,10.,kFALSE)) isOKT=kFALSE;
+ if (!AliTracker::PropagateTrackToBxByBz(&rtrack0, 0,mass,1.,kFALSE)) isOKT=kFALSE;
+ if (!AliTracker::PropagateTrackToBxByBz(&rtrack1, 0,mass,1.,kFALSE)) isOKT=kFALSE;
+ const Double_t *param0=rtrack0.GetParameter();
+ const Double_t *param1=rtrack1.GetParameter();
+ for (Int_t ipar=0; ipar<4; ipar++){
+ if (TMath::Abs(param1[ipar]-param0[ipar])>kMaxDelta[ipar]) isOK=kFALSE;
+ }
+ tracks0.AddAt(rtrack0.Clone(), icorr+1);
+ tracks1.AddAt(rtrack1.Clone(), icorr+1);
+ AliExternalTrackParam out0=*(ftrack0->GetTPCOut());
+ AliExternalTrackParam out1=*(ftrack1->GetTPCOut());
+ Int_t nentries=TMath::Min(ncl0,ncl1);
+
+ if (icorr<0) {
+ (*pcstream)<<"cosmic"<<
+ "isOK="<<isOK<< // correct all propagation update and also residuals
+ "isOKT="<<isOKT<< // correct all propagation update
+ "isOKTrigger="<<isOKTrigger<< // correct? estimate of trigger offset
+ "id="<<cosmicType<<
+ //
+ //
+ "cross="<<crossCounter<<
+ "vDT.="<<&vectorDT<<
+ //
+ "dTime="<<deltaTime<< // delta time using the A-c side cross
+ "dTimeCross="<<deltaTimeCross<< // delta time using missing clusters
+ //
+ "dEdx0Max="<<dEdx0Max<<
+ "dEdx0Tot="<<dEdx0Tot<<
+ "dEdx1Max="<<dEdx1Max<<
+ "dEdx1Tot="<<dEdx1Tot<<
+ //
+ "track0.="<<track0<< // original track 0
+ "track1.="<<track1<< // original track 1
+ "out0.="<<&out0<< // outer track 0
+ "out1.="<<&out1<< // outer track 1
+ "rtrack0.="<<&rtrack0<< // refitted track with current transform
+ "rtrack1.="<<&rtrack1<< //
+ "ncl0="<<ncl0<<
+ "ncl1="<<ncl1<<
+ "entries="<<nentries<< // number of clusters
+ "\n";
+ }
+ }
+ //
+
+ if (isOK){
+ Int_t nentries=TMath::Min(ncl0,ncl1);
+ for (Int_t ipar=0; ipar<5; ipar++){
+ for (Int_t icorr=-1; icorr<ncorr; icorr++){
+ AliTPCCorrection *corr = 0;
+ if (icorr>=0) corr = (AliTPCCorrection*)corrArray->At(icorr);
+ //
+ AliExternalTrackParam *param0=(AliExternalTrackParam *) tracks0.At(icorr+1);
+ AliExternalTrackParam *param1=(AliExternalTrackParam *) tracks1.At(icorr+1);
+ distortions[icorr+1]=param1->GetParameter()[ipar]-param0->GetParameter()[ipar];
+ if (icorr>=0){
+ distortions[icorr+1]-=distortions[0];
+ }
+ //
+ if (icorr<0){
+ Double_t bz=AliTrackerBase::GetBz();
+ Double_t gxyz[3];
+ param0->GetXYZ(gxyz);
+ Int_t dtype=20;
+ Double_t theta=param0->GetParameter()[3];
+ Double_t phi = alpha;
+ Double_t snp = track0->GetInnerParam()->GetSnp();
+ Double_t mean= distortions[0];
+ Int_t index = param0->GetIndex(ipar,ipar);
+ Double_t rms=TMath::Sqrt(param1->GetCovariance()[index]+param1->GetCovariance()[index]);
+ if (crossCounter<1) rms*=1;
+ Double_t sector=9*phi/TMath::Pi();
+ Double_t dsec = sector-TMath::Nint(sector+0.5);
+ Double_t gx=gxyz[0],gy=gxyz[1],gz=gxyz[2];
+ Double_t refX=TMath::Sqrt(gx*gx+gy*gy);
+ Double_t dRrec=0;
+ // Double_t pt=(param0->GetSigned1Pt()+param1->GetSigned1Pt())*0.5;
+ Double_t pt=(param0->GetSigned1Pt()+param1->GetSigned1Pt())*0.5;
+
+ (*pcstream)<<"fit"<< // dump valus for fit
+ "run="<<run<< //run number
+ "bz="<<bz<< // magnetic filed used
+ "dtype="<<dtype<< // detector match type 20
+ "ptype="<<ipar<< // parameter type
+ "theta="<<theta<< // theta
+ "phi="<<phi<< // phi
+ "snp="<<snp<< // snp
+ "mean="<<mean<< // mean dist value
+ "rms="<<rms<< // rms
+ "sector="<<sector<<
+ "dsec="<<dsec<<
+ //
+ "refX="<<refX<< // reference radius
+ "gx="<<gx<< // global position
+ "gy="<<gy<< // global position
+ "gz="<<gz<< // global position
+ "dRrec="<<dRrec<< // delta Radius in reconstruction
+ "pt="<<pt<< //1/pt
+ "id="<<cosmicType<< //type of the cosmic used
+ "entries="<<nentries;// number of entries in bin
+ (*pcstream)<<"cosmicDebug"<<
+ "p0.="<<param0<< // dump distorted track 0
+ "p1.="<<param1; // dump distorted track 1
+ }
+ if (icorr>=0){
+ (*pcstream)<<"fit"<<
+ Form("%s=",corr->GetName())<<distortions[icorr+1]; // dump correction value
+ (*pcstream)<<"cosmicDebug"<<
+ Form("%s=",corr->GetName())<<distortions[icorr+1]<< // dump correction value
+ Form("%sp0.=",corr->GetName())<<param0<< // dump distorted track 0
+ Form("%sp1.=",corr->GetName())<<param1; // dump distorted track 1
+ }
+ } //loop corrections
+ (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n";
+ (*pcstream)<<"cosmicDebug"<<"isOK="<<isOK<<"\n";
+ } //loop over parameters
+ } // dump results
+ }//loop tracks
+}
+
+
+
+Double_t AliTPCcalibCosmic::GetDeltaTime(Double_t rmin0, Double_t rmax0, Double_t rmin1, Double_t rmax1, Double_t tmin0, Double_t tmax0, Double_t tmin1, Double_t tmax1, Double_t dcaR, TVectorD &vectorDT)
+{
+ //
+ // Estimate trigger offset between random cosmic event and "physics" trigger
+ // Efficiency about 50 % of cases:
+ // Cases:
+ // 0. Tracks crossing A side and C side - no match in z - 30 % of cases
+ // 1. Track only on one side and dissapear at small or at high radiuses - 50 % of cases
+ // 1.a) rmax<Rc && tmax<Tcmax && tmax>tmin => deltaT=Tcmax-tmax
+ // 1.b) rmin>Rcmin && tmin<Tcmax && tmin>tmax => deltaT=Tcmax-tmin
+ // // case the z matching gives proper time
+ // 1.c) rmax<Rc && tmax>Tcmin && tmax<tmin => deltaT=-tmax
+ //
+ // check algorithm:
+ // TCut cutStop = "(min(rmax0,rmax1)<235||abs(rmin0-rmin1)>10)"; // tracks not registered for full time
+
+ // Combinations:
+ // 0-1 - forbidden
+ // 0-2 - forbidden
+ // 0-3 - occur - wrong correlation
+ // 1-2 - occur - wrong correlation
+ // 1-3 - forbidden
+ // 2-3 - occur - small number of outlyers -20%
+ // Frequency:
+ // 0 - 106
+ // 1 - 265
+ // 2 - 206
+ // 3 - 367
+ //
+ const Double_t kMaxRCut=235; // max radius
+ const Double_t kMinRCut=TMath::Max(dcaR,90.); // min radius
+ const Double_t kMaxDCut=30; // max distance for minimal radius
+ const Double_t kMinTime=110;
+ const Double_t kMaxTime=950;
+ Double_t deltaT=0;
+ Int_t counter=0;
+ vectorDT[6]=TMath::Min(TMath::Min(tmin0,tmin1),TMath::Min(tmax0,tmax1));
+ vectorDT[7]=TMath::Max(TMath::Max(tmin0,tmin1),TMath::Max(tmax0,tmax1));
+ if (TMath::Min(rmax0,rmax1)<kMaxRCut){
+ // max cross - deltaT>0
+ if (rmax0<kMaxRCut && tmax0 <kMaxTime && tmax0>tmin0) vectorDT[0]=kMaxTime-tmax0; // disapear at CE
+ if (rmax1<kMaxRCut && tmax1 <kMaxTime && tmax1>tmin1) vectorDT[1]=kMaxTime-tmax1; // disapear at CE
+ // min cross - deltaT<0 - OK they are correlated
+ if (rmax0<kMaxRCut && tmax0 >kMinTime && tmax0<tmin0) vectorDT[2]=-tmax0; // disapear at ROC
+ if (rmax1<kMaxRCut && tmax1 >kMinTime && tmax1<tmin1) vectorDT[3]=-tmax1; // disapear at ROC
+ }
+ if (rmin0> kMinRCut+kMaxDCut && tmin0 <kMaxTime && tmin0>tmax0) vectorDT[4]=kMaxTime-tmin0;
+ if (rmin1> kMinRCut+kMaxDCut && tmin1 <kMaxTime && tmin1>tmax1) vectorDT[5]=kMaxTime-tmin1;
+ Bool_t isOK=kTRUE;
+ for (Int_t i=0; i<6;i++) {
+ if (TMath::Abs(vectorDT[i])>0) {
+ counter++;
+ if (vectorDT[i]+vectorDT[6]<0) isOK=kFALSE;
+ if (vectorDT[i]+vectorDT[7]>kMaxTime) isOK=kFALSE;
+ if (isOK) deltaT=vectorDT[i];
+ }
+ }
+ return deltaT;
+}
//
void Init();
void FindPairs(AliESDEvent *event);
+ void FindCosmicPairs(AliESDEvent * event);
+
Bool_t IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1);
static void CalculateBetheParams(TH2F *hist, Double_t * initialParam);
static Double_t CalculateMIPvalue(TH1F * hist);
void UpdateTrack(AliExternalTrackParam &track0, const AliExternalTrackParam &track1);
//
- void FillHistoPerformance(const AliExternalTrackParam *par0, const AliExternalTrackParam *par1, const AliExternalTrackParam *inner0, const AliExternalTrackParam *inner1, AliTPCseed *seed0, AliTPCseed *seed1, const AliExternalTrackParam *param0Combined);
+ void FillHistoPerformance(const AliExternalTrackParam *par0, const AliExternalTrackParam *par1, const AliExternalTrackParam *inner0, const AliExternalTrackParam *inner1, AliTPCseed *seed0, AliTPCseed *seed1, const AliExternalTrackParam *param0Combined, Int_t cross);
void MaterialBudgetDump(AliExternalTrackParam *const par0, AliExternalTrackParam *const par1, const AliExternalTrackParam *inner0, const AliExternalTrackParam *inner1, AliTPCseed *const seed0, AliTPCseed *const seed1, AliExternalTrackParam *const param0Combined, AliExternalTrackParam *const param1Combined);
-
-
+ static void MakeFitTree(TTree * treeInput, TTreeSRedirector *pcstream, const TObjArray * corrArray, Int_t step, Int_t run);
+ TTree * GetCosmicTree() const {return fCosmicTree;}
//
TH1F * GetHistNTracks() const {return fHistNTracks;};
TH1F * GetHistClusters() const {return fClusters;};
void Process(AliESDtrack *const track, Int_t runNo=-1) {AliTPCcalibBase::Process(track,runNo);};
void Process(AliTPCseed *const track) {return AliTPCcalibBase::Process(track);}
-
+ virtual void Terminate();
+ static Double_t GetDeltaTime(Double_t rmin0, Double_t rmax0, Double_t rmin1, Double_t rmax1, Double_t tmin0, Double_t tmax0, Double_t tmin1, Double_t tmax1, Double_t dcaR, TVectorD& vectorDT);
public:
//
// Performance histograms
THnSparse *fHistoPull[6]; // histograms of tracking performance pull
THnSparse *fHistodEdxMax[4]; // histograms of dEdx perfomance - max charge
THnSparse *fHistodEdxTot[4]; // histograms of dEdx perfomance - tot charge
-
+ static void AddTree(TTree* treeOutput, TTree * treeInput);
private:
void FillAcordeHist(AliESDtrack *upperTrack);
Float_t fCutTheta; // maximal distance in theta ditection
Float_t fCutMinDir; // direction vector products
+ TTree *fCosmicTree; // tree with the cosmic tracks
AliTPCcalibCosmic(const AliTPCcalibCosmic&);
AliTPCcalibCosmic& operator=(const AliTPCcalibCosmic&);
- ClassDef(AliTPCcalibCosmic, 2);
+ ClassDef(AliTPCcalibCosmic, 3);
};
#endif
#include "AliDCSSensorArray.h"
#include "AliDCSSensor.h"
#include "AliGRPObject.h"
+#include "AliTPCROC.h"
using namespace std;
fSignals(336),
//
fHisLaser(0), // N dim histogram of laser
+ fHisLaserPad(0), // N dim histogram of laser
+ fHisLaserTime(0), // N dim histogram of laser
fHisNclIn(0), //->Number of clusters inner
fHisNclOut(0), //->Number of clusters outer
fHisNclIO(0), //->Number of cluster inner outer
//
//
fHisLaser(0), // N dim histogram of laser
+ fHisLaserPad(0), // N dim histogram of laser
+ fHisLaserTime(0), // N dim histogram of laser
+
fHisNclIn(0), //->Number of clusters inner
fHisNclOut(0), //->Number of clusters outer
fHisNclIO(0), //->Number of cluster inner outer
//
//
fHisLaser(0), // N dim histogram of laser
+ fHisLaserPad(0), // N dim histogram of laser
+ fHisLaserTime(0), // N dim histogram of laser
+
fHisNclIn(new TH2F(*(calibLaser.fHisNclIn))), //->Number of clusters inner
fHisNclOut(new TH2F(*(calibLaser.fHisNclOut))), //->Number of clusters outer
fHisNclIO(new TH2F(*(calibLaser.fHisNclIO))), //->Number of cluster inner outer
//
if ( fHisNclIn){
delete fHisLaser; //->
+ delete fHisLaserPad; //->
+ delete fHisLaserTime; //->
+
delete fHisNclIn; //->Number of clusters inner
delete fHisNclOut; //->Number of clusters outer
delete fHisNclIO; //->Number of cluster inner outer
xhis[11]= TMath::Sqrt(TMath::Abs(track->GetTPCsignal()));
// }
fHisLaser->Fill(xhis);
+ //
+
}
void AliTPCcalibLaser::FitDriftV(){
}
}
+ //
+ // Fill raw THnSparses
+ //
+ for (Int_t irow=0;irow<159;irow++) {
+ AliTPCclusterMI *c=track->GetClusterPointer(irow);
+ if (!c) continue;
+ if (c->GetMax()>800) continue; // saturation cut
+ //if (TMath::Sqrt(TMath::Abs(c->GetSigmaY2()))>1) continue; // saturation cut
+ //
+ Double_t deltaY=c->GetY()-(*ltrp->GetVecLY())[irow];
+ Double_t deltaZ=c->GetZ()-(*ltrp->GetVecLZ())[irow];
+ //TString axisName[6]={"Delta","bin", "rms shape", "Q", "row","trackID"}
+ Double_t xyz[6]={0, 0, 0,TMath::Sqrt(c->GetMax()),irow,id};
+ xyz[0]=deltaY;
+ xyz[1]=c->GetPad();
+ xyz[2]=TMath::Sqrt(TMath::Abs(c->GetSigmaY2()));
+ fHisLaserPad->Fill(xyz);
+ xyz[0]=deltaZ;
+ xyz[1]=c->GetTimeBin();
+ xyz[2]=TMath::Sqrt(TMath::Abs(c->GetSigmaZ2()));
+ fHisLaserTime->Fill(xyz);
+ }
}
for(Int_t bin = 1; bin<=159; bin++) {
- if(yprof->GetBinEntries(bin)<10&&
- zprof->GetBinEntries(bin)<10) {
+ if(yprof->GetBinEntries(bin)<5&&
+ zprof->GetBinEntries(bin)<5) {
continue;
}
vecX[bin-1] = x;
vecDY[bin-1] = yprof->GetBinContent(bin);
vecDZ[bin-1] = zprof->GetBinContent(bin);
+ if (bin>0&&bin<159){
+ //
+ //truncated mean - skip first and the last pad row
+ //
+ Int_t firstBin=TMath::Max(bin-5,0);
+ Int_t lastBin =TMath::Min(bin+5,158);
+ histAbsY->GetXaxis()->SetRangeUser(firstBin,lastBin);
+ histAbsY->GetYaxis()->SetRangeUser(-2,2);
+ vecEy[bin-1]=histAbsY->GetRMS(2);
+ vecDY[bin-1]=histAbsY->GetMean(2);
+ histAbsY->GetXaxis()->SetRangeUser(firstBin+2,lastBin-2);//use+-2 bins
+ histAbsY->GetYaxis()->SetRangeUser(vecDY[bin-1]-4*vecEy[bin-1],vecDY[bin-1]+4*vecEy[bin-1]);
+ if (yprof->GetBinEntries(bin-1)>0) vecEy[bin-1]=histAbsY->GetRMS(2)/TMath::Sqrt(yprof->GetBinEntries(bin-1));
+ vecDY[bin-1]=histAbsY->GetMean(2);
+ }
+
if(!isOuter) { // inner
vecPhi[bin-1]=lasTanPhiLocIn;
// Calculate local y from residual and database
// merge fDeltaZ histograms
hm = (TH1F*)cal->fDeltaZ.At(id);
h = (TH1F*)fDeltaZ.At(id);
- if (!h) {
- h=new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
+ if (!h &&hm &&hm->GetEntries()>0) {
+ h=(TH1F*)hm->Clone();
h->SetDirectory(0);
fDeltaZ.AddAt(h,id);
}
// merge fP3 histograms
hm = (TH1F*)cal->fDeltaP3.At(id);
h = (TH1F*)fDeltaP3.At(id);
- if (!h) {
- h=new TH1F(Form("hisPar3v%d",id),Form("hisPar3v%d",id),400,-0.06,0.06);
+ if (!h&&hm &&hm->GetEntries()>0) {
+ h=(TH1F*)hm->Clone();
h->SetDirectory(0);
fDeltaP3.AddAt(h,id);
}
// merge fP4 histograms
hm = (TH1F*)cal->fDeltaP4.At(id);
h = (TH1F*)fDeltaP4.At(id);
- if (!h) {
- h=new TH1F(Form("hisPar4v%d",id),Form("hisPar4v%d",id),200,-0.06,0.06);
+ if (!h &&hm &&hm->GetEntries()>0) {
+ h=(TH1F*)hm->Clone();
h->SetDirectory(0);
fDeltaP4.AddAt(h,id);
}
// merge fDeltaPhi histograms
hm = (TH1F*)cal->fDeltaPhi.At(id);
h = (TH1F*)fDeltaPhi.At(id);
- if (!h) {
- h= new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
+ if (!h &&hm &&hm->GetEntries()>0) {
+ h= (TH1F*)hm->Clone();
h->SetDirectory(0);
fDeltaPhi.AddAt(h,id);
}
// merge fDeltaPhiP histograms
hm = (TH1F*)cal->fDeltaPhiP.At(id);
h = (TH1F*)fDeltaPhiP.At(id);
- if (!h) {
- h=new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01);
+ if (!h&&hm &&hm->GetEntries()>0) {
+ h=(TH1F*)hm->Clone();
h->SetDirectory(0);
fDeltaPhiP.AddAt(h,id);
}
// merge fSignals histograms
hm = (TH1F*)cal->fSignals.At(id);
h = (TH1F*)fSignals.At(id);
- if (!h) {
- h=new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),100,0,300);
+ if (!h&&hm &&hm->GetEntries()>0) {
+ h=(TH1F*)hm->Clone();
h->SetDirectory(0);
fSignals.AddAt(h,id);
}
h2m = (TH2F*)cal->fDeltaYresAbs.At(id);
h2 = (TH2F*)fDeltaYresAbs.At(id);
if (h2m&&h2) h2->Add(h2m);
-
+ if (h2m&&!h2) { h2=(TH2F*)h2m->Clone(); h2->SetDirectory(0); fDeltaYresAbs.AddAt(h2,id);}
h2m = (TH2F*)cal->fDeltaZresAbs.At(id);
h2 = (TH2F*)fDeltaZresAbs.At(id);
if (h2m&&h2) h2->Add(h2m);
+ if (h2m&&!h2) { h2=(TH2F*)h2m->Clone(); h2->SetDirectory(0); fDeltaZresAbs.AddAt(h2,id);}
// merge ProfileY histograms - 3
//h2m = (TH2F*)cal->fDeltaYres3.At(id);
//h2 = (TH2F*)fDeltaYres3.At(id);
fHisLaser->GetAxis(iaxis)->SetName(nameLaser[iaxis]);
fHisLaser->GetAxis(iaxis)->SetTitle(titleLaser[iaxis]);
}
+ //
+ // Delta Time bin
+ // Pad SigmaShape Q charge pad row trackID
+ Int_t binsRow[6]={200, 10000, 20, 30, 159, 336};
+ Double_t axisMin[6]={-1, 0, 0, 1, 0 , 0};
+ Double_t axisMax[6]={ 1, 1000, 1, 30, 159, 336};
+ TString axisName[6]={"Delta","bin", "rms shape", "sqrt(Q)", "row","trackID"};
+
+ binsRow[1]=2000;
+ axisMin[1]=0;
+ axisMax[1]=200;
+ fHisLaserPad = new THnSparseS("laserPad","#Delta_{Laser}", 6, binsRow,axisMin, axisMax);
+ //
+ binsRow[0]=1000;
+ axisMin[0]=-20;
+ axisMax[0]=20;
+ binsRow[1]=10000;
+ axisMin[1]=0;
+ axisMax[1]=1000;
+ //
+ fHisLaserTime= new THnSparseS("laserTime","#Delta_{Laser}", 6, binsRow,axisMin, axisMax);
+ //
+ for (Int_t iaxis=0; iaxis<6; iaxis++){
+ fHisLaserPad->GetAxis(iaxis)->SetName(axisName[iaxis]);
+ fHisLaserTime->GetAxis(iaxis)->SetTitle(axisName[iaxis]);
+ }
}
void AliTPCcalibLaser::MergeFitHistos(AliTPCcalibLaser * laser){
//
// Only first histogram is checked - all other should be the same
if (fHisLaser &&laser->fHisLaser) fHisLaser->Add(laser->fHisLaser);
-
- if (!laser->fHisNclIn) return; // empty histograms
+ if (fHisLaserPad &&laser->fHisLaserPad) fHisLaserPad->Add(laser->fHisLaserPad);
+ if (!fHisLaserPad &&laser->fHisLaserPad) fHisLaserPad=(THnSparseS*)laser->fHisLaserPad->Clone();
+ if (fHisLaserTime &&laser->fHisLaserTime) fHisLaserTime->Add(laser->fHisLaserTime);
+ if (!fHisLaserTime &&laser->fHisLaserTime) fHisLaserTime=(THnSparseS*)laser->fHisLaserTime->Clone();
+
+ if (!laser->fHisNclIn) laser->MakeFitHistos(); // empty histograms
if (!fHisNclIn) MakeFitHistos();
+ if (fHisNclIn->GetEntries()<1) MakeFitHistos();
//
fHisNclIn->Add(laser->fHisNclIn ); //->Number of clusters inner
fHisNclOut->Add(laser->fHisNclOut ); //->Number of clusters outer
}
- /*
- TFile f("vscan.root");
- */
-
- /*
- pad binning effect
- chain->Draw("Cl[].fPad-int(Cl[].fPad)",cutA+cutCl+"Cl[].fZ>0&&Cl[].fPad>1","",10000);
- //
- chain->Draw("Cl[].fY-TrYpol1.fElements:Cl[].fPad-int(Cl[].fPad)",cutA+cutCl+"Cl[].fZ>0&&Cl[].fPad>1","prof",10000);
- //
-
-chain->Draw("Cl.fY-TrYpol1.fElements-AliTPCClusterParam::SPosCorrection(0,1,Cl[].fPad,Cl[].fTimeBin,Cl[].fZ,Cl[].fSigmaY2,Cl[].fSigmaZ2,Cl[].fMax):Cl[].fPad-int(Cl[].fPad)",cutA+cutCl+"Cl[].fZ>0&&Cl[].fPad>1","prof",10000);
-
-
-chain->Draw("Cl[].fZ-TrZpol1.fElements-0*AliTPCClusterParam::SPosCorrection(1,1,Cl[].fPad,Cl[].fTimeBin,Cl[].fZ,Cl[].fSigmaY2,Cl[].fSigmaZ2,Cl[].fMax):Cl[].fTimeBin-int(Cl[].fTimeBin)",cutA+cutCl+"Cl[].fZ>0","prof",10000)
-
- */
-
-
-
-
-
- /*
- // check edge effects
- chain->Draw("Cl.fY-TrYpol1.fElements:Cl.fY/Cl.fX",""+cutA+cutCl,"prof",10000)
- //
- chain->Draw("Cl.fY-TrYpol2.fElements:Cl.fPad-int(Cl.fPad)","Cl.fZ>0"+cutA+cutCl+cutE,"prof",100000)
-
- chain->Draw("Cl.fY-TrYpol2.fElements:Cl.fPad-int(Cl.fPad)","Cl.fX>80&&Cl.fZ>0&&Cl.fDetector>35"+cutA+cutCl+cutE,"prof",100000)
-
-
-
- chainFit->Draw("yInOut.fElements[4]:LTr.fP[2]","LTr.fP[1]<0"+cutA,"prof",1000);
- chainFit->Draw("yPol2In.fElements[2]*90*90/4.:LTr.fP[2]","nclO>40&<r.fP[1]<0"+cutA+cutD,"prof")
-
-*/
-
-
-
-
-
- /*
- Edge y effect
-
- dedge = sign(Cl.fY)*(Cl.fX*tan(pi/18)-abs(Cl.fY))
-
-
- chain->Draw("sign(Cl.fY)*(Cl.fY-TrYpol1.fElements):pi/18-abs(Cl.fY/Cl.fX)>>hisYdphi(100,0,0.03)",""+cutA+cutCl,"prof",10000)
-
- chain->Draw("sign(Cl.fY)*(Cl.fY-TrYpol1.fElements):Cl.fX*(pi/18-abs(Cl.fY/Cl.fX))>>hisYdy(100,0,5)",""+cutA+cutCl,"prof",10000)
-
-
-
-
-
- chain->Draw("sign(Cl.fY)*(Cl.fY-TrYpol1.fElements):Cl.fX*(pi/18-abs(Cl.fY/Cl.fX))>>hisYdyIROC(100,0,5)","Cl.fDetector<36"+cutA+cutCl,"prof",100000)
-
- chain->Draw("sign(Cl.fY)*(Cl.fY-TrYpol1.fElements):Cl.fX*(pi/18-abs(Cl.fY/Cl.fX))>>hisYdyOROC(100,0,5)","Cl.fDetector>36"+cutA+cutCl,"prof",100000)
-
-
-
- chain->Draw("Cl.fY-TrYpol1.fElements:sign(Cl.fY)*(Cl.fX*tan(pi/18)-abs(Cl.fY))>>his(100,-5,5)",""+cutA+cutCl,"prof",100000)
-
- chain->Draw("Cl.fY-TrYpol1.fElements:sign(Cl.fY)*(Cl.fX*tan(pi/18)-abs(Cl.fY))>>hisdyInner(100,-5,5)","Cl.fDetector<36"+cutA+cutCl,"prof",100000)
-
-
-
-*/
-
-
-/*
-
-chainFit->Draw("yPol2Out.fElements[2]*90*90/4.:LTr.fP[2]","nclO>40&<r.fP[1]<0"+cutA+cutDY,"prof")
-
-chainFit->Draw("yPol2In.fElements[2]*64*64/4.:LTr.fP[2]","nclI>20&<r.fP[1]<0"+cutA+cutDY,"prof")
-
-
-
-chainFit->Draw("LTr.fId","nclI>10",100000)
-
-chainFit->Draw("yPol2In.fElements[2]:LTr.fId>>his(350,0,350,100,-0.002,0.002)","nclI>20","")
-
-chainFit->Draw("yPol2In.fElements[2]:LTr.fId>>hisPy2In0(350,0,350,100,-0.002,0.002)","nclI>20","");
-
-TH2 * phisPy2In = (TH2*) gROOT->FindObject("hisPy2In0")
-
-*/
-
-
-
-
-
-
-/*
- gSystem->Load("libSTAT.so")
- TStatToolkit toolkit;
- Double_t chi2;
- TVectorD fitParam;
- TMatrixD covMatrix;
- Int_t npoints;
-
- TCut cutA("entries>2&&pphi2<3&&abs(gphiP1-pphiP0)<0.003&&abs(gz1)<6");
-
-
-TString fstring="";
-//
-fstring+="(abs(LTr.fP[1]/250)^3-1)*bz++"; //1
-fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*LTr.fP[2]++"; //2
-fstring+="(abs(LTr.fP[1]/250)^1-1)*bz++"; //3
-fstring+="(abs(LTr.fP[1]/250)-1)*bz*LTr.fP[2]++"; //4
-//
-fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*sin(atan2(lx1,lx0))++" //5
-fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++" //6
-fstring+="(abs(LTr.fP[1]/250)-1)*bz*sin(atan2(lx1,lx0))++" //7
-fstring+="(abs(LTr.fP[1]/250)-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++" //8
-//
-fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*cos(atan2(lx1,lx0))++" //9
-fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++" //10
-fstring+="(abs(LTr.fP[1]/250)-1)*bz*cos(atan2(lx1,lx0))++" //11
-fstring+="(abs(LTr.fP[1]/250)-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++" //12
-
-
-
-
- TString *strq0 = toolkit.FitPlane(treeT,"gphi1-pphi0",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
-
- treeT->SetAlias("fit",strq0->Data());
-
-
- TString *strqP = toolkit.FitPlane(treeT,"1000*(gphiP1-pphiP0)",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
-
- treeT->SetAlias("fitP",strqP->Data());
-
-
- TString *strqDrift = toolkit.FitPlane(treeT,"gz1","LTr.fP[1]++(1-2*(fSide==1))++lx1", cutA, chi2,npoints,fitParam,covMatrix);
- treeT->SetAlias("fitD",strqDrift->Data());
-
-
-treeT->Draw("fit:LTr.fP[1]","abs(bz+0.4)<0.05"+cutA,"");
-{
-for (Int_t i=0; i<6;i++){
-treeT->SetLineColor(i+2);
-treeT->SetMarkerSize(1);
-treeT->SetMarkerStyle(22+i);
-treeT->SetMarkerColor(i+2);
-
-treeT->Draw("fit:LTr.fP[1]",Form("abs(bz+0.4)<0.05&fRod==%d",i)+cutA,"same");
-}
-}
- */
-
-
-
-/*
- TTree * tree = (TTree*)f.Get("FitModels");
-
- TEventList listLFit0("listLFit0","listLFit0");
- TEventList listLFit1("listLFit1","listLFit1");
- tree->Draw(">>listLFit0","seed.fdEdx<200&&seed.fdEdx>40");
- tree->SetEventList(&listLFit0);
-
-
-
-
- gSystem->Load("libSTAT.so")
- TStatToolkit toolkit;
- Double_t chi2;
- TVectorD fitParam;
- TMatrixD covMatrix;
- Int_t npoints;
-
- chain->SetAlias("dp","((Cl.fPad-int(Cl.fPad))*pi)");
- chain->SetAlias("dt","((Cl.fTimeBin-int(Cl.fTimeBin))*pi)");
-
-
- TString fstring="";
- fstring+="cos(dp)++";
- fstring+="sin(dp)++";
- fstring+="cos(dt)++";
- fstring+="sin(dt)++";
-
- TString *str = toolkit.FitPlane(chain,"Cl.fZ-TrZInOut.fElements",fstring->Data(), "Cl.fDetector>35", chi2,npoints,fitParam,covMatrix,-1,0,200);
-
-
-
-*/
-
-
-/*
- Edge effects
+void AliTPCcalibLaser::DumpLaser(const char *finput, Int_t run){
//
//
+ //input="TPCLaserObjects.root"
//
- gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
- gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
- AliXRDPROOFtoolkit tool;
- TChain * chainTrack = tool.MakeChain("laser.txt","Track",0,10200);
- chainTrack->Lookup();
- chainTrack->SetProof(kTRUE);
+ // 0. OBJ: TAxis Delta
+ // 1. OBJ: TAxis bin
+ // 2. OBJ: TAxis rms shape
+ // 3. OBJ: TAxis sqrt(Q)
+ // 4. OBJ: TAxis row
+ // 5. OBJ: TAxis trackID
- TChain * chain = tool.MakeChain("laser.txt","Residuals",0,10200);
- chain->Lookup();
- TChain * chainFit = tool.MakeChain("laser.txt","FitModels",0,10200);
- chainFit->Lookup();
- chainFit->SetProof(kTRUE);
- chain->SetProof(kTRUE);
- //
- // Fit cuts
- //
- TCut cutChi2YOut("sqrt(chi2y2Out*dEdx)<10");
- TCut cutChi2ZOut("sqrt(chi2z2Out*dEdx)<10");
- TCut cutChi2YIn("sqrt(chi2y2In*dEdx)<10");
- TCut cutChi2ZIn("sqrt(chi2z2In*dEdx)<10");
- //
- TCut cutdEdx("sqrt(dEdx)<30&&sqrt(dEdx)>3");
- TCut cutDY("abs(yPol2In.fElements[2]*nclO*nclO/4.)<3")
- TCut cutN("nclO>20&&nclI>20");
- TCut cutA = cutChi2YOut+cutChi2ZOut+cutChi2YIn+cutChi2ZIn+cutN+cutdEdx;
+ const Double_t kSigma=4.;
+ TFile f(finput);
+ AliTPCcalibLaser *laserTPC = (AliTPCcalibLaser*) f.Get("laserTPC");
+ THnSparse * hisPadInput = laserTPC->fHisLaserPad;
+ THnSparse * hisTimeInput = laserTPC->fHisLaserTime;
+ TTreeSRedirector *pcstream= new TTreeSRedirector("hisLasers.root");
+ TVectorD meanY(159), sigmaY(159);
+ TVectorD meanZ(159), sigmaZ(159);
+ TVectorD meanPad(159), sigmaPad(159);
+ TVectorD meanTime(159), sigmaTime(159);
+ TVectorD meanDPad(159), sigmaDPad(159);
+ TVectorD meanDTime(159), sigmaDTime(159);
+ TVectorD meandEdx(159), sigmadEdx(159);
+ TVectorD meanSTime(159), sigmaSTime(159);
+ TVectorD meanSPad(159), sigmaSPad(159);
+ TVectorD entries(159);
//
- // Cluster cuts
+ Int_t indexes[10]={0,1,2,3,4,5,6};
+ TH1 *his=0;
+ AliTPCLaserTrack::LoadTracks();
//
- TCut cutClY("abs(Cl.fY-TrYpol2.fElements)<0.2");
- TCut cutClZ("abs(Cl.fZ-TrZpol2.fElements)<0.4");
- TCut cutClX("abs(Cl.fX)>10");
- TCut cutE("abs(Cl.fY/Cl.fX)<0.14");
- TCut cutCl=cutClY+cutClZ+cutClX;
+ for (Int_t id=0; id<336; id++){ // llop over laser beams
+ printf("id=\t%d\n",id);
+ //
+ AliTPCLaserTrack *ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
+ //
+ hisPadInput->GetAxis(5)->SetRange(id+1,id+1);
+ hisTimeInput->GetAxis(5)->SetRange(id+1,id+1);
+ //
+ his=hisTimeInput->Projection(3);
+ Int_t firstBindEdx=his->FindFirstBinAbove(0);
+ Int_t lastBindEdx=his->FindLastBinAbove(0);
+ hisPadInput->GetAxis(3)->SetRange(firstBindEdx, lastBindEdx);
+ hisTimeInput->GetAxis(3)->SetRange(firstBindEdx, lastBindEdx);
+ delete his;
+ //
+ his=hisTimeInput->Projection(1);
+ Int_t firstBinTime=his->FindFirstBinAbove(0);
+ Int_t lastBinTime=his->FindLastBinAbove(0);
+ //hisTimeInput->GetAxis(1)->SetRange(firstBinTime, lastBinTime);
+ delete his;
+ //
+ //
+ his=hisTimeInput->Projection(2);
+ Int_t firstBinZ=his->FindFirstBinAbove(0);
+ Int_t lastBinZ=his->FindLastBinAbove(0);
+ //hisTimeInput->GetAxis(2)->SetRange(firstBinZ, lastBinZ);
+ delete his;
+ //
+ his=hisPadInput->Projection(2);
+ Int_t firstBinY=his->FindFirstBinAbove(0);
+ Int_t lastBinY=his->FindLastBinAbove(0);
+ //hisPadInput->GetAxis(2)->SetRange(firstBinY, lastBinY);
+ delete his;
+ //
+ //
+ //
+ THnSparse *hisPad0 = hisPadInput->Projection(5,indexes);
+ THnSparse *hisTime0 = hisTimeInput->Projection(5,indexes);
+ //
+ //
+ for (Int_t irow=0; irow<159; irow++){
+ entries[irow]=0;
+ if ((*(ltrp->GetVecSec()))[irow] <0) continue;
+ if ((*(ltrp->GetVecLX()))[irow] <80) continue;
+
+ hisPad0->GetAxis(4)->SetRange(irow+1,irow+1);
+ hisTime0->GetAxis(4)->SetRange(irow+1,irow+1);
+ //THnSparse *hisPad = hisPad0->Projection(4,indexes);
+ //THnSparse *hisTime = hisTime0->Projection(4,indexes);
+ THnSparse *hisPad = hisPad0;
+ THnSparse *hisTime = hisTime0;
+ //
+ // Get mean value of QA variables
+ //
+ // dEdx
+ his=hisTime->Projection(3);
+ his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()+kSigma*his->GetRMS());
+ meandEdx[irow] =his->GetMean();
+ sigmadEdx[irow]=his->GetRMS();
+ Int_t bindedx0= his->FindBin(meandEdx[irow]-kSigma*sigmadEdx[irow]);
+ Int_t bindedx1= his->FindBin(meandEdx[irow]+kSigma*sigmadEdx[irow]);
+ // hisPad->GetAxis(3)->SetRange(bindedx0,bindedx1);
+ //hisTime->GetAxis(3)->SetRange(bindedx0,bindedx1 );
+ delete his;
+ //
+ // sigma Time
+ //
+ his=hisTime->Projection(2);
+ his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()-kSigma*his->GetRMS());
+ meanSTime[irow] =his->GetMean();
+ sigmaSTime[irow]=his->GetRMS();
+ Int_t binSTime0= his->FindBin(his->GetMean()-kSigma*his->GetRMS());
+ Int_t binSTime1= his->FindBin(his->GetMean()+kSigma*his->GetRMS());
+ // hisTime->GetAxis(2)->SetRange(binSTime0, binSTime1);
+ delete his;
+ //
+ // sigma Pad
+ his=hisPad->Projection(2);
+ his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()+kSigma*his->GetRMS());
+ meanSPad[irow] =his->GetMean();
+ sigmaSPad[irow]=his->GetRMS();
+ Int_t binSPad0= his->FindBin(his->GetMean()-kSigma*his->GetRMS());
+ Int_t binSPad1= his->FindBin(his->GetMean()+kSigma*his->GetRMS());
+ // hisPad->GetAxis(2)->SetRange(binSPad0, binSPad1);
+ delete his;
+ //
+ // apply selection on QA variables
+ //
+ //
+ //
+ // Y
+ his=hisPad->Projection(0);
+ entries[irow]=his->GetEntries();
+ his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()+kSigma*his->GetRMS());
+ meanY[irow] =his->GetMean();
+ sigmaY[irow]=his->GetRMS();
+ delete his;
+ // Z
+ his=hisTime->Projection(0);
+ his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()+kSigma*his->GetRMS());
+ meanZ[irow] =his->GetMean();
+ sigmaZ[irow]=his->GetRMS();
+ delete his;
+ // Pad
+ his=hisPad->Projection(1);
+ his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()+kSigma*his->GetRMS());
+ meanPad[irow] =his->GetMean();
+ meanDPad[irow] =his->GetMean()-Int_t(his->GetMean());
+ sigmaPad[irow]=his->GetRMS();
+ delete his;
+ // Time
+ his=hisTime->Projection(1);
+ his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()+kSigma*his->GetRMS());
+ meanTime[irow] = his->GetMean();
+ meanDTime[irow] = his->GetMean()-Int_t(his->GetMean());
+ sigmaTime[irow]=his->GetRMS();
+ delete his;
+ //
+ //delete hisTime;
+ //delete hisPad;
+ }
+ //
+ //
+ //
+ (*pcstream)<<"laserClusters"<<
+ "id="<<id<<
+ "run="<<run<<
+ "LTr.="<<ltrp<<
+ //
+ "entries.="<<&entries<<
+ "my.="<<&meanY<< //mean delta y
+ "rmsy.="<<&sigmaY<< //rms deltay
+ "mz.="<<&meanZ<< //mean deltaz
+ "rmsz.="<<&sigmaZ<< //rms z
+ //
+ "mPad.="<<&meanPad<< // mean pad
+ "mDPad.="<<&meanDPad<< // mead dpad
+ "rmsPad.="<<&sigmaPad<< // rms pad
+ "mTime.="<<&meanTime<<
+ "mDTime.="<<&meanTime<<
+ "rmsTime.="<<&sigmaTime<<
+ //
+ "mdEdx.="<<&meandEdx<< //mean dedx
+ "rmsdEdx.="<<&sigmadEdx<< //rms dedx
+ "mSPad.="<<&meanSPad<< //mean sigma pad
+ "rmsSPad.="<<&sigmaSPad<< //rms sigma pad
+ "mSTime.="<<&meanSTime<< //mean sigma time
+ "rmsSTime.="<<&sigmaSTime<<
+ "\n";
+ //
+ delete hisPad0;
+ delete hisTime0;
+ }
+ delete pcstream;
+ /*
+
+ */
+}
- // check edge effects
- chain->Draw("Cl.fY-TrYpol1.fElements:Cl.fY/Cl.fX",""+cutA+cutCl,"prof",10000)
- //
- chain->Draw("Cl.fY-TrYpol2.fElements:Cl.fPad-int(Cl.fPad)","Cl.fZ>0"+cutA+cutCl+cutE,"prof",100000)
+void AliTPCcalibLaser::FitLaserClusters(Int_t run){
+ //
+ //
+ //input="TPCLaserObjects.root"
+ //Algorithm:
+ // 1. Select cluster candidates, remove outlyers
+ // edge clusters
+ // clusters with atypical spread (e.g due track overlaps)
+ // small amount of entries clusters (absolute minimal cut + raltive -to mean cut)
+ // 2. Fit the tracklets -per sector - in pad and time coordinate frame
+ // Remove outlyers
+ // Store info distance of track to pad, time center
+ // Fit the correction for distance to the center (sin,cos)
+ // 3. Do local fit
+ const Double_t kEpsilon=0.000001;
+ const Int_t kMinClusters=20;
+ const Double_t kEdgeCut=3;
+ const Double_t kDistCut=1.5; // cut distance to the ideal track
+ const Double_t kDistCutFit=0.5;
+ const Double_t kDistCutFitPad=0.25;
+ const Double_t kDistCutFitTime=0.25;
+ const Int_t kSmoothRow=5.;
+ TFile f("hisLasers.root"); // Input file
+ TTree * treeInput=(TTree*)f.Get("laserClusters");
+ TTreeSRedirector *pcstream=new TTreeSRedirector("fitLasers.root");
+ TVectorD *vecN=0;
+ TVectorD *vecMY=0;
+ TVectorD *vecMZ=0;
+ TVectorD *vecPad=0;
+ TVectorD *vecTime=0;
+ TVectorD *vecSY=0;
+ TVectorD *vecSZ=0;
+ TVectorD *meandEdx=0;
+ TVectorD isOK(159);
+ TVectorD fitPad(159);
+ TVectorD fitTime(159);
+ TVectorD fitPadLocal(159);
+ TVectorD fitTimeLocal(159);
+ TVectorD fitDPad(159);
+ TVectorD fitDTime(159);
+ TVectorD fitIPad(159);
+ TVectorD fitITime(159);
+ Double_t chi2PadIROC=0;
+ Double_t chi2PadOROC=0;
+ //
+ treeInput->SetBranchAddress("my.",&vecMY);
+ treeInput->SetBranchAddress("mz.",&vecMZ);
+ treeInput->SetBranchAddress("mPad.",&vecPad);
+ treeInput->SetBranchAddress("mTime.",&vecTime);
+ treeInput->SetBranchAddress("rmsy.",&vecSY);
+ treeInput->SetBranchAddress("rmsz.",&vecSZ);
+ treeInput->SetBranchAddress("entries.",&vecN);
+ treeInput->SetBranchAddress("mdEdx.",&meandEdx);
+
+ AliTPCLaserTrack::LoadTracks();
+ //
+ //
+ TVectorD fitPadIROC(3), fitPadOROC(3);
+ TVectorD fitPadIROCSin(3), fitPadOROCSin(3);
+ TVectorD fitTimeIROC(3), fitTimeOROC(3);
+ TVectorD fitTimeIROCSin(3), fitTimeOROCSin(3);
+ //
+ AliTPCROC * roc = AliTPCROC::Instance();
+ Double_t refX=roc->GetPadRowRadii(0,roc->GetNRows(0)-1);
- chain->Draw("Cl.fY-TrYpol2.fElements:Cl.fPad-int(Cl.fPad)","Cl.fX>80&&Cl.fZ>0&&Cl.fDetector>35"+cutA+cutCl+cutE,"prof",100000)
-
-
-
- chainFit->Draw("yInOut.fElements[4]:LTr.fP[2]","LTr.fP[1]<0"+cutA,"prof",1000);
- chainFit->Draw("yPol2In.fElements[2]*90*90/4.:LTr.fP[2]","nclO>40&<r.fP[1]<0"+cutA+cutD,"prof")
-
-*/
-
+ //
+ for (Int_t id=0; id<336; id++){
+ //
+ treeInput->GetEntry(id);
+ AliTPCLaserTrack *ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
+ Int_t medianEntries = TMath::Nint(TMath::Median(159,vecN->GetMatrixArray()));
+ Double_t medianRMSY = TMath::Median(159,vecSY->GetMatrixArray());
+ Double_t rmsRMSY = TMath::RMS(159,vecSY->GetMatrixArray());
+ Double_t medianRMSZ = TMath::Median(159,vecSZ->GetMatrixArray());
+ Double_t rmsRMSZ = TMath::RMS(159,vecSZ->GetMatrixArray());
+ Double_t mdEdx = TMath::Median(159,meandEdx->GetMatrixArray());
+ Int_t sectorInner= TMath::Nint(ltrp->GetVecSec()->GetMatrixArray()[63/2]);
+ Int_t sectorOuter= TMath::Nint(ltrp->GetVecSec()->GetMatrixArray()[64+96/2]);
+ TLinearFitter fitterY(2,"pol1");
+ TLinearFitter fitterZ(2,"pol1");
+ TLinearFitter fitterPad(2,"pol1");
+ TLinearFitter fitterTime(2,"pol1");
+ TLinearFitter fitterPadSin(2,"hyp1");
+ TLinearFitter fitterTimeSin(3,"hyp2");
+ //
+ //
+ for (UInt_t irow=0; irow<159; irow++){
+ fitPad[irow]=0; fitIPad[irow]=0; fitDPad[irow]=0;
+ fitTime[irow]=0; fitITime[irow]=0; fitDTime[irow]=0;
+ Double_t sign=(ltrp->GetZ()>0) ? 1.:-1.;
+ isOK[irow]=kFALSE;
+ fitPad[irow]=0;
+ fitTime[irow]=0;
+ Int_t sector=(irow<roc->GetNRows(0))? sectorInner:sectorOuter;
+ Int_t npads=(irow<roc->GetNRows(0))? roc->GetNPads(sector,irow):roc->GetNPads(sector,irow-roc->GetNRows(0));
+ (*vecPad)[irow]-=npads*0.5;
+ //
+ if ((irow<roc->GetNRows(0)) &&TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sectorInner)>0.1) continue;
+ if ((irow>=roc->GetNRows(0)) &&TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sectorOuter)>0.1) continue;
+ //
+ if (TMath::Abs((*vecMY)[irow])<kEpsilon) continue; //not determined position
+ if (TMath::Abs((*vecMZ)[irow])<kEpsilon) continue; //not determined position
+ if (TMath::Abs((*vecPad)[irow])<kEpsilon) continue; //not determined position
+ if (TMath::Abs((*vecTime)[irow])<kEpsilon) continue; //not determined position
+ if ((*vecN)[irow]<0.5*medianEntries) continue; //small amount of clusters
+ if ((*vecSY)[irow]>medianRMSY+3*rmsRMSY) continue; //big sigma
+ if ((*vecSZ)[irow]>medianRMSZ+3*rmsRMSZ) continue; //big sigma
+ Double_t dEdge= TMath::Abs((*(ltrp->GetVecLY()))[irow])-(*(ltrp->GetVecLX()))[irow]*TMath::Tan(TMath::Pi()/18.); //edge cut
+ if (TMath::Abs(dEdge)<kEdgeCut) continue;
+ if (irow<roc->GetNRows(0)){
+ if (TMath::Abs(((*ltrp->GetVecLY())[irow])-sign*(*vecPad)[irow]*0.4)>kDistCut) continue;
+ }
+ if (irow>roc->GetNRows(0)){
+ if (TMath::Abs(((*ltrp->GetVecLY())[irow])-sign*(*vecPad)[irow]*0.6)>kDistCut) continue;
+ }
+
+ isOK[irow]=kTRUE;
+ }
+ //
+ //fit OROC - get delta pad and delta time
+ //
+ fitterPad.ClearPoints();
+ fitterTime.ClearPoints();
+ fitterPadSin.ClearPoints();
+ fitterTimeSin.ClearPoints();
+ {for (Int_t irow=2; irow<157; irow++){
+ if (isOK[irow]<0.5) continue;
+ if (TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sectorOuter)>0.1) continue;
+ if (TMath::Abs(ltrp->GetVecLX()->GetMatrixArray()[irow])<80) continue;
+ Double_t y=(*vecPad)[irow];
+ Double_t z=(*vecTime)[irow];
+ Double_t x=ltrp->GetVecLX()->GetMatrixArray()[irow]-refX;
+ fitterPad.AddPoint(&x,y);
+ fitterTime.AddPoint(&x,z);
+ }}
+ chi2PadOROC=0;
+ if (fitterPad.GetNpoints()>kMinClusters&&fitterTime.GetNpoints()>kMinClusters){
+ fitterPad.Eval();
+ fitterTime.Eval();
+ chi2PadOROC=TMath::Sqrt(fitterPad.GetChisquare()/fitterPad.GetNpoints());
+ for (Int_t irow=2; irow<157; irow++){
+ if (isOK[irow]<0.5) continue;
+ if (TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sectorOuter)>0.1) continue;
+ if (TMath::Abs(ltrp->GetVecLX()->GetMatrixArray()[irow])<80) continue;
+ Double_t y=(*vecPad)[irow];
+ Double_t z=(*vecTime)[irow];
+ Double_t x=ltrp->GetVecLX()->GetMatrixArray()[irow]-refX;
+ Double_t fitP=fitterPad.GetParameter(0)+fitterPad.GetParameter(1)*x;
+ Double_t fitT=fitterTime.GetParameter(0)+fitterTime.GetParameter(1)*x;
+ fitPad[irow]=fitterPad.GetParameter(0)+fitterPad.GetParameter(1)*x;
+ fitTime[irow]=fitterTime.GetParameter(0)+fitterTime.GetParameter(1)*x;
+ fitDPad[irow]=y-(fitterPad.GetParameter(0)+fitterPad.GetParameter(1)*x);
+ fitDTime[irow]=z-(fitterTime.GetParameter(0)+fitterTime.GetParameter(1)*x);
+ fitIPad[irow]=fitP-TMath::Nint(fitP-0.5);
+ fitITime[irow]=fitT-TMath::Nint(fitT-0.5);
+ if (fitDPad[irow]>kDistCutFit) isOK[irow]=kFALSE;
+ if (fitDTime[irow]>kDistCutFit) isOK[irow]=kFALSE;
+ if (isOK[irow]>0){
+ Double_t xxxPad[2]={TMath::Sin(2*TMath::Pi()*fitIPad[irow])};
+ Double_t xxxTime[3]={TMath::Sin(2*TMath::Pi()*fitITime[irow]),
+ TMath::Cos(2*TMath::Pi()*fitITime[irow])};
+ fitterPadSin.AddPoint(xxxPad,fitDPad[irow]);
+ fitterTimeSin.AddPoint(xxxTime,fitDTime[irow]);
+ }
+ }
+ fitterPadSin.Eval();
+ fitterTimeSin.Eval();
+ fitterPadSin.FixParameter(0,0);
+ fitterTimeSin.FixParameter(0,0);
+ fitterPadSin.Eval();
+ fitterTimeSin.Eval();
+ //
+ fitterPad.GetParameters(fitPadOROC);
+ fitterTime.GetParameters(fitTimeOROC);
+ fitterPadSin.GetParameters(fitPadOROCSin);
+ fitterTimeSin.GetParameters(fitTimeOROCSin);
+ }
+ //
+ //
+ //fit IROC
+ //
+ fitterPad.ClearPoints();
+ fitterTime.ClearPoints();
+ fitterPadSin.ClearPoints();
+ fitterTimeSin.ClearPoints();
+ for (Int_t irow=2; irow<157; irow++){
+ if (isOK[irow]<0.5) continue;
+ if (TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sectorInner)>0.1) continue;
+ if (TMath::Abs(ltrp->GetVecLX()->GetMatrixArray()[irow])<80) continue;
+ Double_t y=(*vecPad)[irow];
+ Double_t z=(*vecTime)[irow];
+ Double_t x=ltrp->GetVecLX()->GetMatrixArray()[irow]-refX;
+ fitterPad.AddPoint(&x,y);
+ fitterTime.AddPoint(&x,z);
+ }
+ chi2PadIROC=0;
+ if (fitterPad.GetNpoints()>kMinClusters&&fitterTime.GetNpoints()>kMinClusters){
+ fitterPad.Eval();
+ fitterTime.Eval();
+ chi2PadIROC=TMath::Sqrt(fitterPad.GetChisquare()/fitterPad.GetNpoints());
+ for (Int_t irow=2; irow<157; irow++){
+ if (isOK[irow]<0.5) continue;
+ if (TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sectorInner)>0.1) continue;
+ if (TMath::Abs(ltrp->GetVecLX()->GetMatrixArray()[irow])<80) continue;
+ Double_t y=(*vecPad)[irow];
+ Double_t z=(*vecTime)[irow];
+ Double_t x=ltrp->GetVecLX()->GetMatrixArray()[irow]-refX;
+ Double_t fitP=fitterPad.GetParameter(0)+fitterPad.GetParameter(1)*x;
+ Double_t fitT=fitterTime.GetParameter(0)+fitterTime.GetParameter(1)*x;
+ fitPad[irow]=fitterPad.GetParameter(0)+fitterPad.GetParameter(1)*x;
+ fitTime[irow]=fitterTime.GetParameter(0)+fitterTime.GetParameter(1)*x;
+ fitDPad[irow]=y-(fitterPad.GetParameter(0)+fitterPad.GetParameter(1)*x);
+ fitDTime[irow]=z-(fitterTime.GetParameter(0)+fitterTime.GetParameter(1)*x);
+ fitIPad[irow]=fitP-TMath::Nint(fitP-0.5);
+ fitITime[irow]=fitT-TMath::Nint(fitT-0.5);
+ if (fitDPad[irow]>kDistCutFit) isOK[irow]=kFALSE;
+ if (fitDTime[irow]>kDistCutFit) isOK[irow]=kFALSE;
+ if (isOK[irow]>0.5){
+ Double_t xxxPad[3]={TMath::Sin(2*TMath::Pi()*fitIPad[irow]),
+ TMath::Cos(2*TMath::Pi()*fitIPad[irow])};
+ Double_t xxxTime[3]={TMath::Sin(2*TMath::Pi()*fitITime[irow]),
+ TMath::Cos(2*TMath::Pi()*fitITime[irow])};
+ fitterPadSin.AddPoint(xxxPad,fitDPad[irow]);
+ fitterTimeSin.AddPoint(xxxTime,fitDTime[irow]);
+ }
+ }
+ fitterPadSin.Eval();
+ fitterTimeSin.Eval();
+ fitterPadSin.FixParameter(0,0);
+ fitterTimeSin.FixParameter(0,0);
+ fitterPadSin.Eval();
+ fitterTimeSin.Eval();
+ fitterPad.GetParameters(fitPadIROC);
+ fitterTime.GetParameters(fitTimeIROC);
+ fitterPadSin.GetParameters(fitPadIROCSin);
+ fitterTimeSin.GetParameters(fitTimeIROCSin);
+ }
+ for (Int_t irow=0; irow<159; irow++){
+ if (TMath::Abs(fitDPad[irow])<kEpsilon) isOK[irow]=kFALSE;
+ if (TMath::Abs(fitDTime[irow])<kEpsilon) isOK[irow]=kFALSE;
+ if (TMath::Abs(fitDPad[irow])>kDistCutFitPad) isOK[irow]=kFALSE;
+ if (TMath::Abs(fitDTime[irow])>kDistCutFitTime) isOK[irow]=kFALSE;
+ }
+ for (Int_t irow=kSmoothRow/2; irow<159-kSmoothRow/2; irow++){
+ fitPadLocal[irow]=0;
+ fitTimeLocal[irow]=0;
+ if (isOK[irow]<0.5) continue;
+ Int_t sector=(irow<roc->GetNRows(0))? sectorInner:sectorOuter;
+ if (TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sector)>0.1) continue;
+ //
+ TLinearFitter fitterPadLocal(2,"pol1");
+ TLinearFitter fitterTimeLocal(2,"pol1");
+ Double_t xref=ltrp->GetVecLX()->GetMatrixArray()[irow];
+ for (Int_t delta=-kSmoothRow; delta<=kSmoothRow; delta++){
+ Int_t jrow=irow+delta;
+ if (jrow<0) jrow=0;
+ if (jrow>159) jrow=159;
+ if (isOK[jrow]<0.5) continue;
+ if (TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[jrow]-sector)>0.1) continue;
+ Double_t y=(*vecPad)[jrow];
+ Double_t z=(*vecTime)[jrow];
+ Double_t x=ltrp->GetVecLX()->GetMatrixArray()[jrow]-xref;
+ fitterPadLocal.AddPoint(&x,y);
+ fitterTimeLocal.AddPoint(&x,z);
+ }
+ if (fitterPadLocal.GetNpoints()<kSmoothRow) continue;
+ fitterPadLocal.Eval();
+ fitterTimeLocal.Eval();
+ fitPadLocal[irow]=fitterPadLocal.GetParameter(0);
+ fitTimeLocal[irow]=fitterTimeLocal.GetParameter(0);
+ }
+ //
+ //
+ (*pcstream)<<"fit"<<
+ "run="<<run<<
+ "id="<<id<<
+ "chi2PadIROC="<<chi2PadIROC<<
+ "chi2PadOROC="<<chi2PadOROC<<
+ "mdEdx="<<mdEdx<<
+ "LTr.="<<ltrp<<
+ "isOK.="<<&isOK<<
+ // mean measured-ideal values
+ "mY.="<<vecMY<<
+ "mZ.="<<vecMZ<<
+ // local coordinate fit
+ "mPad.="<<vecPad<<
+ "mTime.="<<vecTime<<
+ "fitPad.="<<&fitPad<<
+ "fitTime.="<<&fitTime<<
+ "fitPadLocal.="<<&fitPadLocal<<
+ "fitTimeLocal.="<<&fitTimeLocal<<
+ "fitDPad.="<<&fitDPad<<
+ "fitDTime.="<<&fitDTime<<
+ "fitIPad.="<<&fitIPad<<
+ "fitITime.="<<&fitITime<<
+ //
+ "fitPadIROC.="<<&fitPadIROC<< // pad fit linear IROC
+ "fitPadIROCSin.="<<&fitPadIROCSin<< // pad fit linear+ pad correction
+ "fitPadOROC.="<<&fitPadOROC<<
+ "fitPadOROCSin.="<<&fitPadOROCSin<<
+ //
+ "fitTimeIROC.="<<&fitTimeIROC<<
+ "fitTimeIROCSin.="<<&fitTimeIROCSin<<
+ "fitTimeOROC.="<<&fitTimeOROC<<
+ "fitTimeOROCSin.="<<&fitTimeOROCSin<<
+ "\n";
+ }
+ delete pcstream;
+}
virtual void Process(AliESDEvent *event);
Int_t GetNtracks(){return fNtracks;}
virtual void Analyze();
+ static void DumpLaser(const char *finput, Int_t run);
+ static void FitLaserClusters(Int_t run);
virtual Long64_t Merge(TCollection *li);
virtual void DumpMeanInfo(Int_t run=-1);
static void DumpScanInfo(TTree * tree, const char * cutUser="entries>300&&(gz2<0.15&&gphi2<0.1&&gp42<0.02&&abs(gp41)<0.03)");
// Refit residuals histogram
//
THnSparseS *fHisLaser; // N dim histogram of laser
+ THnSparseS *fHisLaserPad; // N dim histogram of laser
+ THnSparseS *fHisLaserTime; // N dim histogram of laser
+ //
TH2F *fHisNclIn; //->Number of clusters inner
TH2F *fHisNclOut; //->Number of clusters outer
TH2F *fHisNclIO; //->Number of cluster inner outer
Double_t fFixedFitCside1; // Fixed drift v constant 1 - C side
//
private:
- ClassDef(AliTPCcalibLaser,5)
+ ClassDef(AliTPCcalibLaser,6)
};
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff