X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=TPC%2FAliTPCcalibAlign.cxx;h=ef84bc9b14c63e1b987f5011b9d676a8773002a4;hp=7a26695b6c7ecee5755ecdea8fd55ee038906ffb;hb=9f02b62b3f5da006ee8f2d7e1de20f0754fd9765;hpb=e81dc11276498f9f0adc95ec4f6a0943a552eae5 diff --git a/TPC/AliTPCcalibAlign.cxx b/TPC/AliTPCcalibAlign.cxx index 7a26695b6c7..ef84bc9b14c 100644 --- a/TPC/AliTPCcalibAlign.cxx +++ b/TPC/AliTPCcalibAlign.cxx @@ -18,40 +18,183 @@ // // // Class to make a internal alignemnt of TPC chambers // // +// Requierements - Warnings: +// 1. Before using this componenent the magnetic filed has to be set properly // +// 2. The systematic effects - unlinearities has to be understood +// +// If systematic and unlinearities are not under control +// the alignment is just effective alignment. Not second order corrction +// are calculated. +// +// The histograming of the edge effects and unlineratities integral part +// of the component (currently only in debug stream) +// +// 3 general type of linear transformation investigated (see bellow) +// +// By default only 6 parameter alignment to be used - other just for QA purposes + // Different linear tranformation investigated // 12 parameters - arbitrary linear transformation +// a00 a01 a02 a03 p[0] p[1] p[2] p[9] +// a10 a11 a12 a13 ==> p[3] p[4] p[5] p[10] +// a20 a21 a22 a23 p[6] p[7] p[8] p[11] +// // 9 parameters - scaling fixed to 1 +// a00 a01 a02 a03 1 p[0] p[1] p[6] +// a10 a11 a12 a13 ==> p[2] 1 p[3] p[7] +// a20 a21 a22 a23 p[4] p[5] 1 p[8] +// // 6 parameters - x-y rotation x-z, y-z tiliting -//// -//// Only the 12 parameter version works so far... -//// +// a00 a01 a02 a03 1 -p[0] 0 p[3] +// a10 a11 a12 a13 ==> p[0] 1 0 p[4] +// a20 a21 a22 a23 p[1] p[2] 1 p[5] +// +// +// Debug stream supported +// 0. Align - The main output of the Alignment component +// - Used for visualization of the misalignment between sectors +// - Results of the missalignment fit and the mean and sigmas of histograms +// stored there +// 1. Tracklet - StreamLevel >1 +// - Dump all information about tracklet match from sector1 to sector 2 +// - Default histogram residulas created in parallel +// - Check this streamer in case of suspicious content of these histograms +// 2. Track - StreamLevel>5 +// - For debugging of the edge effects +// - All information - extrapolation inside of one sectors +// - Created in order to distinguish between unlinearities inside of o +// sector and missalignment + +// +// +/* + gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros"); + gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+") + AliXRDPROOFtoolkit tool; + TChain * chain = tool.MakeChain("align.txt","Track",0,10200); + chain->Lookup(); + TCut cutA("abs(tp1.fP[1]-tp2.fP[1])<0.3&&abs(tp1.fP[0]-tp2.fP[0])<0.15&&abs(tp1.fP[3]-tp2.fP[3])<0.01&&abs(tp1.fP[2]-tp2.fP[2])<0.01"); + TCut cutS("s1%36==s2%36"); + + .x ~/UliStyle.C + .x $ALICE_ROOT/macros/loadlibsREC.C + + gSystem->Load("$ROOTSYS/lib/libXrdClient.so"); + gSystem->Load("libProof"); + gSystem->Load("libANALYSIS"); + gSystem->Load("libSTAT"); + gSystem->Load("libTPCcalib"); + // + // compare reference + TFile fcalib("CalibObjects.root"); + TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib"); + + AliTPCcalibAlign * align = ( AliTPCcalibAlign *)array->FindObject("alignTPC"); + // + // + align->EvalFitters(); + align->MakeTree("alignTree.root"); + TFile falignTree("alignTree.root"); + TTree * treeAlign = (TTree*)falignTree.Get("Align"); + + +*/ + //// //// #include "TLinearFitter.h" #include "AliTPCcalibAlign.h" +#include "AliTPCROC.h" +#include "AliTPCPointCorrection.h" +#include "AliTrackPointArray.h" + #include "AliExternalTrackParam.h" +#include "AliESDEvent.h" +#include "AliESDfriend.h" +#include "AliESDtrack.h" + #include "AliTPCTracklet.h" #include "TH1D.h" +#include "TH2F.h" +#include "THnSparse.h" #include "TVectorD.h" #include "TTreeStream.h" #include "TFile.h" +#include "TTree.h" #include "TF1.h" +#include "TGraphErrors.h" +#include "AliTPCclusterMI.h" +#include "AliTPCseed.h" +#include "AliTracker.h" +#include "TClonesArray.h" +#include "AliExternalComparison.h" +#include "AliLog.h" +#include "TFile.h" +#include "TProfile.h" +#include "TCanvas.h" +#include "TDatabasePDG.h" -#include +#include "TTreeStream.h" +#include "Riostream.h" #include using namespace std; +AliTPCcalibAlign* AliTPCcalibAlign::fgInstance = 0; ClassImp(AliTPCcalibAlign) + + + +AliTPCcalibAlign* AliTPCcalibAlign::Instance() +{ + // + // Singleton implementation + // Returns an instance of this class, it is created if neccessary + // + if (fgInstance == 0){ + fgInstance = new AliTPCcalibAlign(); + } + return fgInstance; +} + + + + AliTPCcalibAlign::AliTPCcalibAlign() - : fDphiHistArray(72*72), + : AliTPCcalibBase(), + fDphiHistArray(72*72), fDthetaHistArray(72*72), fDyHistArray(72*72), fDzHistArray(72*72), + // + fDyPhiHistArray(72*72), // array of residual histograms y -kYPhi + fDzThetaHistArray(72*72), // array of residual histograms z-z -kZTheta + fDphiZHistArray(72*72), // array of residual histograms phi -kPhiz + fDthetaZHistArray(72*72), // array of residual histograms theta -kThetaz + fDyZHistArray(72*72), // array of residual histograms y -kYz + fDzZHistArray(72*72), // array of residual histograms z -kZz fFitterArray12(72*72), fFitterArray9(72*72), - fFitterArray6(72*72) + fFitterArray6(72*72), + fMatrixArray12(72*72), + fMatrixArray9(72*72), + fMatrixArray6(72*72), + fCombinedMatrixArray6(72), + fNoField(kFALSE), + fXIO(0), + fXmiddle(0), + fXquadrant(0), + fArraySectorIntParam(36), // array of sector alignment parameters + fArraySectorIntCovar(36), // array of sector alignment covariances + // + // Kalman filter for global alignment + // + fSectorParamA(0), // Kalman parameter for A side + fSectorCovarA(0), // Kalman covariance for A side + fSectorParamC(0), // Kalman parameter for A side + fSectorCovarC(0), // Kalman covariance for A side + fUseInnerOuter(kTRUE)// flag- use Inner Outer sector for left righ alignment { // // Constructor @@ -59,6 +202,18 @@ AliTPCcalibAlign::AliTPCcalibAlign() for (Int_t i=0;i<72*72;++i) { fPoints[i]=0; } + AliTPCROC * roc = AliTPCROC::Instance(); + fXquadrant = roc->GetPadRowRadii(36,53); + fXmiddle = ( roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5; + fXIO = ( roc->GetPadRowRadii(0,roc->GetNRows(0)-1)+roc->GetPadRowRadii(36,0))*0.5; + fClusterDelta[0]=0; // cluster residuals - Y + fClusterDelta[1]=0; // cluster residuals - Z + + + fTrackletDelta[0]=0; // tracklet residuals + fTrackletDelta[1]=0; // tracklet residuals + fTrackletDelta[2]=0; // tracklet residuals + fTrackletDelta[3]=0; // tracklet residuals } AliTPCcalibAlign::AliTPCcalibAlign(const Text_t *name, const Text_t *title) @@ -67,9 +222,34 @@ AliTPCcalibAlign::AliTPCcalibAlign(const Text_t *name, const Text_t *title) fDthetaHistArray(72*72), fDyHistArray(72*72), fDzHistArray(72*72), + fDyPhiHistArray(72*72), // array of residual histograms y -kYPhi + fDzThetaHistArray(72*72), // array of residual histograms z-z -kZTheta + fDphiZHistArray(72*72), // array of residual histograms phi -kPhiz + fDthetaZHistArray(72*72), // array of residual histograms theta -kThetaz + fDyZHistArray(72*72), // array of residual histograms y -kYz + fDzZHistArray(72*72), // array of residual histograms z -kZz // fFitterArray12(72*72), fFitterArray9(72*72), - fFitterArray6(72*72) + fFitterArray6(72*72), + fMatrixArray12(72*72), + fMatrixArray9(72*72), + fMatrixArray6(72*72), + fCombinedMatrixArray6(72), + fNoField(kFALSE), + fXIO(0), + fXmiddle(0), + fXquadrant(0), + fArraySectorIntParam(36), // array of sector alignment parameters + fArraySectorIntCovar(36), // array of sector alignment covariances + // + // Kalman filter for global alignment + // + fSectorParamA(0), // Kalman parameter for A side + fSectorCovarA(0), // Kalman covariance for A side + fSectorParamC(0), // Kalman parameter for A side + fSectorCovarC(0), // Kalman covariance for A side + fUseInnerOuter(kTRUE)// flag- use Inner Outer sector for left righ alignment + { // // Constructor @@ -79,46 +259,556 @@ AliTPCcalibAlign::AliTPCcalibAlign(const Text_t *name, const Text_t *title) for (Int_t i=0;i<72*72;++i) { fPoints[i]=0; } + AliTPCROC * roc = AliTPCROC::Instance(); + fXquadrant = roc->GetPadRowRadii(36,53); + fXmiddle = ( roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5; + fXIO = ( roc->GetPadRowRadii(0,roc->GetNRows(0)-1)+roc->GetPadRowRadii(36,0))*0.5; + fClusterDelta[0]=0; // cluster residuals + fClusterDelta[1]=0; // cluster residuals + + fTrackletDelta[0]=0; // tracklet residuals + fTrackletDelta[1]=0; // tracklet residuals + fTrackletDelta[2]=0; // tracklet residuals + fTrackletDelta[3]=0; // tracklet residuals +} + + +AliTPCcalibAlign::AliTPCcalibAlign(const AliTPCcalibAlign &align) + :AliTPCcalibBase(align), + fDphiHistArray(align.fDphiHistArray), + fDthetaHistArray(align.fDthetaHistArray), + fDyHistArray(align.fDyHistArray), + fDzHistArray(align.fDzHistArray), + fDyPhiHistArray(align.fDyPhiHistArray), // array of residual histograms y -kYPhi + fDzThetaHistArray(align.fDzThetaHistArray), // array of residual histograms z-z -kZTheta + fDphiZHistArray(align.fDphiZHistArray), // array of residual histograms phi -kPhiz + fDthetaZHistArray(align.fDthetaZHistArray), // array of residual histograms theta -kThetaz + fDyZHistArray(align.fDyZHistArray), // array of residual histograms y -kYz + fDzZHistArray(align.fDzZHistArray), // array of residual histograms z -kZz + // + fFitterArray12(align.fFitterArray12), + fFitterArray9(align.fFitterArray9), + fFitterArray6(align.fFitterArray6), + + fMatrixArray12(align.fMatrixArray12), + fMatrixArray9(align.fMatrixArray9), + fMatrixArray6(align.fMatrixArray6), + fCombinedMatrixArray6(align.fCombinedMatrixArray6), + fNoField(align.fNoField), + fXIO(align.fXIO), + fXmiddle(align.fXmiddle), + fXquadrant(align.fXquadrant), + fArraySectorIntParam(align.fArraySectorIntParam), // array of sector alignment parameters + fArraySectorIntCovar(align.fArraySectorIntCovar), // array of sector alignment covariances + fSectorParamA(0), // Kalman parameter for A side + fSectorCovarA(0), // Kalman covariance for A side + fSectorParamC(0), // Kalman parameter for A side + fSectorCovarC(0), // Kalman covariance for A side + fUseInnerOuter(kTRUE)// flag- use Inner Outer sector for left righ alignment + +{ + // + // copy constructor - copy also the content + // + TH1 * his = 0; + TObjArray * arr0=0; + const TObjArray *arr1=0; + for (Int_t index =0; index<72*72; index++){ + for (Int_t iarray=0;iarray<10; iarray++){ + if (iarray==kY){ + arr0 = &fDyHistArray; + arr1 = &align.fDyHistArray; + } + if (iarray==kZ){ + arr0 = &fDzHistArray; + arr1 = &align.fDzHistArray; + } + if (iarray==kPhi){ + arr0 = &fDphiHistArray; + arr1 = &align.fDphiHistArray; + } + if (iarray==kTheta){ + arr0 = &fDthetaHistArray; + arr1 = &align.fDthetaHistArray; + } + if (iarray==kYz){ + arr0 = &fDyZHistArray; + arr1 = &align.fDyZHistArray; + } + if (iarray==kZz){ + arr0 = &fDzZHistArray; + arr1 = &align.fDzZHistArray; + } + if (iarray==kPhiZ){ + arr0 = &fDphiZHistArray; + arr1 = &align.fDphiZHistArray; + } + if (iarray==kThetaZ){ + arr0 = &fDthetaZHistArray; + arr1 = &align.fDthetaZHistArray; + } + + if (iarray==kYPhi){ + arr0 = &fDyPhiHistArray; + arr1 = &align.fDyPhiHistArray; + } + if (iarray==kZTheta){ + arr0 = &fDzThetaHistArray; + arr1 = &align.fDzThetaHistArray; + } + + if (arr1->At(index)) { + his = (TH1*)arr1->At(index)->Clone(); + his->SetDirectory(0); + arr0->AddAt(his,index); + } + } + } + // + // + // + if (align.fSectorParamA){ + fSectorParamA = (TMatrixD*)align.fSectorParamA->Clone(); + fSectorParamA = (TMatrixD*)align.fSectorCovarA->Clone(); + fSectorParamC = (TMatrixD*)align.fSectorParamA->Clone(); + fSectorParamC = (TMatrixD*)align.fSectorCovarA->Clone(); + } + fClusterDelta[0]=0; // cluster residuals + fClusterDelta[1]=0; // cluster residuals + + fTrackletDelta[0]=0; // tracklet residuals + fTrackletDelta[1]=0; // tracklet residuals + fTrackletDelta[2]=0; // tracklet residuals + fTrackletDelta[3]=0; // tracklet residuals } + AliTPCcalibAlign::~AliTPCcalibAlign() { // // destructor // + fDphiHistArray.SetOwner(kTRUE); // array of residual histograms phi -kPhi + fDthetaHistArray.SetOwner(kTRUE); // array of residual histograms theta -kTheta + fDyHistArray.SetOwner(kTRUE); // array of residual histograms y -kY + fDzHistArray.SetOwner(kTRUE); // array of residual histograms z -kZ + // + fDyPhiHistArray.SetOwner(kTRUE); // array of residual histograms y -kYPhi + fDzThetaHistArray.SetOwner(kTRUE); // array of residual histograms z-z -kZTheta + // + fDphiZHistArray.SetOwner(kTRUE); // array of residual histograms phi -kPhiz + fDthetaZHistArray.SetOwner(kTRUE); // array of residual histograms theta -kThetaz + fDyZHistArray.SetOwner(kTRUE); // array of residual histograms y -kYz + fDzZHistArray.SetOwner(kTRUE); // array of residual histograms z -kZz + + fDphiHistArray.Delete(); // array of residual histograms phi -kPhi + fDthetaHistArray.Delete(); // array of residual histograms theta -kTheta + fDyHistArray.Delete(); // array of residual histograms y -kY + fDzHistArray.Delete(); // array of residual histograms z -kZ + // + fDyPhiHistArray.Delete(); // array of residual histograms y -kYPhi + fDzThetaHistArray.Delete(); // array of residual histograms z-z -kZTheta + // + fDphiZHistArray.Delete(); // array of residual histograms phi -kPhiz + fDthetaZHistArray.Delete(); // array of residual histograms theta -kThetaz + fDyZHistArray.Delete(); // array of residual histograms y -kYz + fDzZHistArray.Delete(); // array of residual histograms z -kZz + + fFitterArray12.SetOwner(kTRUE); // array of fitters + fFitterArray9.SetOwner(kTRUE); // array of fitters + fFitterArray6.SetOwner(kTRUE); // array of fitters + // + fMatrixArray12.SetOwner(kTRUE); // array of transnformtation matrix + fMatrixArray9.SetOwner(kTRUE); // array of transnformtation matrix + fMatrixArray6.SetOwner(kTRUE); // array of transnformtation matrix + // + fFitterArray12.Delete(); // array of fitters + fFitterArray9.Delete(); // array of fitters + fFitterArray6.Delete(); // array of fitters + // + fMatrixArray12.Delete(); // array of transnformtation matrix + fMatrixArray9.Delete(); // array of transnformtation matrix + fMatrixArray6.Delete(); // array of transnformtation matrix + + + fArraySectorIntParam.SetOwner(kTRUE); // array of sector alignment parameters + fArraySectorIntCovar.SetOwner(kTRUE); // array of sector alignment covariances + fArraySectorIntParam.Delete(); // array of sector alignment parameters + fArraySectorIntCovar.Delete(); // array of sector alignment covariances + for (Int_t i=0; i<2; i++){ + delete fClusterDelta[i]; // cluster residuals + } + + for (Int_t i=0; i<4; i++){ + delete fTrackletDelta[i]; // tracklet residuals + } + + +} + +void AliTPCcalibAlign::Process(AliESDEvent *event) { + // + // Process pairs of cosmic tracks + // + if (!fClusterDelta[0]) MakeResidualHistos(); + if (!fTrackletDelta[0]) MakeResidualHistosTracklet(); + // + fCurrentEvent=event; + ExportTrackPoints(event); // export track points for external calibration + const Int_t kMaxTracks =6; + const Int_t kminCl = 40; + AliESDfriend *eESDfriend=static_cast(event->FindListObject("AliESDfriend")); + if (!eESDfriend) return; + Int_t ntracks=event->GetNumberOfTracks(); + Float_t dca0[2]; + Float_t dca1[2]; + // + // + // process seeds + // + for (Int_t i0=0;i0GetTrack(i0); + AliESDfriendTrack *friendTrack = 0; + TObject *calibObject=0; + AliTPCseed *seed0 = 0; + // + friendTrack = (AliESDfriendTrack *)eESDfriend->GetTrack(i0);; + if (!friendTrack) continue; + for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) { + if ((seed0=dynamic_cast(calibObject))) break; + } + if (!seed0) continue; + fCurrentTrack=track0; + fCurrentFriendTrack=friendTrack; + fCurrentSeed=seed0; + fCurrentEvent=event; + ProcessSeed(seed0); + } + // + // process cosmic pairs + // + if (ntracks>kMaxTracks) return; + // + //select pairs - for alignment + for (Int_t i0=0;i0GetTrack(i0); + // if (track0->GetTPCNcls()GetImpactParameters(dca0[0],dca0[1]); + // if (TMath::Abs(dca0[0])>30) continue; + // + for (Int_t i1=0;i1GetTrack(i1); + // if (track1->GetTPCNcls()GetImpactParameters(dca1[0],dca1[1]); + // fast cuts on dca and theta + // if (TMath::Abs(dca1[0]+dca0[0])>15) continue; + // if (TMath::Abs(dca1[1]-dca0[1])>15) continue; + if (TMath::Abs(track0->GetParameter()[3]+track1->GetParameter()[3])>0.1) continue; + // + AliESDfriendTrack *friendTrack = 0; + TObject *calibObject=0; + AliTPCseed *seed0 = 0,*seed1=0; + // + friendTrack = (AliESDfriendTrack *)eESDfriend->GetTrack(i0);; + if (!friendTrack) continue; + for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) { + if ((seed0=dynamic_cast(calibObject))) break; + } + friendTrack = (AliESDfriendTrack *)eESDfriend->GetTrack(i1);; + if (!friendTrack) continue; + for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) { + if ((seed1=dynamic_cast(calibObject))) break; + } + if (!seed0) continue; + // + // + if (!seed1) continue; + Int_t nclsectors0[72], nclsectors1[72]; + for (Int_t isec=0;isec<72;isec++){ + nclsectors0[isec]=0; + nclsectors1[isec]=0; + } + for (Int_t i=0;i<160;i++){ + AliTPCclusterMI *c0=seed0->GetClusterPointer(i); + AliTPCclusterMI *c1=seed1->GetClusterPointer(i); + if (c0) nclsectors0[c0->GetDetector()]+=1; + if (c1) nclsectors1[c1->GetDetector()]+=1; + } + + for (Int_t isec0=0; isec0<72;isec0++){ + if (nclsectors0[isec0]1.0) isOK=kFALSE; + if (TMath::Abs(dp1)>0.02) isOK=kFALSE; + if (TMath::Abs(dp3)>0.02) isOK=kFALSE; + if (TMath::Abs(pp0)>6) isOK=kFALSE; + if (TMath::Abs(pp1)>6) isOK=kFALSE; + if (TMath::Abs(pp3)>6) isOK=kFALSE; + // + if (isOK){ + FillHisto(parLine0,parLine1,s0,s1); + ProcessAlign(parLine0,parLine1,s0,s1); + UpdateKalman(s0,s1,par0, cov0, par1, cov1); + } + if (fStreamLevel>0){ + TTreeSRedirector *cstream = GetDebugStreamer(); + if (cstream){ + (*cstream)<<"cosmic"<< + "isOK="<(event->FindListObject("AliESDfriend")); + if (!eESDfriend) return; + Int_t ntracks=event->GetNumberOfTracks(); + Int_t kMaxTracks=4; // maximal number of tracks for cosmic pairs + Int_t kMinVertexTracks=5; // maximal number of tracks for vertex mesurement + + //cuts + const Int_t kminCl = 60; + const Int_t kminClSum = 120; + //const Double_t kDistY = 5; + // const Double_t kDistZ = 40; + const Double_t kDistTh = 0.05; + const Double_t kDistThS = 0.002; + const Double_t kDist1Pt = 0.1; + const Double_t kMaxD0 =3; // max distance to the primary vertex + const Double_t kMaxD1 =5; // max distance to the primary vertex + const AliESDVertex *tpcVertex = 0; + // get the primary vertex TPC + if (ntracks>kMinVertexTracks) { + tpcVertex = event->GetPrimaryVertexSPD(); + if (tpcVertex->GetNContributors()GetTrack(i0); + if (!track0) continue; + if ((track0->GetStatus() & AliESDtrack::kTPCrefit)==0) continue; + if (track0->GetOuterParam()==0) continue; + if (track0->GetInnerParam()==0) continue; + if (TMath::Abs(track0->GetInnerParam()->GetSigned1Pt()-track0->GetOuterParam()->GetSigned1Pt())>kDist1Pt) continue; + if (TMath::Abs(track0->GetInnerParam()->GetSigned1Pt())>kDist1Pt) continue; + if (TMath::Abs(track0->GetInnerParam()->GetTgl()-track0->GetOuterParam()->GetTgl())>kDistThS) continue; + AliESDtrack *track1P = 0; + if (track0->GetTPCNcls()GetImpactParameters(dca0[0],dca0[1]); + index0=i0; + index1=-1; + // + if (ntracksGetTrack(i1); + if (!track1) continue; + if ((track1->GetStatus() & AliESDtrack::kTPCrefit)==0) continue; + if (track1->GetOuterParam()==0) continue; + if (track1->GetInnerParam()==0) continue; + if (track1->GetTPCNcls()GetInnerParam()->GetSigned1Pt()-track1->GetOuterParam()->GetSigned1Pt())>kDist1Pt) continue; + if (TMath::Abs(track1->GetInnerParam()->GetTgl()-track1->GetOuterParam()->GetTgl())>kDistThS) continue; + if (TMath::Abs(track1->GetInnerParam()->GetSigned1Pt())>kDist1Pt) continue; + //track1->GetImpactParameters(dca1[0],dca1[1]); + //if (TMath::Abs(dca1[0]-dca0[0])>kDistY) continue; + //if (TMath::Abs(dca1[1]-dca0[1])>kDistZ) continue; + if (TMath::Abs(track0->GetTgl()+track1->GetTgl())>kDistTh) continue; + if (TMath::Abs(track0->GetSigned1Pt()+track1->GetSigned1Pt())>kDist1Pt) continue; + track1P = track1; + index1=i1; + } + AliESDfriendTrack *friendTrack = 0; + TObject *calibObject=0; + AliTPCseed *seed0 = 0,*seed1=0; + // + friendTrack = (AliESDfriendTrack *)eESDfriend->GetTrack(index0);; + if (!friendTrack) continue; + for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) { + if ((seed0=dynamic_cast(calibObject))) break; + } + if (index1>0){ + friendTrack = (AliESDfriendTrack *)eESDfriend->GetTrack(index1);; + if (!friendTrack) continue; + for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) { + if ((seed1=dynamic_cast(calibObject))) break; + } + } + // + Int_t npoints=0, ncont=0; + if (seed0) {npoints+=seed0->GetNumberOfClusters(); ncont++;} + if (seed1) {npoints+=seed1->GetNumberOfClusters(); ncont++;} + if (npointsGetXYZ(dxyz); + tpcVertex->GetCovarianceMatrix(dc); + Float_t xyz[3]={dxyz[0],dxyz[1],dxyz[2]}; + Float_t cov[6]={dc[0]+1,dc[1],dc[2]+1,dc[3],dc[4],dc[5]+100.}; + AliTrackPoint point(xyz,cov,73); // add point to not existing volume + Float_t dtpc[2],dcov[3]; + track0->GetImpactParametersTPC(dtpc,dcov); + if (TMath::Abs(dtpc[0])>kMaxD0) continue; + if (TMath::Abs(dtpc[1])>kMaxD1) continue; + } + + if (seed0) for (Int_t icl = 0; icl<160; icl++){ + AliTPCclusterMI *cluster=seed0->GetClusterPointer(icl); + if (!cluster) continue; + Float_t xyz[3]; + Float_t cov[6]; + cluster->GetGlobalXYZ(xyz); + cluster->GetGlobalCov(cov); + AliTrackPoint point(xyz,cov,cluster->GetDetector()); + array.AddPoint(npoints, &point); + if (cpoint>=npoints) continue; //shoul not happen + array.AddPoint(cpoint, &point); + cpoint++; + } + if (seed1) for (Int_t icl = 0; icl<160; icl++){ + AliTPCclusterMI *cluster=seed1->GetClusterPointer(icl); + if (!cluster) continue; + Float_t xyz[3]; + Float_t cov[6]; + cluster->GetGlobalXYZ(xyz); + cluster->GetGlobalCov(cov); + AliTrackPoint point(xyz,cov,cluster->GetDetector()); + array.AddPoint(npoints, &point); + if (cpoint>=npoints) continue; //shoul not happen + array.AddPoint(cpoint, &point); + cpoint++; + } + // + // + // + TTreeSRedirector *cstream = GetDebugStreamer(); + if (cstream){ + Bool_t isVertex=(tpcVertex)? kTRUE:kFALSE; + Double_t tof0=track0->GetTOFsignal(); + Double_t tof1=(track1P) ? track1P->GetTOFsignal(): 0; + static AliExternalTrackParam dummy; + AliExternalTrackParam *p0In = &dummy; + AliExternalTrackParam *p1In = &dummy; + AliExternalTrackParam *p0Out = &dummy; + AliExternalTrackParam *p1Out = &dummy; + AliESDVertex vdummy; + AliESDVertex *pvertex= (tpcVertex)? (AliESDVertex *)tpcVertex: &vdummy; + if (track0) { + p0In= new AliExternalTrackParam(*track0); + p0Out=new AliExternalTrackParam(*(track0->GetOuterParam())); + } + if (track1P) { + p1In= new AliExternalTrackParam(*track1P); + p1Out=new AliExternalTrackParam(*(track1P->GetOuterParam())); + } + + (*cstream)<<"trackPoints"<< + "run="<(tracklets[0]); - AliTPCTracklet *t2=static_cast(tracklets[1]); - if (t1->GetSector()>t2->GetSector()) { - AliTPCTracklet* tmp=t1; - t1=t2; - t2=tmp; - } - AliExternalTrackParam *common1=0,*common2=0; - if (AliTPCTracklet::PropagateToMeanX(*t1,*t2,common1,common2)) - ProcessTracklets(*common1,*common2,t1->GetSector(),t2->GetSector()); - delete common1; - delete common2; - } - + Int_t ntracklets = tracklets.GetEntries(); + if (ntracklets<2) return; + // + // + for (Int_t i1=0;i1(tracklets[i1]); + AliTPCTracklet *t2=static_cast(tracklets[i2]); + AliExternalTrackParam *common1=0,*common2=0; + if (AliTPCTracklet::PropagateToMeanX(*t1,*t2,common1,common2)){ + ProcessTracklets(*common1,*common2,seed, t1->GetSector(),t2->GetSector()); + UpdateAlignSector(seed,t1->GetSector()); + } + delete common1; + delete common2; + } } void AliTPCcalibAlign::Analyze(){ @@ -134,91 +824,452 @@ void AliTPCcalibAlign::Terminate(){ // Terminate function // call base terminate + Eval of fitters // - if (GetDebugLevel()>0) Info("AliTPCcalibAlign","Teminate"); + Info("AliTPCcalibAlign","Terminate"); EvalFitters(); AliTPCcalibBase::Terminate(); } +void AliTPCcalibAlign::UpdatePointCorrection(AliTPCPointCorrection * correction){ + // + // Update point correction with alignment coefficients + // + for (Int_t isec=0;isec<36;isec++){ + TMatrixD * matCorr = (TMatrixD*)(correction->fArraySectorIntParam.At(isec)); + TMatrixD * matAlign = (TMatrixD*)(fArraySectorIntParam.At(isec)); + TMatrixD * matAlignCovar = (TMatrixD*)(fArraySectorIntCovar.At(isec)); + if (!matAlign) continue; + if (!matCorr) { + correction->fArraySectorIntParam.AddAt(matAlign->Clone(),isec); + correction->fArraySectorIntCovar.AddAt(matAlignCovar->Clone(),isec); + continue; + } + (*matCorr)+=(*matAlign); + correction->fArraySectorIntCovar.AddAt(matAlignCovar->Clone(),isec); + } + // + +} void AliTPCcalibAlign::ProcessTracklets(const AliExternalTrackParam &tp1, const AliExternalTrackParam &tp2, + const AliTPCseed * seed, Int_t s1,Int_t s2) { - - // - // - // - FillHisto(tp1,tp2,s1,s2); // // Process function to fill fitters // - Double_t t1[5],t2[5]; - Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4]; - Double_t &x2=t2[0], &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4]; + Double_t t1[10],t2[10]; + Double_t &x1=t1[0], &y1=t1[1], &z1=t1[3], &dydx1=t1[2], &dzdx1=t1[4]; + Double_t &x2=t2[0], &y2=t2[1], &z2=t2[3], &dydx2=t2[2], &dzdx2=t2[4]; x1 =tp1.GetX(); y1 =tp1.GetY(); z1 =tp1.GetZ(); Double_t snp1=tp1.GetSnp(); - dydx1=snp1/TMath::Sqrt(1.-snp1*snp1); + dydx1=snp1/TMath::Sqrt((1.-snp1)*(1.+snp1)); Double_t tgl1=tp1.GetTgl(); // dz/dx = 1/(cos(theta)*cos(phi)) - dzdx1=tgl1/TMath::Sqrt(1.-snp1*snp1); + dzdx1=tgl1/TMath::Sqrt((1.-snp1)*(1.+snp1)); x2 =tp2.GetX(); y2 =tp2.GetY(); z2 =tp2.GetZ(); Double_t snp2=tp2.GetSnp(); - dydx2=snp2/TMath::Sqrt(1.-snp2*snp2); + dydx2=snp2/TMath::Sqrt((1.-snp2)*(1.+snp2)); Double_t tgl2=tp2.GetTgl(); - dzdx2=tgl2/TMath::Sqrt(1.-snp2*snp2); + dzdx2=tgl2/TMath::Sqrt((1.-snp2)*(1.+snp2)); + + // + // Kalman parameters + // + t1[0]-=fXIO; + t2[0]-=fXIO; + // errors + t1[5]=0; t2[5]=0; + t1[6]=TMath::Sqrt(tp1.GetSigmaY2()); + t1[7]=TMath::Sqrt(tp1.GetSigmaSnp2()); + t1[8]=TMath::Sqrt(tp1.GetSigmaZ2()); + t1[9]=TMath::Sqrt(tp1.GetSigmaTgl2()); + + t2[6]=TMath::Sqrt(tp2.GetSigmaY2()); + t2[7]=TMath::Sqrt(tp2.GetSigmaSnp2()); + t2[8]=TMath::Sqrt(tp2.GetSigmaZ2()); + t2[9]=TMath::Sqrt(tp2.GetSigmaTgl2()); // + // linear parameters // + Double_t parLine1[10]; + Double_t parLine2[10]; + TMatrixD par1(4,1),cov1(4,4),par2(4,1),cov2(4,4); + Bool_t useInnerOuter = kFALSE; + if (s1%36!=s2%36) useInnerOuter = fUseInnerOuter; // for left - right alignment bot sectors refit can be used if specified + Int_t nl1 = RefitLinear(seed,s1, parLine1, s1,par1,cov1,tp1.GetX(), useInnerOuter); + Int_t nl2 = RefitLinear(seed,s2, parLine2, s1,par2,cov2,tp1.GetX(), useInnerOuter); + parLine1[0]=tp1.GetX()-fXIO; // parameters in IROC-OROC boundary + parLine2[0]=tp1.GetX()-fXIO; // parameters in IROC-OROC boundary // - if (fStreamLevel>1){ + // + // + Int_t accept = AcceptTracklet(tp1,tp2); + Int_t acceptLinear = AcceptTracklet(parLine1,parLine2); + + + if (fStreamLevel>1 && seed){ TTreeSRedirector *cstream = GetDebugStreamer(); if (cstream){ static TVectorD vec1(5); static TVectorD vec2(5); + static TVectorD vecL1(9); + static TVectorD vecL2(9); vec1.SetElements(t1); vec2.SetElements(t2); + vecL1.SetElements(parLine1); + vecL2.SetElements(parLine2); AliExternalTrackParam *p1 = &((AliExternalTrackParam&)tp1); AliExternalTrackParam *p2 = &((AliExternalTrackParam&)tp2); (*cstream)<<"Tracklet"<< + "accept="<10 && nl2>10 &&(acceptLinear==0)){ + if (seed) ProcessDiff(tp1,tp2, seed,s1,s2); + 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); + //UpdateKalman(s1,s2,par1, cov1, par2, cov2); - OBSOLETE to be removed - 50 % of time here + } + } + } + if (accept>0) return; + // + // fill resolution histograms - previous cut included + if (TMath::Abs(fMagF)>0.005){ + // + // use Kalman if mag field + // + if (seed) { + ProcessDiff(tp1,tp2, seed,s1,s2); + FillHisto((AliExternalTrackParam*)&tp1,(AliExternalTrackParam*)&tp2,s1,s2); + FillHisto((AliExternalTrackParam*)&tp2,(AliExternalTrackParam*)&tp1,s2,s1); + } + FillHisto(t1,t2,s1,s2); + ProcessAlign(t1,t2,s1,s2); + } +} + +void AliTPCcalibAlign::ProcessAlign(Double_t * t1, + Double_t * t2, + Int_t s1,Int_t s2){ + // + // Do intersector alignment + // + //Process12(t1,t2,GetOrMakeFitter12(s1,s2)); + //Process9(t1,t2,GetOrMakeFitter9(s1,s2)); Process6(t1,t2,GetOrMakeFitter6(s1,s2)); ++fPoints[GetIndex(s1,s2)]; } +void AliTPCcalibAlign::ProcessTree(TTree * chainTracklet, AliExternalComparison *comp){ + // + // Process the debug streamer tree + // Possible to modify selection criteria + // Used with entry list + // + TTreeSRedirector * cstream = new TTreeSRedirector("aligndump.root"); + + AliTPCcalibAlign *align = this; + // + TVectorD * vec1 = 0; + TVectorD * vec2 = 0; + AliExternalTrackParam * tp1 = 0; + AliExternalTrackParam * tp2 = 0; + Int_t s1 = 0; + Int_t s2 = 0; + Int_t npoints =0; + { + Int_t entries=chainTracklet->GetEntries(); + for (Int_t i=0; i< entries; i++){ + chainTracklet->GetBranch("tp1.")->SetAddress(&tp1); + chainTracklet->GetBranch("tp2.")->SetAddress(&tp2); + chainTracklet->GetBranch("v1.")->SetAddress(&vec1); + chainTracklet->GetBranch("v2.")->SetAddress(&vec2); + chainTracklet->GetBranch("s1")->SetAddress(&s1); + chainTracklet->GetBranch("s2")->SetAddress(&s2); + chainTracklet->GetEntry(i); + if (!vec1) continue; + if (!vec2) continue; + if (!tp1) continue; + if (!tp2) continue; + if (!vec1->GetMatrixArray()) continue; + if (!vec2->GetMatrixArray()) continue; + // make a local copy + AliExternalTrackParam par1(*tp1); + AliExternalTrackParam par2(*tp2); + TVectorD svec1(*vec1); + TVectorD svec2(*vec2); + // + if (s1==s2) continue; + if (i%100==0) printf("%d\t%d\t%d\t%d\t\n",i, npoints,s1,s2); + AliExternalTrackParam cpar1(par1); + AliExternalTrackParam cpar2(par2); + Constrain1Pt(cpar1,par2,fNoField); + Constrain1Pt(cpar2,par1,fNoField); + Bool_t acceptComp = kFALSE; + if (comp) acceptComp=comp->AcceptPair(&par1,&par2); + if (comp) acceptComp&=comp->AcceptPair(&cpar1,&cpar2); + // + Int_t reject = align->AcceptTracklet(par1,par2); + Int_t rejectC =align->AcceptTracklet(cpar1,cpar2); + + if (1||fStreamLevel>0){ + (*cstream)<<"Tracklet"<< + "s1="<Process(&cpar1,&cpar2); + // + if (reject>0 || rejectC>0) continue; + npoints++; + align->ProcessTracklets(cpar1,cpar2,0,s1,s2); + align->ProcessTracklets(cpar2,cpar1,0,s2,s1); + } + } + delete cstream; +} + + +Int_t AliTPCcalibAlign::AcceptTracklet(const AliExternalTrackParam &p1, + const AliExternalTrackParam &p2) const +{ + + // + // Accept pair of tracklets? + // + /* + // resolution cuts + TCut cutS0("sqrt(tp2.fC[0]+tp1.fC[0])<0.2"); + TCut cutS1("sqrt(tp2.fC[2]+tp1.fC[2])<0.2"); + TCut cutS2("sqrt(tp2.fC[5]+tp1.fC[5])<0.01"); + TCut cutS3("sqrt(tp2.fC[9]+tp1.fC[9])<0.01"); + TCut cutS4("sqrt(tp2.fC[14]+tp1.fC[14])<0.25"); + TCut cutS=cutS0+cutS1+cutS2+cutS3+cutS4; + // + // parameters matching cuts + TCut cutP0("abs(tp1.fP[0]-tp2.fP[0])<0.6"); + TCut cutP1("abs(tp1.fP[1]-tp2.fP[1])<0.6"); + TCut cutP2("abs(tp1.fP[2]-tp2.fP[2])<0.03"); + TCut cutP3("abs(tp1.fP[3]-tp2.fP[3])<0.03"); + TCut cutP4("abs(tp1.fP[4]-tp2.fP[4])<0.5"); + TCut cutPP4("abs(tp1.fP[4]-tp2.fP[4])/sqrt(tp2.fC[14]+tp1.fC[14])<3"); + TCut cutP=cutP0+cutP1+cutP2+cutP3+cutP4+cutPP4; + */ + // + // resolution cuts + Int_t reject=0; + const Double_t *cp1 = p1.GetCovariance(); + const Double_t *cp2 = p2.GetCovariance(); + if (TMath::Sqrt(cp1[0]+cp2[0])>0.2) reject|=1;; + if (TMath::Sqrt(cp1[2]+cp2[2])>0.2) reject|=2; + if (TMath::Sqrt(cp1[5]+cp2[5])>0.01) reject|=4; + if (TMath::Sqrt(cp1[9]+cp2[9])>0.01) reject|=8; + if (TMath::Sqrt(cp1[14]+cp2[14])>0.2) reject|=16; + + //parameters difference + const Double_t *tp1 = p1.GetParameter(); + const Double_t *tp2 = p2.GetParameter(); + if (TMath::Abs(tp1[0]-tp2[0])>0.6) reject|=32; + if (TMath::Abs(tp1[1]-tp2[1])>0.6) reject|=64; + if (TMath::Abs(tp1[2]-tp2[2])>0.03) reject|=128; + if (TMath::Abs(tp1[3]-tp2[3])>0.03) reject|=526; + if (TMath::Abs(tp1[4]-tp2[4])>0.4) reject|=1024; + if (TMath::Abs(tp1[4]-tp2[4])/TMath::Sqrt(cp1[14]+cp2[14])>4) reject|=2048; + + // + if (TMath::Abs(tp2[1])>235) reject|=2*4096; + + if (fNoField){ + + } + + return reject; +} + + +Int_t AliTPCcalibAlign::AcceptTracklet(const Double_t *t1, const Double_t *t2) const +{ + // + // accept tracklet - + // dist cut + 6 sigma cut + // + Double_t dy = t2[1]-t1[1]; + Double_t dphi = t2[2]-t1[2]; + Double_t dz = t2[3]-t1[3]; + Double_t dtheta = t2[4]-t1[4]; + // + Double_t sy = TMath::Sqrt(t1[6]*t1[6]+t2[6]*t2[6]+0.05*0.05); + Double_t sdydx = TMath::Sqrt(t1[7]*t1[7]+t2[7]*t2[7]+0.001*0.001); + Double_t sz = TMath::Sqrt(t1[8]*t1[8]+t2[8]*t2[8]+0.05*0.05); + Double_t sdzdx = TMath::Sqrt(t1[9]*t1[9]+t2[9]*t2[9]+0.001*0.001); + // + Int_t reject=0; + if (TMath::Abs(dy)>1.) reject|=2; + if (TMath::Abs(dphi)>0.1) reject|=4; + if (TMath::Abs(dz)>1.) reject|=8; + if (TMath::Abs(dtheta)>0.1) reject|=16; + // + if (TMath::Abs(dy/sy)>6) reject|=32; + if (TMath::Abs(dphi/sdydx)>6) reject|=64; + if (TMath::Abs(dz/sz)>6) reject|=128; + if (TMath::Abs(dtheta/sdzdx)>6) reject|=256; + return reject; +} + + +void AliTPCcalibAlign::ProcessDiff(const AliExternalTrackParam &t1, + const AliExternalTrackParam &t2, + const AliTPCseed *seed, + Int_t s1,Int_t s2) +{ + // + // Process local residuals function + // + TVectorD vecX(160); + TVectorD vecY(160); + TVectorD vecZ(160); + TVectorD vecClY(160); + TVectorD vecClZ(160); + TClonesArray arrCl("AliTPCclusterMI",160); + arrCl.ExpandCreateFast(160); + Int_t count1=0, count2=0; + + for (Int_t i=0;i<160;++i) { + AliTPCclusterMI *c=seed->GetClusterPointer(i); + vecX[i]=0; + vecY[i]=0; + vecZ[i]=0; + if (!c) continue; + AliTPCclusterMI & cl = (AliTPCclusterMI&) (*arrCl[i]); + if (c->GetDetector()!=s1 && c->GetDetector()!=s2) continue; + vecClY[i] = c->GetY(); + vecClZ[i] = c->GetZ(); + cl=*c; + const AliExternalTrackParam *par = (c->GetDetector()==s1)? &t1:&t2; + if (c->GetDetector()==s1) ++count1; + if (c->GetDetector()==s2) ++count2; + Double_t gxyz[3],xyz[3]; + t1.GetXYZ(gxyz); + Float_t bz = AliTracker::GetBz(gxyz); + par->GetYAt(c->GetX(), bz, xyz[1]); + par->GetZAt(c->GetX(), bz, xyz[2]); + vecX[i] = c->GetX(); + vecY[i]= xyz[1]; + vecZ[i]= xyz[2]; + } + // + // + if (fStreamLevel>5 && count1>10 && count2>10){ + // + // huge output - cluster residuals to be investigated + // + TTreeSRedirector *cstream = GetDebugStreamer(); + AliExternalTrackParam *p1 = &((AliExternalTrackParam&)t1); + AliExternalTrackParam *p2 = &((AliExternalTrackParam&)t2); + /* + + Track->Draw("Cl[].fY-vtY.fElements:vtY.fElements-vtX.fElements*tan(pi/18.)>>his(100,-10,0)","Cl.fY!=0&&abs(Cl.fY-vtY.fElements)<1","prof"); + + */ + + if (cstream){ + (*cstream)<<"Track"<< + "run="< p[3] p[4] p[5] p[10] - // a20 a21 a22 a23 p[6] p[7] p[8] p[11] - const Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4]; - const Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4]; + // a00 a01 a02 a03 p[0] p[1] p[2] p[9] + // a10 a11 a12 a13 ==> p[3] p[4] p[5] p[10] + // a20 a21 a22 a23 p[6] p[7] p[8] p[11] + - // TODO: - Double_t sy = 0.1; - Double_t sz = 0.1; - Double_t sdydx = 0.001; - Double_t sdzdx = 0.001; + + const Double_t &x1=t1[0], &y1=t1[1], &z1=t1[3], &dydx1=t1[2], &dzdx1=t1[4]; + const Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[3], &dydx2=t2[2], &dzdx2=t2[4]; + + // + Double_t sy = TMath::Sqrt(t1[6]*t1[6]+t2[6]*t2[6]); + Double_t sdydx = TMath::Sqrt(t1[7]*t1[7]+t2[7]*t2[7]); + Double_t sz = TMath::Sqrt(t1[8]*t1[8]+t2[8]*t2[8]); + Double_t sdzdx = TMath::Sqrt(t1[9]*t1[9]+t2[9]*t2[9]); Double_t p[12]; Double_t value; @@ -276,175 +1327,179 @@ void AliTPCcalibAlign::Process12(const Double_t *t1, fitter->AddPoint(p,value,sdzdx); } -void AliTPCcalibAlign::Process9(Double_t *t1, - Double_t *t2, - TLinearFitter *fitter) { +void AliTPCcalibAlign::Process9(const Double_t * const t1, + const Double_t * const t2, + TLinearFitter *fitter) const +{ // x2 = a00*x1 + a01*y1 + a02*z1 + a03 // y2 = a10*x1 + a11*y1 + a12*z1 + a13 // z2 = a20*x1 + a21*y1 + a22*z1 + a23 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1) // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1) // - // a00 a01 a02 a03 p[0] p[1] p[2] p[9] - // a10 a11 a12 a13 ==> p[3] p[4] p[5] p[10] - // a20 a21 a22 a23 p[6] p[7] p[8] p[11] - Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4]; - Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4]; + // a00 a01 a02 a03 1 p[0] p[1] p[6] + // a10 a11 a12 a13 ==> p[2] 1 p[3] p[7] + // a20 a21 a21 a23 p[4] p[5] 1 p[8] - // TODO: - Double_t sy = 1.; - Double_t sz = 1.; - Double_t sdydx = 1.; - Double_t sdzdx = 1.; + const Double_t &x1=t1[0], &y1=t1[1], &z1=t1[3], &dydx1=t1[2], &dzdx1=t1[4]; + const Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[3], &dydx2=t2[2], &dzdx2=t2[4]; + // + Double_t sy = TMath::Sqrt(t1[6]*t1[6]+t2[6]*t2[6]); + Double_t sdydx = TMath::Sqrt(t1[7]*t1[7]+t2[7]*t2[7]); + Double_t sz = TMath::Sqrt(t1[8]*t1[8]+t2[8]*t2[8]); + Double_t sdzdx = TMath::Sqrt(t1[9]*t1[9]+t2[9]*t2[9]); + + // Double_t p[12]; Double_t value; // x2 = a00*x1 + a01*y1 + a02*z1 + a03 // y2 = a10*x1 + a11*y1 + a12*z1 + a13 - // y2' = a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a03)*dydx2 + // y2' = a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2 for (Int_t i=0; i<12;i++) p[i]=0.; - p[3+0] = x1; - p[3+1] = y1; - p[3+2] = z1; - p[9+1] = 1.; - p[0+1] = y1*dydx2; - p[0+2] = z1*dydx2; - p[9+0] = dydx2; - value = y2; + p[2] += x1; // a10 + //p[] +=1; // a11 + p[3] += z1; // a12 + p[7] += 1; // a13 + p[0] += y1*dydx2; // a01 + p[1] += z1*dydx2; // a02 + p[6] += dydx2; // a03 + value = y2-y1; //-a11 fitter->AddPoint(p,value,sy); - + // // x2 = a00*x1 + a01*y1 + a02*z1 + a03 // z2 = a20*x1 + a21*y1 + a22*z1 + a23 - // z2' = a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a03)*dzdx2; + // z2' = a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 + a03)*dzdx2; for (Int_t i=0; i<12;i++) p[i]=0.; - p[6+0] = x1; - p[6+1] = y1; - p[6+2] = z1; - p[9+2] = 1.; - p[0+1] = y1*dzdx2; - p[0+2] = z1*dzdx2; - p[9+0] = dzdx2; - value = z2; + p[4] += x1; // a20 + p[5] += y1; // a21 + //p[] += z1; // a22 + p[8] += 1.; // a23 + p[0] += y1*dzdx2; // a01 + p[1] += z1*dzdx2; // a02 + p[6] += dzdx2; // a03 + value = z2-z1; //-a22 fitter->AddPoint(p,value,sz); // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1) // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0 for (Int_t i=0; i<12;i++) p[i]=0.; - p[3+0] = 1.; - p[3+1] = dydx1; - p[3+2] = dzdx1; - p[0+0] = -dydx2; - p[0+1] = -dydx1*dydx2; - p[0+2] = -dzdx1*dydx2; - value = 0.; + p[2] += 1.; // a10 + //p[] += dydx1; // a11 + p[3] += dzdx1; // a12 + //p[] += -dydx2; // a00 + p[0] += -dydx1*dydx2; // a01 + p[1] += -dzdx1*dydx2; // a02 + value = -dydx1+dydx2; // -a11 + a00 fitter->AddPoint(p,value,sdydx); // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1) // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0 for (Int_t i=0; i<12;i++) p[i]=0.; - p[6+0] = 1; - p[6+1] = dydx1; - p[6+2] = dzdx1; - p[0+0] = -dzdx2; - p[0+1] = -dydx1*dzdx2; - p[0+2] = -dzdx1*dzdx2; - value = 0.; + p[4] += 1; // a20 + p[5] += dydx1; // a21 + //p[] += dzdx1; // a22 + //p[] += -dzdx2; // a00 + p[0] += -dydx1*dzdx2; // a01 + p[1] += -dzdx1*dzdx2; // a02 + value = -dzdx1+dzdx2; // -a22 + a00 fitter->AddPoint(p,value,sdzdx); } -void AliTPCcalibAlign::Process6(Double_t *t1, - Double_t *t2, - TLinearFitter *fitter) { +void AliTPCcalibAlign::Process6(const Double_t *const t1, + const Double_t *const t2, + TLinearFitter *fitter) const +{ // x2 = 1 *x1 +-a01*y1 + 0 +a03 // y2 = a01*x1 + 1 *y1 + 0 +a13 // z2 = a20*x1 + a21*y1 + 1 *z1 +a23 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1) // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1) // - // 1 -a01 0 a03 x -p[0] x p[3] - // a10 1 0 a13 ==> p[0] x x p[4] - // a20 a21 1 a23 p[1] p[2] x p[5] - Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4]; - Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4]; + // a00 a01 a02 a03 1 -p[0] 0 p[3] + // a10 a11 a12 a13 ==> p[0] 1 0 p[4] + // a20 a21 a21 a23 p[1] p[2] 1 p[5] - // TODO: - Double_t sy = 0.1; - Double_t sz = 0.1; - Double_t sdydx = 0.001; - Double_t sdzdx = 0.001; + const Double_t &x1=t1[0], &y1=t1[1], &z1=t1[3], &dydx1=t1[2], &dzdx1=t1[4]; + const Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[3], &dydx2=t2[2], &dzdx2=t2[4]; + // + Double_t sy = TMath::Sqrt(t1[6]*t1[6]+t2[6]*t2[6]); + Double_t sdydx = TMath::Sqrt(t1[7]*t1[7]+t2[7]*t2[7]); + Double_t sz = TMath::Sqrt(t1[8]*t1[8]+t2[8]*t2[8]); + Double_t sdzdx = TMath::Sqrt(t1[9]*t1[9]+t2[9]*t2[9]); + + Double_t p[12]; Double_t value; - - // x2 = 1 *x1 +-a01*y1 + 0 +a03 - // y2 = a01*x1 + 1 *y1 + 0 +a13 + // x2 = a00*x1 + a01*y1 + a02*z1 + a03 + // y2 = a10*x1 + a11*y1 + a12*z1 + a13 // y2' = a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2 for (Int_t i=0; i<12;i++) p[i]=0.; - p[3+0] = x1; - p[3+1] = y1; - p[3+2] = z1; - p[9+1] = 1.; - p[0+1] = y1*dydx2; - p[0+2] = z1*dydx2; - p[9+0] = dydx2; - value = y2; + p[0] += x1; // a10 + //p[] +=1; // a11 + //p[] += z1; // a12 + p[4] += 1; // a13 + p[0] += -y1*dydx2; // a01 + //p[] += z1*dydx2; // a02 + p[3] += dydx2; // a03 + value = y2-y1; //-a11 fitter->AddPoint(p,value,sy); - - // x2 = 1 *x1 +-a01*y1 + 0 + a03 - // z2 = a20*x1 + a21*y1 + 1 *z1 + a23 - // z2' = a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 +a03)*dzdx2; - for (Int_t i=0; i<12;i++) p[i]=0.; - p[6+0] = x1; - p[6+1] = y1; - p[6+2] = z1; - p[9+2] = 1.; - p[0+1] = y1*dzdx2; - p[0+2] = z1*dzdx2; - p[9+0] = dzdx2; - value = z2; + // + // x2 = a00*x1 + a01*y1 + a02*z1 + a03 + // z2 = a20*x1 + a21*y1 + a22*z1 + a23 + // z2' = a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 + a03)*dzdx2; + for (Int_t i=0; i<12;i++) p[i]=0.; + p[1] += x1; // a20 + p[2] += y1; // a21 + //p[] += z1; // a22 + p[5] += 1.; // a23 + p[0] += -y1*dzdx2; // a01 + //p[] += z1*dzdx2; // a02 + p[3] += dzdx2; // a03 + value = z2-z1; //-a22 fitter->AddPoint(p,value,sz); // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1) // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0 for (Int_t i=0; i<12;i++) p[i]=0.; - p[3+0] = 1.; - p[3+1] = dydx1; - p[3+2] = dzdx1; - p[0+0] = -dydx2; - p[0+1] = -dydx1*dydx2; - p[0+2] = -dzdx1*dydx2; - value = 0.; + p[0] += 1.; // a10 + //p[] += dydx1; // a11 + //p[] += dzdx1; // a12 + //p[] += -dydx2; // a00 + //p[0] += dydx1*dydx2; // a01 FIXME- 0912 MI + //p[] += -dzdx1*dydx2; // a02 + value = -dydx1+dydx2; // -a11 + a00 fitter->AddPoint(p,value,sdydx); // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1) // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0 for (Int_t i=0; i<12;i++) p[i]=0.; - p[6+0] = 1; - p[6+1] = dydx1; - p[6+2] = dzdx1; - p[0+0] = -dzdx2; - p[0+1] = -dydx1*dzdx2; - p[0+2] = -dzdx1*dzdx2; - value = 0.; + p[1] += 1; // a20 + // p[2] += dydx1; // a21 FIXME- 0912 MI + //p[] += dzdx1; // a22 + //p[] += -dzdx2; // a00 + //p[0] += dydx1*dzdx2; // a01 FIXME- 0912 MI + //p[] += -dzdx1*dzdx2; // a02 + value = -dzdx1+dzdx2; // -a22 + a00 fitter->AddPoint(p,value,sdzdx); } -void AliTPCcalibAlign::EvalFitters() { +void AliTPCcalibAlign::EvalFitters(Int_t minPoints) { // // Analyze function // // Perform the fitting using linear fitters // - Int_t kMinPoints =50; TLinearFitter *f; TFile fff("alignDebug.root","recreate"); for (Int_t s1=0;s1<72;++s1) for (Int_t s2=0;s2<72;++s2){ - if ((f=GetFitter12(s1,s2))&&fPoints[GetIndex(s1,s2)]>kMinPoints) { + if ((f=GetFitter12(s1,s2))&&fPoints[GetIndex(s1,s2)]>minPoints) { // cerr<Eval()!=0) { cerr<<"Evaluation failed for "<Write(Form("f12_%d_%d",s1,s2)); } } - if ((f=GetFitter9(s1,s2))&&fPoints[GetIndex(s1,s2)]>kMinPoints) { + if ((f=GetFitter9(s1,s2))&&fPoints[GetIndex(s1,s2)]>minPoints) { // cerr<Eval()!=0) { cerr<<"Evaluation failed for "<Write(Form("f9_%d_%d",s1,s2)); } } - if ((f=GetFitter6(s1,s2))&&fPoints[GetIndex(s1,s2)]>kMinPoints) { + if ((f=GetFitter6(s1,s2))&&fPoints[GetIndex(s1,s2)]>minPoints) { // cerr<Eval()!=0) { cerr<<"Evaluation failed for "<Write("align"); - /* - - fitter->Eval(); - fitter->Eval(); - chi212 = align->GetChisquare()/(4.*entries); - - TMatrixD mat(13,13); - TVectorD par(13); - align->GetParameters(par); - align->GetCovarianceMatrix(mat); - - // - // - for (Int_t i=0; i<12;i++){ - palign12(i)= par(i+1); - for (Int_t j=0; j<12;j++){ - pcovar12(i,j) = mat(i+1,j+1); - pcovar12(i,j) *= chi212; - } - } - // - for (Int_t i=0; i<12;i++){ - psigma12(i) = TMath::Sqrt(pcovar12(i,i)); - palignR12(i) = palign12(i)/TMath::Sqrt(pcovar12(i,i)); - for (Int_t j=0; j<12;j++){ - pcovarN12(i,j) = pcovar12(i,j)/TMath::Sqrt(pcovar12(i,i)*pcovar12(j,j)); + TMatrixD mat(4,4); + for (Int_t s1=0;s1<72;++s1) + for (Int_t s2=0;s2<72;++s2){ + if (GetTransformation12(s1,s2,mat)){ + fMatrixArray12.AddAt(mat.Clone(), GetIndex(s1,s2)); + } + if (GetTransformation9(s1,s2,mat)){ + fMatrixArray9.AddAt(mat.Clone(), GetIndex(s1,s2)); + } + if (GetTransformation6(s1,s2,mat)){ + fMatrixArray6.AddAt(mat.Clone(), GetIndex(s1,s2)); + } } - } - */ + //this->Write("align"); + } TLinearFitter* AliTPCcalibAlign::GetOrMakeFitter12(Int_t s1,Int_t s2) { @@ -515,7 +1555,7 @@ TLinearFitter* AliTPCcalibAlign::GetOrMakeFitter12(Int_t s1,Int_t s2) { fitter->StoreData(kFALSE); fFitterArray12.AddAt(fitter,GetIndex(s1,s2)); counter12++; - if (GetDebugLevel()>0) cerr<<"Creating fitter12 "<0) cerr<<"Creating fitter12 "<StoreData(kFALSE); fFitterArray9.AddAt(fitter,GetIndex(s1,s2)); counter9++; - if (GetDebugLevel()>0) cerr<<"Creating fitter12 "<0) cerr<<"Creating fitter12 "<StoreData(kFALSE); fFitterArray6.AddAt(fitter,GetIndex(s1,s2)); counter6++; - if (GetDebugLevel()>0) cerr<<"Creating fitter6 "<0) cerr<<"Creating fitter6 "<GetParameters(p); a.ResizeTo(4,4); - a[0][0]=p[0]; a[0][1]=p[1]; a[0][2]=p[2]; a[0][3]=p[9]; - a[1][0]=p[3]; a[1][1]=p[4]; a[1][2]=p[5]; a[1][3]=p[10]; - a[2][0]=p[6]; a[2][1]=p[7]; a[2][2]=p[8]; a[2][3]=p[11]; + a[0][0]=1; a[0][1]=p[0]; a[0][2]=p[1]; a[0][3]=p[6]; + a[1][0]=p[2]; a[1][1]=1; a[1][2]=p[3]; a[1][3]=p[7]; + a[2][0]=p[4]; a[2][1]=p[5]; a[2][2]=1; a[2][3]=p[8]; a[3][0]=0.; a[3][1]=0.; a[3][2]=0.; a[3][3]=1.; return true; } @@ -608,26 +1648,194 @@ Bool_t AliTPCcalibAlign::GetTransformation6(Int_t s1,Int_t s2,TMatrixD &a) { TVectorD p(6); GetFitter6(s1,s2)->GetParameters(p); a.ResizeTo(4,4); - a[0][0]=p[0]; a[0][1]=p[1]; a[0][2]=p[2]; a[0][3]=p[9]; - a[1][0]=p[3]; a[1][1]=p[4]; a[1][2]=p[5]; a[1][3]=p[10]; - a[2][0]=p[6]; a[2][1]=p[7]; a[2][2]=p[8]; a[2][3]=p[11]; + a[0][0]=1; a[0][1]=-p[0];a[0][2]=0; a[0][3]=p[3]; + a[1][0]=p[0]; a[1][1]=1; a[1][2]=0; a[1][3]=p[4]; + a[2][0]=p[1]; a[2][1]=p[2]; a[2][2]=1; a[2][3]=p[5]; a[3][0]=0.; a[3][1]=0.; a[3][2]=0.; a[3][3]=1.; return true; } } -void AliTPCcalibAlign::FillHisto(const AliExternalTrackParam &tp1, - const AliExternalTrackParam &tp2, - Int_t s1,Int_t s2) { +void AliTPCcalibAlign::MakeResidualHistos(){ + // + // Make cluster residual histograms + // + Double_t xminTrack[9], xmaxTrack[9]; + Int_t binsTrack[9]; + TString axisName[9],axisTitle[9]; + // + // 0 - delta of interest + // 1 - global phi in sector number as float + // 2 - local x + // 3 - local ky + // 4 - local kz + // + axisName[0]="delta"; axisTitle[0]="#Delta (cm)"; + 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; + // + axisName[2]="R"; axisTitle[2]="r (cm)"; + binsTrack[2]=53; xminTrack[2]=85.; xmaxTrack[2]=245.; + // + // + axisName[3]="kZ"; axisTitle[3]="dz/dx"; + binsTrack[3]=36; xminTrack[3]=-1.8; xmaxTrack[3]=1.8; + // + fClusterDelta[0] = new THnSparseS("#Delta_{Y} (cm)","#Delta_{Y} (cm)", 4, binsTrack,xminTrack, xmaxTrack); + fClusterDelta[1] = new THnSparseS("#Delta_{Z} (cm)","#Delta_{Z} (cm)", 4, binsTrack,xminTrack, xmaxTrack); + // + // + // + for (Int_t ivar=0;ivar<2;ivar++){ + for (Int_t ivar2=0;ivar2<4;ivar2++){ + fClusterDelta[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data()); + fClusterDelta[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data()); + } + } + +} + + +void AliTPCcalibAlign::MakeResidualHistosTracklet(){ + // + // Make tracklet residual histograms + // + Double_t xminTrack[9], xmaxTrack[9]; + Int_t binsTrack[9]; + TString axisName[9],axisTitle[9]; + // + // 0 - delta of interest + // 1 - global phi in sector number as float + // 2 - local x + // 3 - local ky + // 4 - local kz + // 5 - sector 1 + // 5 - sector 0 + + axisName[0]="delta"; axisTitle[0]="#Delta (cm)"; + binsTrack[0]=60; xminTrack[0]=-0.6; xmaxTrack[0]=0.6; + // + axisName[1]="phi"; axisTitle[1]="#phi"; + binsTrack[1]=180; xminTrack[1]=-TMath::Pi(); xmaxTrack[1]=TMath::Pi(); + // + axisName[2]="localX"; axisTitle[2]="x (cm)"; + binsTrack[2]=10; xminTrack[2]=120.; xmaxTrack[2]=200.; + // + axisName[3]="kY"; axisTitle[3]="dy/dx"; + binsTrack[3]=10; xminTrack[3]=-0.5; xmaxTrack[3]=0.5; + // + axisName[4]="kZ"; axisTitle[4]="dz/dx"; + binsTrack[4]=22; xminTrack[4]=-1.1; xmaxTrack[4]=1.1; + // + axisName[5]="is1"; axisTitle[5]="is1"; + binsTrack[5]=72; xminTrack[5]=0; xmaxTrack[5]=72; + // + axisName[6]="is0"; axisTitle[6]="is0"; + binsTrack[6]=72; xminTrack[6]=0; xmaxTrack[6]=72; + + // + 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.5; xmaxTrack[0]=0.5; + fTrackletDelta[1] = new THnSparseF("#Delta_{Z} (cm)","#Delta_{Z} (cm)", 7, 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); + // + // + // + for (Int_t ivar=0;ivar<4;ivar++){ + for (Int_t ivar2=0;ivar2<7;ivar2++){ + fTrackletDelta[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data()); + fTrackletDelta[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data()); + } + } + +} + + + +void AliTPCcalibAlign::FillHisto(const Double_t *t1, + const Double_t *t2, + Int_t s1,Int_t s2) { // // Fill residual histograms + // Track2-Track1 + // Innner-Outer + // Left right - x-y + // A-C side + if (1) { + Double_t dy = t2[1]-t1[1]; + Double_t dphi = t2[2]-t1[2]; + Double_t dz = t2[3]-t1[3]; + Double_t dtheta = t2[4]-t1[4]; + Double_t zmean = (t2[3]+t1[3])*0.5; + // + GetHisto(kPhi,s1,s2,kTRUE)->Fill(dphi); + GetHisto(kTheta,s1,s2,kTRUE)->Fill(dtheta); + GetHisto(kY,s1,s2,kTRUE)->Fill(dy); + GetHisto(kZ,s1,s2,kTRUE)->Fill(dz); + // + GetHisto(kPhiZ,s1,s2,kTRUE)->Fill(zmean,dphi); + GetHisto(kThetaZ,s1,s2,kTRUE)->Fill(zmean,dtheta); + GetHisto(kYz,s1,s2,kTRUE)->Fill(zmean,dy); + GetHisto(kZz,s1,s2,kTRUE)->Fill(zmean,dz); + // + GetHisto(kYPhi,s1,s2,kTRUE)->Fill(t2[2],dy); + GetHisto(kZTheta,s1,s2,kTRUE)->Fill(t2[4],dz); + } +} + + +void AliTPCcalibAlign::FillHisto(AliExternalTrackParam *tp1, + AliExternalTrackParam *tp2, + Int_t s1,Int_t s2) { // - if (s2-s1==36) {//only inner-outer - GetHisto(kPhi,s1,s2,kTRUE)->Fill(TMath::ASin(tp1.GetSnp())-TMath::ASin(tp2.GetSnp())); - GetHisto(kTheta,s1,s2,kTRUE)->Fill(TMath::ATan(tp1.GetTgl())-TMath::ATan(tp2.GetTgl())); - GetHisto(kY,s1,s2,kTRUE)->Fill(tp1.GetY()-tp2.GetY()); - GetHisto(kZ,s1,s2,kTRUE)->Fill(tp1.GetZ()-tp2.GetZ()); - } + // Fill residual histograms + // Track2-Track1 + const Double_t kEpsilon=0.001; + Double_t x[8]={0,0,0,0,0,0,0,0}; + AliExternalTrackParam p1(*tp1); + AliExternalTrackParam p2(*tp2); + if (s1%18==s2%18) { + // inner outer - match at the IROC-OROC boundary + if (!p1.PropagateTo(fXIO, AliTrackerBase::GetBz())) return; + } + if (!p2.Rotate(p1.GetAlpha())) return; + if (!p2.PropagateTo(p1.GetX(),AliTrackerBase::GetBz())) return; + if (TMath::Abs(p1.GetX()-p2.GetX())>kEpsilon) return; + Double_t xyz[3]; + p1.GetXYZ(xyz); + x[1]=TMath::ATan2(xyz[1],xyz[0]); + x[2]=p1.GetX(); + x[3]=0.5*(p1.GetSnp()+p2.GetSnp()); // mean snp + x[4]=0.5*(p1.GetTgl()+p2.GetTgl()); // mean tgl + x[5]=s2; + x[6]=s1; + + x[0]=p2.GetY()-p1.GetY(); + fTrackletDelta[0]->Fill(x); + x[0]=p2.GetZ()-p1.GetZ(); + fTrackletDelta[1]->Fill(x); + x[0]=p2.GetSnp()-p1.GetSnp(); + fTrackletDelta[2]->Fill(x); + x[0]=p2.GetTgl()-p1.GetTgl(); + fTrackletDelta[3]->Fill(x); + TTreeSRedirector *cstream = GetDebugStreamer(); + if (cstream){ + (*cstream)<<"trackletMatch"<< + "tp1.="<At(GetIndex(s1,s2)); if (histo) return histo; @@ -661,25 +1881,1939 @@ TH1 * AliTPCcalibAlign::GetHisto(HistoType type, Int_t s1, Int_t s2, Bool_t forc case kY: name<<"hist_y_"<SetDirectory(0); histoArray->AddAt(histo,GetIndex(s1,s2)); return histo; } + +TGraphErrors * AliTPCcalibAlign::MakeGraph(Int_t sec0, Int_t sec1, Int_t dsec, + Int_t i0, Int_t i1, FitType type) +{ + // + // + // + TMatrixD mat; + //TObjArray *fitArray=0; + Double_t xsec[1000]; + Double_t ysec[1000]; + Int_t npoints=0; + for (Int_t isec = sec0; isec<=sec1; isec++){ + Int_t isec2 = (isec+dsec)%72; + switch (type) { + case k6: + GetTransformation6(isec,isec2,mat);break; + case k9: + GetTransformation9(isec,isec2,mat);break; + case k12: + GetTransformation12(isec,isec2,mat);break; + } + xsec[npoints]=isec; + ysec[npoints]=mat(i0,i1); + ++npoints; + } + TGraphErrors *gr = new TGraphErrors(npoints,xsec,ysec,0,0); + Char_t name[1000]; + sprintf(name,"Mat[%d,%d] Type=%d",i0,i1,type); + gr->SetName(name); + return gr; +} + +void AliTPCcalibAlign::MakeTree(const char *fname, Int_t minPoints){ + // + // make tree with alignment cosntant - + // For QA visualization + // + /* + TFile fcalib("CalibObjects.root"); + TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib"); + AliTPCcalibAlign * align = ( AliTPCcalibAlign *)array->FindObject("alignTPC"); + align->EvalFitters(); + align->MakeTree("alignTree.root"); + TFile falignTree("alignTree.root"); + TTree * treeAlign = (TTree*)falignTree.Get("Align"); + */ + TTreeSRedirector cstream(fname); + for (Int_t s1=0;s1<72;++s1) + for (Int_t s2=0;s2<72;++s2){ + TMatrixD m6; + TMatrixD m6FX; + TMatrixD m9; + TMatrixD m12; + TVectorD param6Diff; // align parameters diff + TVectorD param6s1(6); // align parameters sector1 + TVectorD param6s2(6); // align parameters sector2 + + // + // + if (fSectorParamA){ + TMatrixD * kpar = fSectorParamA; + TMatrixD * kcov = fSectorCovarA; + if (s1%36>=18){ + kpar = fSectorParamC; + kcov = fSectorCovarC; + } + for (Int_t ipar=0;ipar<6;ipar++){ + Int_t isec1 = s1%18; + Int_t isec2 = s2%18; + if (s1>35) isec1+=18; + if (s2>35) isec2+=18; + param6s1(ipar)=(*kpar)(6*isec1+ipar,0); + param6s2(ipar)=(*kpar)(6*isec2+ipar,0); + } + } + + Double_t dy=0, dz=0, dphi=0,dtheta=0; + Double_t sy=0, sz=0, sphi=0,stheta=0; + Double_t ny=0, nz=0, nphi=0,ntheta=0; + Double_t chi2v12=0, chi2v9=0, chi2v6=0; + // Int_t npoints=0; + // TLinearFitter * fitter = 0; + if (fPoints[GetIndex(s1,s2)]>minPoints){ + // + // + // +// fitter = GetFitter12(s1,s2); +// npoints = fitter->GetNpoints(); +// chi2v12 = TMath::Sqrt(fitter->GetChisquare()/npoints); + +// // +// fitter = GetFitter9(s1,s2); +// npoints = fitter->GetNpoints(); +// chi2v9 = TMath::Sqrt(fitter->GetChisquare()/npoints); +// // +// fitter = GetFitter6(s1,s2); +// npoints = fitter->GetNpoints(); +// chi2v6 = TMath::Sqrt(fitter->GetChisquare()/npoints); +// fitter->GetParameters(param6Diff); +// // +// GetTransformation6(s1,s2,m6); +// GetTransformation9(s1,s2,m9); +// GetTransformation12(s1,s2,m12); +// // +// fitter = GetFitter6(s1,s2); +// //fitter->FixParameter(3,0); +// //fitter->Eval(); +// GetTransformation6(s1,s2,m6FX); + // + TH1 * his=0; + his = GetHisto(kY,s1,s2); + if (his) { dy = his->GetMean(); sy = his->GetRMS(); ny = his->GetEntries();} + his = GetHisto(kZ,s1,s2); + if (his) { dz = his->GetMean(); sz = his->GetRMS(); nz = his->GetEntries();} + his = GetHisto(kPhi,s1,s2); + if (his) { dphi = his->GetMean(); sphi = his->GetRMS(); nphi = his->GetEntries();} + his = GetHisto(kTheta,s1,s2); + if (his) { dtheta = his->GetMean(); stheta = his->GetRMS(); ntheta = his->GetEntries();} + // + + } + AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField(); + if (!magF) AliError("Magneticd field - not initialized"); + Double_t bz = magF->SolenoidField()/10.; //field in T + + cstream<<"Align"<< + "run="<0) Info("AliTPCcalibAlign","Merge"); + if (!list) + return 0; + if (list->IsEmpty()) + return 1; + + TIterator* iter = list->MakeIterator(); + TObject* obj = 0; + iter->Reset(); + Int_t count=0; + // + TString str1(GetName()); + while((obj = iter->Next()) != 0) + { + AliTPCcalibAlign* entry = dynamic_cast(obj); + if (entry == 0) continue; + if (str1.CompareTo(entry->GetName())!=0) continue; + Add(entry); + count++; + } + return count; +} + + +void AliTPCcalibAlign::Add(AliTPCcalibAlign * align){ + // + // Add entry - used for merging of compoents + // + for (Int_t i=0; i<72;i++){ + for (Int_t j=0; j<72;j++){ + if (align->fPoints[GetIndex(i,j)]<1) continue; + fPoints[GetIndex(i,j)]+=align->fPoints[GetIndex(i,j)]; + // + // + // + for (Int_t itype=0; itype<10; itype++){ + TH1 * his0=0, *his1=0; + his0 = GetHisto((HistoType)itype,i,j); + his1 = align->GetHisto((HistoType)itype,i,j); + if (his1){ + if (his0) his0->Add(his1); + else { + his0 = GetHisto((HistoType)itype,i,j,kTRUE); + his0->Add(his1); + } + } + } + } + } + TLinearFitter *f0=0; + TLinearFitter *f1=0; + for (Int_t i=0; i<72;i++){ + for (Int_t j=0; j<72;j++){ + if (align->fPoints[GetIndex(i,j)]<1) continue; + // + // + // fitter12 + f0 = GetFitter12(i,j); + f1 = align->GetFitter12(i,j); + if (f1){ + if (f0) f0->Add(f1); + else { + fFitterArray12.AddAt(f1->Clone(),GetIndex(i,j)); + } + } + // + // fitter9 + f0 = GetFitter9(i,j); + f1 = align->GetFitter9(i,j); + if (f1){ + if (f0) f0->Add(f1); + else { + fFitterArray9.AddAt(f1->Clone(),GetIndex(i,j)); + } + } + f0 = GetFitter6(i,j); + f1 = align->GetFitter6(i,j); + if (f1){ + if (f0) f0->Add(f1); + else { + fFitterArray6.AddAt(f1->Clone(),GetIndex(i,j)); + } + } + } + } + // + // Add Kalman filter + // + for (Int_t i=0;i<36;i++){ + TMatrixD *par0 = (TMatrixD*)fArraySectorIntParam.At(i); + if (!par0){ + MakeSectorKalman(); + par0 = (TMatrixD*)fArraySectorIntParam.At(i); + } + TMatrixD *par1 = (TMatrixD*)align->fArraySectorIntParam.At(i); + if (!par1) continue; + // + TMatrixD *cov0 = (TMatrixD*)fArraySectorIntCovar.At(i); + TMatrixD *cov1 = (TMatrixD*)align->fArraySectorIntCovar.At(i); + UpdateSectorKalman(*par0,*cov0,*par1,*cov1); + } + if (!fSectorParamA){ + MakeKalman(); + } + if (align->fSectorParamA){ + UpdateKalman(*fSectorParamA,*fSectorCovarA,*align->fSectorParamA,*align->fSectorCovarA); + UpdateKalman(*fSectorParamC,*fSectorCovarC,*align->fSectorParamC,*align->fSectorCovarC); + } + if (!fClusterDelta[0]) MakeResidualHistos(); + + for (Int_t i=0; i<2; i++){ + if (align->fClusterDelta[i]) fClusterDelta[i]->Add(align->fClusterDelta[i]); + } + + 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]); + } + +} + +Double_t AliTPCcalibAlign::Correct(Int_t type, Int_t value, Int_t s1, Int_t s2, Double_t x1, Double_t y1, Double_t z1, Double_t dydx1,Double_t dzdx1){ + // + // GetTransformed value + // + // + // x2 = a00*x1 + a01*y1 + a02*z1 + a03 + // y2 = a10*x1 + a11*y1 + a12*z1 + a13 + // z2 = a20*x1 + a21*y1 + a22*z1 + a23 + // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1) + // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1) + + + const TMatrixD * mat = GetTransformation(s1,s2,type); + if (!mat) { + if (value==0) return x1; + if (value==1) return y1; + if (value==2) return z1; + if (value==3) return dydx1; + if (value==4) return dzdx1; + // + if (value==5) return dydx1; + if (value==6) return dzdx1; + return 0; + } + Double_t valT=0; + + if (value==0){ + valT = (*mat)(0,0)*x1+(*mat)(0,1)*y1+(*mat)(0,2)*z1+(*mat)(0,3); + } + + if (value==1){ + valT = (*mat)(1,0)*x1+(*mat)(1,1)*y1+(*mat)(1,2)*z1+(*mat)(1,3); + } + if (value==2){ + valT = (*mat)(2,0)*x1+(*mat)(2,1)*y1+(*mat)(2,2)*z1+(*mat)(2,3); + } + if (value==3){ + // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1) + valT = (*mat)(1,0) +(*mat)(1,1)*dydx1 +(*mat)(1,2)*dzdx1; + valT/= ((*mat)(0,0) +(*mat)(0,1)*dydx1 +(*mat)(0,2)*dzdx1); + } + + if (value==4){ + // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1) + valT = (*mat)(2,0) +(*mat)(2,1)*dydx1 +(*mat)(2,2)*dzdx1; + valT/= ((*mat)(0,0) +(*mat)(0,1)*dydx1 +(*mat)(0,2)*dzdx1); + } + // + if (value==5){ + // onlys shift in angle + // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1) + valT = (*mat)(1,0) +(*mat)(1,1)*dydx1; + } + + if (value==6){ + // only shift in angle + // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1) + valT = (*mat)(2,0) +(*mat)(2,1)*dydx1; + } + // + return valT; +} + + +void AliTPCcalibAlign::Constrain1Pt(AliExternalTrackParam &track1, const AliExternalTrackParam &track2, Bool_t noField){ + // + // Update track parameters t1 + // + TMatrixD vecXk(5,1); // X vector + TMatrixD covXk(5,5); // X covariance + TMatrixD matHk(1,5); // vector to mesurement + TMatrixD measR(1,1); // measurement error + //TMatrixD matQk(5,5); // prediction noise vector + TMatrixD vecZk(1,1); // measurement + // + TMatrixD vecYk(1,1); // Innovation or measurement residual + TMatrixD matHkT(5,1); + TMatrixD matSk(1,1); // Innovation (or residual) covariance + TMatrixD matKk(5,1); // Optimal Kalman gain + TMatrixD mat1(5,5); // update covariance matrix + TMatrixD covXk2(5,5); // + TMatrixD covOut(5,5); + // + Double_t *param1=(Double_t*) track1.GetParameter(); + Double_t *covar1=(Double_t*) track1.GetCovariance(); + + // + // copy data to the matrix + for (Int_t ipar=0; ipar<5; ipar++){ + vecXk(ipar,0) = param1[ipar]; + for (Int_t jpar=0; jpar<5; jpar++){ + covXk(ipar,jpar) = covar1[track1.GetIndex(ipar, jpar)]; + } + } + // + // + // + vecZk(0,0) = track2.GetParameter()[4]; // 1/pt measurement from track 2 + measR(0,0) = track2.GetCovariance()[14]; // 1/pt measurement error + if (noField) { + measR(0,0)*=0.000000001; + vecZk(0,0)=0.; + } + // + matHk(0,0)=0; matHk(0,1)= 0; matHk(0,2)= 0; + matHk(0,3)= 0; matHk(0,4)= 1; // vector to measurement + // + // + // + vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual + matHkT=matHk.T(); matHk.T(); + matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance + matSk.Invert(); + matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain + vecXk += matKk*vecYk; // updated vector + mat1(0,0)=1; mat1(1,1)=1; mat1(2,2)=1; mat1(3,3)=1; mat1(4,4)=1; + covXk2 = (mat1-(matKk*matHk)); + covOut = covXk2*covXk; + // + // + // + // copy from matrix to parameters + if (0) { + covOut.Print(); + vecXk.Print(); + covXk.Print(); + track1.Print(); + track2.Print(); + } + + for (Int_t ipar=0; ipar<5; ipar++){ + param1[ipar]= vecXk(ipar,0) ; + for (Int_t jpar=0; jpar<5; jpar++){ + covar1[track1.GetIndex(ipar, jpar)]=covOut(ipar,jpar); + } + } + +} + +void AliTPCcalibAlign::GlobalAlign6(Int_t minPoints, Float_t sysError, Int_t niter){ + // + // Global Align -combine the partial alignment of pair of sectors + // minPoints - minimal number of points - don't use sector alignment wit smaller number + // sysError - error added to the alignemnt error + // + AliTPCcalibAlign * align = this; + TMatrixD * arrayAlign[72]; + TMatrixD * arrayAlignDiff[72]; + // + for (Int_t i=0;i<72; i++) { + TMatrixD * mat = new TMatrixD(4,4); + mat->UnitMatrix(); + arrayAlign[i]=mat; + arrayAlignDiff[i]=(TMatrixD*)(mat->Clone()); + } + + TTreeSRedirector *cstream = new TTreeSRedirector("galign6.root"); + for (Int_t iter=0; iterGetTransformation(is0,is1,0); + if (mat){ + npoints = align->GetFitter6(is0,is1)->GetNpoints(); + if (npoints>minPoints){ + align->GetFitter6(is0,is1)->GetCovarianceMatrix(covar); + align->GetFitter6(is0,is1)->GetErrors(errors); + } + } + else{ + invers=kTRUE; + mat = align->GetTransformation(is1,is0,0); + if (mat) { + npoints = align->GetFitter6(is1,is0)->GetNpoints(); + if (npoints>minPoints){ + align->GetFitter6(is1,is0)->GetCovarianceMatrix(covar); + align->GetFitter6(is1,is0)->GetErrors(errors); + } + } + } + if (!mat) continue; + if (npointsis0/36) weight*=2./3.; //IROC-OROC + if (is1/360){ + matDiff*=1/sumw; + matDiff(0,0)=0; + matDiff(1,1)=0; + matDiff(1,1)=0; + matDiff(1,1)=0; + (*arrayAlignDiff[is0]) = matDiff; + } + } + for (Int_t is0=0;is0<72; is0++) { + if (is0<36) (*arrayAlign[is0]) += 0.4*(*arrayAlignDiff[is0]); + if (is0>=36) (*arrayAlign[is0]) += 0.2*(*arrayAlignDiff[is0]); + // + (*cstream)<<"GAlign"<< + "iter="<GetClusterPointer(irow); + if (!c) continue; + // + if (c->GetDetector()%36!=(isec%36)) continue; + if (!both && c->GetDetector()!=isec) continue; + + if (c->GetRow()GetX()-sa*c->GetY(); + Double_t lyR = +sa*c->GetX()+ca*c->GetY(); + Double_t lzR = c->GetZ(); + + Double_t dx = lxR -xRef; // distance to reference X + Double_t x[2]={dx, dx*dx}; + + Double_t yfitR = pyf[0]+pyf[1]*dx; // fit value Y in ref frame + Double_t zfitR = pzf[0]+pzf[1]*dx; // fit value Z in ref frame + // + Double_t yfit = -sa*lxR + ca*yfitR; // fit value Y in local frame + // + if (iter==0 &&c->GetType()<0) continue; + if (iter>0){ + if (TMath::Abs(lyR-yfitR)>kMaxDist*erry) continue; + if (TMath::Abs(lzR-zfitR)>kMaxDist*errz) continue; + Double_t dedge = c->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(yfit); + if (isec<36 && dedge35 && dedgeRPhiCOGCorrection(isec,c->GetRow(), c->GetPad(), + c->GetY(),yfit, c->GetZ(), pyf[1], c->GetMax(),2.5); + Double_t corrclY = + corr->RPhiCOGCorrection(isec,c->GetRow(), c->GetPad(), + c->GetY(),c->GetY(), c->GetZ(), pyf[1], c->GetMax(),2.5); + if (TMath::Abs((corrtrY+corrclY)*0.5)>kMaxCorrY) continue; + if (TMath::Abs(corrtrY)>kMaxCorrY) continue; + } + fyf.AddPoint(x,lyR,erry); + fzf.AddPoint(x,lzR,errz); + } + nf = fyf.GetNpoints(); + if (nfGetImpactParameters(vertexXY,vertexZ); + if (TMath::Abs(vertexXY)>kVertexCut) return; + if (TMath::Abs(vertexZ)>kVertexCut) return; + if (TMath::Abs(seed->Pt())GetNumberOfClusters()GetSnp())>kSnpCut) return; + if (!fClusterDelta[0]) MakeResidualHistos(); + + Int_t detector=-1; + // + // + AliExternalTrackParam trackIn = *(fCurrentTrack->GetInnerParam()); + AliExternalTrackParam trackOut = *(fCurrentFriendTrack->GetTPCOut()); + static Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass(); + // + Int_t ncl=0; + for (Int_t irow=0; irow<160; irow++){ + AliTPCclusterMI *cl=seed->GetClusterPointer(irow); + if (!cl) continue; + if (cl->GetX()<80) continue; + if (detector<0) detector=cl->GetDetector()%36; + if (detector!=cl->GetDetector()%36) return; // cluster from different sectors + // skip such tracks + ncl++; + } + if (nclGetClusterPointer(irow); + if (!cl) continue; + if (cl->GetX()<80) continue; + if (detector<0) detector=cl->GetDetector()%36; + Int_t sector = cl->GetDetector(); + Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackOut.GetAlpha(); + 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 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 + cov[0]*=cov[0]; + cov[2]*=cov[2]; + if (!AliTracker::PropagateTrackToBxByBz(&trackOut, r[0],mass,1.,kFALSE)) continue; + + if (TMath::Abs(dedge)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); + } + // + // Refit in - store residual maps + // + for (Int_t irow=159; irow>=0; irow--){ + AliTPCclusterMI *cl=seed->GetClusterPointer(irow); + if (!cl) continue; + if (cl->GetX()<80) continue; + if (detector<0) detector=cl->GetDetector()%36; + Int_t sector = cl->GetDetector(); + Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackIn.GetAlpha(); + 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 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 +- + cov[0]*=cov[0]; + cov[2]*=cov[2]; + if (!AliTracker::PropagateTrackToBxByBz(&trackIn, r[0],mass,1.,kFALSE)) continue; + if (TMath::Abs(dedge)GetY()-trackIn.GetY())GetDetector()%36!=detector) continue; + // + // fill residual histogram + // + Double_t resVector[5]; + trackIn.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(); + fClusterDelta[0]->Fill(resVector); + resVector[0]= cl->GetZ()-trackIn.GetZ(); + fClusterDelta[1]->Fill(resVector); + } + + +} + + +void AliTPCcalibAlign::UpdateAlignSector(const AliTPCseed * track,Int_t isec){ + // + // Update Kalman filter of Alignment - only setup without filed + // IROC - OROC quadrants + // + if (TMath::Abs(AliTracker::GetBz())>0.5) return; + if (!fClusterDelta[0]) MakeResidualHistos(); + const Int_t kMinClusterF=40; + const Int_t kMinClusterFit=10; + const Int_t kMinClusterQ=10; + // + const Int_t kdrow1Fit =5; // rows to skip from fit at the end + const Int_t kdrow0Fit =10; // rows to skip from fit at beginning + const Float_t kedgey=3.0; + const Float_t kMaxDist=1; + const Float_t kMaxCorrY=0.05; + const Float_t kPRFWidth = 0.6; //cut 2 sigma of PRF + isec = isec%36; // use the hardware numbering + // + // + AliTPCPointCorrection * corr = AliTPCPointCorrection::Instance(); + // + // full track fit parameters + // + static TLinearFitter fyf(2,"pol1"); // make it static - too much time for comiling of formula + static TLinearFitter fzf(2,"pol1"); // calgrind recomendation + TVectorD pyf(2), peyf(2),pzf(2), pezf(2); + TVectorD pyfc(2),pzfc(2); + Int_t nf=0; + // + // make full fit as reference + // + for (Int_t iter=0; iter<2; iter++){ + fyf.ClearPoints(); + fzf.ClearPoints(); + for (Int_t irow=kdrow0Fit;irow<159-kdrow1Fit;irow++) { + AliTPCclusterMI *c=track->GetClusterPointer(irow); + if (!c) continue; + if ((c->GetDetector()%36)!=isec) continue; + if (c->GetRow()GetX()-fXmiddle; + Double_t x[2]={dx, dx*dx}; + if (iter==0 &&c->GetType()<0) continue; + if (iter==1){ + Double_t yfit = pyf[0]+pyf[1]*dx; + Double_t zfit = pzf[0]+pzf[1]*dx; + Double_t dedge = c->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(yfit); + if (TMath::Abs(c->GetY()-yfit)>kMaxDist) continue; + if (TMath::Abs(c->GetZ()-zfit)>kMaxDist) continue; + if (dedgeRPhiCOGCorrection(c->GetDetector(),c->GetRow(), c->GetPad(), + c->GetY(),yfit, c->GetZ(), pyf[1], c->GetMax(),2.5); + if (TMath::Abs(corrtrY)>kMaxCorrY) continue; + } + 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); + } + nf = fyf.GetNpoints(); + if (fyf.GetNpoints()GetTOFsignal(); // tof signal + TVectorF tpcPosG(3); // global position of track at the middle of fXmiddle + Double_t lphi = TMath::ATan2(pyf[0],fXmiddle); // expected local phi angle - if vertex at 0 + Double_t gphi = 2.*TMath::Pi()*(isec%18+0.5)/18.+lphi; // expected global phi if vertex at 0 + Double_t ky = pyf[0]/fXmiddle; + Double_t kyE =0, kzE=0; // ky and kz expected + Double_t alpha =2.*TMath::Pi()*(isec%18+0.5)/18.; + Double_t scos=TMath::Cos(alpha); + Double_t ssin=TMath::Sin(alpha); + const AliESDVertex* vertex = fCurrentEvent->GetPrimaryVertexTracks(); + vertex->GetXYZ(vPosG.GetMatrixArray()); + fCurrentTrack->GetImpactParameters(vImpact[0],vImpact[1]); // track impact parameters + // + tpcPosG[0]= scos*fXmiddle-ssin*pyf[0]; + tpcPosG[1]=+ssin*fXmiddle+scos*pyf[0]; + tpcPosG[2]=pzf[0]; + vPosL[0]= scos*vPosG[0]+ssin*vPosG[1]; + vPosL[1]=-ssin*vPosG[0]+scos*vPosG[1]; + vPosL[2]= vPosG[2]; + kyE = (pyf[0]-vPosL[1])/(fXmiddle-vPosL[0]); + kzE = (pzf[0]-vPosL[2])/(fXmiddle-vPosL[0]); + // + // get constrained parameters + // + Double_t xvertex=vPosL[0]-fXmiddle; + fyf.AddPoint(&xvertex,vPosL[1], 0.00001); + fzf.AddPoint(&xvertex,vPosL[2], 2.); + fyf.Eval(); + fyf.GetParameters(pyfc); + fzf.Eval(); + fzf.GetParameters(pzfc); + // + // + // Make Fitters and params for 5 fitters + // 1-4 OROC quadrants + // 0 IROC + // + static TLinearFitter *fittersY[5]={0,0,0,0,0}; // calgrind recomendation - fater to clear points + static TLinearFitter *fittersZ[5]={0,0,0,0,0}; // than create the fitter + if (fittersY[0]==0){ + for (Int_t i=0;i<5;i++) { + fittersY[i] = new TLinearFitter(2,"pol1"); + fittersZ[i] = new TLinearFitter(2,"pol1"); + } + } + // + Int_t npoints[5]; + TVectorD paramsY[5]; + TVectorD errorsY[5]; + TMatrixD covY[5]; + Double_t chi2CY[5]; + TVectorD paramsZ[5]; + TVectorD errorsZ[5]; + TMatrixD covZ[5]; + Double_t chi2CZ[5]; + for (Int_t i=0;i<5;i++) { + npoints[i]=0; + paramsY[i].ResizeTo(2); + errorsY[i].ResizeTo(2); + covY[i].ResizeTo(2,2); + paramsZ[i].ResizeTo(2); + errorsZ[i].ResizeTo(2); + covZ[i].ResizeTo(2,2); + fittersY[i]->ClearPoints(); + fittersZ[i]->ClearPoints(); + } + // + // Update fitters + // + Int_t countRes=0; + for (Int_t irow=0;irow<159;irow++) { + AliTPCclusterMI *c=track->GetClusterPointer(irow); + if (!c) continue; + if ((c->GetDetector()%36)!=isec) continue; + Double_t dx = c->GetX()-fXmiddle; + Double_t x[2]={dx, dx*dx}; + Double_t yfit = pyf[0]+pyf[1]*dx; + Double_t zfit = pzf[0]+pzf[1]*dx; + Double_t yfitC = pyfc[0]+pyfc[1]*dx; + Double_t zfitC = pzfc[0]+pzfc[1]*dx; + Double_t dedge = c->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(yfit); + if (TMath::Abs(c->GetY()-yfit)>kMaxDist) continue; + if (TMath::Abs(c->GetZ()-zfit)>kMaxDist) continue; + if (dedgeRPhiCOGCorrection(c->GetDetector(),c->GetRow(), c->GetPad(), + c->GetY(),yfit, c->GetZ(), pyf[1], c->GetMax(),2.5); + if (TMath::Abs(corrtrY)>kMaxCorrY) continue; + // + vecX[countRes]=c->GetX(); + vecY[countRes]=c->GetY()-yfit; + vecZ[countRes]=c->GetZ()-zfit; + countRes++; + // + // Fill THnSparse cluster residuals + // use only primary candidates with ITS signal + if (fCurrentTrack->IsOn(0x4)&&TMath::Abs(vImpact[0])<1&&TMath::Abs(vImpact[1])<1){ + Double_t resVector[5]; + resVector[1]= 9.*gphi/TMath::Pi(); + resVector[2]= TMath::Sqrt(c->GetX()*c->GetX()+c->GetY()*c->GetY()); + resVector[3]= c->GetZ()/resVector[2]; + // + // + resVector[0]= c->GetY()-yfitC; + fClusterDelta[0]->Fill(resVector); + resVector[0]= c->GetZ()-zfitC; + fClusterDelta[1]->Fill(resVector); + } + if (c->GetRow()GetRow()>159-kdrow1Fit) continue; + // + + if (c->GetDetector()>35){ + if (c->GetX()AddPoint(x,c->GetY(),0.1); + if (yfit<-kPRFWidth) fittersZ[1]->AddPoint(x,c->GetZ(),0.1); + if (yfit>kPRFWidth) fittersY[2]->AddPoint(x,c->GetY(),0.1); + if (yfit>kPRFWidth) fittersZ[2]->AddPoint(x,c->GetZ(),0.1); + } + if (c->GetX()>fXquadrant){ + if (yfit<-kPRFWidth) fittersY[3]->AddPoint(x,c->GetY(),0.1); + if (yfit<-kPRFWidth) fittersZ[3]->AddPoint(x,c->GetZ(),0.1); + if (yfit>kPRFWidth) fittersY[4]->AddPoint(x,c->GetY(),0.1); + if (yfit>kPRFWidth) fittersZ[4]->AddPoint(x,c->GetZ(),0.1); + } + } + if (c->GetDetector()<36){ + fittersY[0]->AddPoint(x,c->GetY(),0.1); + fittersZ[0]->AddPoint(x,c->GetZ(),0.1); + } + } + // + // + // + for (Int_t i=0;i<5;i++) { + npoints[i] = fittersY[i]->GetNpoints(); + if (npoints[i]>=kMinClusterQ){ + // Y fit + fittersY[i]->Eval(); + Double_t chi2FacY = TMath::Sqrt(fittersY[i]->GetChisquare()/(fittersY[i]->GetNpoints()-2)); + chi2CY[i]=chi2FacY; + fittersY[i]->GetParameters(paramsY[i]); + fittersY[i]->GetErrors(errorsY[i]); + fittersY[i]->GetCovarianceMatrix(covY[i]); + // renormalize errors + errorsY[i][0]*=chi2FacY; + errorsY[i][1]*=chi2FacY; + covY[i](0,0)*=chi2FacY*chi2FacY; + covY[i](0,1)*=chi2FacY*chi2FacY; + covY[i](1,0)*=chi2FacY*chi2FacY; + covY[i](1,1)*=chi2FacY*chi2FacY; + // Z fit + fittersZ[i]->Eval(); + Double_t chi2FacZ = TMath::Sqrt(fittersZ[i]->GetChisquare()/(fittersZ[i]->GetNpoints()-2)); + chi2CZ[i]=chi2FacZ; + fittersZ[i]->GetParameters(paramsZ[i]); + fittersZ[i]->GetErrors(errorsZ[i]); + fittersZ[i]->GetCovarianceMatrix(covZ[i]); + // renormalize errors + errorsZ[i][0]*=chi2FacZ; + errorsZ[i][1]*=chi2FacZ; + covZ[i](0,0)*=chi2FacZ*chi2FacZ; + covZ[i](0,1)*=chi2FacZ*chi2FacZ; + covZ[i](1,0)*=chi2FacZ*chi2FacZ; + covZ[i](1,1)*=chi2FacZ*chi2FacZ; + } + } + // + // void UpdateSectorKalman + // + for (Int_t i0=0;i0<5;i0++){ + for (Int_t i1=i0+1;i1<5;i1++){ + if(npoints[i0]0){ + // get vertex position + // + TTreeSRedirector *cstream = GetDebugStreamer(); + + + if (cstream){ + for (Int_t i0=0;i0<5;i0++){ + for (Int_t i1=i0;i1<5;i1++){ + if (i0==i1) continue; + if(npoints[i0]fXIO){ + if (lx0) quadrant=2; + } + if (lx>fXquadrant) { + if (ly<0) quadrant=3; + if (ly>0) quadrant=4; + } + } + Double_t a10 = (*param)(quadrant*6+0,0); + Double_t a20 = (*param)(quadrant*6+1,0); + Double_t a21 = (*param)(quadrant*6+2,0); + Double_t dx = (*param)(quadrant*6+3,0); + Double_t dy = (*param)(quadrant*6+4,0); + Double_t dz = (*param)(quadrant*6+5,0); + Double_t deltaX = lx-fXIO; + if (coord==0) return dx; + if (coord==1) return dy+deltaX*a10; + if (coord==2) return dz+deltaX*a20+ly*a21; + return 0; +} + +Double_t AliTPCcalibAlign::SGetCorrectionSector(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t lz){ + // + // + // + if (!Instance()) return 0; + return Instance()->GetCorrectionSector(coord,sector,lx,ly,lz); +} + +void AliTPCcalibAlign::MakeKalman(){ + // + // Make a initial Kalman paramaters for sector Alignemnt + // + fSectorParamA = new TMatrixD(6*36+6,1); + fSectorParamC = new TMatrixD(6*36+6,1); + fSectorCovarA = new TMatrixD(6*36+6,6*36+6); + fSectorCovarC = new TMatrixD(6*36+6,6*36+6); + // + // set starting parameters at 0 + // + for (Int_t isec=0;isec<37;isec++) + for (Int_t ipar=0;ipar<6;ipar++){ + (*fSectorParamA)(isec*6+ipar,0) =0; + (*fSectorParamC)(isec*6+ipar,0) =0; + } + // + // set starting covariance + // + for (Int_t isec=0;isec<36;isec++) + for (Int_t ipar=0;ipar<6;ipar++){ + if (ipar<3){ + (*fSectorCovarA)(isec*6+ipar,isec*6+ipar) =0.002*0.002; // 2 mrad + (*fSectorCovarC)(isec*6+ipar,isec*6+ipar) =0.002*0.002; + } + if (ipar>=3){ + (*fSectorCovarA)(isec*6+ipar,isec*6+ipar) =0.02*0.02; // 0.2 mm + (*fSectorCovarC)(isec*6+ipar,isec*6+ipar) =0.02*0.02; + } + } + (*fSectorCovarA)(36*6+0,36*6+0) =0.04; // common shift y up-up + (*fSectorCovarA)(36*6+1,36*6+1) =0.04; // common shift y down-down + (*fSectorCovarA)(36*6+2,36*6+2) =0.04; // common shift y up-down + (*fSectorCovarA)(36*6+3,36*6+3) =0.004; // common shift phi up-up + (*fSectorCovarA)(36*6+4,36*6+4) =0.004; // common shift phi down-down + (*fSectorCovarA)(36*6+5,36*6+5) =0.004; // common shift phi up-down + // + (*fSectorCovarC)(36*6+0,36*6+0) =0.04; // common shift y up-up + (*fSectorCovarC)(36*6+1,36*6+1) =0.04; // common shift y down-down + (*fSectorCovarC)(36*6+2,36*6+2) =0.04; // common shift y up-down + (*fSectorCovarC)(36*6+3,36*6+3) =0.004; // common shift phi up-up + (*fSectorCovarC)(36*6+4,36*6+4) =0.004; // common shift phi down-down + (*fSectorCovarC)(36*6+5,36*6+5) =0.004; // common shift phi up-down +} + +void AliTPCcalibAlign::UpdateKalman(Int_t sector0, Int_t sector1, TMatrixD &p0, TMatrixD &c0, TMatrixD &p1, TMatrixD &c1){ + // + // Update Kalman parameters + // Note numbering from 0..36 0..17 IROC 18..35 OROC + // + // + if (fSectorParamA==NULL) MakeKalman(); + if (CheckCovariance(c0)>0) return; + if (CheckCovariance(c1)>0) return; + const Int_t nelem = 6*36+6; + // + // + static TMatrixD matHk(4,nelem); // vector to mesurement + static TMatrixD measR(4,4); // measurement error + static TMatrixD vecZk(4,1); // measurement + // + static TMatrixD vecYk(4,1); // Innovation or measurement residual + static TMatrixD matHkT(nelem,4); // helper matrix Hk transpose + static TMatrixD matSk(4,4); // Innovation (or residual) covariance + static TMatrixD matKk(nelem,4); // Optimal Kalman gain + static TMatrixD mat1(nelem,nelem); // update covariance matrix + static TMatrixD covXk2(nelem,nelem); // helper matrix + // + TMatrixD *vOrig = 0; + TMatrixD *cOrig = 0; + vOrig = (sector0%36>=18) ? fSectorParamA:fSectorParamC; + cOrig = (sector0%36>=18) ? fSectorCovarA:fSectorCovarC; + // + Int_t sec0= sector0%18; + Int_t sec1= sector1%18; + if (sector0>35) sec0+=18; + if (sector1>35) sec1+=18; + // + TMatrixD vecXk(*vOrig); // X vector + TMatrixD covXk(*cOrig); // X covariance + // + //Unit matrix + // + for (Int_t i=0;i0) dir= 1.; + if (dsec<0) dir=-1.; + if (sector0>35&§or1>35){ + matHk(0,36*6+0)=dir; + matHk(1,36*6+3+0)=dir; + } + if (sector0<36&§or1<36){ + matHk(0,36*6+1)=dir; + matHk(1,36*6+3+1)=dir; + } + if (sector0<36&§or1>35){ + matHk(0,36*6+2)=dir; + matHk(1,36*6+3+2)=dir; + } + if (sector0>35&§or1<36){ + matHk(0,36*6+2)=-dir; + matHk(1,36*6+3+2)=-dir; + } + } + // + // + vecZk =(p1)-(p0); // measurement + measR =(c1)+(c0); // error of measurement + vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual + // + // + matHkT=matHk.T(); matHk.T(); + matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance + matSk.Invert(); + matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain + vecXk += matKk*vecYk; // updated vector + covXk2= (mat1-(matKk*matHk)); + covXk = covXk2*covXk; + + if (CheckCovariance(covXk)>0) return; + + // + // + (*cOrig)=covXk; + (*vOrig)=vecXk; +} + + +void AliTPCcalibAlign::UpdateKalman(TMatrixD &par0, TMatrixD &cov0, TMatrixD &par1, TMatrixD &cov1){ + // + // Update kalman vector para0 with vector par1 + // Used for merging + // + Int_t nelem = 6*36+6; + static TMatrixD matHk(nelem,nelem); // vector to mesurement + static TMatrixD measR(nelem,nelem); // measurement error + static TMatrixD vecZk(nelem,1); // measurement + // + static TMatrixD vecYk(nelem,1); // Innovation or measurement residual + static TMatrixD matHkT(nelem,nelem); // helper matrix Hk transpose + static TMatrixD matSk(nelem,nelem); // Innovation (or residual) covariance + static TMatrixD matKk(nelem,nelem); // Optimal Kalman gain + static TMatrixD mat1(nelem,nelem); // update covariance matrix + static TMatrixD covXk2(nelem,nelem); // helper matrix + // + TMatrixD vecXk(par0); // X vector + TMatrixD covXk(cov0); // X covariance + + // + //Unit matrix + // + for (Int_t i=0;ikEpsilon){ + printf("Error 1 - non symetric matrix %d\t%d\t%f",i0,i1,r1-r0); + nerrors++; + } + if (TMath::Abs(r0)>=1){ + printf("Error 2 - Wrong correlation %d\t%d\t%f\n",i0,i1,r0); + nerrors++; + } + if (TMath::Abs(r1)>=1){ + printf("Error 3 - Wrong correlation %d\t%d\t%f\n",i0,i1,r1); + nerrors++; + } + } + return nerrors; +} + + + +void AliTPCcalibAlign::MakeReportDy(TFile *output){ + // + // Draw histogram of dY + // + Int_t kmicolors[10]={1,2,3,4,6,7,8,9,10,11}; + Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30}; + + AliTPCcalibAlign *align = this; + TVectorD vecOOP(36); + TVectorD vecOOM(36); + TVectorD vecOIP(36); + TVectorD vecOIM(36); + TVectorD vecOIS(36); + TVectorD vecSec(36); + TCanvas * cOROCdy = new TCanvas("OROC dy","OROC dy",900,600); + cOROCdy->Divide(6,6); + TCanvas * cIROCdy = new TCanvas("IROC dy","IROC dy",900,600); + cIROCdy->Divide(6,6); + TCanvas * cDy = new TCanvas("Dy","Dy",600,700); + cDy->Divide(1,2); + for (Int_t isec=0;isec<36;isec++){ + Bool_t isDraw=kFALSE; + vecSec(isec)=isec; + cOROCdy->cd(isec+1); + Int_t secPlus = (isec%18==17)? isec-17:isec+1; + Int_t secMinus= (isec%18==0) ? isec+17:isec-1; + printf("%d\t%d\t%d\n",isec,secPlus,secMinus); + TH1 * hisOOP= align->GetHisto(AliTPCcalibAlign::kY,isec+36,secPlus+36); + TH1 * hisOOM= align->GetHisto(AliTPCcalibAlign::kY,isec+36,secMinus+36); + TH1 * hisOIS= align->GetHisto(AliTPCcalibAlign::kY,isec+36,isec); + + if (hisOIS) { + hisOIS = (TH1F*)(hisOIS->Clone()); + hisOIS->SetDirectory(0); + hisOIS->Scale(1./(hisOIS->GetMaximum()+1)); + hisOIS->SetLineColor(kmicolors[0]); + hisOIS->Draw(); + isDraw = kTRUE; + vecOIS(isec)=10*hisOIS->GetMean(); + } + if (hisOOP) { + hisOOP = (TH1F*)(hisOOP->Clone()); + hisOOP->SetDirectory(0); + hisOOP->Scale(1./(hisOOP->GetMaximum()+1)); + hisOOP->SetLineColor(kmicolors[1]); + if (isDraw) hisOOP->Draw("same"); + if (!isDraw) {hisOOP->Draw(""); isDraw=kTRUE;} + vecOOP(isec)=10*hisOOP->GetMean(); + } + if (hisOOM) { + hisOOM = (TH1F*)(hisOOM->Clone()); + hisOOM->SetDirectory(0); + hisOOM->Scale(1/(hisOOM->GetMaximum()+1)); + hisOOM->SetLineColor(kmicolors[3]); + if (isDraw) hisOOM->Draw("same"); + if (!isDraw) {hisOOM->Draw(""); isDraw=kTRUE;} + vecOOM(isec)=10*hisOOM->GetMean(); + } + } + // + // + for (Int_t isec=0;isec<36;isec++){ + Bool_t isDraw=kFALSE; + cIROCdy->cd(isec+1); + Int_t secPlus = (isec%18==17)? isec-17:isec+1; + Int_t secMinus= (isec%18==0) ? isec+17:isec-1; + printf("%d\t%d\t%d\n",isec,secPlus,secMinus); + TH1 * hisOIP= align->GetHisto(AliTPCcalibAlign::kY,isec+36,secPlus); + TH1 * hisOIM= align->GetHisto(AliTPCcalibAlign::kY,isec+36,secMinus); + TH1 * hisOIS= align->GetHisto(AliTPCcalibAlign::kY,isec+36,isec); + if (hisOIS) { + hisOIS = (TH1F*)(hisOIS->Clone()); + hisOIS->SetDirectory(0); + hisOIS->Scale(1./(hisOIS->GetMaximum()+1)); + hisOIS->SetLineColor(kmicolors[0]); + hisOIS->Draw(); + isDraw = kTRUE; + vecOIS(isec)=10*hisOIS->GetMean(); + } + if (hisOIP) { + hisOIP = (TH1F*)(hisOIP->Clone()); + hisOIP->SetDirectory(0); + hisOIP->Scale(1./(hisOIP->GetMaximum()+1)); + hisOIP->SetLineColor(kmicolors[1]); + if (isDraw) hisOIP->Draw("same"); + if (!isDraw) {hisOIP->Draw(""); isDraw=kTRUE;} + hisOIP->Draw("same"); + vecOIP(isec)=10*hisOIP->GetMean(); + } + if (hisOIM) { + hisOIM = (TH1F*)(hisOIM->Clone()); + hisOIM->SetDirectory(0); + hisOIM->Scale(1/(hisOIM->GetMaximum()+1)); + hisOIM->SetLineColor(kmicolors[3]); + if (isDraw) hisOIM->Draw("same"); + if (!isDraw) {hisOIM->Draw(""); isDraw=kTRUE;} + vecOIM(isec)=10*hisOIM->GetMean(); + } + } + TGraph* grOIM = new TGraph(36,vecSec.GetMatrixArray(),vecOIM.GetMatrixArray()); + TGraph* grOIP = new TGraph(36,vecSec.GetMatrixArray(),vecOIP.GetMatrixArray()); + TGraph* grOIS = new TGraph(36,vecSec.GetMatrixArray(),vecOIS.GetMatrixArray()); + TGraph* grOOM = new TGraph(36,vecSec.GetMatrixArray(),vecOOM.GetMatrixArray()); + TGraph* grOOP = new TGraph(36,vecSec.GetMatrixArray(),vecOOP.GetMatrixArray()); + // + grOIS->SetMarkerStyle(kmimarkers[0]); + grOIP->SetMarkerStyle(kmimarkers[1]); + grOIM->SetMarkerStyle(kmimarkers[3]); + grOOP->SetMarkerStyle(kmimarkers[1]); + grOOM->SetMarkerStyle(kmimarkers[3]); + grOIS->SetMarkerColor(kmicolors[0]); + grOIP->SetMarkerColor(kmicolors[1]); + grOIM->SetMarkerColor(kmicolors[3]); + grOOP->SetMarkerColor(kmicolors[1]); + grOOM->SetMarkerColor(kmicolors[3]); + grOIS->SetLineColor(kmicolors[0]); + grOIP->SetLineColor(kmicolors[1]); + grOIM->SetLineColor(kmicolors[3]); + grOOP->SetLineColor(kmicolors[1]); + grOOM->SetLineColor(kmicolors[3]); + grOIS->SetMaximum(1.5); + grOIS->SetMinimum(-1.5); + grOIS->GetXaxis()->SetTitle("Sector number"); + grOIS->GetYaxis()->SetTitle("#Delta_{y} (mm)"); + + cDy->cd(1); + grOIS->Draw("apl"); + grOIM->Draw("pl"); + grOIP->Draw("pl"); + cDy->cd(2); + grOIS->Draw("apl"); + grOOM->Draw("pl"); + grOOP->Draw("pl"); + cOROCdy->SaveAs("picAlign/OROCOROCdy.eps"); + cOROCdy->SaveAs("picAlign/OROCOROCdy.gif"); + cOROCdy->SaveAs("picAlign/OROCOROCdy.root"); + // + cIROCdy->SaveAs("picAlign/OROCIROCdy.eps"); + cIROCdy->SaveAs("picAlign/OROCIROCdy.gif"); + cIROCdy->SaveAs("picAlign/OROCIROCdy.root"); + // + cDy->SaveAs("picAlign/Sectordy.eps"); + cDy->SaveAs("picAlign/Sectordy.gif"); + cDy->SaveAs("picAlign/Sectordy.root"); + if (output){ + output->cd(); + cOROCdy->Write("OROCOROCDy"); + cIROCdy->Write("OROCIROCDy"); + cDy->Write("SectorDy"); + } +} + +void AliTPCcalibAlign::MakeReportDyPhi(TFile */*output*/){ + // + // + // + Int_t kmicolors[10]={1,2,3,4,6,7,8,9,10,11}; + Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30}; + + AliTPCcalibAlign *align = this; + TCanvas * cOROCdyPhi = new TCanvas("OROC dyphi","OROC dyphi",900,600); + cOROCdyPhi->Divide(6,6); + for (Int_t isec=0;isec<36;isec++){ + cOROCdyPhi->cd(isec+1); + Int_t secPlus = (isec%18==17)? isec-17:isec+1; + Int_t secMinus= (isec%18==0) ? isec+17:isec-1; + //printf("%d\t%d\t%d\n",isec,secPlus,secMinus); + TH2F *htemp=0; + TProfile * profdyphiOOP=0,*profdyphiOOM=0,*profdyphiOOS=0; + htemp = (TH2F*) (align->GetHisto(AliTPCcalibAlign::kYPhi,isec+36,secPlus+36)); + if (htemp) profdyphiOOP= htemp->ProfileX(); + htemp = (TH2F*)(align->GetHisto(AliTPCcalibAlign::kYPhi,isec+36,secMinus+36)); + if (htemp) profdyphiOOM= htemp->ProfileX(); + htemp = (TH2F*)(align->GetHisto(AliTPCcalibAlign::kYPhi,isec+36,isec)); + if (htemp) profdyphiOOS= htemp->ProfileX(); + + if (profdyphiOOS){ + profdyphiOOS->SetLineColor(kmicolors[0]); + profdyphiOOS->SetMarkerStyle(kmimarkers[0]); + profdyphiOOS->SetMarkerSize(0.2); + profdyphiOOS->SetMaximum(0.5); + profdyphiOOS->SetMinimum(-0.5); + profdyphiOOS->SetXTitle("tan(#phi)"); + profdyphiOOS->SetYTitle("#DeltaY (cm)"); + } + if (profdyphiOOP){ + profdyphiOOP->SetLineColor(kmicolors[1]); + profdyphiOOP->SetMarkerStyle(kmimarkers[1]); + profdyphiOOP->SetMarkerSize(0.2); + } + if (profdyphiOOM){ + profdyphiOOM->SetLineColor(kmicolors[3]); + profdyphiOOM->SetMarkerStyle(kmimarkers[3]); + profdyphiOOM->SetMarkerSize(0.2); + } + if (profdyphiOOS){ + profdyphiOOS->Draw(); + }else{ + if (profdyphiOOM) profdyphiOOM->Draw(""); + if (profdyphiOOP) profdyphiOOP->Draw(""); + } + if (profdyphiOOM) profdyphiOOM->Draw("same"); + if (profdyphiOOP) profdyphiOOP->Draw("same"); + + } +} +