PWGLF/totEt/macros/emEt/PlotNeutronDeposits.C

   1 TH1* bayneseffdiv(TH1* numerator, TH1* denominator,Char_t* name);
   2 TH2* bayneseffdiv2D(TH2* numerator, TH2* denominator,Char_t* name) ;
   3 TH2F* flip(TH2F* orig,Char_t* name){
   4   float nbinx  = orig->GetXaxis()->GetNbins();
   5   float xlow = orig->GetXaxis()->GetBinLowEdge(1);
   6   float xhigh = orig->GetXaxis()->GetBinLowEdge(nbinx+1);
   7   float nbiny  = orig->GetYaxis()->GetNbins();
   8   float ylow= orig->GetYaxis()->GetBinLowEdge(1);
   9   float yhigh = orig->GetYaxis()->GetBinLowEdge(nbiny+1);
  10   TH2F *output = new TH2F(name,orig->GetTitle(),nbiny,ylow,yhigh,nbinx,xlow,xhigh);
  11   for(int i = 1;i<=nbinx;i++){
  12     for(int j = 1;j<=nbiny;j++){
  13       output->SetBinContent(j,i,orig->GetBinContent(i,j));
  14     }
  15   }
  16   return output;
  17 }
  18 void SetStyles(TH1 *histo,int marker, int color){
  19   histo->Sumw2();
  20   histo->SetMarkerStyle(marker);
  21   histo->SetMarkerColor(color);
  22   histo->SetLineColor(color);
  23   //histo->GetXaxis()->SetTitle(xtitle);
  24   //histo->GetYaxis()->SetTitle(ytitle);
  25 }
  26 Int_t colors[] = {0,TColor::kRed, TColor::kOrange, TColor::kGreen+3, TColor::kBlue, TColor::kBlack,
  27                     TColor::kRed, TColor::kOrange, TColor::kGreen+3, TColor::kBlue, TColor::kBlack,
  28                     TColor::kRed, TColor::kOrange, TColor::kGreen+3, TColor::kBlue, TColor::kBlack,
  29                     TColor::kRed, TColor::kOrange, TColor::kGreen+3, TColor::kBlue, TColor::kBlack};
  30 Int_t markers[] = {20,21,22,23,33, 24,25,26,32,27, 20,21,22,23,33, 24,25,26,32,27};
  31
  32 void PlotNeutronDeposits(TString filename = "rootFiles/LHC11a10a_bis/Et.ESD.simPbPb.EMCal.LHC11a10a_bis.Run139465.root"){
  33   gROOT->LoadMacro("macros/loadLibraries.C");
  34   loadLibraries();
  35    TFile *f = TFile::Open(filename, "READ");
  36     TList *l = dynamic_cast<TList*>(f->Get("out1"));
  37     TString histoname = "fHistNeutronsEtVsCent";
  38     TH2F *histo = l->FindObject(histoname.Data());
  39     histoname = "fHistNeutronsNumVsCent";
  40     TH2F *histoNum = l->FindObject(histoname.Data());
  41     TString outfilename = "/tmp/";
  42     TString calocorrfilename = "calocorrections.";
  43     TString det = "";
  44     if(filename.Contains("EMC")){
  45       outfilename +="EMCAL";
  46       calocorrfilename +="EMCAL";
  47       det += "Emcal";
  48     }
  49     else{
  50       outfilename +="PHOS";
  51       calocorrfilename +="PHOS";
  52       det += "Phos";
  53     }
  54     calocorrfilename +=".root";
  55     TString textfilename = "Neutrons"+det+".dat";
  56   TH2F  *fHistGammasGeneratedMult = l->FindObject("fHistGammasGeneratedCent");
  57   TH2F  *fHistGammasFoundMult = l->FindObject("fHistGammasFoundCent");
  58   TH2F  *fHistGammasFoundOutOfAccCent = l->FindObject("fHistGammasFoundOutOfAccCent");
  59   fHistGammasFoundMult->Add(fHistGammasFoundOutOfAccCent);
  60   // TH2F *gammaEff2D = (TH2F*) bayneseffdiv2D(fHistGammasFoundMult,fHistGammasGeneratedMult,"gammaEff2D");
  61
  62     TCanvas *c1 = new TCanvas("c1","<E_{T}>",600,400);
  63     c1->SetTopMargin(0.02);
  64     c1->SetRightMargin(0.02);
  65     c1->SetBorderSize(0);
  66     c1->SetFillColor(0);
  67     c1->SetFillColor(0);
  68     c1->SetBorderMode(0);
  69     c1->SetFrameFillColor(0);
  70     c1->SetFrameBorderMode(0);
  71     TH1D *prof = histo->ProfileY();
  72     histo->Draw("colz");
  73     //prof->Draw("same");
  74
  75     TCanvas *c3 = new TCanvas("c3","<E_{T}>",600,400);
  76     c3->SetTopMargin(0.02);
  77     c3->SetRightMargin(0.02);
  78     c3->SetBorderSize(0);
  79     c3->SetFillColor(0);
  80     c3->SetFillColor(0);
  81     c3->SetBorderMode(0);
  82     c3->SetFrameFillColor(0);
  83     c3->SetFrameBorderMode(0);
  84     //for(int i=1;i<=18;i++){
  85
  86     //}
  87     prof->Draw();
  88     prof->GetYaxis()->SetTitle("<E_{T}>");
  89     prof->GetXaxis()->SetTitle("centrality bin");
  90     TString outfilename1 = outfilename+"Et.png";
  91     c3->SaveAs(outfilename1.Data());
  92
  93     TCanvas *c2 = new TCanvas("c2","<Num>",600,400);
  94     c2->SetTopMargin(0.02);
  95     c2->SetRightMargin(0.02);
  96     c2->SetBorderSize(0);
  97     c2->SetFillColor(0);
  98     c2->SetFillColor(0);
  99     c2->SetBorderMode(0);
 100     c2->SetFrameFillColor(0);
 101     c2->SetFrameBorderMode(0);
 102     TH1D *profNum = histoNum->ProfileY();
 103     histoNum->Draw("colz");
 104     //profNum->Draw("same");
 105
 106     TCanvas *c4 = new TCanvas("c4","<Num>",600,400);
 107     c4->SetTopMargin(0.02);
 108     c4->SetRightMargin(0.02);
 109     c4->SetBorderSize(0);
 110     c4->SetFillColor(0);
 111     c4->SetFillColor(0);
 112     c4->SetBorderMode(0);
 113     c4->SetFrameFillColor(0);
 114     c4->SetFrameBorderMode(0);
 115     profNum->Draw("same");
 116     profNum->GetYaxis()->SetTitle("<N_{neutrons}>");
 117     profNum->GetXaxis()->SetTitle("centrality bin");
 118     TString outfilename2 = outfilename+"Num.png";
 119     c4->SaveAs(outfilename2.Data());
 120
 121     TFile *calocorrfile = TFile::Open(calocorrfilename, "READ");
 122     TString recoEffName = "ReCorrections";
 123     recoEffName += det;
 124
 125     AliAnalysisEtRecEffCorrection *recoEff = (AliAnalysisEtRecEffCorrection *) calocorrfile->Get(recoEffName.Data());
 126
 127     cout<<"I am here"<<endl;
 128     ofstream myfile;
 129     myfile.open (textfilename.Data());
 130     for(int i=0;i<20;i++){
 131       Float_t et = prof->GetBinContent(i+1);
 132       Float_t num = profNum->GetBinContent(i+1);
 133       Float_t eff = recoEff->ReconstructionEfficiency(et,i);
 134       float neutroncorr = 0;
 135       if(eff>0) neutroncorr = et*num/eff;
 136       float error = 0;
 137       cout<<"et*num/eff = "<<et<<"*"<<num<<"/"<<eff<<" "<< neutroncorr<<endl;
 138       myfile<<neutroncorr<<" "<<error<<endl;
 139     }
 140     myfile.close();
 141 //     return;
 142 //     TCanvas *c5 = new TCanvas("c5","<Num>",600,400);
 143 //     c5->SetTopMargin(0.02);
 144 //     c5->SetRightMargin(0.02);
 145 //     c5->SetBorderSize(0);
 146 //     c5->SetFillColor(0);
 147 //     c5->SetFillColor(0);
 148 //     c5->SetBorderMode(0);
 149 //     c5->SetFrameFillColor(0);
 150 //     c5->SetFrameBorderMode(0);
 151
 152
 153     TCanvas *c6 = new TCanvas("c6","<Num>",600,400);
 154     c6->SetTopMargin(0.02);
 155     c6->SetRightMargin(0.02);
 156     c6->SetBorderSize(0);
 157     c6->SetFillColor(0);
 158     c6->SetFillColor(0);
 159     c6->SetBorderMode(0);
 160     c6->SetFrameFillColor(0);
 161     c6->SetFrameBorderMode(0);
 162     c6->SetLogz();
 163     TH2F *histoNumAlt = flip(histoNum,"histoNumAlt");
 164     TH1D *profNumAlt = histoNum->ProfileY();
 165     histoNumAlt->Draw("colz");
 166     profNum->Draw("same");
 167
 168 }
 169
 170 TH1* bayneseffdiv(TH1* numerator, TH1* denominator,Char_t* name)
 171 {
 172     if(!numerator){
 173       cerr<<"Error:  numerator does not exist!"<<endl;
 174       return NULL;
 175     }
 176     if(!denominator){
 177       cerr<<"Error:  denominator does not exist!"<<endl;
 178       return NULL;
 179     }
 180     TH1* result = (TH1*)numerator->Clone(name);
 181     Int_t nbins = numerator->GetNbinsX();
 182     for (Int_t ibin=0; ibin<= nbins+1; ++ibin) {
 183       Double_t numeratorVal = numerator->GetBinContent(ibin);
 184       Double_t denominatorVal = denominator->GetBinContent(ibin);
 185       // Check if the errors are right or the thing is scaled
 186       Double_t numeratorValErr = numerator->GetBinError(ibin);
 187       if (!(numeratorValErr==0. || numeratorVal ==0.) ) {
 188         Double_t rescale = numeratorValErr*numeratorValErr/numeratorVal;
 189         numeratorVal /= rescale;
 190       }
 191       Double_t denominatorValErr = denominator->GetBinError(ibin);
 192       if (!(denominatorValErr==0. || denominatorVal==0. )) {
 193         Double_t rescale = denominatorValErr*denominatorValErr/denominatorVal;
 194         denominatorVal /= rescale;
 195       }
 196       Double_t quotient = 0.;
 197       if (denominatorVal!=0.) {
 198         quotient = numeratorVal/denominatorVal;
 199       }
 200       Double_t quotientErr=0;
 201       quotientErr = TMath::Sqrt(
 202                                 (numeratorVal+1.0)/(denominatorVal+2.0)*
 203                                 ((numeratorVal+2.0)/(denominatorVal+3.0)-(numeratorVal+1.0)/(denominatorVal+2.0)));
 204       result->SetBinContent(ibin,quotient);
 205       result->SetBinError(ibin,quotientErr);
 206       //cout<<"Setting bin "<<ibin<<" to "<<quotient<<" "<<numeratorVal<<"/"<<denominatorVal<<endl;
 207     }
 208     return result;
 209 }
 210
 211
 212 TH2* bayneseffdiv2D(TH2* numerator, TH2* denominator,Char_t* name)
 213 {
 214   if(!numerator){
 215     cerr<<"Error:  numerator does not exist!"<<endl;
 216     return NULL;
 217   }
 218   if(!denominator){
 219     cerr<<"Error:  denominator does not exist!"<<endl;
 220     return NULL;
 221   }
 222   TH2* result = (TH2*)numerator->Clone(name);
 223   Int_t nbinsX = numerator->GetNbinsX();
 224   Int_t nbinsY = numerator->GetNbinsY();
 225   for (Int_t ibin=0; ibin<= nbinsX+1; ++ibin) {
 226     for (Int_t jbin=0; jbin<= nbinsY+1; ++jbin) {
 227       Double_t numeratorVal = numerator->GetBinContent(ibin,jbin);
 228       Double_t denominatorVal = denominator->GetBinContent(ibin,jbin);
 229       // Check if the errors are right or the thing is scaled
 230       Double_t numeratorValErr = numerator->GetBinError(ibin,jbin);
 231       if (!(numeratorValErr==0. || numeratorVal ==0.) ) {
 232         Double_t rescale = numeratorValErr*numeratorValErr/numeratorVal;
 233         numeratorVal /= rescale;
 234       }
 235       Double_t denominatorValErr = denominator->GetBinError(ibin,jbin);
 236       if (!(denominatorValErr==0. || denominatorVal==0. )) {
 237         Double_t rescale = denominatorValErr*denominatorValErr/denominatorVal;
 238         denominatorVal /= rescale;
 239       }
 240       Double_t quotient = 0.;
 241       if (denominatorVal!=0.) {
 242         quotient = numeratorVal/denominatorVal;
 243       }
 244       Double_t quotientErr=0;
 245       quotientErr = TMath::Sqrt(
 246                                 (numeratorVal+1.0)/(denominatorVal+2.0)*
 247                                 ((numeratorVal+2.0)/(denominatorVal+3.0)-(numeratorVal+1.0)/(denominatorVal+2.0)));
 248       result->SetBinContent(ibin,jbin,quotient);
 249       result->SetBinError(ibin,jbin,quotientErr);
 250       //cout<<"Setting bin "<<ibin<<" to "<<quotient<<" "<<numeratorVal<<"/"<<denominatorVal<<endl;
 251     }
 252   }
 253   return result;
 254 }