PHOS/macros/pi0Calib/pi0Calib.C

   1 Double_t rpol2(Double_t *x, Double_t *par)
   2 {
   3   Double_t rp2 = par[0] + par[1]*x[0] + par[2]*x[0]*x[0];
   4   return rp2;
   5 }
   6
   7 TF1* InitPol2()
   8 {
   9   TF1* f1 = new TF1("rpol2",rpol2,0.,300.,3);
  10   f1->SetParNames("p0","p1","p2");
  11
  12   f1->SetParameter("p0",1.);
  13   f1->SetParameter("p1",1.);
  14   f1->SetParameter("p2",1.);
  15
  16   return f1;
  17 }
  18
  19 Double_t gausPol2(Double_t *x, Double_t *par)
  20 {
  21   //gaus+second order polynomial.
  22
  23   Double_t gp2 = par[0] * TMath::Exp( -(x[0]-par[1])*(x[0]-par[1]) /
  24                                       (2*par[2]*par[2]) ) /
  25     + par[3] + par[4]*x[0] + par[5]*x[0]*x[0];
  26
  27   return gp2;
  28 }
  29
  30
  31
  32 void pi0Calib(const char* file="Sum_All_Emc.root")
  33 {
  34
  35   AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT");
  36   AliCDBManager::Instance()->SetSpecificStorage("PHOS/*","local://./");
  37
  38   AliPHOSCalibData cdb(0);
  39
  40   Int_t fMinEntr = 50; // min 50 pi0 canditates per cell
  41   Double_t xmin = 100; // fit range lower end
  42   Double_t xmax = 170; // fit range upper end
  43
  44   const Float_t pi0m = 134.98;
  45   Double_t pi0f;
  46   const Int_t fMod=2;
  47   char hname[128];
  48   TH1F* h1=0;
  49   Double_t cc_i, cc_ii;
  50
  51   TH1F* hPi0Mass = new TH1F("pi0_mass","#pi^{0} mass in calibrated cells",100,0.,300.);
  52
  53   TFile f(file);
  54
  55   for(Int_t iMod=0; iMod<5; iMod++) {
  56     for(Int_t iX=0; iX<64; iX++) {
  57       for(Int_t iZ=0; iZ<56; iZ++) {
  58
  59         sprintf(hname,"%d_%d_%d",iMod,iX,iZ);
  60         h1 = (TH1F*)f.Get(hname);
  61         if(!h1) continue;
  62         if(h1->GetEntries()<fMinEntr) continue;
  63
  64         TF1* rp2 = InitPol2();
  65         rp2->SetRange(180,260);
  66 //      rp2->SetRange(50.,110.);
  67         h1->Fit(rp2,"RNQ");
  68         rp2->Print();
  69
  70         h1->GetXaxis()->SetRange(17,210);
  71         Double_t max = h1->GetBinCenter(h1->GetMaximumBin()); // peak
  72
  73         TF1* f1 = new TF1("gpol2",gausPol2,xmin,xmax,6);
  74         f1->SetParNames("Constant","Mean","Sigma","p0","p1","p2");
  75
  76         printf("Trying to fit %s, %d entries, peak at %.3f\n",
  77                h1->GetName(),h1->GetEntries(),max);
  78
  79         f1->SetParameter("Constant",h1->GetMaximum());
  80         f1->SetParameter("Mean",max);
  81         f1->SetParameter("Sigma",1.);
  82         f1->SetParameter("p0",rp2->GetParameter("p0"));
  83         f1->SetParameter("p1",rp2->GetParameter("p1"));
  84         f1->SetParameter("p2",rp2->GetParameter("p2"));
  85
  86         f1->SetLineColor(kBlue);
  87         h1->Fit(f1,"R+");
  88
  89         pi0f = f1->GetParameter("Mean");
  90         hPi0Mass->Fill(pi0f);
  91
  92         //CC from the previous iteration.
  93         cc_i = cdb.GetADCchannelEmc(iMod+1,iZ+1,iX+1);
  94
  95         if(pi0f>100. && pi0f<200.) {
  96           Double_t k = pi0m/pi0f;
  97           cc_ii = cc_i*(1+k**2)/2.;
  98           cdb.SetADCchannelEmc(iMod+1,iZ+1,iX+1,cc_ii);
  99         }
 100
 101         printf("%s: %d entries, mean=%.3f, peak=%.3f, rms= %.3f. pi0 = %.3f,
 102              cc_i=%f, cc_ii=%f\n",
 103                h1->GetName(),h1->GetEntries(),h1->GetMean(),max,h1->GetRMS(),
 104                pi0f,cc_i,cc_ii);
 105
 106       }
 107     }
 108   }
 109
 110   TFile ff("pi0Calib.root","recreate");
 111   hPi0Mass->Write();
 112
 113   //Writing new calibration coefficients (CCs) to OCDB
 114
 115   AliCDBMetaData *md= new AliCDBMetaData();
 116   md->SetResponsible("Boris Polishchuk");
 117
 118   cdb.WriteEmc(0,999999,md);
 119   cdb.WriteCpv(0,999999,md);
 120
 121 }