From 4e093f35fe337cc6a46b1b77945d237dd0f088d9 Mon Sep 17 00:00:00 2001 From: marian Date: Sat, 18 Sep 2010 08:52:06 +0000 Subject: [PATCH] Adding interface to manipulate with TGeoMatrixes. Used in residual minimization, and for visualization Used in the FitAlignCombined.C macros + // + // Alignment manipulation using TGeoMatrix + + void SetAlignGlobal(const TGeoMatrix * matrixGlobal); + void SetAlignSectors(const TObjArray *arraySector); + TGeoMatrix* GetAlignGlobal() const {return fMatrixGlobal;} + TObjArray * GetAlignSectors() const {return fArraySector;} + // + static AliTPCCalibGlobalMisalignment* CreateOCDBAlign(); + static AliTPCCalibGlobalMisalignment* CreateMeanAlign(const AliTPCCalibGlobalMisalignment *alignIn); + static void DumpAlignment( AliTPCCalibGlobalMisalignment* align, TTreeSRedirector *pcstream, const char *name); --- TPC/AliTPCCalibGlobalMisalignment.cxx | 260 +++++- TPC/AliTPCCalibGlobalMisalignment.h | 25 +- TPC/CalibMacros/FitAlignCombined.C | 1232 +++++++++++++++++++++++++ 3 files changed, 1464 insertions(+), 53 deletions(-) create mode 100644 TPC/CalibMacros/FitAlignCombined.C diff --git a/TPC/AliTPCCalibGlobalMisalignment.cxx b/TPC/AliTPCCalibGlobalMisalignment.cxx index a2cb18e91c2..eae48227073 100644 --- a/TPC/AliTPCCalibGlobalMisalignment.cxx +++ b/TPC/AliTPCCalibGlobalMisalignment.cxx @@ -34,7 +34,11 @@ #include "AliTPCcalibDB.h" #include "AliTPCParam.h" #include - +// +#include "AliAlignObjParams.h" +#include "AliGeomManager.h" +#include "AliCDBManager.h" +#include "AliCDBEntry.h" AliTPCCalibGlobalMisalignment::AliTPCCalibGlobalMisalignment() : AliTPCCorrection("mialign","Misalignment"), @@ -42,8 +46,7 @@ AliTPCCalibGlobalMisalignment::AliTPCCalibGlobalMisalignment() fRotPhiA(0.),fRotPhiC(0.), fdRPhiOffsetA(0.), fdRPhiOffsetC(0.), fQuadrantDQ1(0), fQuadrantDQ2(0), fQuadrantQ2(0), fMatrixGlobal(0), - fMatrixASide(0), fMatrixCSide(0), - fUseGeomanager(kFALSE) + fArraySector(0) { // // default constructor @@ -70,13 +73,24 @@ void AliTPCCalibGlobalMisalignment::SetQuadranAlign(const TVectorD *dq1, const T fQuadrantQ2 = new TVectorD(*q2);; //OROC long pads - OROC medium pads } -void AliTPCCalibGlobalMisalignment::SetGlobalAlign(const TGeoMatrix * matrixGlobal, const TGeoMatrix *matrixA, const TGeoMatrix *matrixC ){ +void AliTPCCalibGlobalMisalignment::SetAlignGlobal(const TGeoMatrix * matrixGlobal){ // - // Set global misalignment as TGeoMatrix + // Set global misalignment + // Object is OWNER // + if (fMatrixGlobal) delete fMatrixGlobal; + fMatrixGlobal=0; if (matrixGlobal) fMatrixGlobal = new TGeoHMatrix(*matrixGlobal); - if (matrixA) fMatrixASide = new TGeoHMatrix(*matrixA); - if (matrixC) fMatrixCSide = new TGeoHMatrix(*matrixC); +} + +void AliTPCCalibGlobalMisalignment::SetAlignSectors(const TObjArray *arraySector){ + // + // Set misalignment TObjArray of TGeoMatrices - for each sector + // Object is OWNER + // + if (fArraySector) delete fArraySector; + fArraySector=0; + if (arraySector) fArraySector = (TObjArray*)arraySector->Clone(); } @@ -168,29 +182,7 @@ void AliTPCCalibGlobalMisalignment::GetCorrection(const Float_t x[],const Short_ dx[0] += (posQG[0]-x[0]); dx[1] += (posQG[1]-x[1]); // - // alignment matrix in local frame // - if (fUseGeomanager){ //loading from the OCDB - Double_t posC[3] ={pos[0],pos[1],pos[2]}; - // - //2. correct the point in the local frame - AliTPCParam *param = AliTPCcalibDB::Instance()->GetParameters(); - if (!param){ - //AliFatal("OCDB not initialized"); - } - TGeoHMatrix *mat = param->GetClusterMatrix(isec); - // - if (mat) mat->LocalToMaster(pos,posC); - Double_t posCG[3]={posC[0],posC[1],posC[2]}; - //3. tranform the corrected point to the global frame - posCG[0]= TMath::Cos(alpha)*posC[0]-TMath::Sin(alpha)*posC[1]; - posCG[1]= TMath::Sin(alpha)*posC[0]+TMath::Cos(alpha)*posC[1]; - posCG[2]= posC[2]; - //4. Add delta - dx[0]+=posCG[0]-x[0]; - dx[1]+=posCG[1]-x[1]; - dx[2]+=posCG[2]-x[2]; - } if (fMatrixGlobal){ // apply global alignment matrix Double_t ppos[3]={x[0],x[1],x[2]}; @@ -201,23 +193,17 @@ void AliTPCCalibGlobalMisalignment::GetCorrection(const Float_t x[],const Short_ dx[2]+=pposC[2]-ppos[2]; } - if (fMatrixASide && roc%36<18){ - // apply global alignment matrix - Double_t ppos[3]={x[0],x[1],x[2]}; - Double_t pposC[3]={x[0],x[1],x[2]}; - fMatrixASide->LocalToMaster(ppos,pposC); - dx[0]+=pposC[0]-ppos[0]; - dx[1]+=pposC[1]-ppos[1]; - dx[2]+=pposC[2]-ppos[2]; - } - if (fMatrixCSide && roc%36>=18){ + if (fArraySector){ // apply global alignment matrix - Double_t ppos[3]={x[0],x[1],x[2]}; - Double_t pposC[3]={x[0],x[1],x[2]}; - fMatrixCSide->LocalToMaster(ppos,pposC); - dx[0]+=pposC[0]-ppos[0]; - dx[1]+=pposC[1]-ppos[1]; - dx[2]+=pposC[2]-ppos[2]; + TGeoMatrix *mat = (TGeoMatrix*)fArraySector->At(isec); + if (mat){ + Double_t ppos[3]={x[0],x[1],x[2]}; + Double_t pposC[3]={x[0],x[1],x[2]}; + mat->LocalToMaster(ppos,pposC); + dx[0]+=pposC[0]-ppos[0]; + dx[1]+=pposC[1]-ppos[1]; + dx[2]+=pposC[2]-ppos[2]; + } } } @@ -233,3 +219,187 @@ void AliTPCCalibGlobalMisalignment::Print(Option_t* /*option*/ ) const { } + +void AliTPCCalibGlobalMisalignment::AddAlign(const AliTPCCalibGlobalMisalignment & add){ + // + // Add the alignmnet to current object + // + fXShift+=add.fXShift; // Shift in global X [cm] + fYShift+=add.fYShift; // Shift in global Y [cm] + fZShift+=add.fZShift; // Shift in global Z [cm] + + fRotPhiA+=add.fRotPhiA; // simple rotation of A side read-out plane around the Z axis [rad] + fRotPhiC+=add.fRotPhiC; // simple rotation of C side read-out plane around the Z axis [rad] + fdRPhiOffsetA+=add.fdRPhiOffsetA; // add a constant offset of dRPhi (or local Y) in [cm]: purely for calibration purposes! + fdRPhiOffsetC+=add.fdRPhiOffsetC; // add a constant offset of dRPhi (or local Y) in [cm]: purely for calibration purposes! + // + // Quadrant alignment + // + if (fQuadrantDQ1&&add.fQuadrantDQ1) fQuadrantDQ1->Add(*(add.fQuadrantDQ1)); //OROC medium pads delta ly+ - ly- + if (fQuadrantDQ2&&add.fQuadrantDQ2) fQuadrantDQ2->Add(*(add.fQuadrantDQ2)); //OROC long pads delta ly+ - ly- + if (fQuadrantQ2&&add.fQuadrantQ2) fQuadrantQ2->Add(*(add.fQuadrantQ2)); //OROC long pads - OROC medium pads + if (!fQuadrantDQ1&&add.fQuadrantDQ1) fQuadrantDQ1= new TVectorD(*add.fQuadrantDQ1); //OROC medium pads delta ly+ - ly- + if (!fQuadrantDQ2&&add.fQuadrantDQ2) fQuadrantDQ2 = new TVectorD(*add.fQuadrantDQ2); //OROC long pads delta ly+ - ly- + if (!fQuadrantQ2&&add.fQuadrantQ2) fQuadrantQ2= new TVectorD(*add.fQuadrantQ2); //OROC long pads - OROC medium pads + // + // Global alignment - use native ROOT representation + // + if (add.fMatrixGlobal){ + if (!fMatrixGlobal) fMatrixGlobal = new TGeoHMatrix(*(add.fMatrixGlobal)); + if (fMatrixGlobal) ((TGeoHMatrix*)fMatrixGlobal)->Multiply(add.fMatrixGlobal); + } + if (add.fArraySector){ + if (!fArraySector) {SetAlignSectors(add.fArraySector); + }else{ + for (Int_t isec=0; isec<72; isec++){ + TGeoHMatrix *mat0= (TGeoHMatrix*)fArraySector->At(isec); + TGeoHMatrix *mat1= (TGeoHMatrix*)add.fArraySector->At(isec); + if (mat0&&mat1) mat0->Multiply(mat1); + } + } + } +} + + +AliTPCCalibGlobalMisalignment * AliTPCCalibGlobalMisalignment::CreateOCDBAlign(){ + // + // Create AliTPCCalibGlobalMisalignment from OCDB Alignment entry + // OCDB has to be initialized before in user code + // All storages (defualt and specific) and run number + // + AliCDBEntry * entry = AliCDBManager::Instance()->Get("TPC/Align/Data"); + if (!entry){ + printf("Missing alignmnet entry. OCDB not initialized?\n"); + return 0; + } + TClonesArray * array = (TClonesArray*)entry->GetObject(); + Int_t entries = array->GetEntries(); + TGeoHMatrix matrixGlobal; + TObjArray *alignArrayOCDB= new TObjArray(73); // sector misalignment + global misalignment + // // global is number 72 + // + { for (Int_t i=0;iUncheckedAt(i); + alignP->GetMatrix(matrix); + Int_t imod; + AliGeomManager::ELayerID ilayer; + alignP->GetVolUID(ilayer, imod); + if (ilayer==AliGeomManager::kInvalidLayer) { + alignArrayOCDB->AddAt(matrix.Clone(),72); + alignP->GetMatrix(matrixGlobal); + }else{ + Int_t sector=imod; + if (ilayer==AliGeomManager::kTPC2) sector+=36; + alignArrayOCDB->AddAt(matrix.Clone(),sector); + } + } + } + AliTPCCalibGlobalMisalignment *align = new AliTPCCalibGlobalMisalignment; + align->SetAlignGlobal(&matrixGlobal); + align->SetAlignSectors(alignArrayOCDB); + return align; +} + + +AliTPCCalibGlobalMisalignment * AliTPCCalibGlobalMisalignment::CreateMeanAlign(const AliTPCCalibGlobalMisalignment *alignIn){ + // + // Create new object, disantangle common mean alignmnet and sector alignment + // + // 1. Try to get mean alignment + // 2. Remove mean alignment from sector alignment + // 3. Create new object + // + TObjArray * array = alignIn->GetAlignSectors(); + TObjArray * arrayNew = new TObjArray(72); + // + //Get mean transformation + TGeoHMatrix matrix; + {for (Int_t isec=0; isec<72; isec++){ + const TGeoMatrix* cmatrix=(TGeoMatrix*)array->At(isec); + if (!cmatrix) continue; + matrix.Multiply(cmatrix); + }} + TGeoHMatrix matrixMean(matrix); + matrixMean.SetDx(matrix.GetTranslation()[0]/72.); + matrixMean.SetDy(matrix.GetTranslation()[1]/72.); + matrixMean.SetDz(matrix.GetTranslation()[2]/72.); + Double_t rotation[12]; + {for (Int_t i=0; i<12; i++) { + rotation[i]=1.0; + if (TMath::Abs(matrix.GetRotationMatrix()[i]-1.)>0.1){ + rotation[i]=matrix.GetRotationMatrix()[i]/72.; + } + }} + matrixMean.SetRotation(rotation); + TGeoHMatrix matrixInv = matrixMean.Inverse(); + // + {for (Int_t isec=0; isec<72; isec++){ + TGeoHMatrix* amatrix=(TGeoHMatrix*)(array->At(isec)->Clone()); + if (!amatrix) continue; + amatrix->Multiply(&matrixInv); + arrayNew->AddAt(amatrix,isec); + }} + if (alignIn->GetAlignGlobal()) matrixMean.Multiply((alignIn->GetAlignGlobal())); + AliTPCCalibGlobalMisalignment *alignOut = new AliTPCCalibGlobalMisalignment; + alignOut->SetAlignGlobal(&matrixMean); + alignOut->SetAlignSectors(arrayNew); + /* + Checks transformation: + AliTPCCalibGlobalMisalignment * alignIn = AliTPCCalibGlobalMisalignment::CreateOCDBAlign() + AliTPCCalibGlobalMisalignment * alignOut = AliTPCCalibGlobalMisalignment::CreateMeanAlign(alignIn) + alignOutM= (AliTPCCalibGlobalMisalignment*)alignOut->Clone(); + alignOutS= (AliTPCCalibGlobalMisalignment*)alignOut->Clone(); + alignOutS->SetAlignGlobal(0); + alignOutM->SetAlignSectors(0); + // + AliTPCCorrection::AddVisualCorrection(alignOut,0); + AliTPCCorrection::AddVisualCorrection(alignOutM,1); + AliTPCCorrection::AddVisualCorrection(alignOutS,2); + AliTPCCorrection::AddVisualCorrection(alignIn,3); + + TF1 f0("f0","AliTPCCorrection::GetCorrSector(x,85,0.9,1,0)",0,18); + TF1 f1("f1","AliTPCCorrection::GetCorrSector(x,85,0.9,1,1)",0,18); + TF1 f2("f2","AliTPCCorrection::GetCorrSector(x,85,0.9,1,2)",0,18); + TF1 f3("f3","AliTPCCorrection::GetCorrSector(x,85,0.9,1,3)",0,18); + f0->SetLineColor(1); + f1->SetLineColor(2); + f2->SetLineColor(3); + f3->SetLineColor(4); + f0->Draw(); + f1->Draw("same"); + f2->Draw("same"); + f3->Draw("same"); + + TF2 f2D("f2D","AliTPCCorrection::GetCorrSector(x,y,0.9,1,0)-AliTPCCorrection::GetCorrSector(x,y,0.9,1,3)",0,18,85,245); + */ + return alignOut; +} + + +void AliTPCCalibGlobalMisalignment::DumpAlignment( AliTPCCalibGlobalMisalignment* align, TTreeSRedirector *pcstream, const char *name){ + // + // Dump alignment into tree + // + TObjArray * array = align->GetAlignSectors(); + if (!array) return; + // + //Get mean transformation + TGeoHMatrix matrix; + {for (Int_t isec=0; isec<72; isec++){ + TGeoHMatrix* cmatrix=(TGeoHMatrix*)array->At(isec); + TGeoHMatrix* cmatrixDown=(TGeoHMatrix*)array->At(isec%36); + TGeoHMatrix* cmatrixUp=(TGeoHMatrix*)array->At(isec%36+36); + TGeoHMatrix diff(*cmatrixDown); + diff.Multiply(&(cmatrixUp->Inverse())); + (*pcstream)<Macro("~/rootlogon.C"); + gSystem->Load("libANALYSIS"); + gSystem->Load("libSTAT"); + gSystem->Load("libTPCcalib"); + gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros -I$ALICE_ROOT/TPC/TPC -I$ALICE_ROOT/STAT"); + .L $ALICE_ROOT/TPC/CalibMacros/FitAlignCombined.C+ + .x ConfigCalibTrain.C(119047) + // + // + FitAlignCombinedCorr(); + + +*/ + + + +#if !defined(__CINT__) || defined(__MAKECINT__) +#include "TH1D.h" +#include "TH2F.h" +#include "THnSparse.h" +#include "TVectorD.h" +#include "TTreeStream.h" +#include "TFile.h" +#include "TChain.h" +#include "AliTPCcalibAlign.h" +#include "AliTPCROC.h" +#include "AliTPCCalROC.h" +#include "AliTPCCalPad.h" +#include "TF1.h" +#include "TH2.h" +#include "TH3.h" +#include "TROOT.h" +#include "TProfile.h" +#include "AliTPCPreprocessorOnline.h" +#include "AliTPCcalibDB.h" +#include "AliTPCkalmanAlign.h" +#include "TPostScript.h" +#include "TLegend.h" +#include "TCut.h" +#include "TCanvas.h" +#include "TStyle.h" +#include "AliLog.h" +#include "AliTPCExBEffectiveSector.h" +#include "AliTPCComposedCorrection.h" +#include "AliTPCCalibGlobalMisalignment.h" +// +#include "AliCDBMetaData.h" +#include "AliCDBId.h" +#include "AliCDBManager.h" +#include "AliCDBStorage.h" +#include "AliCDBEntry.h" +#include "TStatToolkit.h" +#include "AliAlignObjParams.h" +#include "AliTPCParam.h" +#endif + + +void DrawFitQA(); + +AliTPCROC *proc = AliTPCROC::Instance(); +AliTPCCalibGlobalMisalignment *combAlignOCDBOld=0; // old misalignment from OCDB - used for reconstruction +AliTPCCalibGlobalMisalignment *combAlignOCDBNew=0; // new misalignment +AliTPCCalibGlobalMisalignment *combAlignGlobal=0; // delta misalignment - global part +AliTPCCalibGlobalMisalignment *combAlignLocal=0; // delta misalignment - sector/local part +// + +TTreeSRedirector *pcstream= 0; // Workspace to store the fit parameters and QA histograms +// +TChain * chain=0; +TChain * chainZ=0; // trees with z measurement +TTree * tree=0; +TTree *itsdy=0; +TTree *itsdyP=0; +TTree *itsdyM=0; +TTree *itsdp=0; +TTree *itsdphiP=0; +TTree *itsdphiM=0; + +TTree *itsdz=0; +TTree *itsdt=0; +TTree *tofdy=0; +TTree *trddy=0; +// +TTree *vdy=0; +TTree *vdyP=0; +TTree *vdyM=0; +TTree *vds=0; +TTree *vdz=0; +TTree *vdt=0; + +TCut cutS="entries>1000&&abs(snp)<0.2&&abs(theta)<1."; + + +void RegisterAlignFunction(){ + // + // Register primitive alignment components. + // Linear conbination of primitev forulas used for fit + // The nominal delta 1 mm in shift and 1 mrad in rotation + // Primitive formulas registeren in AliTPCCoreection::AddvisualCorrection + // 0 - deltaX + // 1 - deltaY + // 2 - deltaZ + // 3 - rot0 (phi) + // 4 - rot1 (theta) + // 5 - rot2 + TGeoHMatrix matrixX; + TGeoHMatrix matrixY; + TGeoHMatrix matrixZ; + TGeoRotation rot0; + TGeoRotation rot1; + TGeoRotation rot2; //transformation matrices + TGeoRotation rot90; //transformation matrices + matrixX.SetDx(0.1); matrixY.SetDy(0.1); matrixZ.SetDz(0.1); //1 mm translation + rot0.SetAngles(0.001*TMath::RadToDeg(),0,0); + rot1.SetAngles(0,0.001*TMath::RadToDeg(),0); + rot2.SetAngles(0,0,0.001*TMath::RadToDeg()); + //how to get rot02 ? + rot90.SetAngles(0,90,0); + rot2.MultiplyBy(&rot90,kTRUE); + rot90.SetAngles(0,-90,0); + rot2.MultiplyBy(&rot90,kFALSE); + AliTPCCalibGlobalMisalignment *alignRot0 =new AliTPCCalibGlobalMisalignment; + alignRot0->SetAlignGlobal(&rot0); + AliTPCCalibGlobalMisalignment *alignRot1=new AliTPCCalibGlobalMisalignment; + alignRot1->SetAlignGlobal(&rot1); + AliTPCCalibGlobalMisalignment *alignRot2=new AliTPCCalibGlobalMisalignment; + alignRot2->SetAlignGlobal(&rot2); + // + AliTPCCalibGlobalMisalignment *alignTrans0 =new AliTPCCalibGlobalMisalignment; + alignTrans0->SetAlignGlobal(&matrixX); + AliTPCCalibGlobalMisalignment *alignTrans1=new AliTPCCalibGlobalMisalignment; + alignTrans1->SetAlignGlobal(&matrixY); + AliTPCCalibGlobalMisalignment *alignTrans2=new AliTPCCalibGlobalMisalignment; + alignTrans2->SetAlignGlobal(&matrixZ); + AliTPCCorrection::AddVisualCorrection(alignTrans0 ,0); + AliTPCCorrection::AddVisualCorrection(alignTrans1 ,1); + AliTPCCorrection::AddVisualCorrection(alignTrans2 ,2); + AliTPCCorrection::AddVisualCorrection(alignRot0 ,3); + AliTPCCorrection::AddVisualCorrection(alignRot1 ,4); + AliTPCCorrection::AddVisualCorrection(alignRot2 ,5); + TObjArray arrAlign(6); + arrAlign.AddAt(alignTrans0->Clone(),0); + arrAlign.AddAt(alignTrans1->Clone(),1); + arrAlign.AddAt(alignTrans2->Clone(),2); + arrAlign.AddAt(alignRot0->Clone(),3); + arrAlign.AddAt(alignRot1->Clone(),4); + arrAlign.AddAt(alignRot2->Clone(),5); + //combAlign.SetCorrections((TObjArray*)arrAlign.Clone()); +} + +AliTPCCalibGlobalMisalignment * MakeAlignFunctionGlobal(TVectorD paramYGlobal){ + // + // Take a fit parameters and make a combined correction + // 1. Take the common part + // 3. Make combined AliTPCCalibGlobalMisalignment - register it + // 4. Compare the aliases with fit values - IT is OK + // + AliTPCCalibGlobalMisalignment *alignGlobal =new AliTPCCalibGlobalMisalignment; + TGeoHMatrix matGlobal; // global parameters + TGeoHMatrix matDelta; // delta A side - C side + // + TGeoHMatrix matrixX; + TGeoHMatrix matrixY; + TGeoHMatrix matrixZ; + TGeoRotation rot0; + TGeoRotation rot1; + TGeoRotation rot2; //transformation matrices + TGeoRotation rot90; //transformation matrices + // + // + matrixX.SetDx(0.1*paramYGlobal[1]); + matrixY.SetDy(0.1*paramYGlobal[2]); + rot0.SetAngles(0.001*TMath::RadToDeg()*paramYGlobal[3],0,0); + rot1.SetAngles(0,0.001*TMath::RadToDeg()*paramYGlobal[4],0); + rot2.SetAngles(0,0,0.001*TMath::RadToDeg()*paramYGlobal[5]); + rot90.SetAngles(0,90,0); + rot2.MultiplyBy(&rot90,kTRUE); + rot90.SetAngles(0,-90,0); + rot2.MultiplyBy(&rot90,kFALSE); + matGlobal.Multiply(&matrixX); + matGlobal.Multiply(&matrixY); + matGlobal.Multiply(&rot0); + matGlobal.Multiply(&rot1); + matGlobal.Multiply(&rot2); + alignGlobal->SetAlignGlobal((TGeoMatrix*)matGlobal.Clone()); + AliTPCCorrection::AddVisualCorrection(alignGlobal ,100); + // + /* + chain->Draw("mean-deltaG:int(sector)",cutS+"type==0&&refX==0&&theta>0","profsame"); + chain->Draw("mean-deltaG:int(sector+18)",cutS+"type==0&&refX==0&&theta<0","profsame"); + + chain->Draw("deltaG:fitYGlobal",cutS+"type==0"); + chain->Draw("deltaG:fitYGlobal",cutS+"type==2"); + + */ + + return alignGlobal; +} + + +AliTPCCalibGlobalMisalignment * MakeAlignFunctionSector(TVectorD paramYLocal){ + // + // Take a fit parameters and make a combined correction: + // Only delta local Y and delta phi are fitted - not sensityvity for other parameters + // Algorithm: + // 1. Loop over sectors + // 2. Make combined AliTPCCalibGlobalMisalignment - register it + // 3. Compare the aliases with fit values - IT is OK + // + AliTPCCalibGlobalMisalignment *alignLocal =new AliTPCCalibGlobalMisalignment; + TGeoHMatrix matrixX; + TGeoHMatrix matrixY; + TGeoHMatrix matrixZ; + TGeoRotation rot0; + TGeoRotation rot1; + TGeoRotation rot2; //transformation matrices + TGeoRotation rot90; //transformation matrices + TObjArray * array = new TObjArray(72); + TVectorD vecLY(72); + TVectorD vecGY(72); + TVectorD vecGX(72); + TVectorD vecPhi(72); + TVectorD vecSec(72); + // + // + {for (Int_t isec=0; isec<72; isec++){ + Int_t sector=isec%18; + Int_t offset=sector*4+1; + if ((isec%36)>=18) offset+=2; + Double_t phi= (Double_t(sector)+0.5)*TMath::Pi()/9.; + // + Double_t dly = paramYLocal[offset]; + Double_t dphi= paramYLocal[offset+1]; + Double_t dgx = TMath::Cos(phi)*0+TMath::Sin(phi)*dly; + Double_t dgy = TMath::Sin(phi)*0-TMath::Cos(phi)*dly; + vecSec[isec]=isec; + vecLY[isec]=dly; + vecGX[isec]=dgx; + vecGY[isec]=dgy; + vecPhi[isec]=dphi; + // + matrixX.SetDx(dgx); + matrixY.SetDy(dgy); + rot0.SetAngles(0.001*TMath::RadToDeg()*dphi,0,0); + TGeoHMatrix matrixSector; // global parameters + matrixSector.Multiply(&matrixX); + matrixSector.Multiply(&matrixY); + matrixSector.Multiply(&rot0); + array->AddAt(matrixSector.Clone(),isec); + }} + alignLocal->SetAlignSectors(array); + AliTPCCorrection::AddVisualCorrection(alignLocal,101); + // + // + TGraph * graphLY = new TGraph(72,vecSec.GetMatrixArray(), vecLY.GetMatrixArray()); + TGraph * graphGX = new TGraph(72,vecSec.GetMatrixArray(), vecGX.GetMatrixArray()); + TGraph * graphGY = new TGraph(72,vecSec.GetMatrixArray(), vecGY.GetMatrixArray()); + TGraph * graphPhi = new TGraph(72,vecSec.GetMatrixArray(), vecPhi.GetMatrixArray()); + graphLY->SetMarkerStyle(25); graphGX->SetMarkerStyle(25); graphGY->SetMarkerStyle(25); graphPhi->SetMarkerStyle(25); + graphLY->SetName("DeltaLY"); graphLY->SetTitle("#Delta_{ly} (mm)"); + graphGX->SetName("DeltaGX"); graphGX->SetTitle("#Delta_{gx} (mm)"); + graphGY->SetName("DeltaGY"); graphGY->SetTitle("#Delta_{gy} (mm)"); + graphPhi->SetName("DeltaPhi"); graphPhi->SetTitle("#Delta_{phi} (mrad)"); + // + graphLY->Write("grDeltaLY"); + graphGX->Write("grDeltaGX"); + graphGY->Write("grDeltaGY"); + graphPhi->Write("grDeltaPhi"); + // Check: + /* + graphLY.Draw("alp") + chain->Draw("mean-deltaG:int(sector)",cutS+"type==0&&refX==0&&theta>0","profsame"); + chain->Draw("mean-deltaG:int(sector+18)",cutS+"type==0&&refX==0&&theta<0","profsame"); + + graphPhi.Draw("alp") + chain->Draw("1000*(mean-deltaG):int(sector)",cutS+"type==2&&refX==0&&theta>0","profsame"); + chain->Draw("1000*(mean-deltaG):int(sector+18)",cutS+"type==2&&refX==0&&theta<0","profsame"); + + // + chain->Draw("deltaL:fitYLocal",cutS+"type==2&&refX==0",""); // phi shift - OK + chain->Draw("deltaL:fitYLocal",cutS+"type==0&&refX==0",""); // OK + chain->Draw("deltaL:fitYLocal",cutS+"type==0",""); // OK + + */ + + return alignLocal; +} + + + + + +void LoadTrees(){ + // + // make sector alignment - using Kalman filter method -AliTPCkalmanAlign + // Combined information is used, mean residuals are minimized: + // + // 1. TPC-ITS alignment + // 2. TPC vertex alignment + // + TFile *f0 = new TFile("../mergeField0/mean.root"); + TFile *fP= new TFile("../mergePlus/mean.root"); + TFile *fM= new TFile("../mergeMinus/mean.root"); + // + itsdy=(TTree*)f0->Get("ITSdy"); + itsdyP=(TTree*)fP->Get("ITSdy"); + itsdyM=(TTree*)fM->Get("ITSdy"); + itsdp=(TTree*)f0->Get("ITSdsnp"); + itsdphiP=(TTree*)fP->Get("ITSdsnp"); + itsdphiM=(TTree*)fM->Get("ITSdsnp"); + + itsdz=(TTree*)f0->Get("ITSdz"); + itsdt=(TTree*)f0->Get("ITSdtheta"); + tofdy=(TTree*)f0->Get("TOFdy"); + trddy=(TTree*)f0->Get("TRDdy"); + // + vdy=(TTree*)f0->Get("Vertexdy"); + vdyP=(TTree*)fP->Get("Vertexdy"); + vdyM=(TTree*)fM->Get("Vertexdy"); + vds=(TTree*)f0->Get("Vertexdsnp"); + vdz=(TTree*)f0->Get("Vertexdz"); + vdt=(TTree*)f0->Get("Vertexdtheta"); + chain = new TChain("mean","mean"); + chainZ = new TChain("mean","mean"); + chain->AddFile("../mergeField0/mean.root",10000000,"ITSdy"); + chain->AddFile("../mergeField0/mean.root",10000000,"ITSdsnp"); + chain->AddFile("../mergeField0/mean.root",10000000,"Vertexdy"); + chain->AddFile("../mergeField0/mean.root",10000000,"Vertexdsnp"); + chain->AddFile("../mergeField0/mean.root",10000000,"TOFdy"); + chain->AddFile("../mergeField0/mean.root",10000000,"TRDdy"); + // + chainZ->AddFile("../mergeField0/mean.root",10000000,"Vertexdz"); + chainZ->AddFile("../mergeField0/mean.root",10000000,"Vertexdtheta"); + chainZ->AddFile("../mergeField0/mean.root",10000000,"ITSdz"); + chainZ->AddFile("../mergeField0/mean.root",10000000,"ITSdtheta"); + + // + itsdy->AddFriend(itsdp,"Phi"); + itsdy->AddFriend(itsdphiP,"PhiP"); + itsdy->AddFriend(itsdphiM,"PhiM"); + itsdy->AddFriend(itsdz,"Z"); + itsdy->AddFriend(itsdt,"T"); + itsdy->AddFriend(itsdyP,"YP"); + itsdy->AddFriend(itsdyM,"YM"); + // + itsdy->AddFriend(vdy,"V"); + itsdy->AddFriend(vdyP,"VP"); + itsdy->AddFriend(vdyP,"VM"); + itsdy->AddFriend(tofdy,"TOF"); + itsdy->AddFriend(trddy,"TRD"); + itsdy->AddFriend(vds,"VPhi"); + itsdy->AddFriend(vdz,"VZ"); + itsdy->AddFriend(vdt,"VT"); + itsdy->AddFriend(tofdy,"TOF."); + itsdy->SetMarkerStyle(25); + itsdy->SetMarkerSize(0.4); + tree=itsdy; + tree->SetAlias("side","(-1+2*(theta>0))"); + chain->SetAlias("side","(-1+2*(theta>0))"); + chain->SetMarkerStyle(25); + chain->SetMarkerSize(0.5); +} + + +void FitAlignCombinedCorr(){ + // + // Fit Global X and globalY shift at vertex and at ITS + // + RegisterAlignFunction(); + LoadTrees(); + combAlignOCDBOld = AliTPCCalibGlobalMisalignment::CreateOCDBAlign(); + // + + pcstream= new TTreeSRedirector("fitAlignCombined.root"); + + TStatToolkit toolkit; + Double_t chi2=0; + Int_t npoints=0; + // Fit vectors + // global fits + TVectorD vecSec(37); + TVectorD paramYGlobal; + TVectorD paramYLocal; + TMatrixD covar; + // make Aliases + {for (Int_t isec=0; isec<18; isec++){ // sectors + tree->SetAlias(Form("sec%d",isec), Form("(abs(sector-%f)<0.5)",isec+0.5)); + chain->SetAlias(Form("sec%d",isec), Form("(abs(sector-%f)<0.5)",isec+0.5)); + }} + for (Int_t isec=-1;isec<36; isec++){ + vecSec[isec+1]=isec; + } + chain->SetAlias("err","((type==0)*0.1+(type==2)*0.001)"); + // delta phi + chain->SetAlias("dpgx","(0+0)"); + chain->SetAlias("dpgy","(0+0)"); + chain->SetAlias("dpgr0","((AliTPCCorrection::GetCorrSector(sector,245,theta,1,3)-AliTPCCorrection::GetCorrSector(sector,85,theta,1,3))/160)"); + chain->SetAlias("dpgr1","((AliTPCCorrection::GetCorrSector(sector,245,theta,1,4)-AliTPCCorrection::GetCorrSector(sector,85,theta,1,4))/160)"); + chain->SetAlias("dpgr2","((AliTPCCorrection::GetCorrSector(sector,245,theta,1,5)-AliTPCCorrection::GetCorrSector(sector,85,theta,1,5))/160)"); + chain->SetAlias("deltaPG","((AliTPCCorrection::GetCorrSector(sector,245,theta,1,100)-AliTPCCorrection::GetCorrSector(sector,85,theta,1,100))/160)"); + chain->SetAlias("deltaPL","((AliTPCCorrection::GetCorrSector(sector,245,theta,1,101)-AliTPCCorrection::GetCorrSector(sector,85,theta,1,101))/160)"); + // delta y at 85 cm + chain->SetAlias("dygxT","(AliTPCCorrection::GetCorrSector(sector,85,theta,1,0)+0)"); + chain->SetAlias("dygyT","(AliTPCCorrection::GetCorrSector(sector,85,theta,1,1)+0)"); + chain->SetAlias("dyr0T","(AliTPCCorrection::GetCorrSector(sector,85,theta,1,3)+0)"); + chain->SetAlias("dyr1T","(AliTPCCorrection::GetCorrSector(sector,85,theta,1,4)+0)"); + chain->SetAlias("dyr2T","(AliTPCCorrection::GetCorrSector(sector,85,theta,1,5)+0)"); + chain->SetAlias("deltaYGT","(AliTPCCorrection::GetCorrSector(sector,85,theta,1,100)+0)"); + chain->SetAlias("deltaYLT","(AliTPCCorrection::GetCorrSector(sector,85,theta,1,101)+0)"); + // + // delta y at reference X + chain->SetAlias("dygx","(dygxT+0)"); // due global x shift + chain->SetAlias("dygy","(dygyT+0)"); // due global y shift + chain->SetAlias("dyr0","(dyr0T+dpgr0*(refX-85.))"); // due rotation 0 + chain->SetAlias("dyr1","(dyr1T+dpgr1*(refX-85.))"); // due rotation 1 + chain->SetAlias("dyr2","(dyr2T+dpgr2*(refX-85.))"); // due rotation 2 + chain->SetAlias("deltaYG","(deltaYGT+deltaPG*(refX-85.))"); // due global distortion + chain->SetAlias("deltaYL","(deltaYLT+deltaPL*(refX-85.))"); // due local distortion + + chain->SetAlias("deltaG","(type==0)*deltaYG+(type==2)*deltaPG"); //alias to global function + chain->SetAlias("deltaL","(type==0)*deltaYL+(type==2)*deltaPL"); //alias to local function + // + TString fstringGlobal=""; + // global part + fstringGlobal+="((type==0)*dygx+0)++"; + fstringGlobal+="((type==0)*dygy+0)++"; + fstringGlobal+="((type==0)*dyr0+(type==2)*dpgr0)++"; + fstringGlobal+="((type==0)*dyr1+(type==2)*dpgr1)++"; + fstringGlobal+="((type==0)*dyr2+(type==2)*dpgr2)++"; + // + // Make global fits + // + TString *strFitYGlobal = TStatToolkit::FitPlane(chain,"mean:err", fstringGlobal.Data(),cutS+"", chi2,npoints,paramYGlobal,covar,-1,0, 10000000, kTRUE); + combAlignGlobal =MakeAlignFunctionGlobal(paramYGlobal); + printf("chi2=%f\n",TMath::Sqrt(chi2/npoints)); + chain->SetAlias("fitYGlobal",strFitYGlobal->Data()); + paramYGlobal.Print(); + paramYGlobal.Write("paramYGlobal"); + // + // + // + // + {for (Int_t isec=0; isec<18; isec++){ + //A side + chain->SetAlias(Form("glyASec%d",isec),Form("((type==0)*sec%d*(theta>0))",isec)); // ly shift + chain->SetAlias(Form("gxASec%d",isec),Form("((type==0)*dygx*sec%d*(theta>0))",isec)); // gx shift + chain->SetAlias(Form("gyASec%d",isec),Form("((type==0)*dygy*sec%d*(theta>0))",isec)); // gy shift + chain->SetAlias(Form("gpASec%d",isec),Form("(((type==0)*dyr0+(type==2)*dpgr0)*sec%d)*(theta>0)",isec)); // phi rotation + //C side + chain->SetAlias(Form("glyCSec%d",isec),Form("((type==0)*sec%d*(theta<0))",isec)); // ly shift + chain->SetAlias(Form("gxCSec%d",isec),Form("((type==0)*dygx*sec%d*(theta<0))",isec)); + chain->SetAlias(Form("gyCSec%d",isec),Form("((type==0)*dygy*sec%d*(theta<0))",isec)); + chain->SetAlias(Form("gpCSec%d",isec),Form("(((type==0)*dyr0+(type==2)*dpgr0)*sec%d)*(theta<0)",isec)); // phi rotation + }} + + TString fstringLocal=""; + {for (Int_t isec=0; isec<18; isec++){ + //fstringLocal+=Form("((type==0)*dygx*sec%d*(theta>0))++",isec); + // fstringLocal+=Form("(gxASec%d)++",isec); + // fstringLocal+=Form("(gyASec%d)++",isec); + fstringLocal+=Form("(glyASec%d)++",isec); + fstringLocal+=Form("(gpASec%d)++",isec); + fstringLocal+=Form("(glyCSec%d)++",isec); + //fstringLocal+=Form("(gxCSec%d)++",isec); + //fstringLocal+=Form("(gyCSec%d)++",isec); + fstringLocal+=Form("(gpCSec%d)++",isec); + }} + + TString *strFitYLocal = TStatToolkit::FitPlane(chain,"(mean-deltaG):err", fstringLocal.Data(),cutS+"", chi2,npoints,paramYLocal,covar,-1,0, 10000000, kTRUE); + printf("chi2=%f\n",TMath::Sqrt(chi2/npoints)); + chain->SetAlias("fitYLocal",strFitYLocal->Data()); + combAlignLocal =MakeAlignFunctionSector(paramYLocal); + //paramYLocal.Print(); + paramYLocal.Write("paramYLocal"); + // + // scaling - why do we need it? + TString fstringScale=""; + fstringScale+="((type==1)*refX+(type==2))*dsec++"; + fstringScale+="((type==1)*refX+(type==2))*dsec*cos(phi)++"; + fstringScale+="((type==1)*refX+(type==2))*dsec*sin(phi)++"; + fstringScale+="((type==1)*refX+(type==2))*dsec*side++"; + fstringScale+="((type==1)*refX+(type==2))*dsec*side*cos(phi)++"; + fstringScale+="((type==1)*refX+(type==2))*dsec*side*sin(phi)++"; + // + // + // + pcstream->GetFile()->cd(); + DrawFitQA(); + // + AliTPCCalibGlobalMisalignment * combAlignOCDBNew0 =( AliTPCCalibGlobalMisalignment *)combAlignOCDBOld->Clone(); + combAlignOCDBNew0->AddAlign(*combAlignGlobal); + combAlignOCDBNew0->AddAlign(*combAlignLocal); + combAlignOCDBNew= AliTPCCalibGlobalMisalignment::CreateMeanAlign(combAlignOCDBNew0); + // + combAlignGlobal->SetName("fcombAlignGlobal"); + combAlignLocal->SetName("fcombAlignLocal"); + combAlignOCDBOld->SetName("fcombAlignOCDBOld"); + combAlignOCDBNew->SetName("fcombAlignOCDBNew"); + combAlignOCDBNew0->SetName("fcombAlignOCDBNew0"); + // + combAlignGlobal->Write("fcombAlignGlobal"); + combAlignLocal->Write("fcombAlignLocal"); + combAlignOCDBOld->Write("fcombAlignOCDBOld"); + combAlignOCDBNew->Write("fcombAlignOCDBNew"); + combAlignOCDBNew0->Write("fcombAlignOCDBNew0"); + combAlignOCDBNew->Print(); + delete pcstream; +} + +void DrawFitQA(){ + // + // + // + // + // MakeQA plot 1D + // + TCanvas c; + c.SetLeftMargin(0.15); + chain->Draw("1000*(mean-deltaG)>>his(100,-1.5,1.5)",cutS+"type==2&&refX==0",""); + chain->GetHistogram()->SetName("TPC_VertexDphiGlobal1D"); + chain->GetHistogram()->SetTitle("TPC-Vertex #Delta_{#phi} - (Global Fit)"); + chain->GetHistogram()->GetXaxis()->SetTitle("#Delta_{#phi} (mrad)"); + chain->GetHistogram()->Fit("gaus","qnr"); + chain->GetHistogram()->Write(); + // + chain->Draw("1000*(mean-deltaG)>>his(100,-1.5,1.5)",cutS+"type==2&&abs(refX-85)<2",""); + chain->GetHistogram()->SetName("TPC_ITSDphiGlobal1D"); + chain->GetHistogram()->SetTitle("TPC-ITS #Delta_{#phi} - (Global Fit)"); + chain->GetHistogram()->GetXaxis()->SetTitle("#Delta_{#phi} (mrad)"); + chain->GetHistogram()->Fit("gaus","qnr"); + chain->GetHistogram()->Write(); + // + chain->Draw("10*(mean-deltaG)>>his(100,-1,1)",cutS+"type==0&&abs(refX-85)<2",""); + chain->GetHistogram()->SetName("TPC_ITSDRPhiGlobal1D"); + chain->GetHistogram()->SetTitle("TPC-ITS #Delta_{r#phi} - (Global Fit)"); + chain->GetHistogram()->GetXaxis()->SetTitle("#Delta_{r#phi} (mm)"); + chain->GetHistogram()->Fit("gaus","qnr"); + chain->GetHistogram()->Write(); + // + chain->Draw("10*(mean-deltaG)>>his(100,-1,1)",cutS+"type==0&&abs(refX-0)<2",""); + chain->GetHistogram()->SetName("TPC-VertexDRPhiGlobal1D"); + chain->GetHistogram()->SetTitle("TPC-Vertex #Delta_{r#phi} - (Global Fit)"); + chain->GetHistogram()->GetXaxis()->SetTitle("#Delta_{r#phi} (mm)"); + chain->GetHistogram()->Fit("gaus","qnr"); + chain->GetHistogram()->Write(); + // + // Make QA plot 3D + // + chain->Draw("1000*(mean-deltaG):sector:abs(theta)",cutS+"theta>0&&type==2&&refX==0","colz"); + chain->GetHistogram()->SetName("TPC_VertexDphi3DGlobal3D"); + chain->GetHistogram()->SetTitle("TPC-Vertex #Delta_{#phi} - (Global Fit)"); + chain->GetHistogram()->GetYaxis()->SetTitle("#Delta_{#phi} (mrad)"); + chain->GetHistogram()->GetXaxis()->SetTitle("sector"); + chain->GetHistogram()->GetZaxis()->SetTitle("tan(#Theta)"); + chain->GetHistogram()->Draw("colz"); + chain->GetHistogram()->Write(); + // + chain->Draw("1000*(mean-deltaG):sector:abs(theta)",cutS+"theta>0&&type==2&&abs(refX-85)<2","colz"); + chain->GetHistogram()->SetName("TPC_ITSDphi3DGlobal3D"); + chain->GetHistogram()->SetTitle("TPC-ITS #Delta_{#phi} - (Global Fit)"); + chain->GetHistogram()->GetYaxis()->SetTitle("#Delta_{#phi} (mrad)"); + chain->GetHistogram()->GetXaxis()->SetTitle("sector"); + chain->GetHistogram()->GetZaxis()->SetTitle("tan(#Theta)"); + chain->GetHistogram()->Draw("colz"); + chain->GetHistogram()->Write(); + // + chain->Draw("10*(mean-deltaG):sector:abs(theta)",cutS+"theta>0&&type==0&&abs(refX-85)<2","colz"); + chain->GetHistogram()->SetName("TPC_ITSDRphi3DGlobal3D"); + chain->GetHistogram()->SetTitle("TPC-ITS #Delta_{r#phi} - (Global Fit)"); + chain->GetHistogram()->GetYaxis()->SetTitle("#Delta_{r#phi} (mm)"); + chain->GetHistogram()->GetXaxis()->SetTitle("sector"); + chain->GetHistogram()->GetZaxis()->SetTitle("tan(#Theta)"); + chain->GetHistogram()->Draw("colz"); + chain->GetHistogram()->Write(); + // + chain->Draw("10*(mean-deltaG):sector:abs(theta)",cutS+"theta>0&&type==0&&abs(refX-0)<2","colz"); + chain->GetHistogram()->SetName("TPC_VertexDRphi3DGlobal3D"); + chain->GetHistogram()->SetTitle("TPC-Vertex #Delta_{r#phi} - (Global Fit)"); + chain->GetHistogram()->GetYaxis()->SetTitle("#Delta_{r#phi} (mm)"); + chain->GetHistogram()->GetXaxis()->SetTitle("sector"); + chain->GetHistogram()->GetZaxis()->SetTitle("tan(#Theta)"); + chain->GetHistogram()->Draw("colz"); + chain->GetHistogram()->Write(); + // + // + // + chain->Draw("1000*(mean-deltaG-deltaL):sector:abs(theta)",cutS+"theta>0&&type==2&&refX==0","colz"); + chain->GetHistogram()->SetName("TPC_VertexDphi3DLocal3D"); + chain->GetHistogram()->SetTitle("TPC-Vertex #Delta_{#phi} - (Local Fit)"); + chain->GetHistogram()->GetYaxis()->SetTitle("#Delta_{#phi} (mrad)"); + chain->GetHistogram()->GetXaxis()->SetTitle("sector"); + chain->GetHistogram()->GetZaxis()->SetTitle("tan(#Theta)"); + chain->GetHistogram()->Draw("colz"); + chain->GetHistogram()->Write(); + // + chain->Draw("1000*(mean-deltaG-deltaL):sector:abs(theta)",cutS+"theta>0&&type==2&&abs(refX-85)<2","colz"); + chain->GetHistogram()->SetName("TPC_ITSDphi3DLocal3D"); + chain->GetHistogram()->SetTitle("TPC-ITS #Delta_{#phi} - (Local Fit)"); + chain->GetHistogram()->GetYaxis()->SetTitle("#Delta_{#phi} (mrad)"); + chain->GetHistogram()->GetXaxis()->SetTitle("sector"); + chain->GetHistogram()->GetZaxis()->SetTitle("tan(#Theta)"); + chain->GetHistogram()->Draw("colz"); + chain->GetHistogram()->Write(); + // + chain->Draw("10*(mean-deltaG-deltaL):sector:abs(theta)",cutS+"theta>0&&type==0&&abs(refX-85)<2","colz"); + chain->GetHistogram()->SetName("TPC_ITSDRphi3DLocal3D"); + chain->GetHistogram()->SetTitle("TPC-ITS #Delta_{r#phi} - (Local Fit)"); + chain->GetHistogram()->GetYaxis()->SetTitle("#Delta_{r#phi} (mm)"); + chain->GetHistogram()->GetXaxis()->SetTitle("sector"); + chain->GetHistogram()->GetZaxis()->SetTitle("tan(#Theta)"); + chain->GetHistogram()->Draw("colz"); + chain->GetHistogram()->Write(); + // + chain->Draw("10*(mean-deltaG-deltaL):sector:abs(theta)",cutS+"theta>0&&type==0&&abs(refX-0)<2","colz"); + chain->GetHistogram()->SetName("TPC_VertexDRphiLocal3D"); + chain->GetHistogram()->SetTitle("TPC-Vertex #Delta_{r#phi} - (Local Fit)"); + chain->GetHistogram()->GetYaxis()->SetTitle("#Delta_{r#phi} (mm)"); + chain->GetHistogram()->GetXaxis()->SetTitle("sector"); + chain->GetHistogram()->GetZaxis()->SetTitle("tan(#Theta)"); + chain->GetHistogram()->Draw("colz"); + chain->GetHistogram()->Write(); +} + + + + +void FitAlignCombined0(){ + // + // Fit Global X and globalY shift at vertex and at ITS + // + + TTreeSRedirector *pcstream= new TTreeSRedirector("fitAlignCombined.root"); + + TStatToolkit toolkit; + Double_t chi2=0; + Int_t npoints=0; + // Fit vectors + // global fits + TVectorD vecSec(37); + TVectorD paramYVGlobal; + TVectorD paramYITSGlobal; + TVectorD paramPhiVGlobal; + TVectorD paramPhiITSGlobal; + // local fits + TVectorD paramYVLocal; + TVectorD paramPhiVLocal; + TVectorD paramYITSLocal; + TVectorD paramPhiITSLocal; + TMatrixD covar; + // make Aliases + {for (Int_t isec=0; isec<18; isec++){ // sectors + tree->SetAlias(Form("sec%d",isec), Form("(abs(sector-%f)<0.5)",isec+0.5)); + }} + for (Int_t isec=-1;isec<36; isec++){ + vecSec[isec+1]=isec; + } + // + TString fstringGlobal=""; + fstringGlobal+="side++"; + fstringGlobal+="theta++"; + fstringGlobal+="cos(phi)*(theta>0)++"; // GX - A side + fstringGlobal+="sin(phi)*(theta>0)++"; // GY - A side + fstringGlobal+="cos(phi)*(theta<0)++"; // GX - C side + fstringGlobal+="sin(phi)*(theta<0)++"; // GY - C side + fstringGlobal+="cos(phi)*(theta)++"; // theta - get rid of rotation + fstringGlobal+="sin(phi)*(theta)++"; // theta - get rid of rotation + // + // Local Fit + // + TString fstringLocal=""; + {for (Int_t sector=0; sector<18; sector++){ + fstringLocal+=Form("((sec%d)*theta>0)++",sector); + fstringLocal+=Form("((sec%d)*theta<0)++",sector); + }} + // + // Make global fits + // + TString *strFitYVGlobal = TStatToolkit::FitPlane(tree,"V.mean", fstringGlobal.Data(),cutS+"", chi2,npoints,paramYVGlobal,covar,-1,0, 10000000, kFALSE); + printf("chi2=%f\n",TMath::Sqrt(chi2/npoints)); + tree->SetAlias("fitYVGlobal",strFitYVGlobal->Data()); + tree->Draw("V.mean:fitYVGlobal",cutS); + // + TString *strFitPhiVGlobal = TStatToolkit::FitPlane(tree,"VPhi.mean", fstringGlobal.Data(),cutS+"", chi2,npoints,paramPhiVGlobal,covar,-1,0, 10000000, kFALSE); + printf("chi2=%f\n",TMath::Sqrt(chi2/npoints)); + tree->SetAlias("fitPhiVGlobal",strFitPhiVGlobal->Data()); + tree->Draw("VPhi.mean:fitPhiVGlobal",cutS); + // + TString *strFitYITSGlobal = TStatToolkit::FitPlane(tree,"mean", fstringGlobal.Data(),cutS+"", chi2,npoints,paramYITSGlobal,covar,-1,0, 10000000, kFALSE); + printf("chi2=%f\n",TMath::Sqrt(chi2/npoints)); + tree->SetAlias("fitYITSGlobal",strFitYITSGlobal->Data()); + tree->Draw("mean:fitYITSGlobal",cutS); + TString *strFitPhiITSGlobal = TStatToolkit::FitPlane(tree,"Phi.mean", fstringGlobal.Data(),cutS+"", chi2,npoints,paramPhiITSGlobal,covar,-1,0, 10000000, kFALSE); + printf("chi2=%f\n",TMath::Sqrt(chi2/npoints)); + tree->SetAlias("fitPhiITSGlobal",strFitPhiITSGlobal->Data()); + tree->Draw("Phi.mean:fitPhiITSGlobal",cutS); + // + // Residuals to the global fit + // + tree->SetAlias("dyVG", "V.mean-fitYVGlobal"); + tree->SetAlias("dphiVG","(VPhi.mean-fitPhiVGlobal)"); + tree->SetAlias("dyITSG","mean-fitYITSGlobal"); + tree->SetAlias("dphiITSG","(Phi.mean-fitPhiITSGlobal)"); + + TString *strFitYVLocal = TStatToolkit::FitPlane(tree,"dyVG", fstringLocal.Data(),cutS+"", chi2,npoints,paramYVLocal,covar,-1,0, 10000000, kFALSE); + printf("chi2=%f\n",TMath::Sqrt(chi2/npoints)); + tree->SetAlias("fitYVLocal",strFitYVLocal->Data()); + tree->Draw("dyVG-fitYVLocal:sector:abs(theta)",cutS+"theta>0","colz"); + // + TString *strFitPhiVLocal = TStatToolkit::FitPlane(tree,"dphiVG", fstringLocal.Data(),cutS+"", chi2,npoints,paramPhiVLocal,covar,-1,0, 10000000, kFALSE); + printf("chi2=%f\n",TMath::Sqrt(chi2/npoints)); + tree->SetAlias("fitPhiVLocal",strFitPhiVLocal->Data()); + tree->Draw("dphiVG-fitPhiVLocal:sector:abs(theta)",cutS+"theta>0","colz"); + // + // + // + TString *strFitYITSLocal = TStatToolkit::FitPlane(tree,"dyITSG", fstringLocal.Data(),cutS+"", chi2,npoints,paramYITSLocal,covar,-1,0, 10000000, kFALSE); + printf("chi2=%f\n",TMath::Sqrt(chi2/npoints)); + tree->SetAlias("fitYITSLocal",strFitYITSLocal->Data()); + tree->Draw("dyITSG-fitYITSLocal:sector:abs(theta)",cutS+"theta>0","colz"); + // + TString *strFitPhiITSLocal = TStatToolkit::FitPlane(tree,"dphiITSG", fstringLocal.Data(),cutS+"", chi2,npoints,paramPhiITSLocal,covar,-1,0, 10000000, kFALSE); + printf("chi2=%f\n",TMath::Sqrt(chi2/npoints)); + tree->SetAlias("fitPhiITSLocal",strFitPhiITSLocal->Data()); + tree->Draw("dphiITSG-fitPhiITSLocal:sector:abs(theta)",cutS+"theta>0","colz"); + + // + TH1D *hisA = 0; + TH1D *hisC = 0; + TVectorD vSec(36); + TVectorD vecDPhi(36); + TVectorD vecDLY(36); + TVectorD vecDGX(36); + TVectorD vecDGY(36); + + // + tree->SetAlias("phiMean","(fitPhiVLocal+fitPhiITSLocal+fitPhiVGlobal+fitPhiITSGlobal)*0.5"); + tree->SetAlias("yMeanV","((fitYITSLocal+fitYITSGlobal-85*phiMean)+(fitYVLocal+fitYVGlobal))*0.5"); + tree->SetAlias("yMeanITS","((fitYITSLocal+fitYITSGlobal)+(fitYVLocal+fitYVGlobal+85*phiMean))*0.5"); + + + tree->Draw("phiMean:int(sector)>>hhhA(18,0.18)",cutS+"abs(theta)<0.15&&theta>0","prof"); + hisA = (TH1D*) tree->GetHistogram()->Clone(); + tree->Draw("phiMean:int(sector)>>hhhC(18,0.18)",cutS+"abs(theta)<0.15&&theta<0","prof"); + hisC = (TH1D*) tree->GetHistogram()->Clone(); + + for (Int_t isec=0; isec<36; isec++){ + vSec[isec]=isec; + if (isec<18) vecDPhi[isec]=hisA->GetBinContent(isec%18+1); + if (isec>=18) vecDPhi[isec]=hisC->GetBinContent(isec%18+1); + } + tree->Draw("yMeanV:int(sector)>>hhhA(18,0.18)",cutS+"abs(theta)<0.15&&theta>0","prof"); + hisA = (TH1D*) tree->GetHistogram()->Clone(); + tree->Draw("yMeanV:int(sector)>>hhhC(18,0.18)",cutS+"abs(theta)<0.15&&theta<0","prof"); + hisC = (TH1D*) tree->GetHistogram()->Clone(); + // + for (Int_t isec=0; isec<36; isec++){ + if (isec<18) vecDLY[isec]=hisA->GetBinContent(isec%18+1); + if (isec>=18) vecDLY[isec]=hisC->GetBinContent(isec%18+1); + vecDGX[isec]=-vecDLY[isec]*TMath::Sin(TMath::Pi()*(isec+0.5)/9.); + vecDGY[isec]= vecDLY[isec]*TMath::Cos(TMath::Pi()*(isec+0.5)/9.); + } + + + + + // + // Store results + // + {(*pcstream)<<"fitLocal"<< + //global fits + "sec.="<<&vSec<< + "pYVGlobal.="<<¶mYVGlobal<< + "pYITSGlobal.="<<¶mYITSGlobal<< + "pPhiVGlobal.="<<¶mPhiVGlobal<< + "pPhiITSGlobal.="<<¶mPhiITSGlobal<< + // local fits + "pYVLocal.="<<¶mYVLocal<< + "pPhiVLocal.="<<¶mPhiVLocal<< + "pYITSLocal.="<<¶mYITSLocal<< + "pPhiITSLocal.="<<¶mPhiITSLocal<< + // + // combined + "vDPhi.="<<&vecDPhi<< + "vDLY.="<<&vecDLY<< + "vDGX.="<<&vecDGX<< + "vDGY.="<<&vecDGY<< + "\n"; + } + paramYVGlobal.Write("paramYVGlobal"); + paramYITSGlobal.Write("paramYITSGlobal"); + paramPhiVGlobal.Write("paramPhiVGlobal"); + paramPhiITSGlobal.Write("paramPhiITSGlobal"); + // local fits + paramYVLocal.Write("paramYVLocal"); + paramPhiVLocal.Write("paramPhiVLocal"); + paramYITSLocal.Write("paramYITSLocal"); + paramPhiITSLocal.Write("paramPhiITSLocal"); + vecDPhi.Write("vecDPhi"); + vecDLY.Write("vecDLY"); + vecDGX.Write("vecDGX"); + vecDGY.Write("vecDGY"); + + + + tree->Draw("dyVG:sector:abs(theta)",cutS+"theta>0","colz"); + tree->GetHistogram()->Write("fitYVGlobal"); + tree->Draw("dphiVG:sector:abs(theta)",cutS+"theta>0","colz"); + tree->GetHistogram()->Write("fitPhiVGlobal"); + tree->Draw("dyITSG:sector:abs(theta)",cutS+"theta>0","colz"); + tree->GetHistogram()->Write("fitYITSGlobal"); + tree->Draw("dphiITSG:sector:abs(theta)",cutS+"theta>0","colz"); + tree->GetHistogram()->Write("fitPhiITSGlobal"); + // + tree->Draw("dyVG-fitYVLocal:sector:abs(theta)",cutS+"theta>0","colz"); + tree->GetHistogram()->Write("fitYVLocal"); + tree->Draw("dphiVG-fitPhiVLocal:sector:abs(theta)",cutS+"theta>0","colz"); + tree->GetHistogram()->Write("fitPhiVLocal"); + tree->Draw("dyITSG-fitYITSLocal:sector:abs(theta)",cutS+"theta>0","colz"); + tree->GetHistogram()->Write("fitYITSLocal"); + tree->Draw("dphiITSG-fitPhiITSLocal:sector:abs(theta)",cutS+"theta>0","colz"); + tree->GetHistogram()->Write("fitPhiITSLocal"); + delete pcstream; +} + +void FitAlignCombined(){ + // + // + // make sector alignment - using Kalman filter method -AliTPCkalmanAlign + // Combined information is used, mean residuals are minimized: + // + // 1. TPC-TPC sector alignment + // 2. TPC-ITS alignment + // 3. TPC vertex alignment + TFile fcalib("../mergeField0/TPCAlignObjects.root"); + AliTPCcalibAlign * align = ( AliTPCcalibAlign *)fcalib.Get("alignTPC"); + + TCut cutQ="entries>1000&&abs(theta)<0.5&&abs(snp)<0.2&&YP.entries>50&&YM.entries>50"; + TH1F his1("hdeltaY1","hdeltaY1",100,-0.5,0.5); + TMatrixD vecAlign(72,1); + TMatrixD covAlign(72,72); + TMatrixD vecAlignY(72,1); + TMatrixD covAlignY(72,72); + TMatrixD vecAlignTheta(72,1); + TMatrixD covAlignTheta(72,72); + TMatrixD vecAlignZ(72,1); + TMatrixD covAlignZ(72,72); + AliTPCkalmanAlign::BookAlign1D(vecAlign,covAlign,0,0.002); + AliTPCkalmanAlign::BookAlign1D(vecAlignY,covAlignY,0,0.1); + AliTPCkalmanAlign::BookAlign1D(vecAlignTheta,covAlignTheta,0,0.002); + AliTPCkalmanAlign::BookAlign1D(vecAlignZ,covAlignZ,0,0.1); + TVectorD vecITSY(72); + TVectorD vecITSYPM(72); + TVectorD vecITSPhi(72); + TVectorD vecITSPhiPM(72); + TVectorD vecVY(72); + TVectorD vecVS(72); + TVectorD vecITSTheta(72); + TVectorD vecVTheta(72); + {for (Int_t isec0=0; isec0<36; isec0++){ + Double_t phi0=(isec0%18+0.5)*TMath::Pi()/9.; + if (phi0>TMath::Pi()) phi0-=TMath::TwoPi(); + Int_t iside0=(isec0%36<18)? 0:1; + TCut cutSector=Form("abs(%f-phi)<0.14",phi0); + TCut cutSide = (iside0==0)? "theta>0":"theta<0"; + // + itsdy->Draw("mean",cutQ+cutSector+cutSide); + Double_t meanITSY=itsdy->GetHistogram()->GetMean(); + vecITSY[isec0]=meanITSY; + vecITSY[isec0+36]=meanITSY; + // + itsdy->Draw("(YP.mean+YM.mean)*0.5",cutQ+cutSector+cutSide); + Double_t meanITSYPM=itsdy->GetHistogram()->GetMean(); + vecITSYPM[isec0]=meanITSYPM; + vecITSYPM[isec0+36]=meanITSYPM; + // + itsdy->Draw("Phi.mean",cutQ+cutSector+cutSide); + Double_t meanITSPhi=itsdy->GetHistogram()->GetMean(); + vecITSPhi[isec0]=meanITSPhi; + vecITSPhi[isec0+36]=meanITSPhi; + // + itsdy->Draw("(PhiP.mean+PhiM.mean)*0.5",cutQ+cutSector+cutSide); + Double_t meanITSPhiPM=itsdy->GetHistogram()->GetMean(); + vecITSPhiPM[isec0]=meanITSPhiPM; + vecITSPhiPM[isec0+36]=meanITSPhiPM; + // + // + itsdy->Draw("VPhi.mean",cutQ+cutSector+cutSide); + Double_t meanVS=itsdy->GetHistogram()->GetMean(); + vecVS[isec0]=meanVS; + vecVS[isec0+36]=meanVS; + // + itsdy->Draw("V.mean",cutQ+cutSector+cutSide); + Double_t meanVY=itsdy->GetHistogram()->GetMean(); + vecVY[isec0]=meanVY; + vecVY[isec0+36]=meanVY; + // + itsdy->Draw("T.mean",cutQ+cutSector+cutSide); + Double_t meanITST=itsdy->GetHistogram()->GetMean(); + vecITSTheta[isec0]=meanITST; + vecITSTheta[isec0+36]=meanITST; + // + itsdy->Draw("VT.mean",cutQ+cutSector+cutSide); + Double_t meanVT=itsdy->GetHistogram()->GetMean(); + vecVTheta[isec0]=meanVT; + vecVTheta[isec0+36]=meanVT; + } + } + AliTPCROC * roc = AliTPCROC::Instance(); + Double_t fXIROC = (roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(0,roc->GetNRows(0)-1))*0.5; + Double_t fXOROC = (roc->GetPadRowRadii(36,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5; + Double_t fXIO = ( roc->GetPadRowRadii(0,roc->GetNRows(0)-1)+roc->GetPadRowRadii(36,0))*0.5; + + TTreeSRedirector *pcstream=new TTreeSRedirector("combAlign.root"); + + { + for (Int_t iter=0; iter<5; iter++){ + for (Int_t isec0=0; isec0<72; isec0++){ + for (Int_t isec1=0; isec1<72; isec1++){ + //if (isec0%36!=isec0%36) continue; + if (iter==0 && isec0%36!=isec1%36) continue; + if (iter==1 && isec0%18!=isec1%18) continue; + TH1 * his = align->GetHisto(AliTPCcalibAlign::kY,isec0,isec1); + TH1 * hisPhi = align->GetHisto(AliTPCcalibAlign::kPhi,isec0,isec1); + TH1 * hisTheta = align->GetHisto(AliTPCcalibAlign::kTheta,isec0,isec1); + TH1 * hisZ = align->GetHisto(AliTPCcalibAlign::kZ,isec0,isec1); + Int_t side0=(isec0%36)<18 ? 1:-1; + Int_t side1=(isec1%36)<18 ? 1:-1; + Double_t weightTPC= 0.002; + if (isec0%18==isec1%18) weightTPC=0.0005; + if (isec0%36==isec1%36) weightTPC=0.00005; + if (side0!=side1) weightTPC*=2.; + Double_t weightTPCT= 0.001; + if (isec0%18==isec1%18) weightTPC=0.0005; + if (isec0%36==isec1%36) weightTPC=0.00005; + if (TMath::Abs(isec0%18-isec1%18)==1) weightTPC = 0.0005; + if (TMath::Abs(isec0%36-isec1%36)==1) weightTPC = 0.0001; + + // + Double_t meanITS0=vecITSY[isec0]; + Double_t meanITSPM0=vecITSYPM[isec0]; + Double_t meanITSPhi0=vecITSPhi[isec0]; + Double_t meanITSPhiPM0=vecITSPhiPM[isec0]; + Double_t meanVY0=vecVY[isec0]; + Double_t meanVPhi0=vecVS[isec0]; + Double_t meanITSTheta0=vecITSTheta[isec0]; + Double_t meanVTheta0=vecVTheta[isec0]; + // + Double_t meanITS1=vecITSY[isec1]; + Double_t meanITSPM1=vecITSYPM[isec1]; + Double_t meanITSPhi1=vecITSPhi[isec1]; + Double_t meanITSPhiPM1=vecITSPhiPM[isec1]; + Double_t meanVY1=vecVY[isec1]; + Double_t meanVPhi1=vecVS[isec1]; + Double_t meanITSTheta1=vecITSTheta[isec1]; + Double_t meanVTheta1=vecVTheta[isec1]; + // + Double_t meanPhi0 = (meanITSPhi0+meanVPhi0+2*meanITS0/83.)/4.; + Double_t meanPhi1 = (meanITSPhi1+meanVPhi1+2*meanITS1/83.)/4.; + // + // + if (iter>2 &&isec0==isec1){ + AliTPCkalmanAlign::UpdateAlign1D(-meanITS0/83, 0.001, isec0, vecAlign,covAlign); + AliTPCkalmanAlign::UpdateAlign1D(-meanITSPhi0, 0.001, isec0, vecAlign,covAlign); + AliTPCkalmanAlign::UpdateAlign1D(-meanVPhi0, 0.001, isec0, vecAlign,covAlign); + // + AliTPCkalmanAlign::UpdateAlign1D(-meanPhi0, 0.001, isec0, vecAlign,covAlign); + AliTPCkalmanAlign::UpdateAlign1D(-meanVTheta0, 0.001, isec0, vecAlignTheta,covAlignTheta); + AliTPCkalmanAlign::UpdateAlign1D(-meanITSTheta0, 0.001, isec0, vecAlignTheta,covAlignTheta); + } + if (iter>2){ + AliTPCkalmanAlign::UpdateAlign1D(meanPhi1-meanPhi0, 0.001, isec0,isec1, vecAlign,covAlign); + AliTPCkalmanAlign::UpdateAlign1D(meanITSPhiPM1-meanITSPhiPM0, 0.001, isec0,isec1, vecAlign,covAlign); + AliTPCkalmanAlign::UpdateAlign1D((meanITSPM1-meanITSPM0)/83., 0.001, isec0,isec1, vecAlign,covAlign); + AliTPCkalmanAlign::UpdateAlign1D(meanVTheta1-meanVTheta0, 0.001, isec0,isec1, vecAlignTheta,covAlignTheta); + AliTPCkalmanAlign::UpdateAlign1D(meanITSTheta1-meanITSTheta0, 0.001, isec0,isec1, vecAlignTheta,covAlignTheta); + } + // + if (!his) continue; + if (!hisPhi) continue; + if (!hisTheta) continue; + if (his->GetEntries()<100) continue; + Double_t xref=fXIO; + if (isec0<36&&isec1<36) xref=fXIROC; + if (isec0>=36&&isec1>=36) xref=fXOROC; + Double_t meanTPC=his->GetMean(); + Double_t meanTPCPhi=hisPhi->GetMean(); + Double_t meanTPCTheta=hisTheta->GetMean(); + Double_t meanTPCZ=hisZ->GetMean(); + // + // + Double_t kalman0= vecAlign(isec0,0); + Double_t kalman1= vecAlign(isec1,0); + Double_t kalmanY0= vecAlignY(isec0,0); + Double_t kalmanY1= vecAlignY(isec1,0); + Double_t kalmanTheta0= vecAlignTheta(isec0,0); + Double_t kalmanTheta1= vecAlignTheta(isec1,0); + Double_t kalmanZ0= vecAlignZ(isec0,0); + Double_t kalmanZ1= vecAlignZ(isec1,0); + // + // + AliTPCkalmanAlign::UpdateAlign1D((meanTPC)/xref, weightTPC, isec0,isec1, vecAlign,covAlign); + AliTPCkalmanAlign::UpdateAlign1D(meanTPCPhi, weightTPC*10,isec0,isec1, vecAlign,covAlign); + // + AliTPCkalmanAlign::UpdateAlign1D(meanTPCTheta, weightTPCT, isec0,isec1, vecAlignTheta,covAlignTheta); + if (side0==side1) AliTPCkalmanAlign::UpdateAlign1D(meanTPCZ, weightTPCT*100., isec0,isec1, vecAlignZ,covAlignZ); + //printf("isec0\t%d\tisec1\t%d\t%f\t%f\t%f\n",isec0,isec1, meanTPC, meanITS0,meanITS1); + + if (iter>=0) (*pcstream)<<"align"<< + "iter="<GetFile()->cd(); + vecAlign.Write("alignPhiMean"); + covAlign.Write("alingPhiCovar"); + vecAlignTheta.Write("alignThetaMean"); + covAlignTheta.Write("alingThetaCovar"); + vecAlignZ.Write("alignZMean"); + covAlignZ.Write("alingZCovar"); + delete pcstream; + TFile f("combAlign.root"); + TTree * treeA = (TTree*)f.Get("align"); + treeA->SetMarkerStyle(25); + treeA->SetMarkerSize(0.5); +} + + + + + +void UpdateOCDBAlign(){ + // + // Store resulting OCDB entry + // 0. Setup OCDB to get necccessary old entries - not done here + // 1. Get old OCDB entry + // 2. Get delta alignment + // 3. Add delta alignment + // 4. Store new alignment in + + AliCDBEntry * entry = AliCDBManager::Instance()->Get("TPC/Align/Data"); + TClonesArray * array = (TClonesArray*)entry->GetObject(); + Int_t entries = array->GetEntries(); + TClonesArray * arrayNew=(TClonesArray*)array->Clone(); + TFile f("combAlign.root"); + TMatrixD *matPhi=(TMatrixD*)f.Get("alignPhiMean"); + // + // + { for (Int_t i=0;iUncheckedAt(i); + AliAlignObjParams *alignPNew = (AliAlignObjParams*)arrayNew->UncheckedAt(i); + Int_t imod; + AliGeomManager::ELayerID ilayer; + alignP->GetVolUID(ilayer, imod); + if (ilayer==AliGeomManager::kInvalidLayer) continue; + Int_t sector=imod; + if (ilayer==AliGeomManager::kTPC2) sector+=36; + Double_t transOld[3], rotOld[3]; + alignP->GetTranslation(transOld); // in cm + alignP->GetAngles(rotOld); // in degrees + printf("%d\t%d\t%d\t",ilayer, imod,sector); + printf("%f mrad \t%f mrad\n",1000*rotOld[2]*TMath::DegToRad(), 1000*((*matPhi)(sector,0))); + alignPNew->SetRotation(rotOld[0],rotOld[1], rotOld[2]+(*matPhi)(sector,0)*TMath::RadToDeg()); + alignPNew->Print("kokot"); + alignP->Print("kokot"); + } + } + + + + TString ocdbStorage="local://"+gSystem->GetFromPipe("pwd")+"/OCDBUpdate"; + TString userName=gSystem->GetFromPipe("echo $USER"); + TString date=gSystem->GetFromPipe("date"); + + AliCDBMetaData *metaData= new AliCDBMetaData(); + metaData->SetObjectClassName("TObjArray"); + metaData->SetResponsible("Marian Ivanov"); + metaData->SetBeamPeriod(1); + metaData->SetAliRootVersion("05-26-02"); //root version + metaData->SetComment(Form("Correction calibration. User: %s\n Data%s",userName.Data(),date.Data())); + + AliCDBId* id1=NULL; + id1=new AliCDBId("TPC/Align/Data", 0, AliCDBRunRange::Infinity()); + AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(ocdbStorage); + gStorage->Put(arrayNew, (*id1), metaData); + + TFile falign("falign.root","recreate"); + arrayNew->Write("new"); + array->Write("old"); + falign.Close(); +} + + + +void UpdateOCDBAlign0(){ + // + // Store resulting OCDB entry + // 0. Setup OCDB to get necccessary old entries - not done here + // 1. Get old OCDB entry + // 2. Get delta alignment + // 3. Add delta alignment + // 4. Store new alignment in + + AliCDBEntry * entry = AliCDBManager::Instance()->Get("TPC/Align/Data"); + TClonesArray * array = (TClonesArray*)entry->GetObject(); + Int_t entries = array->GetEntries(); + TClonesArray * arrayNew=(TClonesArray*)array->Clone(); + TFile f("fitAlignCombined.root"); + TVectorD *vecDPhi = (TVectorD*) f.Get("vecDPhi"); + TVectorD *vecDGX = (TVectorD*) f.Get("vecDGX"); + TVectorD *vecDGY = (TVectorD*) f.Get("vecDGY"); + + // + // + { for (Int_t i=0;iUncheckedAt(i); + AliAlignObjParams *alignPNew = (AliAlignObjParams*)arrayNew->UncheckedAt(i); + Int_t imod; + AliGeomManager::ELayerID ilayer; + alignP->GetVolUID(ilayer, imod); + if (ilayer==AliGeomManager::kInvalidLayer) continue; + Int_t sector=imod; + if (ilayer==AliGeomManager::kTPC2) sector+=36; + Double_t transOld[3], rotOld[3]; + alignP->GetTranslation(transOld); // in cm + alignP->GetAngles(rotOld); // in degrees + printf("%d\t%d\t%d\t",ilayer, imod,sector); + printf("%f mrad \t%f mrad\n",1000*rotOld[2]*TMath::DegToRad(), 1000*((*vecDPhi)[sector%36])); + alignPNew->SetRotation(rotOld[0],rotOld[1], rotOld[2]-(*vecDPhi)[sector%36]*TMath::RadToDeg()); + alignPNew->SetTranslation(transOld[0]-(*vecDGX)[sector%36],transOld[1]-(*vecDGY)[sector%36], transOld[2]); + alignPNew->Print("kokot"); + alignP->Print("kokot"); + } + } + + + + TString ocdbStorage="local://"+gSystem->GetFromPipe("pwd")+"/OCDBUpdate"; + TString userName=gSystem->GetFromPipe("echo $USER"); + TString date=gSystem->GetFromPipe("date"); + + AliCDBMetaData *metaData= new AliCDBMetaData(); + metaData->SetObjectClassName("TObjArray"); + metaData->SetResponsible("Marian Ivanov"); + metaData->SetBeamPeriod(1); + metaData->SetAliRootVersion("05-26-02"); //root version + metaData->SetComment(Form("Correction calibration. User: %s\n Data%s",userName.Data(),date.Data())); + + AliCDBId* id1=NULL; + id1=new AliCDBId("TPC/Align/Data", 0, AliCDBRunRange::Infinity()); + AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(ocdbStorage); + gStorage->Put(arrayNew, (*id1), metaData); + + TFile falign("falign.root","recreate"); + arrayNew->Write("new"); + array->Write("old"); + falign.Close(); +} -- 2.43.0