From 76c58ee2097ab24b0f7d908aeea238d4f369c09c Mon Sep 17 00:00:00 2001 From: marian Date: Tue, 15 Feb 2011 14:28:15 +0000 Subject: [PATCH] AliTPCcalibAlign.cxx - changes in filling lookup table 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 --- TPC/AliTPCCalibGlobalMisalignment.cxx | 29 +- TPC/AliTPCCalibGlobalMisalignment.h | 5 +- TPC/AliTPCcalibAlign.cxx | 211 ++++--- TPC/AliTPCcalibCalib.cxx | 14 + TPC/AliTPCcalibCosmic.cxx | 713 ++++++++++++++++++++-- TPC/AliTPCcalibCosmic.h | 16 +- TPC/AliTPCcalibLaser.cxx | 844 ++++++++++++++++++-------- TPC/AliTPCcalibLaser.h | 7 +- 8 files changed, 1469 insertions(+), 370 deletions(-) diff --git a/TPC/AliTPCCalibGlobalMisalignment.cxx b/TPC/AliTPCCalibGlobalMisalignment.cxx index ce5be288085..d7744a5b4d6 100644 --- a/TPC/AliTPCCalibGlobalMisalignment.cxx +++ b/TPC/AliTPCCalibGlobalMisalignment.cxx @@ -53,6 +53,7 @@ AliTPCCalibGlobalMisalignment::AliTPCCalibGlobalMisalignment() 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 { // @@ -71,6 +72,7 @@ AliTPCCalibGlobalMisalignment::~AliTPCCalibGlobalMisalignment() { 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 } @@ -100,6 +102,16 @@ void AliTPCCalibGlobalMisalignment::SetAlignGlobal(const TGeoMatrix * matrixGlob 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 @@ -137,7 +149,8 @@ void AliTPCCalibGlobalMisalignment::GetCorrection(const Float_t x[],const Short_ // 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]); @@ -162,8 +175,8 @@ void AliTPCCalibGlobalMisalignment::GetCorrection(const Float_t x[],const Short_ 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]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 @@ -214,6 +227,16 @@ void AliTPCCalibGlobalMisalignment::GetCorrection(const Float_t x[],const Short_ 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 diff --git a/TPC/AliTPCCalibGlobalMisalignment.h b/TPC/AliTPCCalibGlobalMisalignment.h index 1dec213f733..4e849f8d1e8 100644 --- a/TPC/AliTPCCalibGlobalMisalignment.h +++ b/TPC/AliTPCCalibGlobalMisalignment.h @@ -53,8 +53,10 @@ public: // 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(); @@ -87,11 +89,12 @@ private: // 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 diff --git a/TPC/AliTPCcalibAlign.cxx b/TPC/AliTPCcalibAlign.cxx index ef84bc9b14c..a0bddbde76c 100644 --- a/TPC/AliTPCcalibAlign.cxx +++ b/TPC/AliTPCcalibAlign.cxx @@ -958,6 +958,8 @@ void AliTPCcalibAlign::ProcessTracklets(const AliExternalTrackParam &tp1, 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 } } @@ -1671,7 +1673,11 @@ void AliTPCcalibAlign::MakeResidualHistos(){ // 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; @@ -1712,7 +1718,8 @@ void AliTPCcalibAlign::MakeResidualHistosTracklet(){ // 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; @@ -1734,21 +1741,27 @@ void AliTPCcalibAlign::MakeResidualHistosTracklet(){ // 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()); } @@ -1796,6 +1809,7 @@ void AliTPCcalibAlign::FillHisto(AliExternalTrackParam *tp1, // // Fill residual histograms // Track2-Track1 + if (s2Fill(x); x[0]=p2.GetZ()-p1.GetZ(); @@ -2217,15 +2235,43 @@ void AliTPCcalibAlign::Add(AliTPCcalibAlign * align){ 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; + } } } @@ -2638,9 +2684,8 @@ void AliTPCcalibAlign::UpdateClusterDeltaField(const AliTPCseed * seed){ // // 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 @@ -2648,6 +2693,10 @@ void AliTPCcalibAlign::UpdateClusterDeltaField(const AliTPCseed * seed){ 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; @@ -2658,12 +2707,31 @@ void AliTPCcalibAlign::UpdateClusterDeltaField(const AliTPCseed * seed){ if (seed->GetNumberOfClusters()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; @@ -2677,7 +2745,6 @@ void AliTPCcalibAlign::UpdateClusterDeltaField(const AliTPCseed * seed){ ncl++; } if (nclGetDetector()%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)GetX()-xIFC)GetX()-xOFC)GetX()-fXIO)GetY()-trackOut.GetY())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); @@ -2740,48 +2800,63 @@ void AliTPCcalibAlign::UpdateClusterDeltaField(const AliTPCseed * seed){ 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)GetX()-xIFC)GetX()-xOFC)GetX()-fXIO)GetY()-trackIn.GetY())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); } - } @@ -2792,7 +2867,7 @@ void AliTPCcalibAlign::UpdateAlignSector(const AliTPCseed * track,Int_t isec){ // 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; // @@ -2842,12 +2917,12 @@ void AliTPCcalibAlign::UpdateAlignSector(const AliTPCseed * track,Int_t isec){ } 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()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; @@ -314,6 +321,8 @@ Bool_t AliTPCcalibCalib::RefitTrack(AliESDtrack * track, AliTPCseed *seed, Floa 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); @@ -353,6 +362,8 @@ Bool_t AliTPCcalibCalib::RefitTrack(AliESDtrack * track, AliTPCseed *seed, Floa 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; @@ -391,6 +402,9 @@ Bool_t AliTPCcalibCalib::RefitTrack(AliESDtrack * track, AliTPCseed *seed, Floa 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); diff --git a/TPC/AliTPCcalibCosmic.cxx b/TPC/AliTPCcalibCosmic.cxx index aea04cfe3d0..373cc739222 100644 --- a/TPC/AliTPCcalibCosmic.cxx +++ b/TPC/AliTPCcalibCosmic.cxx @@ -1,3 +1,5 @@ + + /************************************************************************** * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * * @@ -75,15 +77,20 @@ #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) @@ -99,7 +106,8 @@ AliTPCcalibCosmic::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++){ @@ -125,7 +133,8 @@ AliTPCcalibCosmic::AliTPCcalibCosmic(const Text_t *name, const Text_t *title) 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); @@ -180,10 +189,11 @@ void AliTPCcalibCosmic::Init(){ // 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"; @@ -212,50 +222,54 @@ void AliTPCcalibCosmic::Init(){ 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++){ @@ -265,16 +279,16 @@ void AliTPCcalibCosmic::Init(){ 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;ivar2GetAxis(ivar2)->SetName(axisName[ivar2].Data()); fHistoDelta[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data()); fHistoPull[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data()); @@ -314,6 +328,13 @@ void AliTPCcalibCosmic::Add(const AliTPCcalibCosmic* cosmic){ fHistodEdxTot[ivar]->Add(cosmic->fHistodEdxTot[ivar]); } } + if (cosmic->fCosmicTree){ + if (!fCosmicTree) { + fCosmicTree = new TTree("pairs","pairs"); + fCosmicTree->SetDirectory(0); + } + AliTPCcalibCosmic::AddTree(fCosmicTree,cosmic->fCosmicTree); + } } @@ -333,24 +354,34 @@ void AliTPCcalibCosmic::Process(AliESDEvent *event) { 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 @@ -360,9 +391,8 @@ void AliTPCcalibCosmic::FillHistoPerformance(const AliExternalTrackParam *par0, 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); @@ -380,6 +410,7 @@ void AliTPCcalibCosmic::FillHistoPerformance(const AliExternalTrackParam *par0, x[6] = param0Combined->GetSigned1Pt(); x[7] = param0Combined->Pt(); x[8] = alpha; + x[9] = cross; // deltas Double_t delta[6]; Double_t sigma[6]; @@ -406,8 +437,8 @@ void AliTPCcalibCosmic::FillHistoPerformance(const AliExternalTrackParam *par0, // 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]; @@ -651,13 +682,17 @@ void AliTPCcalibCosmic::FindPairs(AliESDEvent *event) { 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(); @@ -673,7 +708,10 @@ void AliTPCcalibCosmic::FindPairs(AliESDEvent *event) { // // // - 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" << @@ -1027,4 +1065,597 @@ void AliTPCcalibCosmic::UpdateTrack(AliExternalTrackParam &track1, const AliExte +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(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;itrack0GetTrack(itrack0); + if (!track0) continue; + if (!track0->IsOn(AliESDtrack::kTPCrefit)) continue; + + if (TMath::Abs(AliTracker::GetBz())>1&&track0->Pt()GetTPCncls()GetY())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])GetInnerParam(); + for (Int_t itrack1=itrack0+1;itrack1GetTrack(itrack1); + if (!track1) continue; + if (!track1->IsOn(AliESDtrack::kTPCrefit)) continue; + if (track1->GetKinkIndex(0)>0) continue; + if (TMath::Abs(AliTracker::GetBz())>1&&track1->Pt()GetTPCncls()1&&TMath::Max(track1->Pt(), track0->Pt())GetY())GetImpactParametersTPC(dcaTPC,covTPC); + // if (TMath::Abs(dcaTPC[0])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(calibObject))) break; + } + for (Int_t l=0;(calibObject=friendTrack1->GetCalibObject(l));++l) { + if ((seed1=dynamic_cast(calibObject))) break; + } + // + if (pcstream){ + (*pcstream)<<"pairs"<< + "run="<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; iGetEntry(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; iGetEntry(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(calibObject))) break; + } + for (Int_t l=0;(calibObject=ftrack1->GetCalibObject(l));++l) { + if ((seed1=dynamic_cast(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()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()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 ((nclA0nclC0) && (nclA1nclC1)) cosmicType=3; // CA side + //if ((nclA0>nclC0) && (nclA1nclC0) && (nclA11)&&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=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="<=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="<tmin => deltaT=Tcmax-tmax + // 1.b) rmin>Rcmin && tmintmax => deltaT=Tcmax-tmin + // // case the z matching gives proper time + // 1.c) rmaxTcmin && tmax 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)0 + if (rmax0tmin0) vectorDT[0]=kMaxTime-tmax0; // disapear at CE + if (rmax1tmin1) vectorDT[1]=kMaxTime-tmax1; // disapear at CE + // min cross - deltaT<0 - OK they are correlated + if (rmax0kMinTime && tmax0kMinTime && tmax1 kMinRCut+kMaxDCut && tmin0 tmax0) vectorDT[4]=kMaxTime-tmin0; + if (rmin1> kMinRCut+kMaxDCut && 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; +} diff --git a/TPC/AliTPCcalibCosmic.h b/TPC/AliTPCcalibCosmic.h index ffee9e5e13e..9f00377d0eb 100644 --- a/TPC/AliTPCcalibCosmic.h +++ b/TPC/AliTPCcalibCosmic.h @@ -26,6 +26,8 @@ public: // 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); @@ -34,10 +36,10 @@ public: 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;}; @@ -53,7 +55,8 @@ public: 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 @@ -62,7 +65,7 @@ public: 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); @@ -85,10 +88,11 @@ private: 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 diff --git a/TPC/AliTPCcalibLaser.cxx b/TPC/AliTPCcalibLaser.cxx index e44f72c4f6d..83294ad63eb 100644 --- a/TPC/AliTPCcalibLaser.cxx +++ b/TPC/AliTPCcalibLaser.cxx @@ -128,6 +128,7 @@ #include "AliDCSSensorArray.h" #include "AliDCSSensor.h" #include "AliGRPObject.h" +#include "AliTPCROC.h" using namespace std; @@ -151,6 +152,8 @@ AliTPCcalibLaser::AliTPCcalibLaser(): 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 @@ -260,6 +263,9 @@ AliTPCcalibLaser::AliTPCcalibLaser(const Text_t *name, const Text_t *title, Bool // // 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 @@ -373,6 +379,9 @@ AliTPCcalibLaser::AliTPCcalibLaser(const AliTPCcalibLaser& calibLaser): // // 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 @@ -484,6 +493,9 @@ AliTPCcalibLaser::~AliTPCcalibLaser() { // 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 @@ -713,6 +725,8 @@ void AliTPCcalibLaser::MakeDistHisto(Int_t id){ xhis[11]= TMath::Sqrt(TMath::Abs(track->GetTPCsignal())); // } fHisLaser->Fill(xhis); + // + } void AliTPCcalibLaser::FitDriftV(){ @@ -2252,6 +2266,28 @@ void AliTPCcalibLaser::RefitLaserJW(Int_t id){ } } + // + // 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); + } } @@ -2682,8 +2718,8 @@ void AliTPCcalibLaser::DumpMeanInfo(Int_t run){ 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; } @@ -2715,6 +2751,22 @@ void AliTPCcalibLaser::DumpMeanInfo(Int_t run){ 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 @@ -3223,8 +3275,8 @@ Long64_t AliTPCcalibLaser::Merge(TCollection *li) { // 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); } @@ -3232,8 +3284,8 @@ Long64_t AliTPCcalibLaser::Merge(TCollection *li) { // 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); } @@ -3241,8 +3293,8 @@ Long64_t AliTPCcalibLaser::Merge(TCollection *li) { // 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); } @@ -3252,8 +3304,8 @@ Long64_t AliTPCcalibLaser::Merge(TCollection *li) { // 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); } @@ -3261,8 +3313,8 @@ Long64_t AliTPCcalibLaser::Merge(TCollection *li) { // 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); } @@ -3270,8 +3322,8 @@ Long64_t AliTPCcalibLaser::Merge(TCollection *li) { // 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); } @@ -3299,10 +3351,11 @@ Long64_t AliTPCcalibLaser::Merge(TCollection *li) { 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); @@ -3617,6 +3670,32 @@ void AliTPCcalibLaser::MakeFitHistos(){ 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){ @@ -3625,9 +3704,14 @@ 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 @@ -4034,246 +4118,506 @@ void AliTPCcalibLaser::SetBeamParameters(TVectorD& meanOffset, } - /* - 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="<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=(irowGetNRows(0))? sectorInner:sectorOuter; + Int_t npads=(irowGetNRows(0))? roc->GetNPads(sector,irow):roc->GetNPads(sector,irow-roc->GetNRows(0)); + (*vecPad)[irow]-=npads*0.5; + // + if ((irowGetNRows(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])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)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])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=(irowGetNRows(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()Number of clusters inner TH2F *fHisNclOut; //->Number of clusters outer TH2F *fHisNclIO; //->Number of cluster inner outer @@ -178,7 +183,7 @@ public: Double_t fFixedFitCside1; // Fixed drift v constant 1 - C side // private: - ClassDef(AliTPCcalibLaser,5) + ClassDef(AliTPCcalibLaser,6) }; -- 2.43.0