09aded15d3162a4251d3000eac8712466214b084
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / scripts / RunCopyELossFit.C
1 Color_t Color(UShort_t d, Char_t r ) const 
2
3   return ((d == 1 ? kRed : (d == 2 ? kGreen : kBlue))
4           + ((r == 'I' || r == 'i') ? 2 : -2));
5 }
6
7 /** 
8  * 
9  * @param sys       Collision system 
10  * @param cms       Center of mass energy per nucleon in GeV
11  * @param field     Magnetic field 
12  * 
13  * @ingroup pwg2_forward_analysis_scripts
14  */
15 void
16 RunCopyELossFit(UShort_t sys, UShort_t cms, Short_t field, bool mc=false)
17 {
18   gROOT->Macro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadLibs.C");
19   gSystem->Load("libPWG2forward.so");
20
21   AliFMDAnaParameters* p = AliFMDAnaParameters::Instance();
22   p->SetRealData(!mc);
23   p->SetEnergy(Float_t(cms));
24   p->SetMagField(Float_t(field));
25   p->SetCollisionSystem(sys == 1 ? "pp" : "PbPb");
26   p->Init(true, AliFMDAnaParameters::kBackgroundCorrection|
27           AliFMDAnaParameters::kEnergyDistributions);
28   
29   Int_t    nVtx   = p->GetNvtxBins();
30   Double_t minVtx = -p->GetVtxCutZ();
31   Double_t maxVtx = -p->GetVtxCutZ();
32   Int_t    nEta   = p->GetNetaBins();
33   Double_t minEta = p->GetEtaMin();
34   Double_t maxEta = p->GetEtaMax();
35   
36   TAxis vtxAxis(nVtx, minVtx, maxVtx);
37   TAxis etaAxis(nEta, minEta, maxEta);
38   AliFMDCorrELossFit* m = new AliFMDCorrELossFit;
39   m->SetEtaAxis(etaAxis);
40   
41   for (UShort_t d = 1; d <= 3; d++) { 
42     UShort_t nQ = (d == 1 ? 1 : 2);
43     for (UShort_t q = 0; q < nQ; q++) { 
44       Char_t r = (q == 0 ? 'I' : 'O');
45       
46       for (UShort_t b = 1; b < nEta; b++) { 
47         Double_t eta = etaAxis.GetBinCenter(b);
48         Info("RunCopyELossFit", "FMD%d%c, bin %3d, eta %+8.4f", d, r, b, eta);
49         if (eta < minEta || eta > maxEta) continue;
50         TH1F* oldhist = p->GetEnergyDistribution(d,r,eta);
51         if (!oldhist) {
52           Warning("RunCopyELossFit",
53                   "Didn't find energy distribution for FMD%d%c, eta bin %3d", 
54                   d, r, b);
55           continue;
56         }
57         TF1*  oldfunc = oldhist->GetFunction("FMDfitFunc");
58         if (!oldfunc) {
59           // Warning("RunCopyELossFit",
60           // "Didn't find energy fit for FMD%d%c, eta bin %3d", 
61           // d, r, b);
62           continue;
63         }
64         Double_t chi2   = oldfunc->GetChisquare();
65         UShort_t nu     = oldfunc->GetNDF();
66         Double_t c      = oldfunc->GetParameter(0);
67         Double_t delta  = oldfunc->GetParameter(1);
68         Double_t xi     = oldfunc->GetParameter(2);
69         Double_t a2 = 0, a3 = 0, ea2 = 0, ea3 = 0;
70         if (oldfunc->GetNpar() >= 4) {
71           a2     = oldfunc->GetParameter(3);
72           ea2    = oldfunc->GetParError(3);
73         }
74         if (oldfunc->GetNpar() >= 5) { 
75           a3     = oldfunc->GetParameter(4);
76           ea3    = oldfunc->GetParError(4);
77         }
78         Int_t    n      = (a3 < 1e-5 ? (a2 < 1e-5 ? 1 : 2) : 3);
79         Double_t ec     = oldfunc->GetParError(0);
80         Double_t edelta = oldfunc->GetParError(1);
81         Double_t exi    = oldfunc->GetParError(2);
82         delta           = delta - xi * -0.22278298;
83         edelta          = TMath::Sqrt(TMath::Power(edelta,2) - 
84                                       TMath::Power(0.22278298*exi,2));
85         Info("RunCopyELossFit", "FMD%d%c etabin=%3d: n=%d\n" 
86              "  C=%f+/-%f, Delta=%f+/-%f, xi=%f+/-%f, a2=%f+/-%f, a3=%f+/-%f",
87              d, r, b, n, c, ec, delta, edelta, xi, exi, a2, ea2, a3, ea3);
88         Double_t a[]  = { a2, a3, -9999 };
89         Double_t ea[] = { ea2, ea3, -9999 };
90
91         AliFMDCorrELossFit::ELossFit* fit = new 
92           AliFMDCorrELossFit::ELossFit(0, n, 
93                                        chi2,   nu, 
94                                        c,      ec, 
95                                        delta,  edelta, 
96                                        xi,     exi, 
97                                        0,      0, 
98                                        0,      0, 
99                                        a,      ea);
100         fit->CalculateQuality();
101         fit->fBin = b;
102         fit->fDet = d;
103         fit->fRing = r;
104         fit->Print();
105         
106         Info("RunCopyELossFit", "Setting fit FMD%d%c etabin=%3d: %p", 
107              d, r, b, fit);
108         m->SetFit(d, r, b, fit);
109       }
110     }
111   }
112   AliForwardCorrectionManager& mgr = AliForwardCorrectionManager::Instance();
113   TString fname(mgr.GetFileName(AliForwardCorrectionManager::kELossFits,
114                                 sys, cms, field, mc));
115   TFile* output = TFile::Open(fname.Data(), "RECREATE");
116   if (!output) { 
117     Warning("RunCopyELossFit", "Failed to open output file %s", fname.Data());
118     return kFALSE;
119   }
120   m->Write(mgr.GetObjectName(AliForwardCorrectionManager::kELossFits));
121   output->Write();
122   output->Close();
123   Info("RunCopyELossFit", 
124        "File %s created.  It should be copied to %s and stored in SVN",
125        fname.Data(), mgr.GetFileDir(AliForwardCorrectionManager::kELossFits));
126   
127 }
128 //
129 // EOF
130 //