PWG3/dielectron/macros/FitCDFLocal.C

   1 TH1F *psproperMCsecJPSI();
   2 void LoadLib();
   3
   4 void FitCDFLocal(){
   5
   6   ///////////////////////////////////////////////////////////////////
   7   //
   8   // Example macro to read local N-Tuples of JPSI
   9   // and bkg candidates and perform log-likelihood
  10   // minimization using the minimization handler class
  11   //
  12   // Origin: C. Di Giglio
  13   //
  14   ///////////////////////////////////////////////////////////////////
  15   LoadLib();
  16   AliDielectronBtoJPSItoEleCDFfitFCN *likely_obj = 0x0;
  17
  18   // Those initial parameters are evaluated from the fits done separately on
  19   // invariant mass and pseudoproper distributions. For the moment are setted by
  20   // hand
  21
  22   Double_t *resParam = new Double_t[6];
  23   resParam[0] = 2586.32; resParam[1]= -0.163616; resParam[2] = 31.026770; resParam[3] = 433.85;
  24   resParam[4] = 4.67; resParam[5] = 73.31;
  25
  26   Double_t *bkgParam  = new Double_t[7];
  27   bkgParam[0] = 7063.8; bkgParam[1] = 23214.45; bkgParam[2] = 24742.65; bkgParam[3] = 15188.57;
  28   bkgParam[4] = 0.009579; bkgParam[5] = 0.010422; bkgParam[6] = 0.000423;
  29
  30   Double_t *invMassParam = new Double_t[9];
  31   invMassParam[0] = 3.092; invMassParam[1] =  0.0306; invMassParam[2] =  0.4939;
  32   invMassParam[3] = 1.36; invMassParam[4] = 41.97; invMassParam[5] = 80.43; invMassParam[6] = 0.975;
  33   invMassParam[7] = 1.61; invMassParam[8] = 2.129;
  34
  35   Double_t Fsig = 0.454;
  36
  37   TH1F* hCsiMCPithya = new TH1F();
  38
  39   hCsiMCPithya = psproperMCsecJPSI();
  40   Double_t integral = 0;
  41   for(int i=1;i<hCsiMCPithya->GetNbinsX()+1; i++) integral += (hCsiMCPithya->GetBinContent(i)*hCsiMCPithya->GetBinWidth(i));
  42   hCsiMCPithya->Scale(1./integral);
  43
  44   Double_t* x=0x0; Double_t* m=0x0; Int_t n=0;
  45   AliDielectronBtoJPSItoEle* aBtoJPSItoEle =new AliDielectronBtoJPSItoEle();
  46   Double_t paramInputValues [20] =
  47                                    /*paramInputValues[0]   ----> fWeightRes;
  48                                      paramInputValues[1]   ----> fPos;
  49                                      paramInputValues[2]   ----> fNeg;
  50                                      paramInputValues[3]   ----> fSym;
  51                                      paramInputValues[4]   ----> fOneOvLamPlus;
  52                                      paramInputValues[5]   ----> fOneOvLamMinus;
  53                                      paramInputValues[6]   ----> fOneOvLamSym;
  54                                      paramInputValues[7]   ----> fB;
  55                                      paramInputValues[8]   ----> fFsig;
  56                                      paramInputValues[9]   ----> fMmean;
  57                                      paramInputValues[10]  ----> fNexp;
  58                                      paramInputValues[11]  ----> fSigma;
  59                                      paramInputValues[12]  ----> fAlpha;
  60                                      paramInputValues[13]  ----> fNorm;
  61                                      paramInputValues[14]  ----> fBkgNorm;
  62                                      paramInputValues[15]  ----> fBkgMean;
  63                                      paramInputValues[16]  ----> fBkgSlope;
  64                                      paramInputValues[17]  ----> fBkgConst;
  65                                      paramInputValues[18]  ----> fGaus1Norm;
  66                                      paramInputValues[19]  ----> fGaus2Norm;*/
  67                                    { bkgParam[9],  // fWeightRes [0]
  68                                      bkgParam[10], // Fpos [1]
  69                                      bkgParam[11], // FNeg [2]
  70                                      bkgParam[12], // FSym [3]
  71                                      bkgParam[6], // LamdaPos [4]
  72                                      bkgParam[7], // LambdaNeg [5]
  73                                      bkgParam[8], // LamdaSym [6]
  74                                      0.10, // fB [7]
  75                                      Fsig, // FSig [8]
  76                                      invMassParam[0], // FMean [9]
  77                                      invMassParam[4], // fNexp [10]
  78                                      invMassParam[1], // fNsigma [11]
  79                                      invMassParam[2], // fAlpha [12]
  80                                      invMassParam[3], //fNorm [13]
  81                                      invMassParam[5], // fBkgNorm [14]
  82                                      invMassParam[6], // fBkgMean [15]
  83                                      invMassParam[7], // fBkgSlope [16]
  84                                      invMassParam[8], // fBkgConst [17]
  85                                      resParam[0], // Gaus1Norm [18]
  86                                      resParam[3] // Gaus2Norm [19]
  87                                       }; // Starting values for parameters
  88
  89   // retrieve the TNtupla from the file result.root
  90   TFile f("result.root");
  91   TList * list = (TList*)f.Get("resultAOD");
  92   TNtuple *nt=(TNtuple*)list->FindObject("fNtupleJPSI");
  93   nt->ls();
  94   aBtoJPSItoEle->ReadCandidates(nt,x,m,n); // read N-Tuples
  95   printf("+++\n+++ Number of total candidates (prim J/psi + secondary J/psi + bkg) ---> %d candidates \n+++\n",n);
  96
  97   aBtoJPSItoEle->SetFitHandler(x,m,n); // Set the fit handler with given values of x, m, # of candidates
  98
  99   aBtoJPSItoEle->CloneMCtemplate(hCsiMCPithya);    // clone MC template and copy internally
 100                                                    // in this way any model can be setted from outside
 101
 102   aBtoJPSItoEle->SetCsiMC(); // Pass the MC template to the CDF fit function
 103
 104   AliDielectronBtoJPSItoEleCDFfitHandler* fithandler = aBtoJPSItoEle->GetCDFFitHandler(); // Get the fit handler
 105
 106   //
 107   // Set some fit options through the handler class
 108   //
 109   fithandler->SetErrorDef(0.5); // tells Minuit that the error interval is the one in which
 110                                 // the function differs from the minimum for less than setted value
 111
 112   fithandler->SetPrintStatus(kTRUE);
 113
 114   fithandler->SetParamStartValues(paramInputValues);
 115
 116   Double_t *resConst = new Double_t[4];
 117
 118   //resolution constants: fixed from MC
 119   resConst[0]= resParam[1]; // mean1
 120   resConst[1]= resParam[2]; // sigma1
 121   resConst[2]= resParam[4]; // mean2
 122   resConst[3]= resParam[5]; // sigma2
 123
 124
 125   fithandler->SetResolutionConstants(resConst);
 126   fithandler->SetCrystalBallFunction(kTRUE);
 127
 128   // fix parameter
 129   fithandler->FixParam(0,kTRUE);
 130   fithandler->FixParam(1,kTRUE); // Fix Bkg  weights
 131   fithandler->FixParam(2,kTRUE);
 132   fithandler->FixParam(3,kTRUE);
 133
 134   fithandler->FixParam(4,kTRUE);
 135   fithandler->FixParam(5,kTRUE); // Fix bkg exponential
 136   fithandler->FixParam(6,kTRUE);
 137   //fithandler->FixParam(8,kTRUE); // Fsig fix
 138
 139   fithandler->FixParam(9,kTRUE);
 140   fithandler->FixParam(10,kTRUE); // Fix CristalBall param
 141   fithandler->FixParam(11,kTRUE);
 142   fithandler->FixParam(12,kTRUE);
 143   fithandler->FixParam(13,kTRUE);
 144
 145   fithandler->FixParam(14,kTRUE);
 146   fithandler->FixParam(15,kTRUE); // Fix Bkg invMass param
 147   fithandler->FixParam(16,kTRUE);
 148   fithandler->FixParam(17,kTRUE);
 149
 150   fithandler->FixParam(18,kTRUE); // Fix resolution weights
 151   fithandler->FixParam(19,kTRUE);
 152
 153   //invariant mass function normalized between invMassMin, invMassMax  -->
 154   Double_t invMassMin = 2.0; // GeV
 155   Double_t invMassMax = 4.0; // GeV
 156   Double_t massLow = (TDatabasePDG::Instance()->GetParticle(443)->Mass()) - invMassMin;
 157   Double_t massHigh = invMassMax - (TDatabasePDG::Instance()->GetParticle(443)->Mass());
 158   fithandler->SetMassWndLow(massLow);
 159   fithandler->SetMassWndHigh(massHigh);
 160
 161   likely_obj = fithandler->LikelihoodPointer();
 162   likely_obj->SetAllParameters(paramInputValues);
 163   likely_obj->ComputeMassIntegral();
 164   likely_obj->PrintStatus();
 165
 166   aBtoJPSItoEle->DoMinimization();
 167   f.Close();
 168   return;
 169 }
 170
 171
 172 TH1F *psproperMCsecJPSI(){
 173    TH1F *hDecayTimeMCjpsifromB = new TH1F("hDecayTimeMCjpsifromB","B --> J/#Psi MC pseudo proper decay length",150,-3000,3000
 174 );
 175    hDecayTimeMCjpsifromB->SetBinContent(67,2);
 176    hDecayTimeMCjpsifromB->SetBinContent(68,2);
 177    hDecayTimeMCjpsifromB->SetBinContent(69,2);
 178    hDecayTimeMCjpsifromB->SetBinContent(70,4);
 179    hDecayTimeMCjpsifromB->SetBinContent(71,5);
 180    hDecayTimeMCjpsifromB->SetBinContent(72,17);
 181    hDecayTimeMCjpsifromB->SetBinContent(73,35);
 182    hDecayTimeMCjpsifromB->SetBinContent(74,96);
 183    hDecayTimeMCjpsifromB->SetBinContent(75,439);
 184    hDecayTimeMCjpsifromB->SetBinContent(76,11428);
 185    hDecayTimeMCjpsifromB->SetBinContent(77,9094);
 186    hDecayTimeMCjpsifromB->SetBinContent(78,7769);
 187    hDecayTimeMCjpsifromB->SetBinContent(79,6475);
 188    hDecayTimeMCjpsifromB->SetBinContent(80,5695);
 189    hDecayTimeMCjpsifromB->SetBinContent(81,4884);
 190    hDecayTimeMCjpsifromB->SetBinContent(82,4294);
 191    hDecayTimeMCjpsifromB->SetBinContent(83,3780);
 192    hDecayTimeMCjpsifromB->SetBinContent(84,3321);
 193    hDecayTimeMCjpsifromB->SetBinContent(85,2827);
 194    hDecayTimeMCjpsifromB->SetBinContent(86,2531);
 195    hDecayTimeMCjpsifromB->SetBinContent(87,2256);
 196    hDecayTimeMCjpsifromB->SetBinContent(88,2099);
 197    hDecayTimeMCjpsifromB->SetBinContent(89,1782);
 198    hDecayTimeMCjpsifromB->SetBinContent(90,1592);
 199    hDecayTimeMCjpsifromB->SetBinContent(91,1478);
 200    hDecayTimeMCjpsifromB->SetBinContent(92,1286);
 201    hDecayTimeMCjpsifromB->SetBinContent(93,1145);
 202    hDecayTimeMCjpsifromB->SetBinContent(94,980);
 203    hDecayTimeMCjpsifromB->SetBinContent(95,933);
 204    hDecayTimeMCjpsifromB->SetBinContent(96,865);
 205    hDecayTimeMCjpsifromB->SetBinContent(97,761);
 206    hDecayTimeMCjpsifromB->SetBinContent(98,654);
 207    hDecayTimeMCjpsifromB->SetBinContent(99,622);
 208    hDecayTimeMCjpsifromB->SetBinContent(100,587);
 209    hDecayTimeMCjpsifromB->SetBinContent(101,572);
 210    hDecayTimeMCjpsifromB->SetBinContent(102,449);
 211    hDecayTimeMCjpsifromB->SetBinContent(103,416);
 212    hDecayTimeMCjpsifromB->SetBinContent(104,346);
 213    hDecayTimeMCjpsifromB->SetBinContent(105,361);
 214    hDecayTimeMCjpsifromB->SetBinContent(105,361);
 215    hDecayTimeMCjpsifromB->SetBinContent(106,315);
 216    hDecayTimeMCjpsifromB->SetBinContent(107,265);
 217    hDecayTimeMCjpsifromB->SetBinContent(108,239);
 218    hDecayTimeMCjpsifromB->SetBinContent(109,247);
 219    hDecayTimeMCjpsifromB->SetBinContent(110,189);
 220    hDecayTimeMCjpsifromB->SetBinContent(111,223);
 221    hDecayTimeMCjpsifromB->SetBinContent(112,174);
 222    hDecayTimeMCjpsifromB->SetBinContent(113,184);
 223    hDecayTimeMCjpsifromB->SetBinContent(114,171);
 224    hDecayTimeMCjpsifromB->SetBinContent(115,153);
 225    hDecayTimeMCjpsifromB->SetBinContent(116,151);
 226    hDecayTimeMCjpsifromB->SetBinContent(117,126);
 227    hDecayTimeMCjpsifromB->SetBinContent(118,105);
 228    hDecayTimeMCjpsifromB->SetBinContent(119,98);
 229    hDecayTimeMCjpsifromB->SetBinContent(120,88);
 230    hDecayTimeMCjpsifromB->SetBinContent(121,92);
 231    hDecayTimeMCjpsifromB->SetBinContent(122,86);
 232    hDecayTimeMCjpsifromB->SetBinContent(123,73);
 233    hDecayTimeMCjpsifromB->SetBinContent(124,78);
 234    hDecayTimeMCjpsifromB->SetBinContent(125,63);
 235    hDecayTimeMCjpsifromB->SetBinContent(126,69);
 236    hDecayTimeMCjpsifromB->SetBinContent(127,55);
 237    hDecayTimeMCjpsifromB->SetBinContent(128,43);
 238    hDecayTimeMCjpsifromB->SetBinContent(129,46);
 239    hDecayTimeMCjpsifromB->SetBinContent(130,38);
 240    hDecayTimeMCjpsifromB->SetBinContent(131,34);
 241    hDecayTimeMCjpsifromB->SetBinContent(132,33);
 242    hDecayTimeMCjpsifromB->SetBinContent(133,51);
 243    hDecayTimeMCjpsifromB->SetBinContent(134,25);
 244    hDecayTimeMCjpsifromB->SetBinContent(135,44);
 245    hDecayTimeMCjpsifromB->SetBinContent(136,24);
 246    hDecayTimeMCjpsifromB->SetBinContent(137,25);
 247    hDecayTimeMCjpsifromB->SetBinContent(138,21);
 248    hDecayTimeMCjpsifromB->SetBinContent(139,23);
 249    hDecayTimeMCjpsifromB->SetBinContent(140,16);
 250    hDecayTimeMCjpsifromB->SetBinContent(141,24);
 251    hDecayTimeMCjpsifromB->SetBinContent(142,14);
 252    hDecayTimeMCjpsifromB->SetBinContent(143,21);
 253    hDecayTimeMCjpsifromB->SetBinContent(144,22);
 254    hDecayTimeMCjpsifromB->SetBinContent(145,22);
 255    hDecayTimeMCjpsifromB->SetBinContent(146,13);
 256    hDecayTimeMCjpsifromB->SetBinContent(147,10);
 257    hDecayTimeMCjpsifromB->SetBinContent(148,15);
 258    hDecayTimeMCjpsifromB->SetBinContent(149,12);
 259    hDecayTimeMCjpsifromB->SetBinContent(150,12);
 260    hDecayTimeMCjpsifromB->SetBinContent(151,231);
 261    return hDecayTimeMCjpsifromB;
 262 }
 263
 264 void LoadLib(){
 265   gSystem->Load("libTree.so");
 266   gSystem->Load("libGeom.so");
 267   gSystem->Load("libPhysics.so");
 268   gSystem->Load("libVMC.so");
 269   gSystem->Load("libMinuit.so");
 270   gSystem->Load("libSTEERBase.so");
 271   gSystem->Load("libESD.so");
 272   gSystem->Load("libAOD.so");
 273   gSystem->Load("libANALYSIS.so");
 274   gSystem->Load("libANALYSISalice.so");
 275   gSystem->Load("libCORRFW.so");
 276   gSystem->Load("libPWG3dielectron.so");
 277  }