PWG0/dNdEta/systematicsAnalysis.C

   1 changeProductionCrossSections() {
   2   //
   3   // calculates dN/dEta from the standard correction and from
   4   // corrections that has been evaluated with different relative cross
   5   // sections (of single diff, double diff and non diff). The ratios
   6   // (standard to changed x-section) of the different dN/deta
   7   // distributions are saved to a file.
   8   //
   9
  10   gSystem->Load("libPWG0base");
  11
  12   // folder names in systematics root file
  13   const Char_t* folderNames[] = {"triggerBiasDD","triggerBiasSD","triggerBiasND","vertexRecoDD","vertexRecoSD","vertexRecoND"};
  14
  15   const Char_t* changes[]  = {"pythia","ddmore","ddless","sdmore","sdless", "dmore", "dless"};
  16
  17   Float_t scalesDD[] = {1.0, 1.5, 0.5, 1.0, 1.0, 1.5, 0.5};
  18   Float_t scalesSD[] = {1.0, 1.0, 1.0, 1.5, 0.5, 1.5, 0.5};
  19
  20   // cross section from Pythia
  21   //Float_t sigmaND = 55.2;
  22   //Float_t sigmaDD = 9.78;
  23   //Float_t sigmaSD = 14.30;
  24
  25   // getting data
  26   TFile* finCorr = TFile::Open("analysis_esd.root");
  27   dNdEtaAnalysis* analysis = new dNdEtaAnalysis("dndeta", "dndeta");
  28   analysis->LoadHistograms();
  29
  30   // getting corrections
  31   AlidNdEtaCorrection* correctionStandard = new AlidNdEtaCorrection("dndeta_correction","dndeta_correction");
  32
  33   TFile* finCorr = TFile::Open("correction_map.root");
  34   correctionStandard->LoadHistograms();
  35
  36   // getting corrections for pure SD, DD and ND
  37   TFile* finSyst = TFile::Open("systematics.root");
  38   AlidNdEtaCorrection* correctionsSyst[7];
  39   for (Int_t i=0; i<6; i++) {
  40     correctionsSyst[i] = new AlidNdEtaCorrection(folderNames[i],folderNames[i]);
  41     correctionsSyst[i]->LoadHistograms();
  42   }
  43
  44   TH1F* hRatios[21];
  45   for (Int_t j=0; j<3; j++) { // j = 0 (change vtx), j = 1 (change trg), j = 2 (change both)
  46
  47     for (Int_t i=0; i<7; i++) { // changes in cross section
  48
  49       // temporary correction with changed cross sections
  50       AlidNdEtaCorrection* correctionTmp = (AlidNdEtaCorrection*)correctionStandard->Clone();;
  51
  52       correctionTmp->Reset();
  53
  54       if (j==0 || j==2) {
  55         //correctionTmp->GetVertexRecoCorrection()->Reset();
  56         //correctionTmp->GetVertexRecoCorrection()->Add(correctionsSyst[3]->GetVertexRecoCorrection(),scalesDD[i]);
  57         //correctionTmp->GetVertexRecoCorrection()->Add(correctionsSyst[4]->GetVertexRecoCorrection(),scalesSD[i]);
  58         //correctionTmp->GetVertexRecoCorrection()->Add(correctionsSyst[5]->GetVertexRecoCorrection(),1);
  59
  60         correctionTmp->Add(correctionsSyst[3], scalesDD[i]);
  61         correctionTmp->Add(correctionsSyst[4], scalesSD[i]);
  62         correctionTmp->Add(correctionsSyst[5], 1);
  63       }
  64       if (j==1 || j==2) {
  65         //correctionTmp->GetTriggerBiasCorrectionINEL()->Reset();
  66         //correctionTmp->GetTriggerBiasCorrectionINEL()->Add(correctionsSyst[0]->GetTriggerBiasCorrectionINEL(),scalesDD[i]);
  67         //correctionTmp->GetTriggerBiasCorrectionINEL()->Add(correctionsSyst[1]->GetTriggerBiasCorrectionINEL(),scalesSD[i]);
  68         //correctionTmp->GetTriggerBiasCorrectionINEL()->Add(correctionsSyst[2]->GetTriggerBiasCorrectionINEL(),1);
  69
  70         correctionTmp->Add(correctionsSyst[0], scalesDD[i]);
  71         correctionTmp->Add(correctionsSyst[1], scalesSD[i]);
  72         correctionTmp->Add(correctionsSyst[2], 1);
  73       }
  74       correctionTmp->Finish();
  75
  76       dNdEtaAnalysis* analysisTmp = (dNdEtaAnalysis*)analysis->Clone();
  77       analysisTmp->Finish(correctionTmp, 0.3, AlidNdEtaCorrection::kINEL);
  78
  79       Int_t counter = i + j*7;
  80
  81       hRatios[counter] = (TH1F*)analysisTmp->GetdNdEtaHistogram(2)->Clone();
  82
  83       TString name("ratio");
  84       if (j==0) name.Append("_vetexReco_");
  85       if (j==1) name.Append("_triggerBias_");
  86       if (j==2) name.Append("_vertexReco_triggerBias_");
  87       name.Append(changes[i]);
  88       hRatios[counter]->SetName(name.Data());
  89       name.Append(Form(" (DD #times %0.1f, SD #times %0.1f)",scalesDD[i],scalesSD[i]));
  90       hRatios[counter]->SetTitle(name.Data());
  91       hRatios[counter]->SetYTitle("ratio (Pythia x-sections)/(changed x-sections)");
  92
  93       delete analysisTmp;
  94       delete correctionTmp;
  95     }
  96   }
  97
  98   for (Int_t i=1; i<21; i++)
  99     hRatios[i]->Divide(hRatios[0],hRatios[i],1,1);
 100   hRatios[0]->Divide(hRatios[0],hRatios[0],1,1);
 101
 102   TFile* fout = new TFile("systematics_xsections.root","RECREATE");
 103
 104   for (Int_t i=0; i<21; i++)
 105     hRatios[i]->Write();
 106
 107   fout->Write();
 108   fout->Close();
 109 }
 110
 111 changeParticleComposition() {
 112   //
 113   // calculates dN/dEta from the standard correction and from
 114   // corrections that has been evaluated with different relative
 115   // abundancies of kaons and protons and save the ratios in a file
 116
 117   gSystem->Load("libPWG0base");
 118
 119   const Char_t* folderNames[] = {"correction_0","correction_1","correction_2","correction_3"};
 120   Float_t  scalesPi[] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
 121   Float_t  scalesK[]  = {0.5, 1.5, 1.0, 1.0, 0.5, 1.5};
 122   Float_t  scalesP[]  = {1.0, 1.0, 0.5, 1.5, 0.5, 1.5};
 123
 124   // getting data
 125   TFile* finCorr = TFile::Open("analysis_esd.root");
 126   dNdEtaAnalysis* analysis = new dNdEtaAnalysis("dndeta", "dndeta");
 127   analysis->LoadHistograms();
 128
 129   // getting corrections for pi, K, p and other particles
 130   TFile* finSyst = TFile::Open("systematics.root");
 131   AlidNdEtaCorrection* correctionsSyst[4];
 132   for (Int_t i=0; i<4; i++) {
 133     correctionsSyst[i] = new AlidNdEtaCorrection(folderNames[i],folderNames[i]);
 134     correctionsSyst[i]->LoadHistograms();
 135   }
 136   AlidNdEtaCorrection* correctionAll = (AlidNdEtaCorrection*) correctionsSyst[0]->Clone();
 137   correctionAll->Reset();
 138   correctionAll->Add(correctionsSyst[0]);
 139   correctionAll->Add(correctionsSyst[1]);
 140   correctionAll->Add(correctionsSyst[2]);
 141   correctionAll->Add(correctionsSyst[3]);
 142   correctionAll->Finish();
 143   analysis->Finish(correctionAll, 0.3, AlidNdEtaCorrection::kINEL);
 144
 145   TH1F* hRatios[6];
 146   for (Int_t i=0; i<6; i++) {
 147     // temporary correction with changed particle composistion
 148     AlidNdEtaCorrection* correctionTmp = (AlidNdEtaCorrection*)correctionAll->Clone();
 149
 150     correctionTmp->Reset();
 151
 152     correctionTmp->Add(correctionsSyst[0],scalesPi[i]);
 153     correctionTmp->Add(correctionsSyst[1],scalesK[i]);
 154     correctionTmp->Add(correctionsSyst[2],scalesP[i]);
 155     correctionTmp->Add(correctionsSyst[3],1);
 156
 157     correctionTmp->Finish();
 158
 159     dNdEtaAnalysis* analysisTmp = (dNdEtaAnalysis*)analysis->Clone();
 160     analysisTmp->Finish(correctionTmp, 0.3, AlidNdEtaCorrection::kINEL);
 161
 162     hRatios[i] = (TH1F*)analysisTmp->GetdNdEtaHistogram(2)->Clone();
 163     hRatios[i]->Divide((TH1F*)analysis->GetdNdEtaHistogram(2)->Clone(),hRatios[i],1,1,"B");
 164
 165     TString name(Form("ratio_%d",i));
 166     hRatios[i]->SetName(name.Data());
 167     name = TString(Form("#pi #times %0.1f, K #times %0.1f, p #times %0.1f",scalesPi[i],scalesK[i],scalesP[i]));
 168     hRatios[i]->SetTitle(name.Data());
 169     hRatios[i]->SetYTitle("ratio (standard/changed compositions)");
 170   }
 171
 172   TFile* fout = new TFile("systematics_composition.root","RECREATE");
 173
 174   for (Int_t i=0; i<6; i++)
 175     hRatios[i]->Write();
 176
 177   fout->Write();
 178   fout->Close();
 179
 180
 181 }