PWGLF/FORWARD/analysis2/corrs/ExtractAcceptance.C

   1 //_____________________________________________________________________
   2 /**
   3  *
   4  *
   5  * @param d      Detector
   6  * @param r      Ring
   7  * @param vz     Z--coordinate of interaction point
   8  * @param nDead  On returm the number of dead strips
   9  * @param deadScript  Output stream for dead strips
  10  *
  11  * @return
  12  *
  13  * @ingroup pwglf_forward_scripts_corr
  14  */
  15 TH2D* MakeOneRing(UShort_t      d,
  16                   Char_t        r,
  17                   Double_t      vz,
  18                   Int_t&        nDead,
  19                   std::ostream* deadScript)
  20 {
  21   AliFMDGeometry*   geom = AliFMDGeometry::Instance();
  22   AliFMDParameters* pars = AliFMDParameters::Instance();
  23
  24   UShort_t nS = (r == 'I' || r == 'i' ?  20 :  40);
  25   UShort_t nT = (r == 'I' || r == 'i' ? 512 : 256);
  26
  27   // Make our two histograms
  28   TH2D* hAll = new TH2D("all","All",200,-4,6,nS,0,2*TMath::Pi());
  29   hAll->SetXTitle("#eta");
  30   hAll->SetYTitle("#phi");
  31   hAll->Sumw2();
  32   hAll->SetDirectory(0);
  33   TH2D* hOK  = static_cast<TH2D*>(hAll->Clone());
  34   hOK->SetDirectory(0);
  35
  36   if (deadScript)
  37     *deadScript << "\n // FMD" << d << r << std::endl;
  38   // Loop over all sectors and strips in this ring
  39   Int_t nOK  = 0;
  40   Int_t nAll = 0;
  41   Int_t nPhi = hAll->GetNbinsY();
  42   for (UShort_t s = 0; s < nS; s++) {
  43     for (UShort_t t = 0; t < nT; t++) {
  44       // Get eta,phi by quering the geometry (first for (x,y,z), then
  45       // correcting for the vertex position, and then calculating
  46       // (eta, phi))
  47       Double_t x, y, z;
  48       geom->Detector2XYZ(d, r, s, t, x, y, z);
  49       z -= vz;
  50       Double_t q, eta, phi, theta;
  51       AliFMDGeometry::XYZ2REtaPhiTheta(x, y, z, q, eta, phi, theta);
  52       if (phi < 0) phi += 2*TMath::Pi();
  53
  54       // Check if this is a dead channel or not
  55       Bool_t isDead = pars->IsDead(d, r, s, t);
  56
  57       // Special check for FMD2i - upper part of sectors 16/17 have
  58       // have anomalous gains and/or noise - common sources are
  59       // power regulartors for bias currents and the like
  60       Int_t VA = t/128;
  61       if(d==2 && r=='I' && VA>1 && (s==16 || s==17)) isDead =true;
  62
  63       // Find the eta bin number and corresponding overflow bin
  64       Int_t etaBin = hAll->GetXaxis()->FindBin(eta);
  65       Int_t ovrBin = hAll->GetBin(etaBin, nPhi+1);
  66
  67       // Increment all histogram
  68       hAll->Fill(eta, phi);
  69       hAll->AddBinContent(ovrBin);
  70       nAll++;
  71
  72       // If not dead, increment OK histogram
  73       if (!isDead) {
  74         hOK->Fill(eta, phi);
  75         hOK->AddBinContent(ovrBin);
  76         nOK++;
  77       }
  78       else {
  79         nDead++;
  80         if (deadScript)
  81           *deadScript << "  filter->AddDead(" << d << ",'" << r << "',"
  82                       << s << ',' << t << ");" << std::endl;
  83       }
  84     }
  85   }
  86   // Divide out the efficiency.
  87   // Note, that the overflow bins along eta now contains the ratio
  88   // nOK/nAll Strips for a given eta bin.
  89   hOK->Divide(hOK,hAll,1,1,"B");
  90
  91   // Invert overflow bin
  92   for (Int_t etaBin = 1; etaBin <= hOK->GetNbinsX(); etaBin++) {
  93     Double_t ovr    = hOK->GetBinContent(etaBin, nPhi+1);
  94     Double_t novr   = (ovr < 1e-12 ? 0 : 1./ovr);
  95     hOK->SetBinContent(etaBin, nPhi+1, novr);
  96 #if 0
  97     if (ovr > 0 && ovr != 1)
  98       Info("", "Setting overflow bin (%3d,%3d) to 1/%f=%f", etaBin, nPhi+1,
  99            ovr, hOK->GetBinContent(etaBin, nPhi+1));
 100 #endif
 101   }
 102   // Clean up
 103   delete hAll;
 104
 105   Printf("=== FMD%d%c at vz=%+5.1f - %d/%d=%3d%% OK (w/overflow)",
 106          d, r, vz, nOK, nAll, (100*nOK)/nAll);
 107
 108   // Return result
 109   return hOK;
 110 }
 111
 112 //_____________________________________________________________________
 113 /**
 114  *
 115  *
 116  * @param runNo     Run number
 117  * @param nVtxBins  Number of @f$IP_{z}@f$ bins
 118  * @param vtxLow    Least @f$IP_{z}@f$
 119  * @param vtxHigh   Largest @f$IP_{z}@f$
 120  *
 121  * @ingroup pwglf_forward_scripts_corr
 122  */
 123 void ExtractAcceptance(Int_t   runNo=121526,
 124                        Int_t   nVtxBins=10,
 125                        Float_t vtxLow=-10,
 126                        Float_t vtxHigh=10)
 127 {
 128   const char* fwd = "$ALICE_ROOT/PWGLF/FORWARD/analysis2";
 129   gSystem->AddIncludePath(Form("-I%s", fwd));
 130   gROOT->Macro(Form("%s/scripts/LoadLibs.C", fwd));
 131
 132   // gSystem->Load("libANALYSIS");
 133   // gSystem->Load("libANALYSISalice");
 134   // gSystem->Load("libPWGLFforward2");
 135   // Float_t delta = (vtxHigh - vtxLow) / (Float_t)nVtxBins;
 136
 137   Bool_t gridOnline = kTRUE;
 138   if(!(TGrid::Connect("alien://",0,0,"t")))
 139     gridOnline = kFALSE;
 140
 141   // --- Figure out the year --------------------------------------
 142   UShort_t year = 0;
 143   if      (runNo <= 99999)  year = 2009;
 144   else if (runNo <= 139667) year = 2010;
 145   else if (runNo <= 170718) year = 2011;
 146   else if (runNo <= 194306) year = 2012;
 147   else if (runNo <= 197709) year = 2013;
 148   if (year <= 0) {
 149     Error("", "Couldn't deduce the year from the run number");
 150     return;
 151   }
 152
 153   // --- Initialisations ------------------------------------------
 154   //Set up CDB manager
 155   Printf("=== Setting up OCDB");
 156
 157   AliCDBManager* cdb = AliCDBManager::Instance();
 158   if(gridOnline) {
 159     cdb->SetDefaultStorageFromRun(runNo);
 160     // cdb->SetDefaultStorage(Form("alien://Folder=/alice/data/%4d/OCDB", year));
 161   }
 162   else
 163     cdb->SetDefaultStorage("local://$(ALICE_ROOT)/OCDB");
 164   cdb->SetRun(runNo);
 165
 166   // Get the geometry
 167   Printf("=== Loading geometry");
 168   AliGeomManager::LoadGeometry();
 169
 170   // Get an initialize parameters
 171   Printf("=== Intialising parameters");
 172   AliFMDParameters* pars = AliFMDParameters::Instance();
 173   pars->Init();
 174
 175   // Get an initialise geometry
 176   Printf("=== Initialising geomtry");
 177   AliFMDGeometry* geom = AliFMDGeometry::Instance();
 178   geom->Init();
 179   geom->InitTransformations();
 180
 181   // --- Get the general run parameters ------------------------------
 182   AliCDBEntry* grpE = cdb->Get("GRP/GRP/Data");
 183   if (!grpE) {
 184     AliWarningF("No GRP entry found for run %d", runNo);
 185     return;
 186   }
 187   AliGRPObject* grp = static_cast<AliGRPObject*>(grpE->GetObject());
 188   if (!grp) {
 189     AliWarningF("No GRP object found for run %d", runNo);
 190     return;
 191   }
 192   Float_t  beamE = grp->GetBeamEnergy();
 193   TString  beamT = grp->GetBeamType();
 194 # if 0
 195   // This isn't really needed as the acceptance map is indifferent to
 196   // the field settings.
 197   Float_t  l3cur = grp->GetL3Current(AliGRPObject::kMean);
 198   Char_t   l3pol = grp->GetL3Polarity();
 199   Bool_t   l3lhc = grp->IsPolarityConventionLHC();
 200   Bool_t   l3uni = grp->IsUniformBMap();
 201   AliMagF* fldM  =
 202     AliMagF::CreateFieldMap(TMath::Abs(l3cur) * (l3pol ? -1:1), 0,
 203                             (l3lhc ? 0 : 1), l3uni, beamE, beamT.Data());
 204   Float_t  l3fld = fldM->SolenoidField();
 205 #endif
 206
 207   UShort_t sys = AliForwardUtil::ParseCollisionSystem(beamT);
 208   UShort_t sNN = AliForwardUtil::ParseCenterOfMassEnergy(sys, 2 * beamE);
 209   Short_t  fld = 0; // AliForwardUtil::ParseMagneticField(l3fld);
 210   Printf("=== Run=%d, year=%d, sys=%d, sNN=%d, fld=%d",
 211          runNo, year, sys, sNN, fld);
 212
 213   // --- Output object -----------------------------------------------
 214   // Make our correction object
 215   AliFMDCorrAcceptance* corr = new AliFMDCorrAcceptance();
 216   corr->SetVertexAxis(nVtxBins, vtxLow, vtxHigh);
 217
 218   // --- Output script -----------------------------------------------
 219   std::ofstream deadScript("deadstrips.C");
 220   deadScript << "// Automatically generaeted by ExtractAcceptance.C\n"
 221              << "// Add additional dead strips to sharing filter\n"
 222              << "// Information taken from OCDB entry for\n"
 223              << "//\n"
 224              << "//    run       = " << runNo << "\n"
 225              << "//    year      = " << year << "\n"
 226              << "//    system    = " << sys << "\n"
 227              << "//    sqrt{sNN} = " << sNN << "GeV\n"
 228              << "//    L3 field  = " << fld << "kG\n"
 229              << "void deadstrips(AliFMDSharingFilter* filter)\n"
 230              << "{" << std::endl;
 231
 232   // --- Loop over verticies and rings -------------------------------
 233   Int_t  nDead   = 0;
 234   Bool_t gotDead = false;
 235   Float_t dV     = (vtxHigh - vtxLow) / nVtxBins;
 236   Printf("=== Looping over vertices: %d bins from %+6.2f to %+6.2f "
 237          "in steps of %+6.2f",
 238          nVtxBins, vtxLow, vtxHigh, dV);
 239   for (Double_t v = vtxLow+dV/2; v < vtxHigh; v += dV) {
 240     for(UShort_t d = 1; d <= 3;d++) {
 241       UShort_t nR = (d == 1 ? 1 : 2);
 242       for (UShort_t q = 0; q < nR; q++) {
 243         Char_t   r  = (q == 0 ? 'I' : 'O');
 244
 245         // Delegate to other function
 246         Int_t nLocal = 0;
 247         TH2D* ratio = MakeOneRing(d, r, v, nLocal,
 248                                   !gotDead ? &deadScript : 0);
 249         nDead += nLocal;
 250         if (!ratio) {
 251           Warning("ExtractAcceptance", "Didn't get correction from "
 252                   "FMD%d%c @ vz=%+6.2fcm", d, r, v);
 253           continue;
 254         }
 255         // Printf("v=%+6.2f FMD%d%c, got %d dead strips", v, d, r, nLocal);
 256
 257         // Set the correction
 258         corr->SetCorrection(d, r, v, ratio);
 259       }
 260     }
 261     gotDead = true;
 262   }
 263   corr->SetHasOverflow();
 264   corr->Print();
 265   // corr->ls();
 266
 267   deadScript << "}\n"
 268              << "// EOF" << std::endl;
 269   deadScript.close();
 270
 271   // Write to a file
 272   Printf("=== Writing to disk");
 273   AliForwardCorrectionManager& cm = AliForwardCorrectionManager::Instance();
 274   if (!cm.Store(corr, runNo, sys, sNN, fld, false, false,
 275                 "fmd_corrections.root")) {
 276     Error("", "Failed to store acceptance correction in local file");
 277     return;
 278   }
 279
 280   std::ofstream f("Upload.C");
 281   f << "// Generated by ExtractELoss.C\n"
 282     << "TString MakeDest(const TString& dest, const TString& fname)\n"
 283     << "{\n"
 284     << "  TString tmp(dest);\n"
 285     << "  if (!tmp.IsNull()) {\n"
 286     << "    if (!tmp.EndsWith(\"/\")) tmp.Append(\"/\");\n"
 287     << "    tmp.Append(fname);\n"
 288     << "  }\n"
 289     << "  return tmp;\n"
 290     << "}\n\n"
 291     << "void Upload(const TString& dest=\"\")\n"
 292     << "{\n"
 293     << "  gROOT->Macro(\"" << fwd << "/scripts/LoadLibs.C\");\n"
 294     << "  \n"
 295     << "  const char* fmdFile = \"fmd_corrections.root\";\n"
 296     << "  TString fdest = MakeDest(dest, fmdFile);\n"
 297     << "  \n"
 298     << "  AliForwardCorrectionManager::Instance().Append(fmdFile, fdest);\n"
 299     << "}\n"
 300     << "// EOF\n"
 301     << std::endl;
 302   f.close();
 303
 304   Info("ExtracAcceptance",
 305        "Run generated Upload.C(DEST) script to copy files in place");
 306 }
 307
 308 //
 309 // EOF
 310 //
 311
 312