PWG3/muondep/MuonTrackingEfficiency.C

   1 #include <fstream>
   2 #include <stdio.h>
   3 #include <stdlib.h>
   4 #include <TH1I.h>
   5 #include <TH2I.h>
   6 #include <TH1D.h>
   7 #include <TH2D.h>
   8 #include <TAxis.h>
   9 #include <TLegend.h>
  10 #include <TString.h>
  11 #include <Riostream.h>
  12 #include <TFile.h>
  13 #include <TList.h>
  14 #include <TCanvas.h>
  15 #include <TSystem.h>
  16 #include <TStyle.h>
  17 #include <TGraphAsymmErrors.h>
  18 #include <vector>
  19 #include <TGaxis.h>
  20
  21 /*
  22
  23 This macro is meant to be used with the output of the Efficiency Task. It computes the uncorrelated efficiencies per chammber, station and globally. It then looks for correlated dead arrays and corrects the previous calculated efficiency.
  24
  25 Macro created by Lizardo Valencia Palomo
  26
  27 lizardo.valencia.palomo@cern.ch
  28
  29 */
  30
  31
  32
  33 TH1D* GraphToHist(Bool_t ToBePlotted, const TGraphAsymmErrors* Graph, Int_t Chamber);
  34 std::vector<Double_t> FillVector(TH1D* Chamber, Bool_t Errors);
  35 std::vector<Double_t> CorrelatedDeadArray(std::vector<Double_t> Chj, std::vector<Double_t> Chi, Int_t NumberOfDEs, Int_t SymmetricGroups);
  36 void Drawing(std::vector<Double_t> Chi, std::vector<Double_t> Chj, std::vector<Double_t> ChiErr, std::vector<Double_t> ChjErr, Int_t Bins, Int_t SymGroups, Int_t Station);
  37 Double_t N00(Double_t Average, Double_t CDA, Int_t SimGroup, Int_t SimGroupsPerCh);
  38
  39
  40 void MuonTrackingEfficiency(TString fileName = "./AnalysisResults.root")
  41 {
  42
  43   TFile *file = new TFile(fileName, "read");
  44   TDirectoryFile *Directory = (TDirectoryFile*)((file->GetDirectory("MUON_Efficiency")));
  45
  46   std::vector<TH1D*> CorrectedHistos;
  47   std::vector<TH1D*> CorrectedChiHisto;
  48   std::vector<TH1D*> CorrectedChjHisto;
  49
  50   TList *listTT = (TList *)Directory->Get("TotalTracksPerChamber;1");
  51   TList *listTD = (TList *)Directory->Get("TracksDetectedPerChamber;1");
  52   TList *listSD = (TList *)Directory->Get("SingleDetectedPerChamber;1");
  53
  54   TH1D *EfficiencyHisto = 0x0;
  55   TH1D *TD = 0x0;
  56   TH1D *TT = 0x0;
  57
  58   for (Int_t Chamber = 1; Chamber < 12; Chamber++) //Loop over all chambers + 1
  59     {
  60
  61       TH2F *TD2Dim = (TH2F *) listTD->At(Chamber - 1);
  62       TH2F *TT2Dim = (TH2F *) listTT->At(Chamber - 1);
  63       TH2F *SD2Dim = (TH2F *) listSD->At(Chamber - 1);
  64
  65       Int_t nDEch = 0;
  66       nDEch = TD2Dim->GetNbinsX();
  67
  68       TH1D *Ncluster = 0x0;
  69       TH1D *Corrected = 0x0;
  70       TH1D *SD = 0x0;
  71
  72       TH1D *EfficiencyHistoCh11 = 0x0;
  73       TH1D *TDch11 = 0x0;
  74       TH1D *TTch11 = 0x0;
  75
  76       if (Chamber == 11)
  77       {
  78
  79         TDch11 = new TH1D(Form("TDch11 %d",Chamber), "", nDEch, - 0.5, nDEch - 0.5);
  80         TTch11 = new TH1D(Form("TTch11 %d",Chamber), "", nDEch, - 0.5, nDEch - 0.5);
  81
  82       }
  83
  84       TD = new TH1D(Form("TD %d",Chamber), "", nDEch, - 0.5, nDEch - 0.5);
  85       TT = new TH1D(Form("TT %d",Chamber), "", nDEch, - 0.5, nDEch - 0.5);
  86       SD = new TH1D(Form("SD %d",Chamber), "", nDEch, - 0.5, nDEch - 0.5);
  87
  88       std::vector<Double_t> TDentries (nDEch, 0.0);
  89       std::vector<Double_t> TTentries (nDEch, 0.0);
  90       std::vector<Double_t> SDentries (nDEch, 0.0);
  91
  92       for (Int_t ix = 1; ix <= nDEch; ix++)
  93         {
  94           for (Int_t iy = 1; iy <= TD2Dim->GetNbinsY(); iy++)
  95             {
  96
  97               TDentries[ix-1] += TD2Dim->GetBinContent(ix,iy);
  98               TTentries[ix-1] += TT2Dim->GetBinContent(ix,iy);
  99               SDentries[ix-1] += SD2Dim->GetBinContent(ix,iy);
 100
 101             }
 102
 103           TD->SetBinContent(ix,TDentries[ix-1]);
 104           TD->SetBinError(ix,sqrt(TDentries[ix-1]));
 105           TT->SetBinContent(ix,TTentries[ix-1]);
 106           TT->SetBinError(ix,sqrt(TTentries[ix-1]));
 107           SD->SetBinContent(ix,SDentries[ix-1]);
 108           SD->SetBinError(ix,sqrt(SDentries[ix-1]));
 109
 110           if (Chamber != 11) continue;
 111
 112           TDch11->SetBinContent(ix,TDentries[ix-1]);
 113           TDch11->SetBinError(ix,sqrt(TDentries[ix-1]));
 114           TTch11->SetBinContent(ix,TTentries[ix-1]);
 115           TTch11->SetBinError(ix,sqrt(TTentries[ix-1]));
 116
 117         }
 118
 119       TGraphAsymmErrors *EfficiencyGraph = new TGraphAsymmErrors(TD, TT);
 120
 121       EfficiencyHisto = GraphToHist(kFALSE,EfficiencyGraph,Chamber);
 122
 123       Ncluster = new TH1D(Form("Ncluster %d",Chamber), "N_{Cluster}", nDEch, - 0.5, nDEch - 0.5);
 124       Ncluster->Sumw2();
 125       Ncluster->Add(TD,SD,1,1);
 126
 127       Corrected = new TH1D(Form("Corrected %d",Chamber), Form("Corrected N_{Cluster} Ch %d",Chamber), nDEch, - 0.5, nDEch - 0.5);
 128       Corrected->Sumw2();
 129       Corrected->Divide(Ncluster,EfficiencyHisto,1,1);
 130
 131       CorrectedHistos.push_back(Corrected);
 132
 133       if (Chamber != 11) continue;
 134
 135       TGraphAsymmErrors *EfficiencyGraphCh11 = new TGraphAsymmErrors(TDch11, TTch11);
 136
 137       EfficiencyHistoCh11 = GraphToHist(kTRUE,EfficiencyGraphCh11,Chamber);
 138
 139       TCanvas *CanvasEffPerCh = new TCanvas("Efficiency per chambers", "Efficiency per chambers");
 140       gStyle->SetOptStat(0);
 141       CanvasEffPerCh->SetFillColor(0);
 142       CanvasEffPerCh->SetBorderSize(1);
 143       CanvasEffPerCh->SetFrameLineWidth(2);
 144
 145       EfficiencyHistoCh11->SetTitle(kFALSE);
 146       EfficiencyHistoCh11->SetMarkerStyle(20);
 147       EfficiencyHistoCh11->SetMarkerColor(1);
 148       EfficiencyHistoCh11->SetMarkerSize(1.3);
 149       EfficiencyHistoCh11->GetXaxis()->SetTitle("Chamber");
 150       EfficiencyHistoCh11->GetXaxis()->SetNdivisions(Chamber);
 151       EfficiencyHistoCh11->GetXaxis()->SetTitleSize(0.05);
 152       EfficiencyHistoCh11->GetXaxis()->SetTitleOffset(0.91);
 153       EfficiencyHistoCh11->GetXaxis()->SetLabelSize(0.05);
 154       EfficiencyHistoCh11->GetYaxis()->SetTitle("Efficiency");
 155       EfficiencyHistoCh11->GetYaxis()->CenterTitle(kTRUE);
 156       EfficiencyHistoCh11->GetYaxis()->SetTitleSize(0.06);
 157       EfficiencyHistoCh11->GetYaxis()->SetTitleOffset(0.79);
 158       EfficiencyHistoCh11->GetYaxis()->SetLabelSize(0.05);
 159       EfficiencyHistoCh11->Draw("e");
 160
 161     }
 162
 163   /*************************Correlated Dead Arrays check****************************/
 164
 165
 166   Double_t MissingTracks [5] = {0.0};//Station one is in slot cero!!!
 167
 168   for (Int_t CDAperSt = 0; CDAperSt < 5; CDAperSt++)
 169     {
 170
 171       std::vector<Double_t> FilledVectorChi = FillVector(CorrectedHistos[CDAperSt*2], kFALSE);
 172       std::vector<Double_t> FilledVectorChj = FillVector(CorrectedHistos[(CDAperSt*2)+1], kFALSE);
 173       std::vector<Double_t> FilledVectorChiErr = FillVector(CorrectedHistos[CDAperSt*2], kTRUE);
 174       std::vector<Double_t> FilledVectorChjErr = FillVector(CorrectedHistos[(CDAperSt*2)+1], kTRUE);
 175
 176       Int_t DEs = FilledVectorChi.size();
 177       Int_t SimGroups = 0;
 178       Double_t ToFillMissingTracks = 0.0;
 179
 180       if ( CDAperSt < 2 ) SimGroups = static_cast<Int_t>(DEs/2);
 181       else SimGroups = static_cast<Int_t>(DEs/2) + 1;
 182
 183       std::vector<Double_t> CDAs = CorrelatedDeadArray(FilledVectorChi, FilledVectorChj, DEs, SimGroups);
 184       Drawing(FilledVectorChi, FilledVectorChj, FilledVectorChiErr, FilledVectorChjErr, DEs, SimGroups, CDAperSt);
 185
 186       for (Int_t DEst = 0; DEst < DEs; DEst++)
 187         {
 188           if ( CDAs[DEst] != 0.0 )
 189             {
 190
 191               printf ("Looks like there is a CDA in DE %d of Station %d: N00 = %.2f \n", DEst, CDAperSt+1, CDAs[DEst]);
 192               ToFillMissingTracks += CDAs[DEst];
 193
 194             }
 195         }
 196
 197       MissingTracks[CDAperSt] = ToFillMissingTracks;
 198
 199     }
 200
 201
 202
 203   /**********************************Efficiency from Efficiency Task************************************/
 204
 205
 206       vector<Double_t> chambersEff;
 207       vector<Double_t> chambersEffErr;
 208       vector<Double_t> StEfficiency;
 209       vector<Double_t> StEfficiencyErr;
 210
 211       chambersEff.push_back(0.0);
 212       chambersEffErr.push_back(0.0);
 213       StEfficiency.push_back(0.0);
 214       StEfficiencyErr.push_back(0.0);
 215
 216       for (Int_t EfCh = 0; EfCh < 10; EfCh++)
 217         {
 218
 219           chambersEff.push_back(EfficiencyHisto->GetBinContent(EfCh+1));
 220           chambersEffErr.push_back(EfficiencyHisto->GetBinError(EfCh+1));
 221
 222           cout << "Efficiency chamber " << EfCh+1 << " = " << chambersEff[EfCh+1] << " +- " << chambersEffErr[EfCh+1] << endl;
 223
 224         }
 225
 226       for (Int_t EffSt = 0; EffSt < 4; EffSt++)
 227         {
 228
 229           if (EffSt < 3)
 230             {
 231
 232               StEfficiency.push_back(1.0 - (1.0-chambersEff[(2*EffSt)+1])*(1.0-chambersEff[(2*EffSt)+2]));
 233               StEfficiencyErr.push_back(StEfficiency[EffSt+1]*(chambersEffErr[(2*EffSt)+1]/chambersEff[(2*EffSt)+1] + chambersEffErr[(2*EffSt)+2]/chambersEff[(2*EffSt)+2]));
 234
 235           cout << "Efficiency station " << EffSt+1 <<" = " <<StEfficiency[EffSt+1] << " +- " <<StEfficiencyErr[EffSt+1]<<endl;
 236
 237             }
 238
 239           else
 240             {
 241
 242               StEfficiency.push_back(chambersEff[(2*EffSt)+1] * chambersEff[(2*EffSt)+2] * chambersEff[(2*EffSt)+3] * chambersEff[(2*EffSt)+4]
 243               + (1.0 - chambersEff[(2*EffSt)+1]) * chambersEff[(2*EffSt)+2] * chambersEff[(2*EffSt)+3] * chambersEff[(2*EffSt)+4]
 244               + chambersEff[(2*EffSt)+1] * (1.0 - chambersEff[(2*EffSt)+2]) * chambersEff[(2*EffSt)+3] * chambersEff[(2*EffSt)+4]
 245               + chambersEff[(2*EffSt)+1] * chambersEff[(2*EffSt)+2] * (1.0 - chambersEff[(2*EffSt)+3]) * chambersEff[(2*EffSt)+4]
 246               + chambersEff[(2*EffSt)+1] * chambersEff[(2*EffSt)+2] * chambersEff[(2*EffSt)+3] * (1.0 - chambersEff[(2*EffSt)+4]));
 247
 248               StEfficiencyErr.push_back(StEfficiency[EffSt+1]*(chambersEffErr[(2*EffSt)+1]/chambersEff[(2*EffSt)+1]
 249               + chambersEffErr[(2*EffSt)+2]/chambersEff[(2*EffSt)+2] + chambersEffErr[(2*EffSt)+3]/chambersEff[(2*EffSt)+3]
 250               + chambersEffErr[(2*EffSt)+4]/chambersEff[(2*EffSt)+4]));
 251
 252               cout << "Efficiency station 45"<<" = " <<StEfficiency[EffSt+1] << " +- " <<StEfficiencyErr[EffSt+1]<<endl;
 253               //Notice efficiency from st 45 is in slot 4 of the vector!!!
 254
 255             }
 256
 257        }
 258
 259       StEfficiency.push_back(StEfficiency[1]*StEfficiency[2]*StEfficiency[3]*StEfficiency[4]);
 260       StEfficiencyErr.push_back(StEfficiencyErr[1]+StEfficiencyErr[2]+StEfficiencyErr[3]+StEfficiencyErr[4]);
 261
 262       cout << "Total efficiency = " <<StEfficiency[5]<< " +- " <<StEfficiencyErr[5]<< endl;
 263
 264
 265
 266  /*************************************Correcting efficiencies********************************************/
 267
 268
 269   TH1D *TTcorrected = new TH1D("TTcorrected", "TDcorrected", 10, -0.5, 9.5);
 270
 271   for (Int_t St = 0; St < 5; St++)
 272     {
 273
 274       Double_t NewBinContentChi = TT->GetBinContent((2*St)+1) + MissingTracks[St];
 275       TTcorrected->SetBinContent((2*St)+1,NewBinContentChi);
 276       TTcorrected->SetBinContent((2*St)+2,TT->GetBinContent((2*St)+2));
 277
 278     }
 279
 280   TGraphAsymmErrors *EfficiencyGraphCorrected = new TGraphAsymmErrors(TD,TTcorrected);
 281
 282   TH1D *EfficiencyHistoCorrected = GraphToHist(1,EfficiencyGraphCorrected,11);//In this case doesnt matter which value is used in the first slot.
 283
 284   vector<Double_t> chambersEffCorrected;
 285   vector<Double_t> chambersEffErrCorrected;
 286
 287   chambersEffCorrected.push_back(0.0);
 288   chambersEffErrCorrected.push_back(0.0);
 289
 290   for (Int_t EfChCorr = 0; EfChCorr < 10; EfChCorr++)
 291     {
 292
 293       chambersEffCorrected.push_back(EfficiencyHistoCorrected->GetBinContent(EfChCorr+1));
 294       chambersEffErrCorrected.push_back(EfficiencyHistoCorrected->GetBinError(EfChCorr+1));
 295
 296       if ( (chambersEffCorrected[EfChCorr+1]/chambersEff[EfChCorr+1]) == 1 ) continue;
 297
 298       cout << "Corrected efficiency chamber " << EfChCorr+1 << " = " << chambersEffCorrected[EfChCorr+1] << " +- " << chambersEffErrCorrected[EfChCorr+1] << endl;
 299
 300     }
 301
 302   vector<Double_t> StEfficiencyCorrected;
 303   vector<Double_t> StEfficiencyErrCorrected;
 304
 305   StEfficiencyCorrected.push_back(0.0);
 306   StEfficiencyErrCorrected.push_back(0.0);
 307
 308   //By default CDA in Station 4 is not corrected.
 309
 310   for (Int_t Station = 0; Station < 4; Station++)
 311     {
 312
 313       if (MissingTracks[Station] == 0.0)
 314         {
 315
 316           StEfficiencyCorrected.push_back(StEfficiency[Station+1]);
 317           StEfficiencyErrCorrected.push_back(StEfficiencyErr[Station+1]);
 318
 319         }
 320
 321       else
 322         {
 323
 324           if (Station < 3)
 325             {
 326
 327           StEfficiencyCorrected.push_back(StEfficiency[Station+1]*(chambersEffCorrected[(2*Station)+1]/chambersEff[(2*Station)+1]));
 328           StEfficiencyErrCorrected.push_back(StEfficiencyCorrected[Station+1]*sqrt((StEfficiencyErr[Station+1]/StEfficiency[Station+1])*(StEfficiencyErr[Station+1]/StEfficiency[Station+1])+(chambersEffErrCorrected[(2*Station)+1]/chambersEffCorrected[(2*Station)+1])*(chambersEffErrCorrected[(2*Station)+1]/chambersEffCorrected[(2*Station)+1]) + (chambersEffErr[(2*Station)+1]/chambersEff[(2*Station)+1])*(chambersEffErr[(2*Station)+1]/chambersEff[(2*Station)+1])));
 329
 330            cout<<"Corrected efficiency station "<<Station+1<<" = "<<StEfficiencyCorrected[Station+1]<<" +- "<<StEfficiencyErrCorrected[Station+1]<<endl;
 331
 332            }
 333
 334          if (Station == 3)
 335            {
 336
 337              StEfficiencyCorrected.push_back(StEfficiency[Station+1]*(chambersEffCorrected[(2*Station)+3]/chambersEff[(2*Station)+3]));
 338              StEfficiencyErrCorrected.push_back(StEfficiencyCorrected[Station+1]*sqrt((StEfficiencyErr[Station+1]/StEfficiency[Station+1])*(StEfficiencyErr[Station+1]/StEfficiency[Station+1])+(chambersEffErrCorrected[(2*Station)+3]/chambersEffCorrected[(2*Station)+3])*(chambersEffErrCorrected[(2*Station)+3]/chambersEffCorrected[(2*Station)+3]) + (chambersEffErr[(2*Station)+3]/chambersEff[(2*Station)+3])*(chambersEffErr[(2*Station)+3]/chambersEff[(2*Station)+3])));
 339
 340              cout<<"Corrected efficiency station 45"<<" = "<<StEfficiencyCorrected[Station+1]<<" +- "<<StEfficiencyErrCorrected[Station+1]<<endl;
 341
 342            }
 343
 344
 345         }
 346
 347     }
 348
 349   StEfficiencyCorrected.push_back(StEfficiencyCorrected[1]*StEfficiencyCorrected[2]*StEfficiencyCorrected[3]*StEfficiencyCorrected[4]);
 350   StEfficiencyErrCorrected.push_back(StEfficiencyErrCorrected[1]+StEfficiencyErrCorrected[2]+StEfficiencyErrCorrected[3]+StEfficiencyErrCorrected[4]);
 351
 352   cout<<"Corrected total efficiency"<<" = "<<StEfficiencyCorrected[5]<< " +- " <<StEfficiencyErrCorrected[5]<< endl;
 353
 354 }
 355
 356
 357 /*************************************Functions************************************/
 358
 359
 360 std::vector<Double_t> FillVector(TH1D* Chamber, Bool_t Errors)
 361 {
 362
 363   Int_t NumOfDEs = Chamber->GetNbinsX();
 364   std::vector<Double_t> ChVector (NumOfDEs, 0);
 365   for (Int_t i = 1; i <= NumOfDEs; i++)
 366         {
 367           if (Errors == kFALSE) ChVector.at(i-1) = Chamber->GetBinContent(i);
 368           if (Errors == kTRUE) ChVector.at(i-1) = Chamber->GetBinError(i);
 369         }
 370
 371   return ChVector;
 372
 373 }
 374
 375 std::vector<Double_t> CorrelatedDeadArray(std::vector<Double_t> Chi, std::vector<Double_t> Chj, Int_t NumberOfDEs, Int_t SymmetricGroups)
 376 {
 377
 378   std::vector<Double_t> Average (SymmetricGroups, 0.0);
 379   std::vector<Double_t> CDA (NumberOfDEs, 0.0);
 380
 381   for (Int_t Sim = 0; Sim < SymmetricGroups/2; Sim++)
 382     {
 383        if (Sim == 0)
 384         {
 385
 386           if ( (Sim == 0) && (SymmetricGroups/2 == 1) )
 387             {
 388
 389               Average[Sim] = (Chi[Sim]+Chi[Sim+1]+Chi[Sim+2]+Chi[Sim+3])/4;
 390               Average[Sim+(SymmetricGroups/2)] = (Chj[Sim]+Chj[Sim+1]+Chj[Sim+2]+Chj[Sim+3])/4;
 391
 392               for (Int_t i = 1; i < 5; i++)
 393                  {
 394
 395                   if ( (Chi[i-1] < (Average[Sim]-(Average[Sim]*0.2))) && (Chj[i-1] < (Average[Sim+(SymmetricGroups/2)]-(Average[Sim+(SymmetricGroups/2)]*0.2))) )
 396                     {
 397                       CDA[i-1] = N00(Average[Sim], Chi[i-1], Sim, SymmetricGroups/2);
 398                     }
 399
 400                  }
 401
 402             }
 403
 404           else
 405             {
 406
 407           Average[Sim] = (Chi[Sim] + Chi[Sim+(NumberOfDEs/2)])/2;
 408           Average[Sim+(SymmetricGroups/2)] = (Chj[Sim] + Chj[Sim+(NumberOfDEs/2)])/2;
 409
 410
 411          if ( (Chi[Sim] < (Average[Sim]-(Average[Sim]*0.2))) && (Chj[Sim] < (Average[Sim+(SymmetricGroups/2)]-(Average[Sim+(SymmetricGroups/2)]*0.2))) )
 412             {
 413               CDA[Sim] = N00(Average[Sim], Chi[Sim], Sim, SymmetricGroups/2);
 414             }
 415          if ( (Chi[Sim+(NumberOfDEs/2)] < (Average[Sim]-(Average[Sim]*0.2))) && (Chj[Sim+(NumberOfDEs/2)] < (Average[Sim+(SymmetricGroups/2)]-(Average[Sim+(SymmetricGroups/2)]*0.2))) )
 416             {
 417               CDA[Sim+(NumberOfDEs/2)] = N00(Average[Sim], Chi[Sim+(NumberOfDEs/2)], Sim, SymmetricGroups/2);
 418             }
 419
 420            }
 421
 422          }
 423
 424       else
 425         {
 426
 427           Average[Sim] = (Chi[Sim] + Chi[(NumberOfDEs/2)-Sim] + Chi[Sim+(NumberOfDEs/2)] + Chi[NumberOfDEs-Sim])/4;
 428           Average[Sim+(SymmetricGroups/2)] = (Chj[Sim] + Chj[(NumberOfDEs/2)-Sim] + Chj[Sim+(NumberOfDEs/2)] + Chj[NumberOfDEs-Sim])/4;
 429
 430           if ( (Chi[Sim] < (Average[Sim]-(Average[Sim]*0.2))) && (Chj[Sim] < (Average[Sim+(SymmetricGroups/2)]-(Average[Sim+(SymmetricGroups/2)]*0.2))) )
 431              {
 432                CDA[Sim] = N00(Average[Sim], Chi[Sim], Sim, SymmetricGroups/2);
 433              }
 434            if ( (Chi[(NumberOfDEs/2)-Sim] < (Average[Sim]-(Average[Sim]*0.2))) && (Chj[(NumberOfDEs/2)-Sim] < (Average[Sim+(SymmetricGroups/2)]-(Average[Sim+(SymmetricGroups/2)]*0.2))) )
 435              {
 436                CDA[(NumberOfDEs/2)-Sim] = N00(Average[Sim], Chi[(NumberOfDEs/2)-Sim], Sim, SymmetricGroups/2);
 437              }
 438            if ( (Chi[Sim+(NumberOfDEs/2)] < (Average[Sim]-(Average[Sim]*0.2))) && (Chj[Sim+(NumberOfDEs/2)] < (Average[Sim+(SymmetricGroups/2)]-(Average[Sim+(SymmetricGroups/2)]*0.2))) )
 439              {
 440                CDA[Sim+(NumberOfDEs/2)] = N00(Average[Sim], Chi[Sim+(NumberOfDEs/2)], Sim, SymmetricGroups/2);
 441              }
 442
 443            if ( (Chi[NumberOfDEs-Sim] < (Average[Sim]-(Average[Sim]*0.2))) && (Chj[NumberOfDEs-Sim] < (Average[Sim+(SymmetricGroups/2)]-(Average[Sim+(SymmetricGroups/2)]*0.2))) )
 444              {
 445                CDA[NumberOfDEs-Sim] = N00(Average[Sim], Chi[NumberOfDEs-Sim], Sim, SymmetricGroups/2);
 446
 447              }
 448
 449         }
 450
 451     }
 452
 453   return CDA;
 454
 455 }
 456
 457
 458 Double_t N00(Double_t Average, Double_t CDA, Int_t SimGroup, Int_t SimGroupsPerCh)
 459 {
 460
 461   Double_t RecomputedAverage = 0.0;
 462   Double_t ToBeAdded = 0.0;
 463
 464   if ( (SimGroup == 0) && (SimGroupsPerCh != 1) ) RecomputedAverage = (2*Average) - CDA;
 465   else RecomputedAverage = ((4*Average) - CDA)/3;
 466
 467   ToBeAdded = RecomputedAverage - CDA;
 468
 469   return ToBeAdded;
 470
 471 }
 472
 473
 474 void Drawing(std::vector<Double_t> Chi, std::vector<Double_t> Chj, std::vector<Double_t> ChiErr, std::vector<Double_t> ChjErr, Int_t Bins, Int_t SymGroups, Int_t Station)
 475 {
 476
 477   TCanvas *Canvas = new TCanvas(Form("Station %d",Station+1), Form("Station %d",Station+1));
 478   gStyle->SetOptStat(0);
 479   Canvas->SetFillColor(0);
 480   if ( (Station+1) < 3) Canvas->Divide(SymGroups/2);
 481   if ( (Station+1) == 3) Canvas->Divide(Station+1,Station);
 482   if ( (Station+1) > 3) Canvas->Divide(3,3);
 483   Canvas->SetBorderSize(1);
 484   Canvas->SetFrameLineWidth(2);
 485
 486   for (Int_t Sim = 0; Sim < SymGroups/2; Sim++)
 487     {
 488
 489       TH1D *Chamberi = new TH1D(Form("Chi SimGroup %d St %d",Sim,Station+1), "", Bins, - 0.5, Bins - 0.5);
 490       TGaxis::SetMaxDigits(2);
 491       Chamberi->SetTitle(kFALSE);
 492       Chamberi->SetMarkerStyle(20);
 493       Chamberi->SetMarkerColor(1);
 494       Chamberi->SetMarkerSize(1.3);
 495       Chamberi->GetXaxis()->SetTitle("Detection Element");
 496       Chamberi->GetXaxis()->SetNdivisions(Bins+1);
 497       Chamberi->GetXaxis()->SetTitleSize(0.05);
 498       Chamberi->GetXaxis()->SetTitleOffset(0.91);
 499       Chamberi->GetXaxis()->SetLabelSize(0.05);
 500       Chamberi->GetYaxis()->SetTitle("Corrected N_{Cluster}");
 501       Chamberi->GetYaxis()->CenterTitle(kTRUE);
 502       Chamberi->GetYaxis()->SetTitleSize(0.06);
 503       Chamberi->GetYaxis()->SetTitleOffset(0.79);
 504       Chamberi->GetYaxis()->SetLabelSize(0.05);
 505       TH1D *Chamberj = new TH1D(Form("Chj SimGroup %d St %d",Sim,Station+1), "", Bins, - 0.5, Bins - 0.5);
 506       Chamberj->SetTitle(kFALSE);
 507       Chamberj->SetMarkerStyle(20);
 508       Chamberj->SetMarkerColor(2);
 509       Chamberj->SetMarkerSize(1.3);
 510
 511       if (Sim == 0)
 512         {
 513
 514           if ( (Sim == 0) && (SymGroups/2 == 1) )
 515             {
 516
 517               for (Int_t i = 1; i < 5; i++)
 518                  {
 519
 520                    Chamberi->SetBinContent(i, Chi[i-1]);
 521                    Chamberi->SetBinError(i, ChiErr[i-1]);
 522                    Chamberj->SetBinContent(i, Chj[i-1]);
 523                    Chamberj->SetBinError(i, ChjErr[i-1]);
 524
 525                  }
 526
 527             }
 528
 529           else
 530             {
 531
 532               Chamberi->SetBinContent(Sim+1, Chi[Sim]);
 533               Chamberi->SetBinContent((Sim+1)+(Bins/2), Chi[Sim+(Bins/2)]);
 534               Chamberi->SetBinError(Sim+1, ChiErr[Sim]);
 535               Chamberi->SetBinError((Sim+1)+(Bins/2), ChiErr[Sim+(Bins/2)]);
 536               Chamberj->SetBinContent(Sim+1, Chj[Sim]);
 537               Chamberj->SetBinContent((Sim+1)+(Bins/2), Chj[Sim+(Bins/2)]);
 538               Chamberj->SetBinError(Sim+1, ChjErr[Sim]);
 539               Chamberj->SetBinError((Sim+1)+(Bins/2), ChjErr[Sim+(Bins/2)]);
 540
 541             }
 542
 543       }
 544
 545       else
 546         {
 547
 548           Chamberi->SetBinContent(Sim+1, Chi[Sim]);
 549           Chamberi->SetBinContent((Bins/2)-(Sim-1), Chi[(Bins/2)-Sim]);
 550           Chamberi->SetBinContent((Sim+1)+(Bins/2), Chi[Sim+(Bins/2)]);
 551           Chamberi->SetBinContent(Bins-(Sim-1), Chi[Bins-Sim]);
 552
 553           Chamberi->SetBinError(Sim+1, ChiErr[Sim]);
 554           Chamberi->SetBinError((Bins/2)-(Sim-1), ChiErr[(Bins/2)-Sim]);
 555           Chamberi->SetBinError((Sim+1)+(Bins/2), ChiErr[Sim+(Bins/2)]);
 556           Chamberi->SetBinError(Bins-(Sim-1), ChiErr[Bins-Sim]);
 557
 558           Chamberj->SetBinContent(Sim+1, Chj[Sim]);
 559           Chamberj->SetBinContent((Bins/2)-(Sim-1), Chj[(Bins/2)-Sim]);
 560           Chamberj->SetBinContent((Sim+1)+(Bins/2), Chj[Sim+(Bins/2)]);
 561           Chamberj->SetBinContent(Bins-(Sim-1), Chj[Bins-Sim]);
 562
 563           Chamberj->SetBinError(Sim+1, ChjErr[Sim]);
 564           Chamberj->SetBinError((Bins/2)-(Sim-1), ChjErr[(Bins/2)-Sim]);
 565           Chamberj->SetBinError((Sim+1)+(Bins/2), ChjErr[Sim+(Bins/2)]);
 566           Chamberj->SetBinError(Bins-(Sim-1), ChjErr[Bins-Sim]);
 567
 568         }
 569
 570       Canvas->cd(Sim+1);
 571       Chamberi->Draw("e");
 572       Chamberj->Draw("esame");
 573
 574       TLegend *Legend = new TLegend (0.75, 0.8, 1.0, 1.0);
 575       Legend->SetBorderSize(0);
 576       Legend->SetFillColor(0);
 577       if ( (Station+1) < 3 ) Legend->SetTextSize(0.05);
 578       if ( (Station+1) > 2 ) Legend->SetTextSize(0.05);
 579       Legend->AddEntry(Chamberi,Form("Chamber %d",2*Station+1),"p");
 580       Legend->AddEntry(Chamberj,Form("Chamber %d",2*Station+2),"p");
 581       Legend->Draw();
 582
 583     }
 584
 585 }
 586
 587
 588 TH1D* GraphToHist(Bool_t ToBePlotted,const TGraphAsymmErrors* Graph,Int_t Chamber)
 589 {
 590
 591   Int_t nBins = 0;
 592   if ( Chamber < 5 ) nBins = 4;
 593   if ( (Chamber == 5) || (Chamber == 6) ) nBins = 18;
 594   if ( (Chamber > 6) && (Chamber < 11) ) nBins = 26;
 595   if ( Chamber == 11 ) nBins = 10;
 596   TArrayF  bins(nBins+1);
 597   TArrayF  y(nBins);
 598   TArrayF  ey(nBins);
 599   Double_t dx = 1;
 600   Double_t xmin = 10000;
 601   Double_t xmax = -10000;
 602   for (Int_t i = 0; i < nBins; i++) {
 603     Double_t x   = Graph->GetX()[i];
 604     Double_t exl = Graph->GetEXlow()[i];
 605     Double_t exh = Graph->GetEXhigh()[i];
 606     xmin             = TMath::Min(x-exl, xmin);
 607     xmax             = TMath::Max(x+exh, xmax);
 608     bins.fArray[i]   = x-exl;
 609     bins.fArray[i+1] = x+exh;
 610     Double_t dxi = exh+exl;
 611     if (dxi == 0 && i != 0) dxi = bins.fArray[i]-bins.fArray[i-1];
 612     if (dx == 0) dx  = dxi;
 613     else if (dxi != dx) dx = 0;
 614
 615     y.fArray[i]  = Graph->GetY()[i];
 616     ey.fArray[i] = TMath::Max(Graph->GetEYlow()[i],Graph->GetEYhigh()[i]);
 617
 618   }
 619   TString name(Graph->GetName());
 620   TString title(Graph->GetTitle());
 621   TH1D* h = 0;
 622   if (ToBePlotted == kTRUE) {
 623     h = new TH1D(name.Data(), title.Data(), nBins, bins[0]-dx/2 + 1.5, bins[nBins]-dx/2 + 1.5);
 624   }
 625   else {
 626     if ( Chamber < 5 ) h = new TH1D(name.Data(), title.Data(), 4, -0.5, 3.5);
 627     if ( (Chamber == 5) || (Chamber == 6) ) h = new TH1D(name.Data(), title.Data(), 18,  -0.5, 17.5);
 628     if ( (Chamber > 6) && (Chamber < 11) ) h = new TH1D(name.Data(), title.Data(), 26, -0.5, 25.5);
 629     if ( Chamber == 11 ) h = new TH1D(name.Data(), title.Data(), 10, -0.5, 9.5);
 630   }
 631   for (Int_t i = 1; i <= nBins; i++) {
 632     h->SetBinContent(i, y.fArray[i-1]);
 633     h->SetBinError(i, ey.fArray[i-1]);
 634   }
 635   h->SetMarkerStyle(Graph->GetMarkerStyle());
 636   h->SetMarkerColor(Graph->GetMarkerColor());
 637   h->SetMarkerSize(Graph->GetMarkerSize());
 638   h->SetDirectory(0);
 639
 640   return h;
 641 }