]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TRD/qaRec/AliTRDresolution.cxx
fix warnings (thanks to Federico)
[u/mrichter/AliRoot.git] / TRD / qaRec / AliTRDresolution.cxx
index 02bd9fa956794b105082a0781cae7bef7a611047..a022e4f7471583799955052cd94da3eeca168004 100644 (file)
@@ -302,9 +302,9 @@ TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track)
   UChar_t s;
   Int_t pdg = fMC->GetPDG(), det=-1;
   Int_t label = fMC->GetLabel();
-  Double_t x, y;
-  Float_t p, pt, x0, y0, z0, dx, dy, dz, dydx, dzdx;
-  Double_t covR[3], cov[3];
+  Double_t x, y, z, pt, dydx, dzdx;
+  Float_t p, pt0, x0, y0, z0, dx, dy, dz, dydx0, dzdx0;
+  Double_t covR[3]/*, cov[3]*/;
 
   if(fDebugLevel>=1){
     Double_t DX[12], DY[12], DZ[12], DPt[12], COV[12][15];
@@ -335,56 +335,55 @@ TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track)
     x0  = fTracklet->GetX0();
     //radial shift with respect to the MC reference (radial position of the pad plane)
     x= fTracklet->GetX();
-    if(!fMC->GetDirections(x0, y0, z0, dydx, dzdx, pt, s)) continue;
+    if(!fMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, s)) continue;
     // MC track position at reference radial position
     dx  = x0 - x;
     if(fDebugLevel>=1){
       (*fDebugStream) << "MC"
         << "det="     << det
         << "pdg="     << pdg
-        << "pt="      << pt
-        << "dx="      << dx
-        << "x0="      << x0
-        << "y0="      << y0
-        << "z0="      << z0
-        << "dydx="    << dydx
-        << "dzdx="    << dzdx
+        << "pt="      << pt0
+        << "x="       << x0
+        << "y="       << y0
+        << "z="       << z0
+        << "dydx="    << dydx0
+        << "dzdx="    << dzdx0
         << "\n";
     }
-    Float_t yt = y0 - dx*dydx;
-    Float_t zt = z0 - dx*dzdx;
-    p = pt*(1.+dzdx*dzdx); // pt -> p
+    Float_t yt = y0 - dx*dydx0;
+    Float_t zt = z0 - dx*dzdx0;
+    p = pt0*(1.+dzdx0*dzdx0); // pt -> p
 
     // add Kalman residuals for y, z and pt
     dx = fTracklet->GetX0() - x;
-    Float_t yr = fTracklet->GetYref(0) - dx*fTracklet->GetYref(1);
-    dy = yt - yr;
-    Float_t zr = fTracklet->GetZref(0) - dx*fTracklet->GetZref(1);
-    dz = zt - zr;
-    Float_t tgl = fTracklet->GetTgl();
-    Float_t dpt = pt - fTracklet->GetMomentum()/(1.+tgl*tgl);
+     = fTracklet->GetYref(0) - dx*fTracklet->GetYref(1);
+    dy = yt - y;
+     = fTracklet->GetZref(0) - dx*fTracklet->GetZref(1);
+    dz = zt - z;
+    dzdx = fTracklet->GetTgl();
+    pt = fTracklet->GetMomentum()/(1.+dzdx*dzdx);
     fTracklet->GetCovRef(covR);
 
-    ((TH2I*)fContainer->At(kMCtrackY))->Fill(dydx, dy);
-    ((TH2I*)fContainer->At(kMCtrackYPull))->Fill(dydx, dy/TMath::Sqrt(covR[0]));
+    ((TH2I*)fContainer->At(kMCtrackY))->Fill(dydx0, dy);
+    ((TH2I*)fContainer->At(kMCtrackYPull))->Fill(dydx0, dy/TMath::Sqrt(covR[0]));
     if(ily==0){
-      ((TH2I*)fContainer->At(kMCtrackZIn))->Fill(dzdx, dz);
-      ((TH2I*)fContainer->At(kMCtrackZInPull))->Fill(dzdx, dz/TMath::Sqrt(covR[2]));
+      ((TH2I*)fContainer->At(kMCtrackZIn))->Fill(dzdx0, dz);
+      ((TH2I*)fContainer->At(kMCtrackZInPull))->Fill(dzdx0, dz/TMath::Sqrt(covR[2]));
     } else if(ily==AliTRDgeometry::kNlayer-1) {
-      ((TH2I*)fContainer->At(kMCtrackZOut))->Fill(dzdx, dz);
-      ((TH2I*)fContainer->At(kMCtrackZOutPull))->Fill(dzdx, dz/TMath::Sqrt(covR[2]));
+      ((TH2I*)fContainer->At(kMCtrackZOut))->Fill(dzdx0, dz);
+      ((TH2I*)fContainer->At(kMCtrackZOutPull))->Fill(dzdx0, dz/TMath::Sqrt(covR[2]));
     }
-    if(pdg!=kElectron && pdg!=kPositron) ((TH2I*)fContainer->At(kMCtrackPt))->Fill(1./pt, dpt);
+    if(pdg!=kElectron && pdg!=kPositron) ((TH2I*)fContainer->At(kMCtrackPt))->Fill(1./pt0, pt-pt0);
     // Fill Debug stream for Kalman track
     if(fDebugLevel>=1){
-      Float_t dydxr = fTracklet->GetYref(1);
+      dydx = fTracklet->GetYref(1);
       (*fDebugStream) << "MCtrack"
+        << "pt="      << pt
         << "x="       << x
-        << "dpt="     << dpt
-        << "dy="      << dy
-        << "dz="      << dz
-        << "dydxr="   << dydxr
-        << "dzdxr="   << tgl
+        << "y="       << y
+        << "z="       << z
+        << "dydx="    << dydx
+        << "dzdx="    << dzdx
         << "s2y="     << covR[0]
         << "s2z="     << covR[2]
         << "\n";
@@ -392,41 +391,45 @@ TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track)
 
     // recalculate tracklet based on the MC info
     AliTRDseedV1 tt(*fTracklet);
-    tt.SetZref(0, z0);
-    tt.SetZref(1, dzdx); 
+    tt.SetZref(0, z0 - (x0-tt.GetX0())*dzdx0);
+    tt.SetZref(1, dzdx0); 
     tt.Fit(kTRUE);
-    x= tt.GetX(); // the true one 
+    x= tt.GetX();y= tt.GetY();z= tt.GetZ();
+    dydx = tt.GetYfit(1);
     dx = x0 - x;
-    yt = y0 - dx*dydx;
+    yt = y0 - dx*dydx0;
+    zt = z0 - dx*dzdx0;
     Bool_t rc = tt.IsRowCross(); 
     
     // add tracklet residuals for y and dydx
-    Float_t yf = tt.GetY();
-    dy = yt - yf; dz = 100.;
-    Float_t dphi   = (tt.GetYfit(1) - dydx);
-    dphi /= 1.- tt.GetYfit(1)*dydx;
-    tt.GetCovAt(x, cov);
     if(!rc){
-      ((TH2I*)fContainer->At(kMCtrackletY))->Fill(dydx, dy);
-      if(cov[0]>0.) ((TH2I*)fContainer->At(kMCtrackletYPull))->Fill(dydx, dy/TMath::Sqrt(cov[0]));
-      ((TH2I*)fContainer->At(kMCtrackletPhi))->Fill(dydx, dphi*TMath::RadToDeg());
+      dy    = yt-y;
+
+      Float_t dphi  = (dydx - dydx0);
+      dphi /= 1.- dydx*dydx0;
+
+      ((TH2I*)fContainer->At(kMCtrackletY))->Fill(dydx0, dy);
+      if(tt.GetS2Y()>0.) ((TH2I*)fContainer->At(kMCtrackletYPull))->Fill(dydx0, dy/TMath::Sqrt(tt.GetS2Y()));
+      ((TH2I*)fContainer->At(kMCtrackletPhi))->Fill(dydx0, dphi*TMath::RadToDeg());
     } else {
       // add tracklet residuals for z
-      dz = tt.GetZ() - (z0 - dx*dzdx) ;
-      ((TH2I*)fContainer->At(kMCtrackletZ))->Fill(dzdx, dz);
-      if(cov[2]>0.) ((TH2I*)fContainer->At(kMCtrackletZPull))->Fill(dzdx, dz/TMath::Sqrt(cov[2]));
+      dz = zt-z;
+      ((TH2I*)fContainer->At(kMCtrackletZ))->Fill(dzdx0, dz);
+      if(tt.GetS2Z()>0.) ((TH2I*)fContainer->At(kMCtrackletZPull))->Fill(dzdx0, dz/TMath::Sqrt(tt.GetS2Z()));
     }
   
     // Fill Debug stream for tracklet
     if(fDebugLevel>=1){
+      Float_t s2y = tt.GetS2Y();
+      Float_t s2z = tt.GetS2Z();
       (*fDebugStream) << "MCtracklet"
         << "rc="    << rc
         << "x="     << x
-        << "dy="    << dy
-        << "dz="    << dz
-        << "dphi="  << dphi
-        << "s2y="   << cov[0]
-        << "s2z="   << cov[2]
+        << "y="     << y
+        << "z="     << z
+        << "dydx="  << dydx
+        << "s2y="   << s2y
+        << "s2z="   << s2z
         << "\n";
     }
 
@@ -441,16 +444,14 @@ TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track)
       Float_t  q = TMath::Abs(c->GetQ());
       AliTRDseedV1::GetClusterXY(c,x,y);
       //x = c->GetX(); y = c->GetY();
-      Float_t xc = x;
-      Float_t yc = y;
-      Float_t zc = c->GetZ();
-      dx = x0 - xc; 
-      yt = y0 - dx*dydx;
-      zt = z0 - dx*dzdx;
-      dy = yt - (yc - tilt*(zc-zt));
+      z = c->GetZ();
+      dx = x0 - x; 
+      yt = y0 - dx*dydx0;
+      zt = z0 - dx*dzdx0;
+      dy = yt - (y - tilt*(z-zt));
 
       // Fill Histograms
-      if(q>20. && q<250.) ((TH2I*)fContainer->At(kMCcluster))->Fill(dydx, dy);
+      if(q>20. && q<250.) ((TH2I*)fContainer->At(kMCcluster))->Fill(dydx0, dy);
 
       // Fill calibration container
       Float_t d = zr0 - zt;
@@ -460,7 +461,7 @@ TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track)
       fMCcl->Add(clInfo);
       clInfo->SetCluster(c);
       clInfo->SetMC(pdg, label);
-      clInfo->SetGlobalPosition(yt, zt, dydx, dzdx);
+      clInfo->SetGlobalPosition(yt, zt, dydx0, dzdx0);
       clInfo->SetResolution(dy);
       clInfo->SetAnisochronity(d);
       clInfo->SetDriftLength(dx-.5*AliTRDgeometry::CamHght());
@@ -533,9 +534,9 @@ Bool_t AliTRDresolution::GetRefFigure(Int_t ifig)
   case kMCcluster:
     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
     ax = g->GetHistogram()->GetYaxis();
-    y[0] = -.05; y[1] = 0.6;
+    y[0] = -50.; y[1] = 600.;
     ax->SetRangeUser(y[0], y[1]);
-    ax->SetTitle("Y_{cluster} #sigma/#mu [mm]");
+    ax->SetTitle("Y_{cluster} #sigma/#mu [#mum]");
     ax = g->GetHistogram()->GetXaxis();
     ax->SetRangeUser(-.3, .3);
     ax->SetTitle("tg(#phi)");
@@ -549,9 +550,9 @@ Bool_t AliTRDresolution::GetRefFigure(Int_t ifig)
   case kMCtrackletY:
     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
     ax = g->GetHistogram()->GetYaxis();
-    y[0] = -.05; y[1] = 0.25;
+    y[0] = -50.; y[1] = 250.;
     ax->SetRangeUser(y[0], y[1]);
-    ax->SetTitle("Y_{tracklet} #sigma/#mu [mm]");
+    ax->SetTitle("Y_{tracklet} #sigma/#mu [#mum]");
     ax = g->GetHistogram()->GetXaxis();
     ax->SetRangeUser(-.2, .2);
     ax->SetTitle("tg(#phi)");
@@ -565,8 +566,8 @@ Bool_t AliTRDresolution::GetRefFigure(Int_t ifig)
   case kMCtrackletZ:
     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
     ax = g->GetHistogram()->GetYaxis();
-    ax->SetRangeUser(-.5, 1.);
-    ax->SetTitle("Z_{tracklet} #sigma/#mu [mm]");
+    ax->SetRangeUser(-50., 700.);
+    ax->SetTitle("Z_{tracklet} #sigma/#mu [#mum]");
     ax = g->GetHistogram()->GetXaxis();
     ax->SetTitle("tg(#theta)");
     g->Draw("apl");
@@ -588,9 +589,9 @@ Bool_t AliTRDresolution::GetRefFigure(Int_t ifig)
   case kMCtrackY:
     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
     ax = g->GetHistogram()->GetYaxis();
-    y[0] = -.05; y[1] = 0.25;
+    y[0] = -50.; y[1] = 200.;
     ax->SetRangeUser(y[0], y[1]);
-    ax->SetTitle("Y_{track} #sigma/#mu [mm]");
+    ax->SetTitle("Y_{track} #sigma/#mu [#mum]");
     ax = g->GetHistogram()->GetXaxis();
     ax->SetRangeUser(-.2, .2);
     ax->SetTitle("tg(#phi)");
@@ -605,8 +606,8 @@ Bool_t AliTRDresolution::GetRefFigure(Int_t ifig)
   case kMCtrackZOut:
     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
     ax = g->GetHistogram()->GetYaxis();
-    ax->SetRangeUser(-.5, 2.);
-    ax->SetTitle(Form("Z_{track}^{%s} #sigma/#mu [mm]", ifig==kMCtrackZIn ? "in" : "out"));
+    ax->SetRangeUser(-500., 2000.);
+    ax->SetTitle(Form("Z_{track}^{%s} #sigma/#mu [#mum]", ifig==kMCtrackZIn ? "in" : "out"));
     ax = g->GetHistogram()->GetXaxis();
     ax->SetTitle("tg(#theta)");
     g->Draw("apl");
@@ -704,11 +705,11 @@ Bool_t AliTRDresolution::PostProcess()
   //PROCESS MC RESIDUAL DISTRIBUTIONS
 
   // cluster y resolution
-  Process(kMCcluster, &f);
+  Process(kMCcluster, &f, 1.e4);
 
   // tracklet resolution
-  Process(kMCtrackletY, &f); // y
-  Process(kMCtrackletZ, &f); // z
+  Process(kMCtrackletY, &f, 1.e4); // y
+  Process(kMCtrackletZ, &f, 1.e4); // z
   Process(kMCtrackletPhi, &f); // phi
 
   // tracklet pulls
@@ -716,9 +717,9 @@ Bool_t AliTRDresolution::PostProcess()
   Process(kMCtrackletZPull, &f); // z
 
   // track resolution
-  Process(kMCtrackY, &f);    // y
-  Process(kMCtrackZIn, &f);  // z towards TPC
-  Process(kMCtrackZOut, &f); // z towards TOF
+  Process(kMCtrackY, &f, 1.e4);    // y
+  Process(kMCtrackZIn, &f, 1.e4);  // z towards TPC
+  Process(kMCtrackZOut, &f, 1.e4); // z towards TOF
 
   // track pulls
   Process(kMCtrackYPull, &f);    // y
@@ -749,11 +750,11 @@ Bool_t AliTRDresolution::PostProcess()
     h->Fit(&fgl, "NQ", "", -5., 20.);
 
     Float_t invpt = ax->GetBinCenter(ip);
-    Int_t ip = gm->GetN();
-    gm->SetPoint(ip, invpt, fgl.GetParameter(1));
-    gm->SetPointError(ip, 0., fgl.GetParError(1));
-    gs->SetPoint(ip, invpt, fgl.GetParameter(2)*invpt);
-    gs->SetPointError(ip, 0., fgl.GetParError(2));
+    Int_t jp = gm->GetN();
+    gm->SetPoint(jp, invpt, fgl.GetParameter(1));
+    gm->SetPointError(jp, 0., fgl.GetParError(1));
+    gs->SetPoint(jp, invpt, fgl.GetParameter(2)*invpt);
+    gs->SetPointError(jp, 0., fgl.GetParError(2));
     // fgl.GetParameter(4) // Landau MPV
     // fgl.GetParameter(5) // Landau Sigma
   }
@@ -966,7 +967,7 @@ TObjArray* AliTRDresolution::Histos()
 
 
 //________________________________________________________
-Bool_t AliTRDresolution::Process(ETRDresolutionPlot ip, TF1 *f)
+Bool_t AliTRDresolution::Process(ETRDresolutionPlot plot, TF1 *f, Float_t k)
 {
   if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
   Bool_t kBUILD = kFALSE;
@@ -976,27 +977,26 @@ Bool_t AliTRDresolution::Process(ETRDresolutionPlot ip, TF1 *f)
   }
 
   TH2I *h2 = 0x0;
-  if(!(h2 = (TH2I *)(fContainer->At(ip)))) return kFALSE;
+  if(!(h2 = (TH2I *)(fContainer->At(plot)))) return kFALSE;
   TGraphErrors *gm = 0x0, *gs = 0x0;
-  if(!(gm=(TGraphErrors*)fGraphM->At(ip))) return kFALSE;
+  if(!(gm=(TGraphErrors*)fGraphM->At(plot))) return kFALSE;
   if(gm->GetN()) for(Int_t ip=gm->GetN(); ip--;) gm->RemovePoint(ip);
-  if(!(gs=(TGraphErrors*)fGraphS->At(ip))) return kFALSE;
+  if(!(gs=(TGraphErrors*)fGraphS->At(plot))) return kFALSE;
   if(gs->GetN()) for(Int_t ip=gs->GetN(); ip--;) gs->RemovePoint(ip);
+  Char_t pn[10]; sprintf(pn, "p%02d", plot);
   for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
     Double_t x = h2->GetXaxis()->GetBinCenter(ibin);
-    TH1D *h = h2->ProjectionY("py", ibin, ibin);
+    TH1D *h = h2->ProjectionY(pn, ibin, ibin);
     if(h->GetEntries()<100) continue;
     AdjustF1(h, f);
 
-//    if(IsVisual()){c->cd(); c->SetLogy();}
-    h->Fit(f, "Q", "", -0.5, 0.5);
-/*    if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}*/
+    h->Fit(f, "QN");
     
     Int_t ip = gm->GetN();
-    gm->SetPoint(ip, x, 10.*f->GetParameter(1));
-    gm->SetPointError(ip, 0., 10.*f->GetParError(1));
-    gs->SetPoint(ip, x, 10.*f->GetParameter(2));
-    gs->SetPointError(ip, 0., 10.*f->GetParError(2));
+    gm->SetPoint(ip, x, k*f->GetParameter(1));
+    gm->SetPointError(ip, 0., k*f->GetParError(1));
+    gs->SetPoint(ip, x, k*f->GetParameter(2));
+    gs->SetPointError(ip, 0., k*f->GetParError(2));
   }
 
   if(kBUILD) delete f;