X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TPC%2FAliTPCPreprocessorOffline.cxx;h=c79ab8017ed5a4b07398877d77ad254bea79a1a5;hb=849102f9ff0bac3c810a30643290d9db9cf7f045;hp=ec859cdfa073632fc6bd5903a44a54ea458fe595;hpb=68d461b2b82187a73ffb41ee8ac5f73e17b21875;p=u%2Fmrichter%2FAliRoot.git diff --git a/TPC/AliTPCPreprocessorOffline.cxx b/TPC/AliTPCPreprocessorOffline.cxx index ec859cdfa07..c79ab8017ed 100644 --- a/TPC/AliTPCPreprocessorOffline.cxx +++ b/TPC/AliTPCPreprocessorOffline.cxx @@ -36,8 +36,10 @@ // default storage ""- data stored at current working directory e.g. - AliTPCPreprocessorOffline process; - process.CalibTimeGain("CalibObjects.root",114000,121040,0); + gSystem->Load("libANALYSIS"); + gSystem->Load("libTPCcalib"); + AliTPCPreprocessorOffline proces; + proces.CalibTimeGain("TPCMultObjects.root",114000,140040,0); TFile oo("OCDB/TPC/Calib/TimeGain/Run114000_121040_v0_s0.root") TObjArray * arr = AliCDBEntry->GetObject() arr->At(4)->Draw("alp") @@ -72,7 +74,15 @@ #include "AliRelAlignerKalman.h" #include "AliTPCParamSR.h" #include "AliTPCcalibTimeGain.h" +#include "AliTPCcalibGainMult.h" #include "AliSplineFit.h" +#include "AliTPCComposedCorrection.h" +#include "AliTPCExBTwist.h" +#include "AliTPCCalibGlobalMisalignment.h" +#include "TStatToolkit.h" +#include "TChain.h" +#include "TCut.h" +#include "AliTrackerBase.h" #include "AliTPCPreprocessorOffline.h" @@ -96,7 +106,12 @@ AliTPCPreprocessorOffline::AliTPCPreprocessorOffline(): fGainArray(new TObjArray), // array to be stored in the OCDB fGainMIP(0), // calibration component for MIP fGainCosmic(0), // calibration component for cosmic - fSwitchOnValidation(kFALSE) // flag to switch on validation of OCDB parameters + fGainMult(0), + fAlignTree(0), // alignment tree + fSwitchOnValidation(kFALSE), // flag to switch on validation of OCDB parameters + fMinGain(2.0), + fMaxGain(3.0), + fMaxVdriftCorr(0.03) { // // default constructor @@ -119,8 +134,8 @@ void AliTPCPreprocessorOffline::GetRunRange(AliTPCcalibTime * const timeDrift){ TObjArray *hisArray =timeDrift->GetHistoDrift(); {for (Int_t i=0; iGetEntriesFast(); i++){ THnSparse* addHist=(THnSparse*)hisArray->UncheckedAt(i); - if (addHist->GetEntries()GetEntries()Projection(3); TH1D* histoTime=addHist->Projection(0); printf("%s\t%f\t%d\t%d\n",histo->GetName(), histo->GetEntries(),histo->FindFirstBinAbove(0),histo->FindLastBinAbove(0)); @@ -203,7 +218,15 @@ void AliTPCPreprocessorOffline::CalibTimeVdrift(const Char_t* file, Int_t ustart Printf("TPC time drift OCDB parameters out of range!"); return; } - + // + //4.b make alignment + // + MakeFitTime(); + TFile * ftime= TFile::Open("fitITSVertex.root"); + if (ftime){ + TObject * alignmentTime=ftime->Get("FitCorrectionTime"); + if (alignmentTime) fVdriftArray->AddLast(alignmentTime); + } // // // 5. update of OCDB @@ -228,15 +251,21 @@ void AliTPCPreprocessorOffline::UpdateOCDBDrift( Int_t ustartRun, Int_t uendRun, gStorage->Put(fVdriftArray, (*id1), metaData); } -Bool_t AliTPCPreprocessorOffline::ValidateTimeGain(Double_t minGain, Double_t maxGain) +Bool_t AliTPCPreprocessorOffline::ValidateTimeGain() { // // Validate time gain corrections // Printf("ValidateTimeGain..." ); + Float_t minGain = fMinGain; + Float_t maxGain = fMaxGain; TGraphErrors *gr = (TGraphErrors*)fGainArray->FindObject("TGRAPHERRORS_MEAN_GAIN_BEAM_ALL"); - if(!gr) return kFALSE; + if (!gr) { + gr = (TGraphErrors*)fGainArray->FindObject("TGRAPHERRORS_MEAN_GAIN_COSMIC_ALL"); + if (!gr) return kFALSE; + Printf("Assuming given run is a cosmic run. Using gain calibration from Fermi-plateau muons."); + } if(gr->GetN()<1) return kFALSE; // check whether gain in the range @@ -250,20 +279,28 @@ return kTRUE; } -Bool_t AliTPCPreprocessorOffline::ValidateTimeDrift(Double_t maxVDriftCorr) +Bool_t AliTPCPreprocessorOffline::ValidateTimeDrift() { // // Validate time drift velocity corrections // Printf("ValidateTimeDrift..." ); + Float_t maxVDriftCorr = fMaxVdriftCorr; + TGraphErrors* gr = (TGraphErrors*)fVdriftArray->FindObject("ALIGN_ITSB_TPC_DRIFTVD"); + Printf("ALIGN_ITSB_TPC_DRIFTVD graph = %p",gr); + if(!gr) return kFALSE; - if(gr->GetN()<1) return kFALSE; + if(gr->GetN()<1) { + Printf("ALIGN_ITSB_TPC_DRIFTVD number of points = %d",gr->GetN()); + return kFALSE; + } // check whether drift velocity corrections in the range for(Int_t iPoint = 0; iPointGetN(); iPoint++) { + Printf("Y value from the graph: %f",TMath::Abs(gr->GetY()[iPoint])); if(TMath::Abs(gr->GetY()[iPoint]) > maxVDriftCorr) return kFALSE; } @@ -359,7 +396,13 @@ TGraphErrors* AliTPCPreprocessorOffline::FilterGraphMedianAbs(TGraphErrors * gra Double_t *erry=new Double_t[npoints0]; // // - if (npoints01) graphOut= new TGraphErrors(npoints,outx,outy,errx,erry); + delete []outx; + delete []outy; + delete []errx; + delete []erry; return graphOut; } @@ -402,6 +449,10 @@ void AliTPCPreprocessorOffline::AddHistoGraphs( TObjArray * vdriftArray, AliTPC THnSparse* newHist=hist->Projection(4,dim); newHist->SetName(name); TGraphErrors* graph=AliTPCcalibBase::FitSlices(newHist,2,0,400,100,0.05,0.95, kTRUE); + if (!graph) { + printf("Graph =%s filtered out\n", name.Data()); + continue; + } printf("name=%s graph=%i, N=%i\n", name.Data(), graph==0, graph->GetN()); Int_t pos=name.Index("_"); name=name(pos,name.Capacity()-pos); @@ -410,19 +461,17 @@ void AliTPCPreprocessorOffline::AddHistoGraphs( TObjArray * vdriftArray, AliTPC graphName.ToUpper(); // graph = FilterGraphDrift(graph, kErrSigmaCut, kMedianCutAbs); - if (!graph) { - printf("Graph =%s filtered out\n", name.Data()); - continue; - } // - graph->SetMarkerStyle(i%8+20); - graph->SetMarkerColor(i%7); - graph->GetXaxis()->SetTitle("Time"); - graph->GetYaxis()->SetTitle("v_{dcor}"); - graph->SetName(graphName); - graph->SetTitle(graphName); - printf("Graph %d\t=\t%s\n", i, graphName.Data()); - vdriftArray->Add(graph); + if (graph){ + graph->SetMarkerStyle(i%8+20); + graph->SetMarkerColor(i%7); + graph->GetXaxis()->SetTitle("Time"); + graph->GetYaxis()->SetTitle("v_{dcor}"); + graph->SetName(graphName); + graph->SetTitle(graphName); + printf("Graph %d\t=\t%s\n", i, graphName.Data()); + vdriftArray->Add(graph); + } } } } @@ -445,9 +494,9 @@ void AliTPCPreprocessorOffline::AddAlignmentGraphs( TObjArray * vdriftArray, Al arrayTRD=timeDrift->GetAlignTRDTPC(); arrayTOF=timeDrift->GetAlignTOFTPC(); - if (arrayITS->GetEntries()>0) mstatITS= AliTPCcalibDButil::MakeStatRelKalman(arrayITS,0.9,50,0.025); - if (arrayTOF->GetEntries()>0) mstatTOF= AliTPCcalibDButil::MakeStatRelKalman(arrayTOF,0.9,1000,0.025); - if (arrayTRD->GetEntries()>0) mstatTRD= AliTPCcalibDButil::MakeStatRelKalman(arrayTRD,0.9,50,0.025); + if (arrayITS->GetEntries()>0) mstatITS= AliTPCcalibDButil::MakeStatRelKalman(arrayITS,0.7,50,fMaxVdriftCorr); + if (arrayTOF->GetEntries()>0) mstatTOF= AliTPCcalibDButil::MakeStatRelKalman(arrayTOF,0.7,1000,fMaxVdriftCorr); + if (arrayTRD->GetEntries()>0) mstatTRD= AliTPCcalibDButil::MakeStatRelKalman(arrayTRD,0.7,50,fMaxVdriftCorr); // TObjArray * arrayITSP= AliTPCcalibDButil::SmoothRelKalman(arrayITS,*mstatITS, 0, 5.); TObjArray * arrayITSM= AliTPCcalibDButil::SmoothRelKalman(arrayITS,*mstatITS, 1, 5.); @@ -518,9 +567,9 @@ void AliTPCPreprocessorOffline::AddLaserGraphs( TObjArray * vdriftArray, AliTPC // add graphs for laser // const Double_t delayL0L1 = 0.071; //this is hack for 1/2 weeks - THnSparse *hisN=0; + //THnSparse *hisN=0; TGraphErrors *grLaser[6]={0,0,0,0,0,0}; - hisN = timeDrift->GetHistVdriftLaserA(0); + //hisN = timeDrift->GetHistVdriftLaserA(0); if (timeDrift->GetHistVdriftLaserA(0)){ grLaser[0]=MakeGraphFilter0(timeDrift->GetHistVdriftLaserA(0),0,2,5,delayL0L1); grLaser[0]->SetName("GRAPH_MEAN_DELAY_LASER_ALL_A"); @@ -703,7 +752,7 @@ void AliTPCPreprocessorOffline::MakeDefaultPlots(TObjArray * const arr, TObjArra //picArray->AddLast(pad); } - if (itstpcP&&itstpcM){ + if (itstpcP&&itstpcM&&itstpcB){ pad = new TCanvas("ITSTPC","ITSTPC"); itstpcP->Draw("alp"); SetPadStyle(pad,mx0,mx1,my0,my1); @@ -720,7 +769,7 @@ void AliTPCPreprocessorOffline::MakeDefaultPlots(TObjArray * const arr, TObjArra //picArray->AddLast(pad); } - if (itstpcB&&laserA){ + if (itstpcB&&laserA&&itstpcP&&itstpcM){ pad = new TCanvas("ITSTPC_LASER","ITSTPC_LASER"); SetPadStyle(pad,mx0,mx1,my0,my1); laserA->Draw("alp"); @@ -787,7 +836,8 @@ void AliTPCPreprocessorOffline::CalibTimeGain(const Char_t* fileName, Int_t star // AnalyzeGain(startRunNumber,endRunNumber, 1000,1.43); AnalyzeAttachment(startRunNumber,endRunNumber); - + AnalyzePadRegionGain(); + AnalyzeGainMultiplicity(); // // 3. Make control plots // @@ -816,9 +866,15 @@ void AliTPCPreprocessorOffline::ReadGainGlobal(const Char_t* fileName){ if (array){ fGainMIP = ( AliTPCcalibTimeGain *)array->FindObject("calibTimeGain"); fGainCosmic = ( AliTPCcalibTimeGain *)array->FindObject("calibTimeGainCosmic"); + fGainMult = ( AliTPCcalibGainMult *)array->FindObject("calibGainMult"); }else{ fGainMIP = ( AliTPCcalibTimeGain *)fcalib.Get("calibTimeGain"); fGainCosmic = ( AliTPCcalibTimeGain *)fcalib.Get("calibTimeGainCosmic"); + fGainMult = ( AliTPCcalibGainMult *)fcalib.Get("calibGainMult"); + } + if (!fGainMult){ + TFile fcalibMult("TPCMultObjects.root"); + fGainMult = ( AliTPCcalibGainMult *)fcalibMult.Get("calibGainMult"); } TH1 * hisT=0; Int_t firstBinA =0, lastBinA=0; @@ -968,6 +1024,126 @@ Bool_t AliTPCPreprocessorOffline::AnalyzeAttachment(Int_t startRunNumber, Int_t } +Bool_t AliTPCPreprocessorOffline::AnalyzePadRegionGain(){ + // + // Analyze gain for different pad regions - produce the calibration graphs 0,1,2 + // + if (fGainMult) + { + TH2D * histQmax = (TH2D*) fGainMult->GetHistPadEqual()->Projection(0,2); + TH2D * histQtot = (TH2D*) fGainMult->GetHistPadEqual()->Projection(1,2); + // + TObjArray arr; + histQmax->FitSlicesY(0,0,-1,0,"QNR",&arr); + Double_t xMax[3] = {0,1,2}; + Double_t yMax[3] = {((TH1D*)arr.At(1))->GetBinContent(1), + ((TH1D*)arr.At(1))->GetBinContent(2), + ((TH1D*)arr.At(1))->GetBinContent(3)}; + Double_t yMaxErr[3] = {((TH1D*)arr.At(1))->GetBinError(1), + ((TH1D*)arr.At(1))->GetBinError(2), + ((TH1D*)arr.At(1))->GetBinError(3)}; + TGraphErrors * fitPadRegionQmax = new TGraphErrors(3, xMax, yMax, 0, yMaxErr); + // + histQtot->FitSlicesY(0,0,-1,0,"QNR",&arr); + Double_t xTot[3] = {0,1,2}; + Double_t yTot[3] = {((TH1D*)arr.At(1))->GetBinContent(1), + ((TH1D*)arr.At(1))->GetBinContent(2), + ((TH1D*)arr.At(1))->GetBinContent(3)}; + Double_t yTotErr[3] = {((TH1D*)arr.At(1))->GetBinError(1), + ((TH1D*)arr.At(1))->GetBinError(2), + ((TH1D*)arr.At(1))->GetBinError(3)}; + TGraphErrors * fitPadRegionQtot = new TGraphErrors(3, xTot, yTot, 0, yTotErr); + // + fitPadRegionQtot->SetName("TGRAPHERRORS_MEANQTOT_PADREGIONGAIN_BEAM_ALL");// set proper names according to naming convention + fitPadRegionQmax->SetName("TGRAPHERRORS_MEANQMAX_PADREGIONGAIN_BEAM_ALL");// set proper names according to naming convention + // + fGainArray->AddLast(fitPadRegionQtot); + fGainArray->AddLast(fitPadRegionQmax); + return kTRUE; + } + return kFALSE; + +} + + +Bool_t AliTPCPreprocessorOffline::AnalyzeGainMultiplicity() { + // + // Analyze gain as a function of multiplicity and produce calibration graphs + // + if (!fGainMult) return kFALSE; + fGainMult->GetHistGainMult()->GetAxis(3)->SetRangeUser(3,3); + TH2D * histMultMax = fGainMult->GetHistGainMult()->Projection(0,4); + TH2D * histMultTot = fGainMult->GetHistGainMult()->Projection(1,4); + histMultMax->RebinX(4); + histMultTot->RebinX(4); + // + TObjArray arrMax; + TObjArray arrTot; + histMultMax->FitSlicesY(0,0,-1,0,"QNR",&arrMax); + histMultTot->FitSlicesY(0,0,-1,0,"QNR",&arrTot); + // + TH1D * meanMax = (TH1D*)arrMax.At(1); + TH1D * meanTot = (TH1D*)arrTot.At(1); + Float_t meanMult = histMultMax->GetMean(); + if(meanMax->GetBinContent(meanMax->FindBin(meanMult))) { + meanMax->Scale(1./meanMax->GetBinContent(meanMax->FindBin(meanMult))); + } + else { + return kFALSE; + } + if(meanTot->GetBinContent(meanTot->FindBin(meanMult))) { + meanTot->Scale(1./meanTot->GetBinContent(meanTot->FindBin(meanMult))); + } + else { + return kFALSE; + } + Float_t xMultMax[50]; + Float_t yMultMax[50]; + Float_t yMultErrMax[50]; + Float_t xMultTot[50]; + Float_t yMultTot[50]; + Float_t yMultErrTot[50]; + // + Int_t nCountMax = 0; + for(Int_t iBin = 1; iBin < meanMax->GetXaxis()->GetNbins(); iBin++) { + Float_t yValMax = meanMax->GetBinContent(iBin); + if (yValMax < 0.7) continue; + if (yValMax > 1.3) continue; + if (meanMax->GetBinError(iBin)/yValMax > 0.01) continue; + xMultMax[nCountMax] = meanMax->GetXaxis()->GetBinCenter(iBin); + yMultMax[nCountMax] = yValMax; + yMultErrMax[nCountMax] = meanMax->GetBinError(iBin); + nCountMax++; + } + // + if (nCountMax < 10) return kFALSE; + TGraphErrors * fitMultMax = new TGraphErrors(nCountMax, xMultMax, yMultMax, 0, yMultErrMax); + fitMultMax->SetName("TGRAPHERRORS_MEANQMAX_MULTIPLICITYDEPENDENCE_BEAM_ALL"); + // + Int_t nCountTot = 0; + for(Int_t iBin = 1; iBin < meanTot->GetXaxis()->GetNbins(); iBin++) { + Float_t yValTot = meanTot->GetBinContent(iBin); + if (yValTot < 0.7) continue; + if (yValTot > 1.3) continue; + if (meanTot->GetBinError(iBin)/yValTot > 0.1) continue; + xMultTot[nCountTot] = meanTot->GetXaxis()->GetBinCenter(iBin); + yMultTot[nCountTot] = yValTot; + yMultErrTot[nCountTot] = meanTot->GetBinError(iBin); + nCountTot++; + } + // + if (nCountTot < 10) return kFALSE; + TGraphErrors * fitMultTot = new TGraphErrors(nCountTot, xMultTot, yMultTot, 0, yMultErrTot); + fitMultTot->SetName("TGRAPHERRORS_MEANQTOT_MULTIPLICITYDEPENDENCE_BEAM_ALL"); + // + fGainArray->AddLast(fitMultMax); + fGainArray->AddLast(fitMultTot); + // + return kTRUE; + +} + + void AliTPCPreprocessorOffline::UpdateOCDBGain(Int_t startRunNumber, Int_t endRunNumber, const Char_t *storagePath){ // // Update OCDB entry @@ -1010,9 +1186,11 @@ void AliTPCPreprocessorOffline::MakeQAPlot(Float_t FPtoMIPratio) { fGraphCosmic->GetY()[i] *= FPtoMIPratio; } } - fGraphCosmic->Draw("lp"); - grfFitCosmic->SetLineColor(2); - grfFitCosmic->Draw("lu"); + fGraphCosmic->Draw("lp"); + if (grfFitCosmic) { + grfFitCosmic->SetLineColor(2); + grfFitCosmic->Draw("lu"); + } fGainArray->AddLast(gainHistoCosmic); //fGainArray->AddLast(canvasCosmic->Clone()); delete canvasCosmic; @@ -1037,5 +1215,363 @@ void AliTPCPreprocessorOffline::MakeQAPlot(Float_t FPtoMIPratio) { } } +void AliTPCPreprocessorOffline::MakeFitTime(){ + // + // mak aligment fit - store results in the file + // + const Int_t kMinEntries=1000; + MakeChainTime(); + MakePrimitivesTime(); + if (!fAlignTree) return; + if (fAlignTree->GetEntries()SetAlias("ptype","type"); + fAlignTree->SetAlias("hasITS","(1+0)"); + fAlignTree->SetAlias("dITS","1-2*(refX<40)"); + fAlignTree->SetAlias("isITS","refX>10"); + fAlignTree->SetAlias("isVertex","refX<10"); + // + Int_t npointsMax=30000000; + TStatToolkit toolkit; + Double_t chi2=0; + Int_t npoints=0; + TVectorD param; + TMatrixD covar; + + TString fstringFast=""; + fstringFast+="FExBTwistX++"; + fstringFast+="FExBTwistY++"; + fstringFast+="FAlignRot0D++"; + fstringFast+="FAlignTrans0D++"; + fstringFast+="FAlignTrans1D++"; + // + fstringFast+="hasITS*FAlignTrans0++"; + fstringFast+="hasITS*FAlignTrans1++"; + fstringFast+="hasITS*FAlignRot0++"; + fstringFast+="hasITS*FAlignRot1++"; + fstringFast+="hasITS*FAlignRot2++"; + // + fstringFast+="dITS*FAlignTrans0++"; + fstringFast+="dITS*FAlignTrans1++"; + fstringFast+="dITS*FAlignRot0++"; + fstringFast+="dITS*FAlignRot1++"; + fstringFast+="dITS*FAlignRot2++"; + + TCut cutFit="entries>10&&abs(mean)>0.00001"; + fAlignTree->SetAlias("err","rms"); + + TString *strDeltaITS = TStatToolkit::FitPlaneConstrain(fAlignTree,"mean:err", fstringFast.Data(),cutFit, chi2,npoints,param,covar,-1,0, npointsMax, 1); + strDeltaITS->Tokenize("++")->Print(); + fAlignTree->SetAlias("fitYFast",strDeltaITS->Data()); + // + TVectorD paramC= param; + TMatrixD covarC= covar; + TStatToolkit::Constrain1D(fstringFast,"Trans0D",paramC,covarC,0, 0.1); + TStatToolkit::Constrain1D(fstringFast,"Trans1D",paramC,covarC,0, 0.1); + TStatToolkit::Constrain1D(fstringFast,"TwistX",paramC,covarC,0, 0.1); + TStatToolkit::Constrain1D(fstringFast,"TwistY",paramC,covarC,0, 0.1); + TString strFitConst=TStatToolkit::MakeFitString(fstringFast, paramC,covar); + fAlignTree->SetAlias("fitYFastC",strFitConst.Data()); + CreateAlignTime(fstringFast,paramC); + + +} + + +void AliTPCPreprocessorOffline::MakeChainTime(){ + // + TFile f("CalibObjects.root"); + // const char *cdtype[7]={"ITS","TRD","Vertex","TOF","TPC","TPC0","TPC1"}; + //const char *cptype[5]={"dy","dz","dsnp","dtheta","d1pt"}; + const char * hname[5]={"dy","dz","dsnp","dtheta","d1pt"}; + Int_t run=0; + AliTPCcalibTime *calibTime = 0; + TObjArray * array = (TObjArray*)f.Get("TPCCalib"); + if (array){ + calibTime = (AliTPCcalibTime *)array->FindObject("calibTime"); + } else { + calibTime = (AliTPCcalibTime*)f.Get("calibTime"); + } + if (!calibTime) return; + TTreeSRedirector *pcstream = new TTreeSRedirector("meanITSVertex.root"); + // + Int_t ihis=0; + THnSparse *his = calibTime->GetResHistoTPCITS(ihis); + if (his){ + his->GetAxis(1)->SetRangeUser(-1.1,1.1); + his->GetAxis(2)->SetRange(0,1000000); + his->GetAxis(3)->SetRangeUser(-0.35,0.35); + AliTPCCorrection::MakeDistortionMap(his,pcstream, Form("ITS%s",hname[ihis]),run,85.,ihis,5); + } + ihis=1; + his = calibTime->GetResHistoTPCITS(ihis); + if (his){ + his->GetAxis(1)->SetRangeUser(-1.1,1.1); + his->GetAxis(2)->SetRange(0,1000000); + his->GetAxis(3)->SetRangeUser(-0.35,0.35); + AliTPCCorrection::MakeDistortionMap(his,pcstream, Form("ITS%s",hname[ihis]),run,85.,ihis,5); + } + ihis=2; + his = calibTime->GetResHistoTPCITS(ihis); + if (his){ + his->GetAxis(1)->SetRangeUser(-1.1,1.1); + his->GetAxis(2)->SetRange(0,1000000); + his->GetAxis(3)->SetRangeUser(-0.35,0.35); + AliTPCCorrection::MakeDistortionMap(his,pcstream, Form("ITS%s",hname[ihis]),run,85.,ihis,5); + } + ihis=0; + his = calibTime->GetResHistoTPCvertex(ihis); + if (his){ + his->GetAxis(1)->SetRangeUser(-1.1,1.1); + his->GetAxis(2)->SetRange(0,1000000); + his->GetAxis(3)->SetRangeUser(-0.35,0.35); + AliTPCCorrection::MakeDistortionMap(his,pcstream, Form("Vertex%s",hname[ihis]),run,0.,ihis,5); + } + ihis=2; + his = calibTime->GetResHistoTPCvertex(ihis); + if (his){ + his->GetAxis(1)->SetRangeUser(-1.1,1.1); + his->GetAxis(2)->SetRange(0,1000000); + his->GetAxis(3)->SetRangeUser(-0.35,0.35); + AliTPCCorrection::MakeDistortionMap(his,pcstream, Form("Vertex%s",hname[ihis]),run,0.,ihis,5); + + } + ihis=1; + his = calibTime->GetResHistoTPCvertex(ihis); + if (his){ + his->GetAxis(1)->SetRangeUser(-1.1,1.1); + his->GetAxis(2)->SetRange(0,1000000); + his->GetAxis(3)->SetRangeUser(-0.35,0.35); + AliTPCCorrection::MakeDistortionMap(his,pcstream, Form("Vertex%s",hname[ihis]),run,0.,ihis,5); + + } + ihis=0; + his = calibTime->GetResHistoTPCTOF(ihis); + if (his){ + his->GetAxis(1)->SetRangeUser(-1.1,1.1); + his->GetAxis(2)->SetRange(0,1000000); + his->GetAxis(3)->SetRangeUser(-0.35,0.35); + AliTPCCorrection::MakeDistortionMap(his,pcstream, Form("TOF%s",hname[ihis]),run,0.,ihis,10); + + } + ihis=0; + his = calibTime->GetResHistoTPCTRD(ihis); + if (his){ + his->GetAxis(1)->SetRangeUser(-1.1,1.1); + his->GetAxis(2)->SetRange(0,1000000); + his->GetAxis(3)->SetRangeUser(-0.35,0.35); + AliTPCCorrection::MakeDistortionMap(his,pcstream, Form("TRD%s",hname[ihis]),run,0.,ihis,10); + + } + delete pcstream; +} + + +Double_t AliTPCPreprocessorOffline::EvalAt(Double_t phi, Double_t refX, Double_t theta, Int_t corr, Int_t ptype){ + // + // + // + Double_t sector = 9*phi/TMath::Pi(); + if (sector<0) sector+=18; + Double_t y85=AliTPCCorrection::GetCorrSector(sector,85,theta,1,corr); + Double_t y245=AliTPCCorrection::GetCorrSector(sector,245,theta,1,corr); + if (ptype==0) return y85+(y245-y85)*(refX-85.)/(245.-85.); + if (ptype==2) return (y245-y85)/(245.-85.); + return 0; +} + + + +void AliTPCPreprocessorOffline::MakePrimitivesTime(){ + // + // Create primitive transformation to fit + // + fAlignTree=new TChain("fit","fit"); + fAlignTree->AddFile("meanITSVertex.root",10000000,"ITSdy"); + fAlignTree->AddFile("meanITSVertex.root",10000000,"ITSdsnp"); + fAlignTree->AddFile("meanITSVertex.root",10000000,"Vertexdy"); + fAlignTree->AddFile("meanITSVertex.root",10000000,"Vertexdsnp"); + // + AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters(); + Double_t bzField=AliTrackerBase::GetBz(); + Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us] // From dataBase: to be updated: per second (ideally) + Double_t ezField = 400; // [V/cm] // to be updated: never (hopefully) + Double_t wtP = -10.0 * (bzField) * vdrift / ezField ; + AliTPCExBTwist *fitExBTwistX= new AliTPCExBTwist; + AliTPCExBTwist *fitExBTwistY= new AliTPCExBTwist; + AliTPCCalibGlobalMisalignment *trans0 =new AliTPCCalibGlobalMisalignment; + AliTPCCalibGlobalMisalignment *trans1 =new AliTPCCalibGlobalMisalignment; + AliTPCCalibGlobalMisalignment *trans0D =new AliTPCCalibGlobalMisalignment; + AliTPCCalibGlobalMisalignment *trans1D =new AliTPCCalibGlobalMisalignment; + AliTPCCalibGlobalMisalignment *rot0 =new AliTPCCalibGlobalMisalignment; + AliTPCCalibGlobalMisalignment *rot1 =new AliTPCCalibGlobalMisalignment; + AliTPCCalibGlobalMisalignment *rot2 =new AliTPCCalibGlobalMisalignment; + AliTPCCalibGlobalMisalignment *rot3 =new AliTPCCalibGlobalMisalignment; + // + // + fitExBTwistX->SetXTwist(0.001); + fitExBTwistX->SetOmegaTauT1T2(wtP,1,1); + // + fitExBTwistY->SetYTwist(0.001); + fitExBTwistY->SetOmegaTauT1T2(wtP,1,1); + // + TGeoHMatrix *matrixRot = new TGeoHMatrix; + TGeoHMatrix *matrixX = new TGeoHMatrix; + TGeoHMatrix *matrixY = new TGeoHMatrix; + matrixX->SetDx(0.1); + matrixY->SetDy(0.1); + Double_t rotAngles0[9]={0}; + Double_t rotAngles1[9]={0}; + Double_t rotAngles2[9]={0}; + // + Double_t rotAngles3[9]={0}; + + rotAngles0[0]=1; rotAngles0[4]=1; rotAngles0[8]=1; + rotAngles1[0]=1; rotAngles1[4]=1; rotAngles1[8]=1; + rotAngles2[0]=1; rotAngles2[4]=1; rotAngles2[8]=1; + rotAngles3[0]=1; rotAngles3[4]=1; rotAngles3[8]=1; + + rotAngles0[1]=-0.001;rotAngles0[3]=0.001; + rotAngles1[5]=-0.001;rotAngles1[7]=0.001; + rotAngles2[2]=0.001;rotAngles2[6]=-0.001; + rotAngles3[1]=0.001;rotAngles3[3]=-0.001; + matrixRot->SetRotation(rotAngles0); + rot0->SetAlignGlobal(matrixRot); + matrixRot->SetRotation(rotAngles1); + rot1->SetAlignGlobal(matrixRot); + matrixRot->SetRotation(rotAngles2); + rot2->SetAlignGlobal(matrixRot); + matrixRot->SetRotation(rotAngles3); + rot3->SetAlignGlobalDelta(matrixRot); + // + trans0->SetAlignGlobal(matrixX); + trans1->SetAlignGlobal(matrixY); + trans0D->SetAlignGlobalDelta(matrixX); + trans1D->SetAlignGlobalDelta(matrixY); + fitExBTwistX->Init(); + fitExBTwistY->Init(); + // + fitExBTwistX->AddVisualCorrection((AliTPCExBTwist*)(fitExBTwistX->Clone()),100); + fitExBTwistY->AddVisualCorrection((AliTPCExBTwist*)(fitExBTwistY->Clone()),101); + // + fitExBTwistX->AddVisualCorrection((AliTPCExBTwist*)(rot0->Clone()),102); + fitExBTwistY->AddVisualCorrection((AliTPCExBTwist*)(rot1->Clone()),103); + fitExBTwistX->AddVisualCorrection((AliTPCExBTwist*)(rot2->Clone()),104); + fitExBTwistX->AddVisualCorrection((AliTPCExBTwist*)(rot3->Clone()),105); + + fitExBTwistX->AddVisualCorrection((AliTPCExBTwist*)(trans0->Clone()),106); + fitExBTwistX->AddVisualCorrection((AliTPCExBTwist*)(trans1->Clone()),107); + fitExBTwistX->AddVisualCorrection((AliTPCExBTwist*)(trans0D->Clone()),108); + fitExBTwistX->AddVisualCorrection((AliTPCExBTwist*)(trans1D->Clone()),109); + // + fAlignTree->SetAlias("FExBTwistX", "AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,100,ptype)+0"); + fAlignTree->SetAlias("FExBTwistY","AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,101,ptype)+0"); + fAlignTree->SetAlias("FAlignRot0","AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,102,ptype)+0"); + fAlignTree->SetAlias("FAlignRot0D","AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,105,ptype)+0"); + fAlignTree->SetAlias("FAlignRot1","AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,103,ptype)+0"); + fAlignTree->SetAlias("FAlignRot2","AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,104,ptype)+0"); + fAlignTree->SetAlias("FAlignTrans0","AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,106,ptype)+0"); + fAlignTree->SetAlias("FAlignTrans1","AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,107,ptype)+0"); + fAlignTree->SetAlias("FAlignTrans0D","AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,108,ptype)+0"); + fAlignTree->SetAlias("FAlignTrans1D","AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,109,ptype)+0"); + // + // test fast function + // +// fAlignTree->Draw("FExBTwistX:ExBTwistX","isITS&&ptype==0&&abs(snp)<0.05",""); +// fAlignTree->Draw("FExBTwistY:ExBTwistY","isITS&&ptype==0&&abs(snp)<0.05",""); +// fAlignTree->Draw("FAlignRot0:alignRot0","isITS&&ptype==0&&abs(snp)<0.05",""); +// fAlignTree->Draw("FAlignRot1:alignRot1","isITS&&ptype==0&&abs(snp)<0.05",""); +// fAlignTree->Draw("FAlignRot2:alignRot2","isITS&&ptype==0&&abs(snp)<0.05",""); +// // +// fAlignTree->Draw("FAlignTrans0:alignTrans0","isITS&&ptype==0&&abs(snp)<0.05",""); +// fAlignTree->Draw("FAlignTrans1:alignTrans1","isITS&&ptype==0&&abs(snp)<0.05",""); + +} + + +void AliTPCPreprocessorOffline::CreateAlignTime(TString fstring, TVectorD paramC){ + // + // + // + // + TGeoHMatrix *matrixDelta = new TGeoHMatrix; + TGeoHMatrix *matrixGlobal = new TGeoHMatrix; + Double_t rAngles[9]; + Int_t index=0; + // + index=TStatToolkit::GetFitIndex(fstring,"FAlignTrans0D"); + if (index>=0) matrixDelta->SetDx(paramC[index+1]*0.1); + index=TStatToolkit::GetFitIndex(fstring,"FAlignTrans1D"); + if (index>=0) matrixDelta->SetDy(paramC[index+1]*0.1); + rAngles[0]=1; rAngles[4]=1; rAngles[8]=1; + index=TStatToolkit::GetFitIndex(fstring,"FAlignRot0D"); + rAngles[1]=-paramC[index+1]*0.001; rAngles[3]=paramC[index+1]*0.001; + rAngles[5]=0; rAngles[7] =0; + rAngles[2]=0; rAngles[6] =0; + matrixDelta->SetRotation(rAngles); + // + // + // + index=TStatToolkit::GetFitIndex(fstring,"hasITS*FAlignTrans0"); + if (index>=0) matrixGlobal->SetDx(paramC[index+1]*0.1); + index=TStatToolkit::GetFitIndex(fstring,"hasITS*FAlignTrans1"); + if (index>=0) matrixGlobal->SetDy(paramC[index+1]*0.1); + rAngles[0]=1; rAngles[4]=1; rAngles[8]=1; + index=TStatToolkit::GetFitIndex(fstring,"hasITS*FAlignRot0"); + rAngles[1]=-paramC[index+1]*0.001; rAngles[3]=paramC[index+1]*0.001; + index=TStatToolkit::GetFitIndex(fstring,"hasITS*FAlignRot1"); + rAngles[5]=-paramC[index+1]*0.001; rAngles[7]=paramC[index+1]*0.001; + index=TStatToolkit::GetFitIndex(fstring,"hasITS*FAlignRot2"); + rAngles[2]=paramC[index+1]*0.001; rAngles[6] =-paramC[index+1]*0.001; + matrixGlobal->SetRotation(rAngles); + // + AliTPCCalibGlobalMisalignment *fitAlignTime =0; + fitAlignTime =new AliTPCCalibGlobalMisalignment; + fitAlignTime->SetName("FitAlignTime"); + fitAlignTime->SetTitle("FitAlignTime"); + fitAlignTime->SetAlignGlobalDelta(matrixDelta); + fitAlignTime->SetAlignGlobal(matrixGlobal); + // + AliTPCExBTwist * fitExBTwist= new AliTPCExBTwist; + Int_t indexX=TStatToolkit::GetFitIndex(fstring,"ExBTwistX"); + Int_t indexY=TStatToolkit::GetFitIndex(fstring,"ExBTwistY"); + fitExBTwist->SetXTwist(0.001*paramC[indexX+1]); // 1 mrad twist in x + fitExBTwist->SetYTwist(0.001*paramC[indexY+1]); // 1 mrad twist in x + fitExBTwist->SetName("FitExBTwistTime"); + fitExBTwist->SetTitle("FitExBTwistTime"); + AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters(); + Double_t bzField=AliTrackerBase::GetBz(); + Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us] // From dataBase: to be updated: per second (ideally) + + Double_t ezField = 400; // [V/cm] // to be updated: never (hopefully) + Double_t wt = -10.0 * (bzField) * vdrift / ezField ; + // + fitExBTwist->SetOmegaTauT1T2(wt,1,1); + fitExBTwist->Init(); + + AliTPCComposedCorrection *corrTime = new AliTPCComposedCorrection; + TObjArray *arr = new TObjArray; + corrTime->SetCorrections(arr); + + corrTime->GetCorrections()->Add(fitExBTwist); + corrTime->GetCorrections()->Add(fitAlignTime); + corrTime->SetName("FitCorrectionTime"); + corrTime->SetTitle("FitCorrectionTime"); + + fitExBTwist->AddVisualCorrection((AliTPCExBTwist*)(fitExBTwist->Clone()),1001); + fitAlignTime->AddVisualCorrection((AliTPCExBTwist*)(fitAlignTime->Clone()),1002); + fitAlignTime->AddVisualCorrection((AliTPCExBTwist*)(corrTime->Clone()),1003); + + + fAlignTree->SetAlias("ExBTwistTime","AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,1001,ptype)+0"); + fAlignTree->SetAlias("AlignTime","AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,1002,ptype)+0"); + fAlignTree->SetAlias("FitCorrectionTime","AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,1003,ptype)+0"); + + + TFile *f = new TFile("fitITSVertex.root","update"); + corrTime->Write("FitCorrectionTime"); + f->Close(); +} + +