New Calorimeter Calibration directory, now only PHOS
[u/mrichter/AliRoot.git] / PWG4 / CaloCalib / macros / 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 }