From b1c27e4ff4ef33fbc95154fceee7fc6868c878b3 Mon Sep 17 00:00:00 2001 From: marian Date: Mon, 8 Mar 2010 12:55:50 +0000 Subject: [PATCH] M AliTPCcalibLaser.cxx - Adding Integral of B field M AliTPCcalibTime.cxx - Protection in case of missing friends In case of filtered ESDs CalibExB.C - Obsolete macro removed CalibLaserExBscan.C - Replaced by another Marian --- TPC/AliTPCcalibLaser.cxx | 45 +++- TPC/AliTPCcalibTime.cxx | 1 + TPC/CalibMacros/CalibExB.C | 277 ----------------------- TPC/CalibMacros/CalibLaserExBscan.C | 329 ++++++++++++++++++++++++++++ 4 files changed, 367 insertions(+), 285 deletions(-) delete mode 100644 TPC/CalibMacros/CalibExB.C create mode 100644 TPC/CalibMacros/CalibLaserExBscan.C diff --git a/TPC/AliTPCcalibLaser.cxx b/TPC/AliTPCcalibLaser.cxx index bc999cc4678..ba72a711c87 100644 --- a/TPC/AliTPCcalibLaser.cxx +++ b/TPC/AliTPCcalibLaser.cxx @@ -611,7 +611,7 @@ void AliTPCcalibLaser::Process(AliESDEvent * event) { if (counter (TGeoGlobalMagField::Instance()->GetField()); - Double_t xyz[3]={90,0,10}; - Double_t bxyz[3]={90,0,10}; + Double_t xyz[3]={90,0,10}; // tmp. global position + Double_t bxyz[3]={90,0,10}; // tmp. mag field integral - cylindrical + Double_t bgxyz[3]={90,0,10}; // tmp. mag field integral - local // AliTPCcalibLaser *laser = this; TTreeSRedirector *pcstream = new TTreeSRedirector("laserMean.root"); @@ -2556,21 +2557,42 @@ void AliTPCcalibLaser::DumpMeanInfo(Int_t run){ // magnetic field integrals TVectorD vecIBR(159); // radial TVectorD vecIBRPhi(159); // r-phi + TVectorD vecIBLX(159); // local x + TVectorD vecIBLY(159); // local y + TVectorD vecIBGX(159); // local x + TVectorD vecIBGY(159); // local y TVectorD vecIBZ(159); // z // for (Int_t irow=0;irow<159;irow++){ vecIBR[irow]=0; vecIBRPhi[irow]=0; + vecIBLX[irow]=0; + vecIBLY[irow]=0; + vecIBGX[irow]=0; + vecIBGY[irow]=0; vecIBZ[irow]=0; - Double_t gx=(*(ltrp->fVecGX))[irow]; - Double_t gy=(*(ltrp->fVecGY))[irow]; + Double_t gx =(*(ltrp->fVecGX))[irow]; + Double_t gy =(*(ltrp->fVecGY))[irow]; + Int_t lsec =TMath::Nint((*(ltrp->fVecSec))[irow]); + Double_t ca =TMath::Cos(TMath::Pi()*(lsec+0.5)/9.); + Double_t sa =TMath::Sin(TMath::Pi()*(lsec+0.5)/9.); xyz[2]=(*(ltrp->fVecGZ))[irow]; xyz[0]=TMath::Sqrt(gx*gx+gy*gy); xyz[1]=TMath::ATan2(gy,gx); + Double_t gxyz[3]={gx,gy,(*(ltrp->fVecGZ))[irow]}; if (magF){ magF->GetTPCIntCyl(xyz,bxyz); + magF->GetTPCInt(gxyz,bgxyz); vecIBR[irow]=bxyz[0]; vecIBRPhi[irow]=bxyz[1]; + // + vecIBGX[irow]=bgxyz[0]; + vecIBGY[irow]=bgxyz[1]; + // + vecIBLX[irow]= bgxyz[0]*ca+bgxyz[1]*sa; + vecIBLY[irow]= -bgxyz[0]*sa+bgxyz[1]*ca; + // + vecIBZ[irow]=bxyz[2]; } } @@ -2914,7 +2936,14 @@ void AliTPCcalibLaser::DumpMeanInfo(Int_t run){ "lasTanPhiLocOut="<Load("libSTAT.so"); - AliTPCExBFirst *exbfirst1 = new AliTPCExBFirst(fieldC1,0.88*2.6400e+04,50,50,50); - AliTPCExB::SetInstance(exbfirst1); - -} - -void FitField(){ - // - // - // - AliTPCExB * exb = AliTPCExB::Instance(); - exb->TestExB("field.root"); - TFile f("field.root"); - TTree * tree = (TTree*)f.Get("positions"); - // - TStatToolkit toolkit; - Double_t chi2; - TVectorD fitParam; - TMatrixD covMatrix; - Int_t npoints; - // - // SetAliases - // - tree->SetAlias("sa","sin(phi+0.0)"); - tree->SetAlias("ca","cos(phi+0.0)"); - tree->SetAlias("sa2","sin(phi*2+0.0)"); - tree->SetAlias("ca2","cos(phi*2+0.0)"); - tree->SetAlias("zn","(x2/250.)"); - tree->SetAlias("rn","(r/250.)"); - - TString fstringSym=""; - // - fstringSym+="zn++"; - fstringSym+="rn++"; - fstringSym+="zn*rn++"; - fstringSym+="zn*zn++"; - fstringSym+="zn*zn*rn++"; - fstringSym+="zn*rn*rn++"; - // - fstringSym+="sa++"; - fstringSym+="ca++"; - fstringSym+="ca2++"; - fstringSym+="sa2++"; - fstringSym+="ca*zn++"; - fstringSym+="sa*zn++"; - fstringSym+="ca2*zn++"; - fstringSym+="sa2*zn++"; - fstringSym+="ca*zn*zn++"; - fstringSym+="sa*zn*zn++"; - fstringSym+="ca*zn*rn++"; - fstringSym+="sa*zn*rn++"; - - // BrBz - TString *strBrBz = toolkit.FitPlane(tree,"br/bz",fstringSym->Data(), "abs(r)<250&&abs(r)>90", chi2,npoints,fitParam,covMatrix); - strBrBz->Tokenize("+")->Print(); - printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints)); - tree->SetAlias("fitBrBz",strBrBz->Data()); - exb->fMatBrBz = new TVectorD(fitParam); - // BrfiBz - TString *strBrfiBz = toolkit.FitPlane(tree,"brfi/bz",fstringSym->Data(), "abs(r)<250&&abs(r)>90", chi2,npoints,fitParam,covMatrix); - strBrfiBz->Tokenize("+")->Print(); - printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints)); - tree->SetAlias("fitBrfiBz",strBrfiBz->Data()); - exb->fMatBrfiBz = new TVectorD(fitParam); - - - // - // BR integral parameterization - // - TString *strBRI0 = toolkit.FitPlane(tree,"bri",fstringSym->Data(), "abs(r)<250&&abs(r)>90&&x2>0&&x2<250", chi2,npoints,fitParam,covMatrix); - strBRI0->Tokenize("+")->Print(); - printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints)); - tree->SetAlias("fitBRI0",strBRI0->Data()); - exb->fMatBrBzI0 = new TVectorD(fitParam); - - TString *strBRI1 = toolkit.FitPlane(tree,"bri",fstringSym->Data(), "abs(r)<250&&abs(r)>90&&x2<-10&&x2>-250", chi2,npoints,fitParam,covMatrix); - strBRI1->Tokenize("+")->Print(); - printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints)); - tree->SetAlias("fitBRI1",strBRI1->Data()); - exb->fMatBrBzI1 = new TVectorD(fitParam); - // - - TString *strBRFII0 = toolkit.FitPlane(tree,"brfii",fstringSym->Data(), "abs(r)<250&&abs(r)>90&&x2>0&&x2<250", chi2,npoints,fitParam,covMatrix); - strBRFII0->Tokenize("+")->Print(); - printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints)); - tree->SetAlias("fitBRFII0",strBRFII0->Data()); - exb->fMatBrfiBzI0 = new TVectorD(fitParam); - - TString *strBRFII1 = toolkit.FitPlane(tree,"brfii",fstringSym->Data(), "abs(r)<250&&abs(r)>90&&x2<-10&&x2>-240", chi2,npoints,fitParam,covMatrix); - strBRFII1->Tokenize("+")->Print(); - printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints)); - tree->SetAlias("fitBRFII1",strBRFII1->Data()); - exb->fMatBrfiBzI1 = new TVectorD(fitParam); - - TFile * fout = new TFile("bfit.root","recreate"); - exb->Write("bfit"); - fout->Close(); - -} - - - - - - -void CalibExB(){ - // - // Macro to test ExB correction versus laser B filed scan data - // - // Before running it be sure the file laserScan.root exist - // - // Create this file from scan - See AliTPCCalibLaser.C - // - TFile fscan("laserScan.root"); - TTree * treeT = (TTree*)fscan.Get("Mean"); - TStatToolkit toolkit; - Double_t chi2; - TVectorD fitParam; - TMatrixD covMatrix; - Int_t npoints; - - TCut cutF0("abs(gphi1-(pphi0+pphi1*bz))<0.05"); - TCut cutF1("abs(gphiP1-(pphiP0+pphiP1*bz))<0.0005"); - TCut cutN("entries>2"); - - TCut cutA = cutF0+cutF1+cutN; - - treeT->SetAlias("side","(-1+(LTr.fP[1]>0)*2)"); // side - treeT->SetAlias("dr","(abs(LTr.fP[1]/250.))"); - treeT->SetAlias("sa","sin(atan2(lx1+0.0,lx0+0.0))"); - treeT->SetAlias("ca","cos(atan2(lx1+0.0,lx0+0.0))"); - treeT->SetAlias("ta","tan(asin(LTr.fP[2]+0.0))"); - - TString fstring=""; - // - fstring+="((dr)^3-1)*bz++"; //1 - fstring+="((dr)^3-1)*bz*sa++"; //3 - fstring+="((dr)^3-1)*bz*ca++"; //5 - // - fstring+="((dr)^1-1)*bz++"; //1 - fstring+="((dr)^1-1)*bz*sa++"; //3 - fstring+="((dr)^1-1)*bz*ca++"; //5 - // - fstring+="side*((dr)^3-1)*bz++"; //1 - fstring+="side*((dr)^3-1)*bz*sa++"; //3 - fstring+="side*((dr)^3-1)*bz*ca++"; //5 - // - fstring+="side*((dr)^1-1)*bz++"; //7 - fstring+="side*((dr)^1-1)*bz*sa++"; //9 - fstring+="side*((dr)^1-1)*bz*ca++"; //11 - // - fstring+="((dr)^3-1)*bz^2*ta++"; //2 - fstring+="((dr)^3-1)*bz^2*sa*ta++"; //4 - fstring+="((dr)^3-1)*bz^2*ca*ta++"; //6 - - fstring+="((dr)^1-1)*bz^2*ta++"; //2 - fstring+="((dr)^1-1)*bz^2*sa*ta++"; //4 - fstring+="((dr)^1-1)*bz^2*ca*ta++"; //6 - // - fstring+="side*((dr)^1-1)*bz^2*ta++"; //8 - fstring+="side*((dr)^1-1)*bz^2*sa*ta++"; //10 - fstring+="side*((dr)^1-1)*bz^2*ca*ta++"; //12 - - - - - TString *strq0 = toolkit.FitPlane(treeT,"gphi1-pphi0",fstring->Data(), cutA, chi2,npoints,fitParam,covMatrix); - strq0->Tokenize("+")->Print(); - treeT->SetAlias("fit",strq0->Data()); - printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints)); - - - - treeT->SetAlias("bir1", "-AliTPCExB::GetBrI(254,atan2(lx1,lx0+0),LTr.fP[1],2)"); - treeT->SetAlias("birfi1","-AliTPCExB::GetBrfiI(254,atan2(lx1,lx0+0),LTr.fP[1],2)"); - - treeT->SetAlias("bir0", "AliTPCExB::GetBrI(254,atan2(lx1,lx0+0),LTr.fP[1],2)"); - treeT->SetAlias("birfi0","AliTPCExB::GetBrfiI(254,atan2(lx1,lx0+0),LTr.fP[1],2)"); - - treeT->SetAlias("fbz00", "(bir0+birfi0*ta)"); - treeT->SetAlias("fbz02", "(birfi0+bir0*ta)"); - treeT->SetAlias("fbz10", "(bir1+birfi1*ta)"); - treeT->SetAlias("fbz12", "(birfi1+bir1*ta)"); - // - treeT->SetAlias("fbz0", "((fSide==0)*fbz00+(fSide==1)*fbz10)"); - treeT->SetAlias("fbz2", "((fSide==0)*fbz02+(fSide==1)*fbz12)"); - - // - TString fstringeb=""; - // - fstringeb+="bz*fbz0++"; //1 - fstringeb+="sign(bz)*bz^2*fbz2++"; //2 - fstringeb+="((dr)-1)*bz*sa++"; //9 - fstringeb+="side*((dr)-1)*bz*sa++"; //9 - - - fstringeb+="((dr)-1)*bz*ca++"; //9 - fstringeb+="side*((dr)^3-1)*bz*sa++"; //9 - fstringeb+="side*((dr)^3-1)*bz++"; //7 - - fstringeb+="((dr)^3-1)*bz++"; //7 - fstringeb+="side*((dr)^3-1)*bz*ca++"; //11 - - // - TString *strExB = toolkit.FitPlane(treeT,"gphi1-pphi0",fstringeb->Data(), "abs(gphi1-pphi0-fit)<0.06&&abs(bz)>0.1"+cutA, chi2,npoints,fitParam,covMatrix); - strExB->Tokenize("+")->Print(); - printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints)); - treeT->SetAlias("fitEB",strExB->Data()); - - - - TString fstringeb=""; - // - fstringeb+="AliTPCExB::GetDrphi(254,atan2(lx1,lx0),LTr.fP[1],bz*10)++"; //1 - fstringeb+="AliTPCExB::GetDr(254,atan2(lx1,lx0),LTr.fP[1],-bz*10)*ta++"; //2 - // - fstringeb+="side*((dr)^1-1)*bz++"; //9 - fstringeb+="side*((dr)^1-1)*bz*sa++"; //9 - // - fstringeb+="side*((dr)^1-1)*bz*ta++"; //9 - fstringeb+="side*((dr)^1-1)*bz*sa*ta++"; //9 - fstringeb+="side*((dr)^1-1)*bz*ca*ta++"; //9 - - TString *strExB = toolkit.FitPlane(treeT,"fit",fstringeb->Data(), "abs(gphi1-pphi0-fit)<0.06&&abs(bz)>0.1"+cutA, chi2,npoints,fitParam,covMatrix); - strExB->Tokenize("+")->Print(); - printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints)); - treeT->SetAlias("fitEB",strExB->Data()); - - - - - - - - { - Float_t otcor0 =0.7; - Float_t omegatau0 =0.35; - Float_t rms0 =10; - for (Float_t otcor=0.6; otcor<1.1;otcor+=0.1) - for (Float_t ometau=0.4; ometau<0.45;ometau+=0.01){ - char hname[100]; - sprintf(hname,"hist_ometau_%f(50,-0.10,0.10)",ometau); - char expr[100]; - sprintf(expr,"gphi1-pphi0-bz*fbz0*%f-sign(bz)*(%f*%f*bz)^2*fbz2>>%s",ometau,ometau,otcor,hname); - treeT->Draw(expr,"abs(gphi1-pphi0-fit)<0.06&&abs(bz)>0.1"+cutA); - printf("Ometau=%f\tCor=%f\tRMS=%f\n",ometau,otcor,treeT->GetHistogram()->GetRMS()); - if (rms0>treeT->GetHistogram()->GetRMS()){ - otcor0=otcor; - omegatau0=ometau; - rms0=treeT->GetHistogram()->GetRMS(); - } - } - } - -} - - - - -void MakePic(){ - // - // - // - -} - diff --git a/TPC/CalibMacros/CalibLaserExBscan.C b/TPC/CalibMacros/CalibLaserExBscan.C new file mode 100644 index 00000000000..abe66fa49f6 --- /dev/null +++ b/TPC/CalibMacros/CalibLaserExBscan.C @@ -0,0 +1,329 @@ +/* +// + .x ~/rootlogon.C + .L $ALICE_ROOT/TPC/CalibMacros/CalibLaserExBscan.C + + + // 0. Make a calibration + // 1. Make a laser scan list + // e.g in TPC workscape + + // 2. Define a reference data e.g: + // for a in `cat laserScan.txt`; do echo `pwd`/../mergerunMag0.list/laserMean.root; done >laserScanRef.txt + + Init(); // init chain + MakeAliases(); // make alaises for variables + + +*/ + + + + +TCut cutE="eY.fElements<0.02&&Rr.eY.fElements<0.02"; // error cut 200 microns +TCut cutN="nCl.fElements>10&&Rr.nCl.fElements>10"; // number of clusters cut +TCut cutEd="(abs(Rr.Y.fElements)-Rr.X.fElements*0.155<-1)"; // edge cut +TCut cutB="LTr.fVecLX.fElements>0"; // local x cut +TCut cutR="(LTr.fVecLX.fElements%3)==1&&((abs(LTr.fVecLZ.fElements)+abs(LTr.fVecLY.fElements))%3)==1"; // cut - skip points +TCut cutA=cutE+cutN+cutEd+cutB; +// +TCut cutAM5 =cutA+"LTr.fP[1]>0&&abs(bz+5)<0.1"; +TCut cutAP2 =cutA+"LTr.fP[1]>0&&abs(bz-2)<0.1"; +TCut cutAP5 =cutA+"LTr.fP[1]>0&&abs(bz-5)<0.1"; +// +TCut cutCM5 =cutA+"LTr.fP[1]<0&&abs(bz+5)<0.1"; +TCut cutCP2 =cutA+"LTr.fP[1]<0&&abs(bz-2)<0.1"; +TCut cutCP5 =cutA+"LTr.fP[1]<0&&abs(bz-5)<0.1"; + + +TChain * chain =0; + +TVectorD fitParamB[6]; // magnetic fit param +TVectorD fitParamBT[6]; // + electric field tilt +TVectorD fitParamBT0[6]; // + electric field close to ROC +TVectorD fitParamA[6]; // + electric field rotation +// +TMatrixD covarB[6]; // covariance matrices +TMatrixD covarBT[6]; // +TMatrixD covarBT0[6]; // +TMatrixD covarA[6]; // +// +TMatrixD *errB[6]; // covariance matrices +TMatrixD *errBT[6]; // +TMatrixD *errBT0[6]; // +TMatrixD *errA[6]; // +// +Double_t chi2B[6]; // chi2 +Double_t chi2BT[6]; // chi2 +Double_t chi2BT0[6]; // +Double_t chi2A[6]; // + + +void Init(){ + // + // + // + gSystem->Load("libANALYSIS"); + gSystem->Load("libTPCcalib"); + gSystem->Load("libSTAT.so"); + gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros"); + gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+"); + AliXRDPROOFtoolkit tool; + chain = tool.MakeChainRandom("laserScan.txt","Mean",0,10200); + chain->Lookup(); + chainRef = tool.MakeChain("laserScanRef.txt","Mean",0,10200); + chainRef->Lookup(); + chain->AddFriend(chainRef,"Rr"); +} + + + +void MakeAliases(){ + // + // shortcuts for variables to fit + // + // bserved distrotions + chain->SetAlias("dy","(dY.fElements-Rr.dY.fElements)"); // y ~ (r-phi) distortion + chain->SetAlias("ctany","(dY.fElements-Rr.dY.fElements)/(250*dr)"); // mean distortion (angle) over drift length + // + chain->SetAlias("dr","(250.-abs(LTr.fP[1]))"); // drift length 0-250 + chain->SetAlias("dr1","(1.-abs(LTr.fP[1]/250.))"); // drift length 0-250 + chain->SetAlias("phiL","(LTr.fVecLY.fElements/LTr.fVecLX.fElements)"); // tann of local phi + chain->SetAlias("Esign","(1.-LTr.fSide*2.)"); // E field sign: +1 for A side : -1 for c side + // + // Br and Brfi + // + chain->SetAlias("ablx","(iblx.fElements/(ibz.fElements))"); + chain->SetAlias("ably","(ibly.fElements/(ibz.fElements))"); + chain->SetAlias("abr","(ibr.fElements/(ibz.fElements))"); + chain->SetAlias("abrf","(ibrphi.fElements/(ibz.fElements))"); + + + chain->SetAlias("r","(R.fElements-165.5)/80.3"); // local X - 0 at middle -1 first +1 last padrow + chain->SetAlias("ky","(kY.fElements)"); // local inclination angle + chain->SetAlias("sec","LTr.fVecSec.fElements"); // sector number + chain->SetAlias("ca","cos(LTr.fVecPhi.fElements+0.)"); // cos of global phi position + chain->SetAlias("sa","sin(LTr.fVecPhi.fElements+0.)"); // sin of global phi position +} + + +TMatrixD * MakeErrVector(TMatrixD & mat){ + // + // get error vector + // + Int_t nrows=mat.GetNrows(); + TMatrixD *err = new TMatrixD(nrows,1); + for (Int_t i=0; iSetAlias(aName+"_FB",strFit->Data()); + cutF=cutR+cutI + Form("abs(dy-%s)<0.2",(aName+"_FB").Data()); + // + // + // + strFit = TStatToolkit::FitPlane(chain,"dy:0.1", fstringB.Data(),cutF, chi2,npoints,fitParamB[i],covarB[i],1,0, ntracks); + chain->SetAlias(aName+"_FB",strFit->Data()); + chi2B[i]=TMath::Sqrt(chi2/npoints); + printf("B fit sqrt(Chi2/npoints)=%f\n",chi2B[i]); + // + strFit = TStatToolkit::FitPlane(chain,"dy:0.1", fstringBT.Data(),cutF, chi2,npoints,fitParamBT[i],covarBT[i],1,0, ntracks); + chain->SetAlias(aName+"_FBT",strFit->Data()); + chi2BT[i]=TMath::Sqrt(chi2/npoints); + printf("BT fit sqrt(Chi2/npoints)=%f\n",chi2BT[i]); + // + strFit = TStatToolkit::FitPlane(chain,"dy:0.1", fstringBT0.Data(),cutF, chi2,npoints,fitParamBT0[i],covarBT0[i],1,0, ntracks); + chain->SetAlias(aName+"_FBT0",strFit->Data()); + chi2BT0[i]=TMath::Sqrt(chi2/npoints); + printf("BT0 Fit: sqrt(Chi2/npoints)=%f\n",chi2BT0[i]); + // + // + strFit = TStatToolkit::FitPlane(chain,"dy:0.1", fstringA.Data(),cutF, chi2,npoints,fitParamA[i],covarA[i],1, 0, ntracks); + chain->SetAlias(aName+"_FA",strFit->Data()); + chi2A[i]=TMath::Sqrt(chi2/npoints); + printf("A fit: sqrt(Chi2/npoints)=%f\n",chi2A[i]); + errB[i] =(MakeErrVector(covarB[i])); + errBT[i] =(MakeErrVector(covarBT[i])); + errBT0[i] =(MakeErrVector(covarBT0[i])); + errA[i] =(MakeErrVector(covarA[i])); +} + + + + + +void MakeFitDy(){ + // + // + // + + + MakeFit(0,cutAM5,"dyAM5"); + MakeFit(1,cutAP2,"dyAP2"); + MakeFit(2,cutAP5,"dyAP5"); + // + MakeFit(3,cutCM5,"dyCM5"); + MakeFit(4,cutCP2,"dyCP2"); + // MakeFit(5,cutCP5,"dyAP5"); + DumpFit(); +} + + +void DrawPhi(){ + // + // + // + chain->Draw("(dy-dyAM5_FA):LTr.fVecPhi.fElements>>hisAM5(60,-3.14,3.14,100,-0.2,0.2)",cutAM5,""); + chain->Draw("(dy-dyAP5_FA):LTr.fVecPhi.fElements>>hisAP5(60,-3.14,3.14,100,-0.2,0.2)",cutAP5,""); + chain->Draw("(dy-dyAP2_FA):LTr.fVecPhi.fElements>>hisAP2(60,-3.14,3.14,100,-0.2,0.2)",cutAP2,""); + hisAM5->FitSlicesY(0,0,-1,20); + hisAP5->FitSlicesY(0,0,-1,20); + hisAP2->FitSlicesY(0,0,-1,20); + hisAM5_1->SetMinimum(-0.2); + hisAM5_1->SetMaximum(0.2); + hisAM5_1->SetMarkerStyle(20); + hisAP5_1->SetMarkerStyle(21); + hisAP2_1->SetMarkerStyle(22); + hisAM5_1->SetMarkerColor(1); + hisAP5_1->SetMarkerColor(2); + hisAP2_1->SetMarkerColor(4); + hisAM5_1->Draw(); + hisAP5_1->Draw("same"); + hisAP2_1->Draw("same"); +} + +void MakeGraphs(){ + Double_t bz[3] ={-5,2,5}; + Double_t p0A[3]={fitParam[0][0],fitParam[1][0],fitParam[2][0]}; + Double_t p1A[3]={fitParam[0][1],fitParam[1][1],fitParam[2][1]}; + Double_t p2A[3]={fitParam[0][2],fitParam[1][2],fitParam[2][2]}; + Double_t p3A[3]={fitParam[0][3],fitParam[1][3],fitParam[2][3]}; + Double_t p4A[3]={fitParam[0][4],fitParam[1][4],fitParam[2][4]}; + Double_t p5A[3]={fitParam[0][5],fitParam[1][5],fitParam[2][5]}; + Double_t p6A[3]={fitParam[0][6],fitParam[1][6],fitParam[2][6]}; + Double_t p7A[3]={fitParam[0][7],fitParam[1][7],fitParam[2][7]}; + Double_t p8A[3]={fitParam[0][8],fitParam[1][8],fitParam[2][8]}; + Double_t p9A[3]={fitParam[0][9],fitParam[1][9],fitParam[2][9]}; + TGraph grA0(3,bz,p0A); + TGraph grA1(3,bz,p1A); + TGraph grA2(3,bz,p2A); + TGraph grA3(3,bz,p3A); + TGraph grA4(3,bz,p4A); + TGraph grA5(3,bz,p5A); + TGraph grA6(3,bz,p6A); + TGraph grA7(3,bz,p7A); + TGraph grA8(3,bz,p8A); + TGraph grA9(3,bz,p9A); + TF1 f1("f1","[0]*x/(1+([0]*x)^2)"); + TF1 f2("f2","([0]*x)^2/(1+([0]*x)^2)"); + TMatrixD matB(5,5); + TMatrixD matBT(6,5); + TMatrixD matBT0(6,5); + for (Int_t i=0; i<5;i++) for (Int_t j=0; j<5;j++) matB[i][j]=fitParamB[j][i]; + for (Int_t i=0; i<5;i++) for (Int_t j=0; j<5;j++) matBT[i][j]=fitParamBT[j][i]; + for (Int_t i=0; i<6;i++) for (Int_t j=0; j<5;j++) matB0[i][j]=fitParamBT0[j][i]; +} + + +void DumpFit(){ + // + // + // + TTreeSRedirector *pcstream = new TTreeSRedirector("exbFits.root"); + + for (Int_t i=0; i<5;i++){ + Double_t bz=0; + Double_t side=0; + if (i==0) { bz=-5; side=1;} + if (i==1) { bz=2; side=1;} + if (i==2) { bz=5; side=1;} + if (i==3) { bz=-5; side=-1;} + if (i==4) { bz=2; side=-1;} + (*pcstream)<<"fit"<< + "bz="<