TFile f("CalibObjects.root")
AliTPCcalibPID *cal = f->Get("TPCCalib")->FindObject("calibPID")
-cal->GetHistLongMediumRatio()->Projection(0)->Fit("gaus")
-cal->GetHistShortMediumRatio()->Projection(0)->Fit("gaus")
+cal->GetHistRatioMaxTot()->Projection(0)->Fit("gaus")
+cal->GetHistRatioTracklet()->Projection(0)->Fit("gaus")
1.b) Update OCDB:
.x $ALICE_ROOT/TPC/macros/ConfigOCDB.C
AliTPCClusterParam * parcl = AliTPCcalibDB::Instance()->GetClusterParam();
-(*parcl->fQpadMnorm)[ipad] = oldvalue*corrFactor
+*parcl->fQpadMnorm)[ipad] = oldvalue*corrFactor
Int_t runNumber = 0;
AliCDBMetaData *metaData= new AliCDBMetaData();
fLandau(0),
fDeDxQmax(0),
fDeDxQtot(0),
- fDeDxRatio(0),
- fDeDxShortMediumRatio(0),
- fDeDxLongMediumRatio(0)
+ fDeDxRatioMaxTot(0),
+ fDeDxRatioQmax(0),
+ fDeDxRatioQtot(0),
+ fDeDxRatioTruncQtot(0),
+ fDeDxRatioTruncQmax(0)
{
//
// Empty default cosntructor
fLandau(0),
fDeDxQmax(0),
fDeDxQtot(0),
- fDeDxRatio(0),
- fDeDxShortMediumRatio(0),
- fDeDxLongMediumRatio(0)
+ fDeDxRatioMaxTot(0),
+ fDeDxRatioQmax(0),
+ fDeDxRatioQtot(0) ,
+ fDeDxRatioTruncQtot(0),
+ fDeDxRatioTruncQmax(0)
{
//
//
fUpperTrunc = 0.6;
//
fUseShapeNorm = kTRUE;
- fUsePosNorm = kFALSE;
+ fUsePosNorm = 0;
fUsePadNorm = 1;
//
fIsCosmic = kTRUE;
//
- // dE/dx, z, phi, theta, p, bg, ncls
- Int_t binsQA[7] = {150, 10, 10, 10, 50, 1, 8};
- Double_t xminQA[7] = {0., 0, 0, 0, 0.01, 0.1, 60};
- Double_t xmaxQA[7] = {10., 1, 0.6, 1.5, 50, 500, 160};
+ // dE/dx, z, phi, theta, p, bg, ncls, tracklet type
+ Int_t binsQA[8] = {150, 10, 10, 10, 50, 50, 8, 5};
+ Double_t xminQA[8] = {0., 0, 0, 0, 0.01, 0.1, 60, 0};
+ Double_t xmaxQA[8] = {10., 1, 0.6, 1.5, 50, 500, 160, 5};
//
- Double_t xminRa[7] = {0.5, 0, 0, 0, 0.01, 0.1, 60};
- Double_t xmaxRa[7] = { 2., 1, 0.6, 1.5, 50, 500, 160};
+ //
+ //
+ // dE/dx, z, phi, theta, dEdx, dEdx*dl, ncls, tracklet type
+ Int_t binsRA[9] = {100, 10, 10, 10, 25, 25, 4, 5};
+ Double_t xminRa[9] = {0.5, 0, 0, 0, 0.2, 0.2, 60, 0};
+ Double_t xmaxRa[9] = {4.0, 1, 0.6, 1.5, 5, 5, 160, 5};
// z,sin(phi),tan(theta),p,betaGamma,ncls
- fDeDxQmax = new THnSparseS("fDeDxQmax","Qmax;(z,sin(phi),tan(theta),p,betaGamma,ncls); TPC signal Qmax (a.u.)",7,binsQA,xminQA,xmaxQA);
- fDeDxQtot = new THnSparseS("fDeDxQtot","Qmax;(z,sin(phi),tan(theta),p,betaGamma,ncls); TPC signal Qmax (a.u.)",7,binsQA,xminQA,xmaxQA);
- fDeDxRatio = new THnSparseS("fDeDxRatio","Qmax;(z,sin(phi),tan(theta),p,betaGamma,ncls); TPC signal Qmax (a.u.)",7,binsQA,xminRa,xmaxRa);
- fDeDxShortMediumRatio = new THnSparseS("fDeDxQmax","Qmax;(z,sin(phi),tan(theta),p,betaGamma,ncls); TPC signal Qmax (a.u.)",7,binsQA,xminRa,xmaxRa);
- fDeDxLongMediumRatio = new THnSparseS("fDeDxQmax","Qmax;(z,sin(phi),tan(theta),p,betaGamma,ncls); TPC signal Qmax (a.u.)",7,binsQA,xminRa,xmaxRa);
+ fDeDxQmax = new THnSparseS("fDeDxQmax","Qmax;(z,sin(phi),tan(theta),p,betaGamma,ncls,type); TPC signal Qmax (a.u.)",8,binsQA,xminQA,xmaxQA);
+ fDeDxQtot = new THnSparseS("fDeDxQtot","Q_{tot};(z,sin(phi),tan(theta),p,betaGamma,ncls,type); TPC signal Qmax (a.u.)",8,binsQA,xminQA,xmaxQA);
+ //
+ // ratio histograms
+ //
+ fDeDxRatioMaxTot = new THnSparseS("fDeDxRatioMaxTot","Q_{max}/Q_{tot};(z,sin(phi),tan(theta),dEdx,dEdx*dl,ncls,type); TPC signal Qmax/Qtot (a.u.)",8,binsRA,xminRa,xmaxRa);
+ fDeDxRatioQmax = new THnSparseS("fDeDxRatioQmax","Q_{mtracklet}/Q_{mtrack};(z,sin(phi),tan(theta),dEdx,dEdx*dl,ncls,type,qtupe); TPC signal Tracklet/Track (a.u.)",8,binsRA,xminRa,xmaxRa);
+ fDeDxRatioQtot = new THnSparseS("fDeDxRatioQtot","Q_{ttracklet}/Q_{ttrack};(z,sin(phi),tan(theta),dEdx,dEdx*dl,ncls,type,qtupe); TPC signal Tracklet/Track (a.u.)",8,binsRA,xminRa,xmaxRa);
+ fDeDxRatioTruncQmax = new THnSparseS("fDeDxRatioTrunQmax","Q_{max}/Q_{maxtrunc};(z,sin(phi),tan(theta),dEdx,dEdx*dl,ncls,type,qtupe); TPC signal Full/Trunc. (a.u.)",8,binsRA,xminRa,xmaxRa);
+ fDeDxRatioTruncQtot = new THnSparseS("fDeDxRatioTruncQtot","Q_{tot}/Q_{tottrunc};(z,sin(phi),tan(theta),dEdx,dEdx*dl,ncls,type,qtupe); TPC signal Full/Trunc (a.u.)",8,binsRA,xminRa,xmaxRa);
+
+
BinLogX(fDeDxQmax,4); BinLogX(fDeDxQmax,5);
BinLogX(fDeDxQtot,4); BinLogX(fDeDxQtot,5);
- BinLogX(fDeDxRatio,4); BinLogX(fDeDxRatio,5);
- BinLogX(fDeDxShortMediumRatio,4); BinLogX(fDeDxShortMediumRatio,5);
- BinLogX(fDeDxLongMediumRatio,4); BinLogX(fDeDxLongMediumRatio,5);
+ //
+ BinLogX(fDeDxRatioMaxTot,4); BinLogX(fDeDxRatioMaxTot,5);
+ BinLogX(fDeDxRatioQmax,4); BinLogX(fDeDxRatioQmax,5);
+ BinLogX(fDeDxRatioQtot,4); BinLogX(fDeDxRatioQtot,5);
+ BinLogX(fDeDxRatioTruncQmax,4); BinLogX(fDeDxRatioTruncQmax,5);
+ BinLogX(fDeDxRatioTruncQtot,4); BinLogX(fDeDxRatioTruncQtot,5);
//
fHistNTracks = new TH1F("ntracks","Number of Tracks per Event; number of tracks per event; number of tracks",1001,-0.5,1000.5);
fClusters = new TH1F("signal","Number of Clusters per track; number of clusters per track n_{cl}; counts",40,0,160);
//
//
//
+
}
Printf("ERROR: ESD not available");
return;
}
-
+ const Int_t kMinClustersTracklet = 25;
Int_t ntracks=event->GetNumberOfTracks();
fHistNTracks->Fill(ntracks);
// calculate necessary track parameters
//Double_t meanP = 0.5*(trackIn->GetP() + trackOut->GetP());
Double_t meanP = trackIn->GetP();
- Double_t d = trackIn->GetLinearD(0,0);
+ //Double_t d = trackIn->GetLinearD(0,0);
Int_t NclsDeDx = track->GetTPCNcls();
//if (meanP > 0.7 || meanP < 0.2) continue;
if (meanP > 30 && TMath::Abs(trackIn->GetSnp()) < 0.2) fClusters->Fill(track->GetTPCNcls());
// calculate dEdx
// (Double_t low, Double_t up, Int_t type, Int_t i1, Int_t i2, Bool_t shapeNorm,Bool_t posNorm, Int_t padNorm, Int_t returnVal)
- Double_t TPCsignalTot = (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,0,0,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
+ //Double_t TPCsignalTot = (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,0,0,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
Double_t TPCsignalMax = (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,0,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
//
- Double_t TPCsignalShort = seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,0,63,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
- Double_t TPCsignalMedium = seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,63,63+64,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
- Double_t TPCsignalLong = seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,63+64,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
+ //Double_t TPCsignalShort = seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,0,63,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
+ //Double_t TPCsignalMedium = seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,63,63+64,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
+ //Double_t TPCsignalLong = seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,63+64,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
//
//Double_t driftMismatch = seed->GetDriftTimeMismatch(0,0.95,0.5);
Double_t driftMismatch = 0;
- Double_t drift = 1 - (TMath::Abs(trackIn->GetZ()) + TMath::Abs(trackOut->GetZ()))/500.;
+ // Double_t drift = 1 - (TMath::Abs(trackIn->GetZ()) + TMath::Abs(trackOut->GetZ()))/500.;
// particle identification
Double_t mass = 0.105658369;// muon
fLandau->Fill(0.1,0.6);
}
//var. of interest, z,sin(phi),tan(theta),p,betaGamma,ncls
- Double_t snp = TMath::Abs(trackIn->GetSnp());
- Double_t tgl = TMath::Abs(trackIn->GetTgl());
+ Double_t snpIn = TMath::Abs(trackIn->GetSnp());
+ Double_t snpOut = TMath::Abs(trackOut->GetSnp());
+ Double_t tglIn = TMath::Abs(trackIn->GetTgl());
+ Double_t tglOut = TMath::Abs(trackOut->GetTgl());
+ Double_t driftIn = 1 - (TMath::Abs(trackIn->GetZ()))/250.;
+ Double_t driftOut = 1 - (TMath::Abs(trackIn->GetZ()))/250.;
//
- Double_t vecQmax[7] = {TPCsignalMax,drift,snp,tgl,meanP,meanP/mass,NclsDeDx};
- Double_t vecQtot[7] = {TPCsignalTot,drift,snp,tgl,meanP,meanP/mass,NclsDeDx};
- Double_t vecRatio[7] = {TPCsignalMax/TPCsignalTot,drift,snp,tgl,meanP,meanP/mass,NclsDeDx};
- Double_t vecShortMediumRatio[7] = {TPCsignalShort/TPCsignalMedium,drift,snp,tgl,meanP,meanP/mass,NclsDeDx};
- Double_t vecLongMediumRatio[7] = {TPCsignalLong/TPCsignalMedium,drift,snp,tgl,meanP,meanP/mass,NclsDeDx};
-
- fDeDxQmax->Fill(vecQmax);
- fDeDxQtot->Fill(vecQtot);
- fDeDxRatio->Fill(vecRatio);
- fDeDxShortMediumRatio->Fill(vecShortMediumRatio);
- fDeDxLongMediumRatio->Fill(vecLongMediumRatio);
-
- }
-
- }
-
-
+ // dEdx in region
+ //
+ Float_t dEdxTot[5],dEdxTotFull[5];
+ Float_t dEdxMax[5],dEdxMaxFull[5];
+ Int_t ncl[5];
+ for (Int_t itype=0; itype<5;itype++){
+ Int_t row0=0, row1 =159;
+ if (itype==1) {row0=0; row1 = 63;};
+ if (itype==2) {row0= 64; row1=63+64;}
+ if (itype==3) {row0= 63+64; row1=159;}
+ dEdxTot[itype]= (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,0,row0,row1,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
+ dEdxMax[itype]= (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,row0,row1,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
+ // non trucated dEdx
+ dEdxTotFull[itype]= (1/fMIP)*seed->CookdEdxNorm(0.0,0.99,0,row0,row1,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
+ dEdxMaxFull[itype]= (1/fMIP)*seed->CookdEdxNorm(0.0,0.99,1,row0,row1,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
+
+
+ ncl[itype]=TMath::Nint(seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,row0,row1,fUseShapeNorm,fUsePosNorm,fUsePadNorm,2));
+ }
+ //
+ //
+ //
+ Float_t wmeanTot=0, wmeanMax=0, sumW=0;
+ Double_t length[3] = {0.75,1,1.5};
+
+ for (Int_t ipad=0; ipad<3; ipad++){
+ if (ncl[1+ipad]<3) continue;
+ Double_t weight = Double_t(ncl[1+ipad])*TMath::Sqrt(length[ipad]);
+ wmeanTot+=dEdxTot[1+ipad]*weight;
+ wmeanMax+=dEdxMax[1+ipad]*weight;
+ sumW+=weight;
+ }
+ if (sumW>0){
+ dEdxTot[4]= wmeanTot/sumW;
+ dEdxMax[4]= wmeanMax/sumW;
+ }
+ for (Int_t itype=0;itype<5;itype++){
+ //
+ Float_t snp=(TMath::Abs(snpIn)+TMath::Abs(snpOut))*0.5;
+ Float_t tgl=(TMath::Abs(tglIn)+TMath::Abs(tglOut))*0.5;
+ Float_t drift = (driftIn+driftOut)*0.5;
+ if (itype==1) {snp = TMath::Abs(snpIn); tgl = TMath::Abs(tglIn); drift= driftIn;};
+ if (itype==3) {snp = TMath::Abs(snpOut); tgl = TMath::Abs(tglOut);drift=driftOut;};
+ if (ncl[itype]<kMinClustersTracklet) continue;
+ Float_t deltaL = TMath::Sqrt(1+snp*snp+tgl*tgl);
+ //
+ Double_t vecQmax[8] = {dEdxMax[itype],drift,snp,tgl,meanP,meanP/mass,NclsDeDx, itype};
+ Double_t vecQtot[8] = {dEdxTot[itype],drift,snp,tgl,meanP,meanP/mass,NclsDeDx, itype};
+ //
+ //
+ //
+ Double_t ratioMaxTot = (dEdxTot[itype]>0) ? dEdxMax[itype]/dEdxTot[itype]:0;
+ Double_t ratioTrackletTot = (dEdxTot[0]>0) ? dEdxTot[itype]/dEdxTot[0]:0;
+ Double_t ratioTrackletMax = (dEdxMax[0]>0) ? dEdxMax[itype]/dEdxMax[0]:0;
+ Double_t ratioTrackletTruncTot = (dEdxTot[0]>0) ? dEdxTotFull[itype]/dEdxTot[itype]:0;
+ Double_t ratioTrackletTruncMax = (dEdxMax[0]>0) ? dEdxMaxFull[itype]/dEdxMax[itype]:0;
+
+ Double_t vecRatioMaxTot[8] = {ratioMaxTot, drift,snp,tgl,dEdxTot[0], dEdxTot[0]*deltaL,NclsDeDx,itype};
+ Double_t vecRatioTrackletTot[8] = {ratioTrackletTot, drift,snp,tgl,dEdxTot[0], dEdxTot[0]*deltaL,NclsDeDx,itype};
+ Double_t vecRatioTrackletMax[8] = {ratioTrackletMax, drift,snp,tgl,dEdxMax[0], dEdxMax[0]*deltaL,NclsDeDx,itype};
+ Double_t vecRatioTrackletTruncTot[8] = {ratioTrackletTruncTot, drift,snp,tgl,dEdxTot[0], dEdxTot[0]*deltaL,NclsDeDx,itype};
+ Double_t vecRatioTrackletTruncMax[8] = {ratioTrackletTruncMax, drift,snp,tgl,dEdxMax[0], dEdxMax[0]*deltaL,NclsDeDx,itype};
+ fDeDxQmax->Fill(vecQmax);
+ fDeDxQtot->Fill(vecQtot);
+ fDeDxRatioMaxTot->Fill(vecRatioMaxTot);
+ fDeDxRatioQmax->Fill(vecRatioTrackletTot);
+ fDeDxRatioQtot->Fill(vecRatioTrackletMax);
+ fDeDxRatioTruncQmax->Fill(vecRatioTrackletTruncTot);
+ fDeDxRatioTruncQtot->Fill(vecRatioTrackletTruncMax);
+ //
+ TTreeSRedirector * cstream = GetDebugStreamer();
+ if (cstream){
+ TVectorD vQT(9,vecQtot);
+ TVectorD vQM(9,vecQmax);
+ TVectorD vQMT(9, vecRatioMaxTot);
+ TVectorD vQRT(9,vecRatioTrackletTot);
+ TVectorD vQRM(9,vecRatioTrackletMax);
+ TVectorD vQRTT(9,vecRatioTrackletTruncTot);
+ TVectorD vQRTM(9,vecRatioTrackletTruncMax);
+ (*cstream) << "dEdx" <<
+ "run="<<fRun<< // run number
+ "event="<<fEvent<< // event number
+ "time="<<fTime<< // time stamp of event
+ "trigger="<<fTrigger<< // trigger
+ "mag="<<fMagF<< // magnetic field
+ "track.="<<seed<< // original seed
+ //
+ "vQT.="<<&vQT<< // trunc mean - total charge
+ "vQM.="<<&vQM<< // trunc mean - max charge
+ //
+ "vQMT.="<<&vQMT<< // ratio max/tot
+ "vQRT.="<<&vQRT<< // ratio tracklet/track - total charge
+ "vQRM.="<<&vQRM<< // ratio tracklet/track - max charge
+ "vQRTT.="<<&vQRTT<< // ratio trunc/full - total charge
+ "vQRTM.="<<&vQRTM<< // ratio trunc/full - max charge
+ "\n";
+ }
+ }
+ }
+ }
}
//
if (cal->GetHistQmax()) fDeDxQmax->Add(cal->GetHistQmax());
if (cal->GetHistQtot()) fDeDxQtot->Add(cal->GetHistQtot());
- if (cal->GetHistRatio()) fDeDxRatio->Add(cal->GetHistRatio());
- if (cal->GetHistShortMediumRatio()) fDeDxShortMediumRatio->Add(cal->GetHistShortMediumRatio());
- if (cal->GetHistLongMediumRatio()) fDeDxLongMediumRatio->Add(cal->GetHistLongMediumRatio());
+ if (cal->GetHistRatioMaxTot()) fDeDxRatioMaxTot->Add(cal->GetHistRatioMaxTot());
+ if (cal->GetHistRatioQmax()) fDeDxRatioQmax->Add(cal->GetHistRatioQmax());
+ if (cal->GetHistRatioQtot()) fDeDxRatioQtot->Add(cal->GetHistRatioQtot());
+ if (cal->GetHistRatioTruncQmax()) fDeDxRatioTruncQmax->Add(cal->GetHistRatioTruncQmax());
+ if (cal->GetHistRatioTruncQtot()) fDeDxRatioTruncQtot->Add(cal->GetHistRatioTruncQtot());
}
return 0;
-
-void AliTPCcalibPID::MakeReport() {
- //
- // 1. standard dEdx picture
- TCanvas *cDeDx = new TCanvas("cDeDx","cDeDx");
- GetHistQmax()->Projection(0,4)->Draw("colz");
- return;
-}
-
-
-
void AliTPCcalibPID::BinLogX(THnSparse *h, Int_t axisDim) {
// Method for the correct logarithmic binning of histograms
}
+
+
+
+
+void AliTPCcalibPID::MakeReport(const char *outputpath) {
+ //
+ // Make a standard QA plots
+ //
+ for (Int_t ipad=0;ipad<4;ipad++){
+ DrawRatioTot(ipad,outputpath);
+ DrawRatioMax(ipad,outputpath);
+ }
+ DrawRatiodEdx(0.5,3,outputpath);
+ DrawResolBGQtot(140,160,outputpath);
+ DrawResolBGQmax(140,160,outputpath);
+ return;
+}
+
+void AliTPCcalibPID::DrawRatioTot(Int_t ipad, const char* outputpath){
+ //
+ // Draw default ratio histogram for given pad type
+ // ipad - 0 - short pads
+ // 1 - medium pads
+ // 2 - long pads
+ //
+ // Int_t kmicolors[10]={1,2,3,4,6,7,8,9,10,11};
+ Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
+ AliTPCcalibPID * pid = this;
+ TCanvas *canvas= new TCanvas(Form("QtotRatio%d)",ipad),Form("QtotRatio%d)",ipad));
+ canvas->Divide(3,2);
+ pid->GetHistRatioQtot()->GetAxis(7)->SetRange(ipad+2,ipad+2);
+ canvas->cd(1);
+ TH1 * his0 = pid->GetHistRatioQtot()->Projection(0);
+ his0->SetXTitle("dEdx_{tracklet}/dEdx_{track}");
+ his0->SetYTitle("");
+ his0->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
+ his0->Draw();
+ //
+ canvas->cd(2);
+ TH1 * hprofz = (TH1*) (pid->GetHistRatioQtot()->Projection(0,1)->ProfileX());
+ hprofz->SetXTitle("drift length");
+ hprofz->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
+ hprofz->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
+ hprofz->SetMarkerStyle(kmimarkers[0]);
+ hprofz->SetMaximum(1.1);
+ hprofz->SetMinimum(0.9);
+ hprofz->Draw();
+ //
+ canvas->cd(3);
+ TH1 * hprofphi = (TH1*) (pid->GetHistRatioQtot()->Projection(0,2)->ProfileX());
+ hprofphi->SetXTitle("sin(#phi)");
+ hprofphi->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
+ hprofphi->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
+ hprofphi->SetMarkerStyle(kmimarkers[0]);
+ hprofz->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
+ hprofphi->SetMaximum(1.1);
+ hprofphi->SetMinimum(0.9);
+ hprofphi->Draw();
+ //
+ canvas->cd(4);
+ TH1 * hproftheta = (TH1*) (pid->GetHistRatioQtot()->Projection(0,3)->ProfileX());
+ hproftheta->SetXTitle("tan(#theta)");
+ hproftheta->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
+ hproftheta->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
+ hproftheta->SetMarkerStyle(kmimarkers[0]);
+ hproftheta->SetMaximum(1.1);
+ hproftheta->SetMinimum(0.9);
+ hproftheta->Draw();
+ //
+ canvas->cd(5);
+ TH1 * hprofdedx = (TH1*) (pid->GetHistRatioQtot()->Projection(0,4)->ProfileX());
+ hprofdedx->SetXTitle("dEdx_{track}");
+ hprofdedx->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
+ hprofdedx->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
+ hprofdedx->SetMarkerStyle(kmimarkers[0]);
+ hprofdedx->SetMaximum(1.1);
+ hprofdedx->SetMinimum(0.9);
+ hprofdedx->Draw();
+
+ canvas->cd(6);
+ TH1 * hprofdedxL = (TH1*) (pid->GetHistRatioQtot()->Projection(0,5)->ProfileX());
+ hprofdedxL->SetXTitle("dEdx_{track}#Delta_{x}");
+ hprofdedxL->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
+ hprofdedxL->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
+ hprofdedxL->SetMarkerStyle(kmimarkers[0]);
+ hprofdedxL->SetMaximum(1.1);
+ hprofdedxL->SetMinimum(0.9);
+ hprofdedxL->Draw();
+
+
+
+ canvas->SaveAs(Form("%s/QtotRatioType%d.eps",outputpath,ipad));
+ canvas->SaveAs(Form("%s/QtotRatioType%d.gif",outputpath,ipad));
+}
+
+void AliTPCcalibPID::DrawRatioMax(Int_t ipad, const char* outputpath){
+ //
+ // Draw default ration histogram for given pad type
+ // ipad - 0 - short pads
+ // 1 - medium pads
+ // 2 - long pads
+ //
+ // Int_t kmicolors[10]={1,2,3,4,6,7,8,9,10,11};
+ Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
+ AliTPCcalibPID * pid = this;
+ TCanvas *canvas= new TCanvas(Form("QmaxRatio%d)",ipad),Form("QmaxRatio%d)",ipad));
+ canvas->Divide(3,2);
+ pid->GetHistRatioQmax()->GetAxis(7)->SetRange(ipad+2,ipad+2);
+ canvas->cd(1);
+ TH1 * his0 = pid->GetHistRatioQmax()->Projection(0);
+ his0->SetXTitle("dEdx_{tracklet}/dEdx_{track}");
+ his0->SetYTitle("");
+ his0->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
+ his0->Draw();
+ //
+ canvas->cd(2);
+ TH1 * hprofz = (TH1*) (pid->GetHistRatioQmax()->Projection(0,1)->ProfileX());
+ hprofz->SetXTitle("drift length");
+ hprofz->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
+ hprofz->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
+ hprofz->SetMarkerStyle(kmimarkers[0]);
+ hprofz->SetMaximum(1.1);
+ hprofz->SetMinimum(0.9);
+ hprofz->Draw();
+ //
+ canvas->cd(3);
+ TH1 * hprofphi = (TH1*) (pid->GetHistRatioQmax()->Projection(0,2)->ProfileX());
+ hprofphi->SetXTitle("sin(#phi)");
+ hprofphi->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
+ hprofphi->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
+ hprofphi->SetMarkerStyle(kmimarkers[0]);
+ hprofphi->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
+ hprofphi->SetMaximum(1.1);
+ hprofphi->SetMinimum(0.9);
+ hprofphi->Draw();
+ //
+ canvas->cd(4);
+ TH1 * hproftheta = (TH1*) (pid->GetHistRatioQmax()->Projection(0,3)->ProfileX());
+ hproftheta->SetXTitle("tan(#theta)");
+ hproftheta->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
+ hproftheta->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
+ hproftheta->SetMarkerStyle(kmimarkers[0]);
+ hproftheta->SetMaximum(1.1);
+ hproftheta->SetMinimum(0.9);
+ hproftheta->Draw();
+
+ canvas->cd(5);
+ TH1 * hprofdedx = (TH1*) (pid->GetHistRatioQmax()->Projection(0,4)->ProfileX());
+ hprofdedx->SetXTitle("dEdx_{track}");
+ hprofdedx->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
+ hprofdedx->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
+ hprofdedx->SetMarkerStyle(kmimarkers[0]);
+ hprofdedx->SetMaximum(1.1);
+ hprofdedx->SetMinimum(0.9);
+ hprofdedx->Draw();
+
+ canvas->cd(6);
+ TH1 * hprofdedxL = (TH1*) (pid->GetHistRatioQmax()->Projection(0,5)->ProfileX());
+ hprofdedxL->SetXTitle("dEdx_{track}#Delta_{x}");
+ hprofdedxL->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
+ hprofdedxL->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
+ hprofdedxL->SetMarkerStyle(kmimarkers[0]);
+ hprofdedxL->SetMaximum(1.1);
+ hprofdedxL->SetMinimum(0.9);
+ hprofdedxL->Draw();
+
+
+ canvas->SaveAs(Form("%s/QmaxRatioType%d.eps",outputpath,ipad));
+ canvas->SaveAs(Form("%s/QmaxRatioType%d.gif",outputpath,ipad));
+}
+
+
+
+void AliTPCcalibPID::DrawRatiodEdx(Float_t demin, Float_t demax, const char* outputpath){
+ //
+ //
+ //
+ //
+ // Int_t kmicolors[10]={1,2,3,4,6,7,8,9,10,11};
+ Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
+ AliTPCcalibPID * pid = this;
+ TCanvas *canvas= new TCanvas("QRatiodEdx","QRatiodEdx");
+ canvas->Divide(2,4);
+ canvas->SetLogx(kTRUE);
+ TH1 * hprofP=0;
+ for (Int_t ipad=0;ipad<4;ipad++){
+ pid->GetHistRatioQmax()->GetAxis(7)->SetRange(ipad+2,ipad+2);
+ pid->GetHistRatioQtot()->GetAxis(7)->SetRange(ipad+2,ipad+2);
+ pid->GetHistRatioQmax()->GetAxis(5)->SetRangeUser(demin,demax);
+ pid->GetHistRatioQtot()->GetAxis(5)->SetRangeUser(demin,demax);
+
+ canvas->cd(ipad*2+1)->SetLogx(kFALSE);
+ hprofP = (TH1*) (pid->GetHistRatioQmax()->Projection(0,5)->ProfileX());
+ hprofP->SetXTitle("dEdx_{track}");
+ hprofP->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
+ hprofP->SetTitle(Form("Q_{max} dEdx_{tracklet}/dEdx_{track} type %d",0));
+ hprofP->SetMarkerStyle(kmimarkers[0]);
+ hprofP->SetMaximum(1.1);
+ hprofP->SetMinimum(0.9);
+ hprofP->Draw();
+ // pad Tot
+ canvas->cd(ipad*2+2)->SetLogx(kFALSE);
+ hprofP = (TH1*) (pid->GetHistRatioQtot()->Projection(0,5)->ProfileX());
+ hprofP->SetXTitle("dEdx_{track}");
+ hprofP->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
+ hprofP->SetTitle(Form("Q_{tot} dEdx_{tracklet}/dEdx_{track} type %d",0));
+ hprofP->SetMarkerStyle(kmimarkers[0]);
+ hprofP->SetMaximum(1.1);
+ hprofP->SetMinimum(0.9);
+ hprofP->Draw();
+ }
+ //
+ //
+ canvas->SaveAs(Form("%s/QratiodEdx%d.eps",outputpath));
+ canvas->SaveAs(Form("%s/QratiodEdx%d.gif",outputpath));
+}
+
+
+
+void AliTPCcalibPID::DrawResolBGQtot(Int_t minClusters, Int_t maxClusters, const char *outputpath){
+ //
+ //
+ //
+ AliTPCcalibPID * pid = this;
+ Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
+
+ TObjArray arr;
+ TH2 * his =0;
+ TH1 * hmean =0;
+ TH1 * hsigma=0;
+ //
+ // set cut
+ pid->GetHistQtot()->GetAxis(6)->SetRangeUser(minClusters,maxClusters);
+ pid->GetHistQtot()->GetAxis(5)->SetRangeUser(1,10000);
+ TCanvas *canvas= new TCanvas("dEdxResolQ_{Tot}","dEdxResolQ_{Tot}");
+ canvas->Divide(2,3);
+ //
+ //
+ //
+ for (Int_t ipad=0;ipad<5;ipad++){
+ canvas->cd(1+ipad)->SetLogx(kTRUE);
+ if (ipad<4) pid->GetHistQtot()->GetAxis(7)->SetRange(ipad+2,ipad+2);
+ if (ipad==4) pid->GetHistQtot()->GetAxis(7)->SetRange(1,1);
+ his = (TH2*)(pid->GetHistQtot()->Projection(0,5));
+ his->FitSlicesY(0,0,-1,0,"QNR",&arr);
+ hmean = (TH1*)arr.At(1);
+ hsigma = (TH1*)arr.At(2)->Clone();
+ hsigma->Divide(hmean);
+ arr.SetOwner(kTRUE);
+ arr.Delete();
+ delete his;
+
+ hsigma->SetXTitle("#beta#gamma");
+ hsigma->SetYTitle("#sigma_{dEdx}/dEdx");
+ hsigma->SetTitle(Form("#sigma_{dEdx}/dEdx_{tot} Pad %d",ipad));
+ hsigma->SetName(Form("#sigma_{dEdx}/dEdx_{tot} Pad %d",ipad));
+ if (ipad==4) {
+ hsigma->SetTitle(Form("#sigma_{dEdx}/dEdx_{tot} Full"));
+ hsigma->SetName(Form("#sigma_{dEdx}/dEdx_{tot} Full"));
+ }
+ hsigma->SetMarkerStyle(kmimarkers[0]);
+ hsigma->SetMaximum(0.15);
+ hsigma->SetMinimum(0.0);
+ hsigma->Draw();
+ }
+ canvas->SaveAs(Form("%s/dEdxResolMax.eps",outputpath));
+ canvas->SaveAs(Form("%s/dEdxResolMax.gif",outputpath));
+}
+
+void AliTPCcalibPID::DrawResolBGQmax(Int_t minClusters, Int_t maxClusters, const char *outputpath){
+ //
+ // Int_t minClusters=140; Int_t maxClusters=200; const char *outputpath="picPID06"
+ //
+ Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
+ AliTPCcalibPID * pid = this;
+ TObjArray arr;
+ TH2 * his =0;
+ TH1 * hmean =0;
+ TH1 * hsigma=0;
+ //
+ // set cut
+ pid->GetHistQmax()->GetAxis(6)->SetRangeUser(minClusters,maxClusters);
+ pid->GetHistQmax()->GetAxis(5)->SetRangeUser(1,10000);
+ TCanvas *canvas= new TCanvas("dEdxResolQ_{Max}","dEdxResolQ_{Max}");
+ canvas->Divide(2,3);
+ //
+ //
+ //
+ for (Int_t ipad=0;ipad<5;ipad++){
+ canvas->cd(1+ipad)->SetLogx(kTRUE);
+ if (ipad<4) pid->GetHistQmax()->GetAxis(7)->SetRange(ipad+2,ipad+2);
+ if (ipad==4) pid->GetHistQmax()->GetAxis(7)->SetRange(1,1);
+ his = (TH2*)(pid->GetHistQmax()->Projection(0,5));
+ his->FitSlicesY(0,0,-1,0,"QNR",&arr);
+ hmean = (TH1*)arr.At(1);
+ hsigma = (TH1*)arr.At(2)->Clone();
+ hsigma->Divide(hmean);
+ arr.SetOwner(kTRUE);
+ arr.Delete();
+ delete his;
+ hsigma->SetXTitle("#beta#gamma");
+ hsigma->SetYTitle("#sigma_{dEdx}/dEdx");
+ hsigma->SetTitle(Form("#sigma_{dEdx}/dEdx_{max} Pad %d",ipad));
+ hsigma->SetName(Form("#sigma_{dEdx}/dEdx_{max} Pad %d",ipad));
+ if (ipad==4){
+ hsigma->SetTitle(Form("#sigma_{dEdx}/dEdx_{max} Full"));
+ hsigma->SetName(Form("#sigma_{dEdx}/dEdx_{max} Full"));
+ }
+ hsigma->SetMarkerStyle(kmimarkers[0]);
+ hsigma->SetMaximum(0.15);
+ hsigma->SetMinimum(0.0);
+ hsigma->Draw();
+ }
+ canvas->SaveAs(Form("%s/dEdxResolMax.eps",outputpath));
+ canvas->SaveAs(Form("%s/dEdxResolMax.gif",outputpath));
+}
+
+
+
+void AliTPCcalibPID::DumpTree(THnSparse * hndim, const char * outname){
+ //
+ // make a tree
+ // the tree will be written to the file - outname
+ //
+ TTreeSRedirector *pcstream = new TTreeSRedirector(outname);
+ //
+ //
+ Double_t position[10];
+ Double_t value;
+ Int_t *bins = new Int_t[10];
+ //
+ //
+ Int_t ndim = hndim->GetNdimensions();
+ //
+ for (Long64_t i = 0; i < hndim->GetNbins(); ++i) {
+ value = hndim->GetBinContent(i, bins);
+ for (Int_t idim = 0; idim < ndim; idim++) {
+ position[idim] = hndim->GetAxis(idim)->GetBinCenter(bins[idim]);
+ }
+ Double_t ty = TMath::Tan(TMath::ASin(position[2]));
+ Double_t tz = position[3]*TMath::Sqrt(1+ty*ty);
+ //
+ (*pcstream)<<"Dump"<<
+ "bincont="<<value<< // bin content
+ "val="<<position[0]<< // parameter difference
+ "dr="<<position[1]<< //drift
+ "ty="<<ty<< //phi
+ "tz="<<tz<< //theta
+ "p2="<<position[2]<< //p2
+ "p3="<<position[3]<< //p3
+ "p="<<position[4]<< //p
+ "bg="<<position[5]<< //bg
+ "ncl="<<position[6]<< //ncl
+ "type="<<position[7]<< //tracklet
+ "tot="<<position[8]<< //is tot
+ "\n";
+ }
+ delete pcstream;
+}
+
+
+void AliTPCcalibPID::DumpTrees(){
+ //
+ // dump the content of histogram to the tree
+ //
+ AliTPCcalibPID *pid = this;
+ DumpTree(pid->GetHistQtot(),"dumpQtot.root");
+ DumpTree(pid->GetHistQmax(),"dumpQmax.root");
+ DumpTree(pid->GetHistRatioQtot(),"dumpRatioQtot.root");
+ DumpTree(pid->GetHistRatioQmax(),"dumpRatioQmax.root");
+}
+
+
+
virtual void Process(AliESDEvent *event);
virtual Long64_t Merge(TCollection *li);
virtual void Analyze();
- void MakeReport();
- //
+ void MakeReport(const char * outputpath);
+ void DrawRatioTot(Int_t ipad, const char* outputpath);
+ void DrawRatioMax(Int_t ipad, const char* outputpath);
+ void DrawRatiodEdx(Float_t demin, Float_t demax, const char* outputpath);
+ void DrawResolBGQtot(Int_t minClusters, Int_t maxClusters, const char *outputpath); //
+ void DrawResolBGQmax(Int_t minClusters, Int_t maxClusters, const char *outputpath);
//
TH1F * GetHistNTracks(){return fHistNTracks;};
TH1F * GetHistClusters(){return fClusters;};
//
THnSparseS * GetHistQmax(){return fDeDxQmax;};
THnSparseS * GetHistQtot(){return fDeDxQtot;};
- THnSparseS * GetHistRatio(){return fDeDxRatio;};
- THnSparseS * GetHistShortMediumRatio(){return fDeDxShortMediumRatio;};
- THnSparseS * GetHistLongMediumRatio(){return fDeDxLongMediumRatio;};
+ THnSparseS * GetHistRatioMaxTot(){return fDeDxRatioMaxTot;};
+ THnSparseS * GetHistRatioQmax(){return fDeDxRatioQmax;};
+ THnSparseS * GetHistRatioQtot(){return fDeDxRatioQtot;};
+ THnSparseS * GetHistRatioTruncQmax(){return fDeDxRatioTruncQmax;};
+ THnSparseS * GetHistRatioTruncQtot(){return fDeDxRatioTruncQtot;};
//
void SetMIPvalue(Float_t mip){fMIP = mip;};
void SetLowerTrunc(Float_t lowerTrunc){fLowerTrunc = lowerTrunc;};
void SetUpperTrunc(Float_t upperTrunc){fUpperTrunc = upperTrunc;};
void SetUseShapeNorm(Bool_t useShapeNorm){fUseShapeNorm = useShapeNorm;};
- void SetUsePosNorm(Bool_t usePosNorm){fUsePosNorm = usePosNorm;};
+ void SetUsePosNorm(Int_t usePosNorm){fUsePosNorm = usePosNorm;};
void SetPadNorm(Int_t padNorm){fUsePadNorm = padNorm;};
void SetIsCosmic(Bool_t isCosmic){fIsCosmic = isCosmic;};
//
//
static void BinLogX(THnSparse * h, Int_t axisDim); // method for correct histogram binning
-
-
+ void DumpTree(THnSparse * hndim, const char * outname);
+ void DumpTrees();
+ void Process(AliESDtrack *track, Int_t runNo=-1){AliTPCcalibBase::Process(track,runNo);};
+ void Process(AliTPCseed *track){return AliTPCcalibBase::Process(track);}
private:
//
// parameter specifications
Float_t fLowerTrunc;
Float_t fUpperTrunc;
Bool_t fUseShapeNorm;
- Bool_t fUsePosNorm;
+ Int_t fUsePosNorm;
Int_t fUsePadNorm;
//
Bool_t fIsCosmic;
//
THnSparseS * fDeDxQmax; // histogram which shows dEdx (Qmax) as a function of z,sin(phi),tan(theta),p,betaGamma
THnSparseS * fDeDxQtot; // histogram which shows dEdx (Qtot) as a function of z,sin(phi),tan(theta),p,betaGamma
- THnSparseS * fDeDxRatio; // histogram which shows dEdx ratio (Qmax/Qtot) as a function of z,sin(phi),tan(theta),p,betaGamma
- THnSparseS * fDeDxShortMediumRatio; // histogram which shows dEdx ratio (QmaxShort/QmaxMedium) as a function of z,sin(phi),tan(theta),p,betaGamma
- THnSparseS * fDeDxLongMediumRatio; // histogram which shows dEdx ratio (QmaxLong/QmaxMedium) as a function of z,sin(phi),tan(theta),p,betaGamma
+ //
+ // ratio histograms
+ //
+ THnSparseS * fDeDxRatioMaxTot; // histogram which shows dEdx ratio (Qmax/Qtot) as a function of z,sin(phi),tan(theta),dEdx,dEdx*dl
+ THnSparseS * fDeDxRatioQmax; // dEdx ratio (tracklet/track) as a function of z,sin(phi),tan(theta),dEdx,dEdx*dl
+ THnSparseS * fDeDxRatioQtot; // dEdx ratio (tracklet/track) as a function of z,sin(phi),tan(theta),dEdx,dEdx*dl
+ THnSparseS * fDeDxRatioTruncQtot; // dEdx ratio (tracklet/track) as a function of z,sin(phi),tan(theta),dEdx,dEdx*dl
+ THnSparseS * fDeDxRatioTruncQmax; // dEdx ratio (tracklet/track) as a function of z,sin(phi),tan(theta),dEdx,dEdx*dl
+
//
AliTPCcalibPID(const AliTPCcalibPID&);
AliTPCcalibPID& operator=(const AliTPCcalibPID&);
dy = +dr(r,z)*tan(phi)
-
- //
+ .x ~/NimStyle.C
+ gSystem->Load("libANALYSIS");
+ gSystem->Load("libTPCcalib");
+ TFile fcalib("CalibObjects.root");
+ TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
+ AliTPCcalibUnlinearity * calibUnlin = ( AliTPCcalibUnlinearity *)array->FindObject("calibUnlinearity");
+ //
+
*/
#include "AliTPCTracklet.h"
#include "TTimeStamp.h"
#include "AliTPCcalibDB.h"
-#include "AliTPCcalibLaser.h"
#include "AliDCSSensorArray.h"
#include "AliDCSSensor.h"
#include "TLinearFitter.h"
#include "AliTPCcalibUnlinearity.h"
+#include "AliTPCPointCorrection.h"
-AliTPCcalibUnlinearity* AliTPCcalibUnlinearity::fgInstance = 0;
ClassImp(AliTPCcalibUnlinearity)
AliTPCcalibUnlinearity::AliTPCcalibUnlinearity():
AliTPCcalibBase(),
- fDiffHistoLine(0),
- fDiffHistoPar(0),
+ fInit(kFALSE),
+ fDiffHistoLineY(0),
+ fDiffHistoParY(0),
+ fDiffHistoLineZ(0),
+ fDiffHistoParZ(0),
fFittersOutR(38),
fFittersOutZ(38),
fParamsOutR(38),
fParamsOutZ(38),
fErrorsOutR(38),
- fErrorsOutZ(38)
+ fErrorsOutZ(38),
+ fDistRPHIPlus(74),
+ fDistRPHIMinus(74),
+ fFitterQuadrantY(16*38),
+ fFitterQuadrantPhi(16*38)
{
//
// Defualt constructor
AliTPCcalibUnlinearity::AliTPCcalibUnlinearity(const Text_t *name, const Text_t *title):
AliTPCcalibBase(name,title),
- fDiffHistoLine(0),
- fDiffHistoPar(0),
+ fInit(kFALSE),
+ fDiffHistoLineY(0),
+ fDiffHistoParY(0),
+ fDiffHistoLineZ(0),
+ fDiffHistoParZ(0),
fFittersOutR(38),
fFittersOutZ(38),
fParamsOutR(38),
fParamsOutZ(38),
fErrorsOutR(38),
- fErrorsOutZ(38)
+ fErrorsOutZ(38),
+ fDistRPHIPlus(74),
+ fDistRPHIMinus(74),
+ fFitterQuadrantY(16*38),
+ fFitterQuadrantPhi(16*38)
{
//
// Non default constructor
//
- MakeFitters();
+ MakeHisto();
}
AliTPCcalibUnlinearity::~AliTPCcalibUnlinearity(){
//
//
//
- if (fDiffHistoLine) delete fDiffHistoLine;
- if (fDiffHistoPar) delete fDiffHistoPar;
+ if (fDiffHistoLineZ) delete fDiffHistoLineY;
+ if (fDiffHistoParZ) delete fDiffHistoParY;
+ if (fDiffHistoLineZ) delete fDiffHistoLineZ;
+ if (fDiffHistoParZ) delete fDiffHistoParZ;
}
-AliTPCcalibUnlinearity* AliTPCcalibUnlinearity::Instance()
-{
- //
- // Singleton implementation
- // Returns an instance of this class, it is created if neccessary
- //
- if (fgInstance == 0){
- fgInstance = new AliTPCcalibUnlinearity();
- }
- return fgInstance;
-}
-
//
//
//
+ if (!fInit) {
+ Init(); // work around for PROOF - initialize firs time used
+ }
const Int_t kdrow=10;
const Int_t kMinCluster=35;
- if (!fDiffHistoLine) MakeHisto();
if (!seed) return;
if (TMath::Abs(fMagF)>0.01) return; // use only data without field
Int_t counter[72];
//
for (Int_t is=0; is<72;is++){
if (counter[is]>kMinCluster) ProcessDiff(seed, is);
+ if (counter[is]>kMinCluster) {
+ AlignOROC(seed, is);
+ }
}
}
//
const Double_t kXmean=134;
const Int_t kdrow=10;
+ const Float_t kSagitaCut = 1;
static TLinearFitter fy1(2,"pol1");
static TLinearFitter fy2(3,"pol2");
static TLinearFitter fz1(2,"pol1");
AliTPCclusterMI *c=track->GetClusterPointer(irow);
if (!c) continue;
if (c->GetDetector()!=isec) continue;
- if (c->GetType()<0) continue;
Double_t dx = c->GetX()-kXmean;
Double_t y1 = py1[0]+py1[1]*dx;
Double_t y2 = py2[0]+py2[1]*dx+py2[2]*dx*dx;
Double_t z1 = pz1[0]+pz1[1]*dx;
Double_t z2 = pz2[0]+pz2[1]*dx+pz2[2]*dx*dx;
//
+ Double_t edgeMinus= -y1+c->GetX()*TMath::Tan(TMath::Pi()/18.);
+ Double_t edgePlus = y1+c->GetX()*TMath::Tan(TMath::Pi()/18.);
+ Int_t npads = AliTPCROC::Instance()->GetNPads(isec,c->GetRow());
+ Float_t dpad = TMath::Min(c->GetPad(), npads-1- c->GetPad());
+ Float_t dy =c->GetY()-y1;
+ //
+ // Corrections
+ //
+ AliTPCPointCorrection * corr = AliTPCPointCorrection::Instance();
+ Double_t corrclY=0, corrtrY=0, corrR=0;
+ corrclY = corr->RPhiCOGCorrection(isec,c->GetRow(), c->GetPad(),
+ c->GetY(),c->GetY(), c->GetZ(), py1[1],c->GetMax(),2.5);
+ corrtrY = corr->RPhiCOGCorrection(isec,c->GetRow(), c->GetPad(),
+ c->GetY(),y1, c->GetZ(), py1[1], c->GetMax(),2.5);
+ corrR = corr->CorrectionOutR0(kFALSE,kFALSE,c->GetX(),c->GetY(),c->GetZ(),isec);
+ //
//
- Double_t x[10];
- x[0]=isec;
- x[1]=irow;
- x[2]=c->GetY();
- x[3]=c->GetZ();
- x[4]=py1[1];
- x[5]=pz1[1];
- x[6]=py2[2]*150*150/4; // sagita 150 cm
//
- x[7]=c->GetY()-y1;
- x[8]=c->GetZ()-z1;
- fDiffHistoLine->Fill(x);
- x[7]=c->GetY()-y2;
- x[8]=c->GetZ()-z2;
- fDiffHistoPar->Fill(x);
-
- if (TMath::Abs(py2[2]*150*150/4)<1 && TMath::Abs(pz2[2]*150*150/4)<1){
- // Apply sagita cut
- AddPoint(isec,irow,c->GetZ()-z1, c->GetY()-y1,
- py1[1],pz1[1],1-TMath::Abs(c->GetZ()/250.),1);
- }
-
if (fStreamLevel>1){
TTreeSRedirector *cstream = GetDebugStreamer();
if (cstream){
"time="<<fTime<< // time stamp of event
"trigger="<<fTrigger<< // trigger
"mag="<<fMagF<< // magnetic field
- "isec="<<isec<<
+ "isec="<<isec<< // current sector
+ "npads="<<npads<< // number of pads at given sector
+ "dpad="<<dpad<< // distance to edge pad
+ //
+ "crR="<<corrR<< // radial correction
+ "cclY="<<corrclY<< // edge effect correction cl
+ "ctrY="<<corrtrY<< // edge effect correction using track
+ //
"Cl.="<<c<<
"y1="<<y1<<
"y2="<<y2<<
"py2.="<<&py2<<
"pz1.="<<&pz1<<
"pz2.="<<&pz2<<
+ "eP="<<edgePlus<<
+ "eM="<<edgeMinus<<
"\n";
}
}
+ if (TMath::Min(edgeMinus,edgePlus)<6){
+ AddPointRPHI(isec,c->GetX(),c->GetY(),c->GetZ(),y1,z1,py1[1],pz1[1],1);
+ TTreeSRedirector *cstream = GetDebugStreamer();
+ if (TMath::Abs(dy)<2 && cstream){
+ (*cstream)<<"rphi"<<
+ "isec="<<isec<< // current sector
+ "npads="<<npads<< // number of pads at given sector
+ "dpad="<<dpad<< // distance to edge pad
+ "eP="<<edgePlus<< // distance to edge
+ "eM="<<edgeMinus<< // distance to edge minus
+ //
+ "crR="<<corrR<< // radial correction
+ "cclY="<<corrclY<< // edge effect correction cl
+ "ctrY="<<corrtrY<< // edge effect correction using track
+ //
+ "dy="<<dy<< // dy
+ "Cl.="<<c<<
+ "y1="<<y1<<
+ "y2="<<y2<<
+ "z1="<<z1<<
+ "z2="<<z2<<
+ "py1.="<<&py1<<
+ "pz1.="<<&pz1<<
+ "\n";
+ }
+ }
+ if (c->GetType()<0) continue; // don't use edge rphi cluster
+ //
+ //
+ Double_t x[10];
+ x[0]=c->GetDetector();
+ x[1]=c->GetRow();
+ x[2]=c->GetY()/c->GetX();
+ x[3]=c->GetZ();
+ //
+ // apply sagita cut
+ //
+ if (TMath::Abs(py2[2]*150*150/4)<kSagitaCut &&
+ TMath::Abs(pz2[2]*150*150/4)<kSagitaCut){
+ //
+ x[4]=py1[1];
+ x[5]=c->GetY()-y1;
+ fDiffHistoLineY->Fill(x);
+ x[5]=c->GetY()-y2;
+ //fDiffHistoParY->Fill(x);
+ //
+ x[4]=pz1[1];
+ x[5]=c->GetY()-z1;
+ fDiffHistoLineZ->Fill(x);
+ x[5]=c->GetY()-z2;
+ //fDiffHistoParZ->Fill(x);
+ // Apply sagita cut
+ AddPoint(isec,c->GetX(),c->GetY(),c->GetZ(),y1,z1,py1[1],pz1[1],1);
+ }
+
}
}
+void AliTPCcalibUnlinearity::AlignOROC(AliTPCseed *track, Int_t isec){
+ //
+ // fit the tracklet in OROC sepatately in Quadrant
+ //
+ //
+ if (isec<36) return;
+ const Int_t kMinClusterF=35;
+ const Int_t kMinClusterQ=10;
+ //
+ const Int_t kdrow1 =8; // rows to skip at the end
+ const Int_t kdrow0 =2; // rows to skip at beginning
+ const Float_t kedgey=3;
+ const Float_t kMaxDist=0.5;
+ const Float_t kMaxCorrY=0.1;
+ const Float_t kPRFWidth = 0.6; //cut 2 sigma of PRF
+ //
+ //
+ AliTPCROC * roc = AliTPCROC::Instance();
+ AliTPCPointCorrection * corr = AliTPCPointCorrection::Instance();
+ const Double_t kXmean=roc->GetPadRowRadii(isec,53);
+ //
+ // full track fit parameters
+ //
+ TLinearFitter fyf(2,"pol1");
+ TLinearFitter fzf(2,"pol1");
+ TVectorD pyf(2), peyf(2),pzf(2), pezf(2);
+ Int_t nf=0;
+ //
+ // make full fit as reference
+ //
+ for (Int_t iter=0; iter<2; iter++){
+ fyf.ClearPoints();
+ for (Int_t irow=kdrow0;irow<159-kdrow1;irow++) {
+ AliTPCclusterMI *c=track->GetClusterPointer(irow);
+ if (!c) continue;
+ if (c->GetDetector()!=isec) continue;
+ if (c->GetRow()<kdrow0) continue;
+ Double_t dx = c->GetX()-kXmean;
+ Double_t x[2]={dx, dx*dx};
+ if (iter==0 &&c->GetType()<0) continue;
+ if (iter==1){
+ Double_t yfit = pyf[0]+pyf[1]*dx;
+ Double_t dedge = c->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(yfit);
+ if (TMath::Abs(c->GetY()-yfit)>kMaxDist) continue;
+ if (dedge<kedgey) continue;
+ Double_t corrtrY =
+ corr->RPhiCOGCorrection(isec,c->GetRow(), c->GetPad(),
+ c->GetY(),yfit, c->GetZ(), pyf[1], c->GetMax(),2.5);
+ if (TMath::Abs(corrtrY)>kMaxCorrY) continue;
+ }
+ fyf.AddPoint(x,c->GetY(),0.1);
+ fzf.AddPoint(x,c->GetZ(),0.1);
+ }
+ nf = fyf.GetNpoints();
+ if (nf<kMinClusterF) return; // not enough points - skip
+ fyf.Eval();
+ fyf.GetParameters(pyf);
+ fyf.GetErrors(peyf);
+ fzf.Eval();
+ fzf.GetParameters(pzf);
+ fzf.GetErrors(pezf);
+ }
+ //
+ // Make Fitters and params for 3 fitters
+ //
+ TLinearFitter *fitters[4];
+ Int_t npoints[4];
+ TVectorD params[4];
+ TVectorD errors[4];
+ Double_t chi2C[4];
+ for (Int_t i=0;i<4;i++) {
+ fitters[i] = new TLinearFitter(2,"pol1");
+ npoints[i]=0;
+ params[i].ResizeTo(2);
+ errors[i].ResizeTo(2);
+ }
+ //
+ // Update fitters
+ //
+ for (Int_t irow=kdrow0;irow<159-kdrow1;irow++) {
+ AliTPCclusterMI *c=track->GetClusterPointer(irow);
+ if (!c) continue;
+ if (c->GetDetector()!=isec) continue;
+ if (c->GetRow()<kdrow0) continue;
+ Double_t dx = c->GetX()-kXmean;
+ Double_t x[2]={dx, dx*dx};
+ Double_t yfit = pyf[0]+pyf[1]*dx;
+ Double_t dedge = c->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(yfit);
+ if (TMath::Abs(c->GetY()-yfit)>kMaxDist) continue;
+ if (dedge<kedgey) continue;
+ Double_t corrtrY =
+ corr->RPhiCOGCorrection(isec,c->GetRow(), c->GetPad(),
+ c->GetY(),yfit, c->GetZ(), pyf[1], c->GetMax(),2.5);
+ if (TMath::Abs(corrtrY)>kMaxCorrY) continue;
+ if (dx<0){
+ if (yfit<-kPRFWidth) fitters[0]->AddPoint(x,c->GetY(),0.1);
+ if (yfit>kPRFWidth) fitters[1]->AddPoint(x,c->GetY(),0.1);
+ }
+ if (dx>0){
+ if (yfit<-kPRFWidth) fitters[2]->AddPoint(x,c->GetY(),0.1);
+ if (yfit>kPRFWidth) fitters[3]->AddPoint(x,c->GetY(),0.1);
+ }
+ }
+
+ for (Int_t i=0;i<4;i++) {
+ npoints[i] = fitters[i]->GetNpoints();
+ if (npoints[i]>=kMinClusterQ){
+ fitters[i]->Eval();
+ Double_t chi2Fac = TMath::Sqrt(fitters[i]->GetChisquare()/(fitters[i]->GetNpoints()-2));
+ chi2C[i]=chi2Fac;
+ fitters[i]->GetParameters(params[i]);
+ fitters[i]->GetErrors(errors[i]);
+ // renormalize errors
+ errors[i][0]*=chi2Fac;
+ errors[i][1]*=chi2Fac;
+ }
+ }
+ //
+ // Fill fitters
+ //
+ for (Int_t i0=0;i0<4;i0++){
+ for (Int_t i1=0;i1<4;i1++){
+ if (i0==i1) continue;
+ if(npoints[i0]<kMinClusterQ) continue;
+ if(npoints[i1]<kMinClusterQ) continue;
+ Int_t index0=i0*4+i1;
+ Double_t xy[1];
+ Double_t xfi[1];
+ xy[0] = pyf[1];
+ xfi[0] = pyf[1];
+ //
+ Double_t dy = params[i1][0]-params[i0][0];
+ Double_t sy = TMath::Sqrt(errors[i1][0]*errors[i1][0]+errors[i0][0]*errors[i0][0]);
+ Double_t dphi = params[i1][1]-params[i0][1];
+ Double_t sphi = TMath::Sqrt(errors[i1][1]*errors[i1][1]+errors[i0][1]*errors[i0][1]);
+ //
+ Int_t sector = isec-36;
+ //
+ if (TMath::Abs(dy/(sy+0.1))>5.) continue; // 5 sigma cut
+ if (TMath::Abs(dphi/(sphi+0.004))>5.) continue; // 5 sigma cut
+ TLinearFitter * fitterY,*fitterPhi;
+ fitterY = (TLinearFitter*) fFitterQuadrantY.At(16*sector+index0);
+ if (TMath::Abs(dy/(sy+0.1))<5.) fitterY->AddPoint(xy,dy,sy);
+ fitterY = (TLinearFitter*) fFitterQuadrantY.At(16*36+index0);
+ if (TMath::Abs(dy/(sy+0.1))<5.) fitterY->AddPoint(xy,dy,sy);
+ //
+ fitterPhi = (TLinearFitter*) fFitterQuadrantPhi.At(16*sector+index0);
+ if (TMath::Abs(dphi/(sphi+0.001))<5.) fitterPhi->AddPoint(xfi,dphi,sphi);
+ fitterPhi = (TLinearFitter*) fFitterQuadrantPhi.At(16*36+index0);
+ if (TMath::Abs(dphi/(sphi+0.001))<5.) fitterPhi->AddPoint(xfi,dphi,sphi);
+ }
+ }
+ //
+ // Dump debug information
+ //
+ if (fStreamLevel>0){
+ TTreeSRedirector *cstream = GetDebugStreamer();
+ Int_t sector = isec-36;
+ if (cstream){
+ for (Int_t i0=0;i0<4;i0++){
+ for (Int_t i1=i0;i1<4;i1++){
+ if (i0==i1) continue;
+ if(npoints[i0]<kMinClusterQ) continue;
+ if(npoints[i1]<kMinClusterQ) continue;
+ (*cstream)<<"fitQ"<<
+ "run="<<fRun<< // run number
+ "event="<<fEvent<< // event number
+ "time="<<fTime<< // time stamp of event
+ "trigger="<<fTrigger<< // trigger
+ "mag="<<fMagF<< // magnetic field
+ "sec="<<sector<< // current sector from 0 to 36
+ "isec="<<isec<< // current sector
+ // full fit
+ "nf="<<nf<< // total number of points
+ "pyf.="<<&pyf<< // full OROC fit y
+ "pzf.="<<&pzf<< // full OROC fit z
+ // quadrant fit
+ "i0="<<i0<< // quadrant number
+ "i1="<<i1<<
+ "n0="<<npoints[i0]<< // number of points
+ "n1="<<npoints[i1]<<
+ "p0.="<<¶ms[i0]<< // parameters
+ "p1.="<<¶ms[i1]<<
+ "e0.="<<&errors[i0]<< // errors
+ "e1.="<<&errors[i1]<<
+ "chi0="<<chi2C[i0]<< // chi2s
+ "chi1="<<chi2C[i1]<<
+
+ "\n";
+ }
+ }
+ }
+ }
+}
+
+
+
+
+
void AliTPCcalibUnlinearity::MakeHisto(){
//
//
//
axisName[0]="sector";
xmin[0]=0; xmax[0]=72; nbins[0]=72;
- //
+ //
axisName[1]="row";
- xmin[1]=0; xmax[1]=159; nbins[1]=159;
- //
- axisName[2]="ly";
- xmin[2]=-50; xmax[2]=50; nbins[2]=10;
- //
- axisName[3]="lz";
- xmin[3]=-250; xmax[3]=250; nbins[3]=50;
- //
- axisName[4]="p2";
- xmin[4]=-0.6; xmax[4]=0.6; nbins[4]=12;
- //
- axisName[5]="p3";
- xmin[5]=-0.6; xmax[5]=0.6; nbins[5]=12;
+ xmin[1]=0; xmax[1]=96; nbins[1]=96;
//
- axisName[6]="p4";
- xmin[6]=-2; xmax[6]=2; nbins[6]=20;
- //
- axisName[7]="dy";
- xmin[7]=-0.5; xmax[7]=0.5; nbins[7]=100;
+ axisName[2]="tphi"; // tan phi - ly/lx
+ xmin[2]=-0.17; xmax[2]=0.17; nbins[2]=5;
+ //
+ axisName[3]="lz"; // global z
+ xmin[3]=-250; xmax[3]=250; nbins[3]=10;
+ //
+ axisName[4]="k"; // tangent - in angle
+ xmin[4]=-0.5; xmax[4]=0.5; nbins[4]=10;
+ //
+ //
+ axisName[5]="delta"; // delta
+ xmin[5]=-0.5; xmax[5]=0.5; nbins[5]=100;
+ //
+ //
+ fDiffHistoLineY = new THnSparseF("hnDistHistoLineY","hnDistHistoLineY",6,nbins,xmin,xmax);
+ fDiffHistoParY = new THnSparseF("hnDistHistoParY","hnDistHistoParY",6,nbins,xmin,xmax);
+ fDiffHistoLineZ = new THnSparseF("hnDistHistoLineZ","hnDistHistoLineZ",6,nbins,xmin,xmax);
+ fDiffHistoParZ = new THnSparseF("hnDistHistoParZ","hnDistHistoParZ",6,nbins,xmin,xmax);
+ // for (Int_t i=0; i<8;i++){
+// fDiffHistoLineY->GetAxis(i)->SetName(axisName[i].Data());
+// fDiffHistoParY->GetAxis(i)->SetName(axisName[i].Data());
+// fDiffHistoLineY->GetAxis(i)->SetTitle(axisName[i].Data());
+// fDiffHistoParY->GetAxis(i)->SetTitle(axisName[i].Data());
+// fDiffHistoLineZ->GetAxis(i)->SetName(axisName[i].Data());
+// fDiffHistoParZ->GetAxis(i)->SetName(axisName[i].Data());
+// fDiffHistoLineZ->GetAxis(i)->SetTitle(axisName[i].Data());
+// fDiffHistoParZ->GetAxis(i)->SetTitle(axisName[i].Data());
+// }
+
//
- axisName[8]="dz";
- xmin[8]=-0.5; xmax[8]=0.5; nbins[8]=100;
+ //
//
- fDiffHistoLine = new THnSparseF("DistHistoLine","DistHistoLine",9,nbins,xmin,xmax);
- fDiffHistoPar = new THnSparseF("DistHistoPar","DistHistoPar",9,nbins,xmin,xmax);
- for (Int_t i=0; i<9;i++){
- fDiffHistoLine->GetAxis(i)->SetName(axisName[i].Data());
- fDiffHistoPar->GetAxis(i)->SetName(axisName[i].Data());
- fDiffHistoLine->GetAxis(i)->SetTitle(axisName[i].Data());
- fDiffHistoPar->GetAxis(i)->SetTitle(axisName[i].Data());
+ char hname[1000];
+ TH2F * his=0;
+ for (Int_t isec=0;isec<74;isec++){
+ sprintf(hname,"DeltaRPhiPlus_Sector%d",isec);
+ his = new TH2F(hname,hname,100,1,4,100,-0.5,0.5);
+ his->SetDirectory(0);
+ fDistRPHIPlus.AddAt(his,isec);
+ sprintf(hname,"DeltaRPhiMinus_Sector%d",isec);
+ his = new TH2F(hname,hname,100,1,4,100,-0.5,0.5);
+ his->SetDirectory(0);
+ fDistRPHIMinus.AddAt(his,isec);
}
}
AliTPCcalibBase::Terminate();
}
+Long64_t AliTPCcalibUnlinearity::Merge(TCollection *list) {
+ //
+ // merge results
+ //
+ TIterator* iter = list->MakeIterator();
+ AliTPCcalibUnlinearity* cal = 0;
+ //
+ while ((cal = (AliTPCcalibUnlinearity*)iter->Next())) {
+ if (!cal->InheritsFrom(AliTPCcalibUnlinearity::Class())) {
+ return -1;
+ }
+ Add(cal);
+ }
+ return 0;
+}
+
+void AliTPCcalibUnlinearity::Add(AliTPCcalibUnlinearity * calib){
+ //
+ //
+ //
+ if (!fInit) Init();
+ if (fDiffHistoLineY && calib->fDiffHistoLineY){
+ fDiffHistoLineY->Add(calib->fDiffHistoLineY);
+ }
+ if (fDiffHistoParY && calib->fDiffHistoParY){
+ fDiffHistoParY->Add(calib->fDiffHistoParY);
+ }
+ if (fDiffHistoLineZ && calib->fDiffHistoLineZ){
+ fDiffHistoLineZ->Add(calib->fDiffHistoLineZ);
+ }
+ if (fDiffHistoParZ && calib->fDiffHistoParZ){
+ fDiffHistoParZ->Add(calib->fDiffHistoParZ);
+ }
+ for (Int_t isec=0;isec<38;isec++){
+ TLinearFitter *f0r = (TLinearFitter*)fFittersOutR.At(isec);
+ TLinearFitter *f1r = (TLinearFitter*)(calib->fFittersOutR.At(isec));
+ if (f0r&&f1r) f0r->Add(f1r);
+ TLinearFitter *f0z = (TLinearFitter*)fFittersOutZ.At(isec);
+ TLinearFitter *f1z = (TLinearFitter*)(calib->fFittersOutZ.At(isec));
+ if (f0z&&f1z) f0z->Add(f1z);
+ }
+
+ for (Int_t isec=0;isec<16*38;isec++){
+ TLinearFitter *f0y = (TLinearFitter*)fFitterQuadrantY.At(isec);
+ TLinearFitter *f1y = (TLinearFitter*)(calib->fFitterQuadrantY.At(isec));
+ if (f0y&&f1y) f0y->Add(f1y);
+ TLinearFitter *f0phi = (TLinearFitter*)fFitterQuadrantPhi.At(isec);
+ TLinearFitter *f1phi = (TLinearFitter*)(calib->fFitterQuadrantPhi.At(isec));
+ if (f0phi&&f1phi) f0phi->Add(f1phi);
+ }
+
+
+ for (Int_t isec=0;isec<74;isec++){
+ TH2F * h0p= (TH2F*)(fDistRPHIPlus.At(isec));
+ TH2F * h1p= (TH2F*)(calib->fDistRPHIPlus.At(isec));
+ if (h0p&&h1p) h0p->Add(h1p);
+ TH2F * h0m= (TH2F*)(fDistRPHIMinus.At(isec));
+ TH2F * h1m= (TH2F*)(calib->fDistRPHIMinus.At(isec));
+ if (h0m&&h1m) h0m->Add(h1m);
+ }
+
+
+}
+
+
-void AliTPCcalibUnlinearity::DumpTree(){
+void AliTPCcalibUnlinearity::DumpTree(const char *fname){
//
//
//
- if (fStreamLevel==0) return;
- TTreeSRedirector *cstream = GetDebugStreamer();
+ TTreeSRedirector *cstream =new TTreeSRedirector(fname);
if (!cstream) return;
//
THnSparse *his=0;
Int_t *bins = new Int_t[10];
//
for (Int_t ihis=0; ihis<2; ihis++){
- his = (ihis==0)? fDiffHistoLine:fDiffHistoPar;
+ his = (ihis==0)? fDiffHistoLineY:fDiffHistoParY;
if (!his) continue;
//
Int_t ndim = his->GetNdimensions();
"\n";
}
}
+ delete cstream;
}
-void AliTPCcalibUnlinearity::AddPoint(Int_t sector, Int_t row, Double_t dz, Double_t dy, Double_t p2, Double_t p3, Double_t dr, Int_t npoints){
+void AliTPCcalibUnlinearity::AddPoint(Int_t sector, Double_t cx, Double_t cy, Double_t cz, Double_t ty, Double_t tz, Double_t ky, Double_t kz, Int_t npoints){
//
//
// Simple distortion fit in outer sectors
// - Decay length 10 and 5 cm
// - Decay amplitude - Z dependent
//
- /*
- chainUnlin->SetAlias("side","(-1+((sector%36)<18)*2)");
- chainUnlin->SetAlias("sa","sin(pi*sector*20/180)");
- chainUnlin->SetAlias("ca","cos(pi*sector*20/180)");
- chainUnlin->SetAlias("dzt","dz*side");
- chainUnlin->SetAlias("dr","(1-abs(lz/250))");
- chainUnlin->SetAlias("dout","(159-row)*1.5");
- chainUnlin->SetAlias("din","row*0.75");
- chainUnlin->SetAlias("eout10","exp(-(dout)/10.)");
- chainUnlin->SetAlias("eout5","exp(-(dout)/5.)");
- */
-
+
Double_t side = (-1+((sector%36)<18)*2); // A side +1 C side -1
- Double_t dzt = dz*side;
- Double_t dout = (159-row)*1.5; // distance to the last pad row - (valid only for long pads)
- if (dout>30) return; // use only edge pad rows
- dr-=0.5; // dr shifted to the middle - reduce numerical instabilities
-
+ Double_t dzt = (cz-tz)*side;
+ Double_t radius = TMath::Sqrt(cx*cx+cy*cy);
+ AliTPCROC * roc = AliTPCROC::Instance();
+ Double_t xout = roc->GetPadRowRadiiUp(roc->GetNRows(37)-1);
+ Double_t dout = xout-radius;
+ if (dout>30) return;
+ //drift
+ Double_t dr = 0.5 - TMath::Abs(cz/250.);
+ Double_t dy = cy-ty;
Double_t eout10 = TMath::Exp(-dout/10.);
Double_t eout5 = TMath::Exp(-dout/5.);
//
xxxZ[4]=dr*dr*eoutp; //fstring+="dr*dr*eoutp++";
xxxZ[5]=dr*dr*eoutm; //fstring+="dr*dr*eoutm++";
//
- xxxR[0]=p2*eoutp; //fstring+="eoutp++";
- xxxR[1]=p2*eoutm; //fstring+="eoutm++";
- xxxR[2]=p2*dr*eoutp; //fstring+="dr*eoutp++";
- xxxR[3]=p2*dr*eoutm; //fstring+="dr*eoutm++";
- xxxR[4]=p2*dr*dr*eoutp; //fstring+="dr*dr*eoutp++";
- xxxR[5]=p2*dr*dr*eoutm; //fstring+="dr*dr*eoutm++";
+ xxxR[0]=ky*eoutp; //fstring+="eoutp++";
+ xxxR[1]=ky*eoutm; //fstring+="eoutm++";
+ xxxR[2]=ky*dr*eoutp; //fstring+="dr*eoutp++";
+ xxxR[3]=ky*dr*eoutm; //fstring+="dr*eoutm++";
+ xxxR[4]=ky*dr*dr*eoutp; //fstring+="dr*dr*eoutp++";
+ xxxR[5]=ky*dr*dr*eoutm; //fstring+="dr*dr*eoutm++";
//
- xxxRZ[0]=p3*eoutp; //fstring+="eoutp++";
- xxxRZ[1]=p3*eoutm; //fstring+="eoutm++";
- xxxRZ[2]=p3*dr*eoutp; //fstring+="dr*eoutp++";
- xxxRZ[3]=p3*dr*eoutm; //fstring+="dr*eoutm++";
- xxxRZ[4]=p3*dr*dr*eoutp; //fstring+="dr*dr*eoutp++";
- xxxRZ[5]=p3*dr*dr*eoutm; //fstring+="dr*dr*eoutm++";
+ xxxRZ[0]=kz*eoutp; //fstring+="eoutp++";
+ xxxRZ[1]=kz*eoutm; //fstring+="eoutm++";
+ xxxRZ[2]=kz*dr*eoutp; //fstring+="dr*eoutp++";
+ xxxRZ[3]=kz*dr*eoutm; //fstring+="dr*eoutm++";
+ xxxRZ[4]=kz*dr*dr*eoutp; //fstring+="dr*dr*eoutp++";
+ xxxRZ[5]=kz*dr*dr*eoutm; //fstring+="dr*dr*eoutm++";
//
TLinearFitter * fitter=0;
Double_t err=0.1;
}
}
}
+void AliTPCcalibUnlinearity::AddPointRPHI(Int_t sector, Double_t cx, Double_t cy, Double_t cz, Double_t ty, Double_t tz, Double_t /*ky*/, Double_t /*kz*/, Int_t npoints){
+ //
+ //
+ //
+ Float_t kMaxDz = 0.5; // cut on z in cm
+ Double_t dy = cy-ty;
+ Double_t dz = cz-tz;
+ if (TMath::Abs(dz)>kMaxDz) return;
+ Double_t edgeMinus= -ty+cx*TMath::Tan(TMath::Pi()/18.);
+ Double_t edgePlus = ty+cx*TMath::Tan(TMath::Pi()/18.);
+ //
+ TH2F * his =0;
+ his = (TH2F*)fDistRPHIPlus.At(sector);
+ his->Fill(edgePlus,dy,npoints);
+ his = (TH2F*)((sector<36)? fDistRPHIPlus.At(72):fDistRPHIPlus.At(73));
+ his->Fill(edgePlus,dy,npoints);
+ his = (TH2F*)fDistRPHIMinus.At(sector);
+ his->Fill(edgeMinus,dy,npoints);
+ his = (TH2F*)((sector<36)? fDistRPHIMinus.At(72):fDistRPHIMinus.At(73));
+ his->Fill(edgeMinus,dy,npoints);
+}
+
-void AliTPCcalibUnlinearity::MakeFitters(){
+void AliTPCcalibUnlinearity::Init(){
//
//
//
+ //
+ MakeHisto();
// Make outer fitters
TLinearFitter * fitter = 0;
for (Int_t ifit=0; ifit<38; ifit++){
fitter->StoreData(kFALSE);
fFittersOutZ.AddAt(fitter,ifit);
}
+ for (Int_t ifit=0; ifit<16*38;ifit++){
+ fitter = new TLinearFitter(2,"hyp1");
+ fitter->StoreData(kFALSE);
+ fFitterQuadrantY.AddAt(fitter,ifit);
+ fitter = new TLinearFitter(2,"hyp1");
+ fitter->StoreData(kFALSE);
+ fFitterQuadrantPhi.AddAt(fitter,ifit);
+ }
+ fInit= kTRUE;
}
void AliTPCcalibUnlinearity::EvalFitters(){
}
}
-Double_t AliTPCcalibUnlinearity::GetDr(Int_t sector, Double_t dout, Double_t dr){
- //
- //
- //
- TVectorD * vec = GetParamOutR(sector);
- if (!vec) return 0;
- dr-=0.5; // dr=0 at the middle drift length
- Double_t eout10 = TMath::Exp(-dout/10.);
- Double_t eout5 = TMath::Exp(-dout/5.);
- Double_t eoutp = (eout10+eout5)*0.5;
- Double_t eoutm = (eout10-eout5)*0.5;
-
- Double_t result=0;
- result+=(*vec)[1]*eoutp;
- result+=(*vec)[2]*eoutm;
- result+=(*vec)[3]*eoutp*dr;
- result+=(*vec)[4]*eoutm*dr;
- result+=(*vec)[5]*eoutp*dr*dr;
- result+=(*vec)[6]*eoutm*dr*dr;
- return result;
-}
-
-
-Double_t AliTPCcalibUnlinearity::GetGDr(Int_t stype,Float_t gx, Float_t gy,Float_t gz){
- //
- //
- //
- Double_t alpha = TMath::ATan2(gy,gx);
- if (alpha<0) alpha+=TMath::Pi()*2;
- Int_t lsector = Int_t(18*alpha/(TMath::Pi()*2));
- alpha = (lsector+0.5)*TMath::Pi()/9.;
- //
- Double_t r[3];
- r[0]=gx; r[1]=gy;
- Double_t cs=TMath::Cos(-alpha), sn=TMath::Sin(-alpha), x=r[0];
- r[0]=x*cs - r[1]*sn; r[1]=x*sn + r[1]*cs;
- //
- AliTPCROC * roc = AliTPCROC::Instance();
- Double_t xout = roc->GetPadRowRadiiUp(roc->GetNRows(37)-1);
- Double_t dout = xout-r[0];
- if (dout<0) return 0;
- Double_t dr = 1-TMath::Abs(gz/250.);
- if (gz<0) lsector+=18;
- if (stype==0) lsector = (gz>0) ? 36:37;
- if (stype<0) return lsector; //
- return GetDr(lsector,dout,dr);
-}
-
-
-Double_t AliTPCcalibUnlinearity::GetDz(Int_t sector, Double_t dout, Double_t dr){
- //
- //
- //
- TVectorD * vec = GetParamOutZ(sector);
- if (!vec) return 0;
- dr-=0.5; // 0 at the middle
- Double_t eout10 = TMath::Exp(-dout/10.);
- Double_t eout5 = TMath::Exp(-dout/5.);
- Double_t eoutp = (eout10+eout5)*0.5;
- Double_t eoutm = (eout10-eout5)*0.5;
-
-
- Double_t result=0;
- result+=(*vec)[1]*eoutp;
- result+=(*vec)[2]*eoutm;
- result+=(*vec)[3]*eoutp*dr;
- result+=(*vec)[4]*eoutm*dr;
- result+=(*vec)[5]*eoutp*dr*dr;
- result+=(*vec)[6]*eoutm*dr*dr;
- return result;
-}
-
-
-Double_t AliTPCcalibUnlinearity::SGetDr(Int_t sector, Double_t dout, Double_t dr){
- //
- //
- //
- return fgInstance->GetDr(sector,dout,dr);
-}
-Double_t AliTPCcalibUnlinearity::SGetGDr(Int_t stype,Float_t gx, Float_t gy,Float_t gz){
- //
- //
- //
- return fgInstance->GetGDr(stype,gx,gy,gz);
-}
-Double_t AliTPCcalibUnlinearity::SGetDz(Int_t sector, Double_t dout, Double_t dr){
- //
- //
- //
- return fgInstance->GetDz(sector,dout,dr);
-}
if (i%10000==0) printf("%d\n",(Int_t)i);
Int_t row= cl->GetRow();
if (cl->GetDetector()>36) row+=64;
+ calib->AddPointRPHI(cl->GetDetector(),cl->GetX(),cl->GetY(),cl->GetZ(),
+ ty,tz,(*vy)[1],(*vz)[1],1);
+
if (cl->GetType()<0) continue;
Double_t dy = cl->GetY()-ty;
Double_t dz = cl->GetZ()-tz;
- Double_t dr = 1-TMath::Abs(cl->GetZ()/250);
//
//
if (TMath::Abs(dy)>0.25) continue;
if (TMath::Abs((*vy2)[2]*150*150/4)>1) continue;
if (TMath::Abs((*vz2)[2]*150*150/4)>1) continue;
// Apply sagita cut
- calib->AddPoint(cl->GetDetector(), row, dz, dy,
- (*vy)[1],(*vz)[1], dr, 1);
+ if (cl->GetType()>=0) {
+ calib->AddPoint(cl->GetDetector(),cl->GetX(),cl->GetY(),cl->GetZ(),
+ ty,tz,(*vy)[1],(*vz)[1],1);
+ }
}
}
}
class AliESDEvent;
class AliESDtrack;
class TLinearFitter;
+class AliTPCClusterParam;
class AliTPCcalibUnlinearity:public AliTPCcalibBase {
virtual void Process(AliTPCseed *track);
virtual void Analyze(){return;}
virtual void Terminate();
- virtual Long64_t Merge(TCollection* list){return 0;}
+ virtual Long64_t Merge(TCollection* list);
+ void Add(AliTPCcalibUnlinearity * calib);
//
void ProcessTree(TTree * tree, Long64_t nmaxPoints);
- void AddPoint(Int_t sector, Int_t row, Double_t dz, Double_t dy, Double_t p2, Double_t p3, Double_t dr, Int_t npoints=1);
+ void AddPoint(Int_t sector, Double_t cx, Double_t cy, Double_t cz, Double_t ty, Double_t tz, Double_t ky, Double_t kz, Int_t npoints=1);
+ void AddPointRPHI(Int_t sector, Double_t cx, Double_t cy, Double_t cz, Double_t ty, Double_t tz, Double_t ky, Double_t kz, Int_t npoints=1);
//
void MakeHisto();
void ProcessDiff(AliTPCseed *track, Int_t isec);
- void DumpTree();
- void MakeFitters();
+ void AlignOROC(AliTPCseed *track, Int_t isec);
+ //
+ void DumpTree(const char *fname="unlinResidual.root");
+ void Init();
void EvalFitters();
TLinearFitter * GetFitterOutR(Int_t sector) {return (TLinearFitter*)fFittersOutR.At(sector);}
TLinearFitter * GetFitterOutZ(Int_t sector) {return (TLinearFitter*)fFittersOutZ.At(sector);}
TVectorD * GetParamOutR(Int_t sector) {return (TVectorD*)fParamsOutR.At(sector);}
TVectorD * GetParamOutZ(Int_t sector) {return (TVectorD*)fParamsOutZ.At(sector);}
//
- Double_t GetDr(Int_t sector, Double_t dout, Double_t dr);
- Double_t GetDz(Int_t sector, Double_t dout, Double_t dr);
- Double_t GetGDr(Int_t stype,Float_t gx, Float_t gy,Float_t gz);
- //
- static Double_t SGetDr(Int_t sector, Double_t dout, Double_t dr);
- static Double_t SGetDz(Int_t sector, Double_t dout, Double_t dr);
- static Double_t SGetGDr(Int_t stype,Float_t gx, Float_t gy,Float_t gz);
-
- static AliTPCcalibUnlinearity* Instance();
- void SetInstance(AliTPCcalibUnlinearity*param){fgInstance = param;}
//TMatrixD * GetNormCovariance(Int_t sector, Int_t type);
//TMatrixD * GetNormCovariance(Int_t sector, Int_t type);
static void MakeQPosNormAll(TTree * chainres, AliTPCClusterParam * param, Int_t maxPoints);
+ void Process(AliESDEvent *event) {AliTPCcalibBase::Process(event);};
+ void Process(AliESDtrack *track, Int_t runNo=-1){AliTPCcalibBase::Process(track,runNo);};
public:
- THnSparse * fDiffHistoLine; // matrix with cluster residuals - linear fit
- THnSparse * fDiffHistoPar; // matrix with cluster residuals - parabolic fit
+ Bool_t fInit; // initialization flag
+ THnSparse * fDiffHistoLineY; // matrix with cluster residuals - linear fit
+ THnSparse * fDiffHistoParY; // matrix with cluster residuals - parabolic fit
+ THnSparse * fDiffHistoLineZ; // matrix with cluster residuals - linear fit
+ THnSparse * fDiffHistoParZ; // matrix with cluster residuals - parabolic fit
//
+ // Outer residula fit
//
TObjArray fFittersOutR; // Unlinearity fitters for radial distortion - outer field cage
TObjArray fFittersOutZ; // Unlinearity fitters for z distortion - outer field cage
TObjArray fErrorsOutR; // Parameters for radial distortion - outer field cage
TObjArray fErrorsOutZ; // Parameters for z distortion - outer field cage
//
-
+ // R-phi residual histogram
+ //
+ TObjArray fDistRPHIPlus; // edge effect histograms - plus direction
+ TObjArray fDistRPHIMinus; // edge effect histograms - minus direction
+ //
+ // Quadrant fitters
+ //
+ TObjArray fFitterQuadrantY; //qudrant misalignemnt fit Y
+ TObjArray fFitterQuadrantPhi; //qudrant misalignemnt fit Phi
private:
AliTPCcalibUnlinearity(const AliTPCcalibUnlinearity&);
AliTPCcalibUnlinearity& operator=(const AliTPCcalibUnlinearity&);
- static AliTPCcalibUnlinearity* fgInstance; //! Instance of this class (singleton implementation)
- ClassDef(AliTPCcalibUnlinearity, 1);
+ ClassDef(AliTPCcalibUnlinearity, 4);
};
#endif
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff