PWGLF/FORWARD/analysis/AliFMDDndeta.cxx

   1 // Code to analyse dN/deta from the forward analysis
   2 // This can plot the results
   3 // Also works for MC data
   4 //
   5 // -- Author: Hans Hjersing Dalsgaard <canute@gmail.com>
   6 #include "AliFMDDndeta.h"
   7 #include "TFile.h"
   8 #include "AliLog.h"
   9 #include "TH1.h"
  10 #include "AliFMDAnaParameters.h"
  11 #include "AliFMDAnaCalibSharingEfficiency.h"
  12 #include "TStyle.h"
  13 //#include "TObjArray.h"
  14 #include "TCanvas.h"
  15 #include "TLine.h"
  16 #include "TGraphErrors.h"
  17 #include "TGraphAsymmErrors.h"
  18 #include "TPad.h"
  19 #include "iostream"
  20 #include "TH3.h"
  21 #include "TMath.h"
  22 #include "TProfile.h"
  23 #include "TProfile2D.h"
  24 #include "TProfile3D.h"
  25 #include "TLegend.h"
  26 #include "TPad.h"
  27 #include "TLatex.h"
  28 #include "TStyle.h"
  29 #include "TF1.h"
  30
  31 #define SMALLNUMBER 0.0001
  32
  33 ClassImp(AliFMDDndeta)
  34 //_____________________________________________________________________
  35
  36 AliFMDDndeta::AliFMDDndeta()
  37 : TObject(),
  38   fList(0),
  39   fMultList(),
  40   fNbinsToCut(2),
  41   fVtxCut1(-10),
  42   fVtxCut2(10),
  43   fIsInit(kFALSE),
  44   fIsGenerated(),
  45   fPrimEvents(),
  46   fEvents(),
  47   fPrimdNdeta(),
  48   fDrawAll(kFALSE)
  49 {
  50   //AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  51   /* fDataObject = new TProfile3D("dataObject","dataObject",
  52                                pars->GetNetaBins(),-6,6,
  53                                20,0,2*TMath::Pi(),
  54                                pars->GetNvtxBins(),-0.5,pars->GetNvtxBins()-0.5);
  55   fDataObject->SetXTitle("#eta");
  56   fDataObject->SetYTitle("#varphi [radians]");
  57   fDataObject->SetZTitle("v_{z} [cm]");*/
  58
  59   fAnalysisNames[0] = "Hits";
  60   fAnalysisNames[1] = "HitsTrVtx";
  61   fAnalysisNames[2] = "dNdeta";
  62   fAnalysisNames[3] = "dNdetaTrVtx";
  63   fAnalysisNames[4] = "dNdetaNSD";
  64
  65   for(Int_t i=0; i<5;i++)
  66     fMultList[i] = new TList();
  67 }
  68 //_____________________________________________________________________
  69 void AliFMDDndeta::SetNames(Analysis what) {
  70   // Set names of histograms from analysis
  71
  72   switch(what) {
  73   case kHits :
  74     fPrimEvents.Form("nMCEventsNoCuts"); //was nMCEvents
  75     fEvents.Form("nEvents");
  76     fPrimdNdeta.Form("hMultvsEtaNoCuts");
  77     break;
  78   case kHitsTrVtx :
  79     fPrimEvents.Form("nMCEvents");
  80     fEvents.Form("nEvents");
  81     fPrimdNdeta.Form("hMultvsEtaNoCuts");
  82     break;
  83   case kMult :
  84     fPrimEvents.Form("nMCEventsNoCuts");
  85     fEvents.Form("nEvents");
  86     fPrimdNdeta.Form("hMultvsEtaNoCuts");
  87     break;
  88   case kMultTrVtx :
  89     fPrimEvents.Form("nMCEvents");
  90     fEvents.Form("nEvents");
  91     fPrimdNdeta.Form("hMultvsEta");
  92     break;
  93   case kMultNSD :
  94     fPrimEvents.Form("nMCEventsNSDNoCuts");
  95     fEvents.Form("nNSDEvents");
  96     fPrimdNdeta.Form("hMultvsEtaNSDNoCuts");
  97     break;
  98
  99   default:
 100     break;
 101   }
 102 }
 103 //_____________________________________________________________________
 104 const char* AliFMDDndeta::GetAnalysisName(Analysis what, UShort_t det, Char_t ring, Int_t vtxbin) {
 105   //Get names of histograms
 106   const char* name = "";
 107
 108   switch(what) {
 109   case kHits :
 110     name = Form("hits_NoCuts_FMD%d%c_vtxbin%d_proj",det,ring,vtxbin); //NoCuts added
 111     break;
 112   case kHitsTrVtx :
 113     name = Form("hits_FMD%d%c_vtxbin%d_proj",det,ring,vtxbin);
 114     break;
 115   case kMult :
 116     name = Form("dNdeta_FMD%d%c_vtxbin%d_proj",det,ring,vtxbin);
 117     break;
 118   case kMultTrVtx :
 119     name = Form("dNdeta_FMD%d%c_TrVtx_vtxbin%d_proj",det,ring,vtxbin);
 120     break;
 121   case kMultNSD :
 122     name = Form("dNdetaNSD_FMD%d%c_vtxbin%d_proj",det,ring,vtxbin);
 123     break;
 124
 125   default:
 126     break;
 127   }
 128   //std::cout<<name.Data()<<std::endl;
 129   return name;
 130 }
 131 //_____________________________________________________________________
 132 const char* AliFMDDndeta::GetPrimName(Analysis what, UShort_t det, Char_t ring, Int_t vtxbin) {
 133   //Get names of primaries
 134   const char* name = "";
 135
 136   switch(what) {
 137   case kHits :
 138     name = Form("hMCHits_nocuts_FMD%d%c_vtxbin%d",det,ring,vtxbin); //nocuts added
 139     break;
 140   case kHitsTrVtx :
 141     name = Form("hMCHits_FMD%d%c_vtxbin%d",det,ring,vtxbin);
 142     break;
 143   case kMult :
 144     name = Form("primmult_vtxbin%d",vtxbin);
 145     break;
 146   case kMultTrVtx :
 147     name = Form("primmult_NoCuts_vtxbin%d",vtxbin);
 148     break;
 149   case kMultNSD :
 150     name = Form("primmult_NSD_vtxbin%d",vtxbin);
 151     break;
 152   default:
 153     break;
 154   }
 155   return name;
 156 }
 157 //_____________________________________________________________________
 158 void AliFMDDndeta::GenerateHits(Analysis what) {
 159   //Generate the hits distributions
 160   AliFMDAnaParameters* pars =  AliFMDAnaParameters::Instance();
 161   TH2F* hTmp         =  pars->GetBackgroundCorrection(1,'I',0);
 162   Int_t nVertexBins  =  pars->GetNvtxBins();
 163   Int_t nEtaBins     =  hTmp->GetNbinsX();
 164   Float_t etaMin     =  hTmp->GetXaxis()->GetXmin();
 165   Float_t etaMax     =  hTmp->GetXaxis()->GetXmax();
 166
 167   for(Int_t i = 0; i<nVertexBins;i++) { //
 168     TH1F* hHits = new TH1F(Form("hMCHits_vtxbin%d_%s",i,fAnalysisNames[what].Data()),Form("hHits_vtxbin%d_%s",i,fAnalysisNames[what].Data()),nEtaBins,etaMin,etaMax);
 169     hHits->Sumw2();
 170     fMultList[what]->Add(hHits);
 171   }
 172
 173   for(Int_t det = 1; det<=3; det++) {
 174     Int_t maxRing = (det == 1 ? 0 : 1);
 175     for(Int_t ring = 0;ring<=maxRing;ring++) {
 176
 177       Char_t ringChar = (ring == 0 ? 'I' : 'O');
 178       for(Int_t v=0; v< nVertexBins; v++) {
 179         TH1F* hits = (TH1F*)fList->FindObject(GetPrimName(what,det,ringChar,v));
 180
 181         Int_t nNonZero = 0, nNonZeroInData = 0;
 182
 183         //removing edges
 184         for(Int_t i =1;i<=hits->GetNbinsX();i++) {
 185           if(hits->GetBinContent(i) !=0)
 186             nNonZero++;
 187         }
 188
 189         TH1F* sumMultHist = (TH1F*)fMultList[what]->FindObject(Form("hMCHits_vtxbin%d_%s",v,fAnalysisNames[what].Data()));
 190         for(Int_t i =1;i<=sumMultHist->GetNbinsX();i++) {
 191           if(hits->GetBinContent(i) == 0 ) continue;
 192
 193           nNonZeroInData++;
 194
 195           if(nNonZeroInData <=fNbinsToCut || nNonZeroInData > (nNonZero - fNbinsToCut)) {
 196             continue;
 197           }
 198           if(sumMultHist->GetBinContent(i) < SMALLNUMBER && TMath::Abs(hits->GetBinContent(i)) > 0){
 199               sumMultHist->SetBinContent(i,hits->GetBinContent(i));
 200               sumMultHist->SetBinError(i,hits->GetBinError(i));
 201             }
 202           else {
 203
 204
 205             Float_t sumofw = (1/TMath::Power(hits->GetBinError(i),2)) + (1/TMath::Power(sumMultHist->GetBinError(i),2));
 206             Float_t wav =  (hits->GetBinContent(i)*(1/TMath::Power(hits->GetBinError(i),2)) + sumMultHist->GetBinContent(i)*(1/TMath::Power(sumMultHist->GetBinError(i),2)))/sumofw;
 207             Float_t error = 1/TMath::Sqrt(sumofw);
 208             sumMultHist->SetBinContent(i, wav);
 209             sumMultHist->SetBinError(i, error);
 210           }
 211         }
 212       }
 213     }
 214   }
 215 }
 216
 217 //_____________________________________________________________________
 218 void AliFMDDndeta::Init(const Char_t* filename) {
 219   //Initialize everything from file
 220   TFile* fin = TFile::Open(filename);
 221
 222   if(!fin) {
 223     AliWarning("No file - exiting !");
 224     return;
 225   }
 226   AliFMDAnaParameters* pars =  AliFMDAnaParameters::Instance();
 227   //pars->Init();
 228
 229   TList* list = (TList*)fin->Get(Form("%s/BackgroundCorrected",pars->GetDndetaAnalysisName()));
 230
 231   if(!list) //an old file ? Perhaps...
 232     list = (TList*)fin->Get("BackgroundCorrected");
 233
 234   Init(list);
 235
 236 }
 237 //_____________________________________________________________________
 238 void AliFMDDndeta::Init(TList* list) {
 239   //Initialize everything from TList
 240
 241   if(!list) {
 242     AliWarning("No list - exiting !");
 243     return;
 244   }
 245
 246   fList = (TList*)list->Clone("inputList");
 247
 248   fIsGenerated[kHits]      = kFALSE;
 249   fIsGenerated[kMult]      = kFALSE;
 250   fIsGenerated[kMultTrVtx] = kFALSE;
 251   fIsGenerated[kHitsTrVtx] = kFALSE;
 252   fIsGenerated[kMultNSD]   = kFALSE;
 253
 254   fIsInit = kTRUE;
 255 }
 256 //_____________________________________________________________________
 257 void AliFMDDndeta::GenerateMult(Analysis what) {
 258   //Generate dNdeta
 259
 260   if(!fIsInit) {
 261     AliWarning("Not initialised - call Init to remedy");
 262     return;
 263   }
 264
 265   if(fIsGenerated[what]) {
 266     AliWarning("Already generated - have a look at the results!");
 267     return;
 268   }
 269   else fIsGenerated[what] = kTRUE;
 270
 271   SetNames(what);
 272
 273   if(what == kHits || what == kHitsTrVtx)
 274     GenerateHits(what);
 275
 276   TH1I* hMCEvents = (TH1I*)fList->FindObject(fPrimEvents.Data());
 277
 278   AliFMDAnaParameters* pars =  AliFMDAnaParameters::Instance();
 279
 280   TH2F* hTmp         =  pars->GetBackgroundCorrection(1,'I',0);
 281   Int_t nVertexBins  =  pars->GetNvtxBins();
 282   Int_t nEtaBins     =  hTmp->GetNbinsX();
 283   Float_t etaMin     =  hTmp->GetXaxis()->GetXmin();
 284   Float_t etaMax     =  hTmp->GetXaxis()->GetXmax();
 285
 286   for(Int_t i = 0; i<nVertexBins;i++) {
 287     TH1F* hMult = new TH1F(Form("hMult_vtxbin%d_%s",i,fAnalysisNames[what].Data()),Form("hMult_vtxbin%d_%s",i,fAnalysisNames[what].Data()),nEtaBins,etaMin,etaMax);
 288     hMult->Sumw2();
 289     fMultList[what]->Add(hMult);
 290   }
 291
 292   for(Int_t det = 1; det<=3;det++) {
 293     Int_t maxRing = (det == 1 ? 0 : 1);
 294     for(Int_t iring = 0; iring<=maxRing; iring++) {
 295       Char_t ringChar = (iring == 0 ? 'I' : 'O');
 296       //Int_t nsec  = (iring == 0 ? 20 : 40);
 297       TH1F* hRingMult= new TH1F(Form("hRingMult_FMD%d%c_%s",det,ringChar,fAnalysisNames[what].Data()),Form("hRingMult_FMD%d%c_%s",det,ringChar,fAnalysisNames[what].Data()),nEtaBins,etaMin,etaMax);
 298       fMultList[what]->Add(hRingMult);
 299       //      TProfile* phiprofile     = new TProfile(Form("dNdphiFMD%d%c",det,ringChar), Form("dNdphiFMD%d%c;#Phi",det,ringChar), nsec , 0, 2*TMath::Pi());
 300       //fMultList[what]->Add(phiprofile);
 301     }
 302   }
 303   TH1I*  hEvents  = (TH1I*)fList->FindObject(fEvents.Data());
 304   TH1F*  hPrimVtx = 0;
 305
 306
 307
 308   for(Int_t det = 1; det<=3; det++) {
 309     Int_t maxRing = (det == 1 ? 0 : 1);
 310     for(Int_t ring = 0;ring<=maxRing;ring++) {
 311
 312       Char_t ringChar = (ring == 0 ? 'I' : 'O');
 313       for(Int_t v=0; v< nVertexBins; v++) {
 314         if(det == 1) {
 315           if(what == kHits || what == kHitsTrVtx)
 316             hPrimVtx = (TH1F*)fMultList[what]->FindObject(Form("hMCHits_vtxbin%d_%s",v,fAnalysisNames[what].Data()));
 317           else
 318             hPrimVtx = (TH1F*)fList->FindObject(GetPrimName(what,det,ringChar,v));
 319
 320
 321           Float_t nPrimevents = hMCEvents->GetBinContent(v+1);
 322           Float_t xb = hPrimVtx->GetNbinsX();
 323           Float_t xr = hPrimVtx->GetXaxis()->GetXmax() - hPrimVtx->GetXaxis()->GetXmin();
 324           hPrimVtx->Scale(xb / xr );
 325           if(nPrimevents > 0)
 326             hPrimVtx->Scale(1/nPrimevents);
 327
 328         }
 329
 330         Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
 331         Float_t vtxZ1   = (delta*v) - pars->GetVtxCutZ();
 332         Float_t vtxZ2   = (delta*(v+1)) - pars->GetVtxCutZ();
 333
 334         if(vtxZ1< fVtxCut1 || vtxZ2 >fVtxCut2)
 335           continue;
 336         Float_t nEvents = (Float_t)hEvents->GetBinContent(v+1);
 337
 338         //TH1F* multhistproj = (TH1F*)fList->FindObject(Form("dNdeta_FMD%d%c_vtxbin%d_proj",det,ringChar,v));
 339
 340
 341
 342
 343
 344
 345         TH1F* multhistproj = (TH1F*)fList->FindObject(GetAnalysisName(what,det,ringChar,v));
 346         if(nEvents)
 347           multhistproj->Scale(1/nEvents);
 348         Float_t xnBins = multhistproj->GetNbinsX();
 349         Float_t xrange = multhistproj->GetXaxis()->GetXmax() - multhistproj->GetXaxis()->GetXmin();
 350         multhistproj->Scale(xnBins / xrange);
 351
 352         //TH2F* multhist = (TH2F*)fList->FindObject(Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,v));
 353
 354
 355         //if(nEvents)
 356         //  multhist->Scale(1/nEvents);
 357         /*      for(Int_t i=1;i<=multhist->GetNbinsX();i++) {
 358           for(Int_t j=1;j<=multhist->GetNbinsY();j++) {
 359
 360             if (multhist->GetBinContent(i,j) <= 0.0001) continue;
 361             //std::cout<<multhist->GetBinContent(i,j)<<std::endl;
 362             fDataObject->Fill(multhist->GetXaxis()->GetBinCenter(i),
 363                               multhist->GetYaxis()->GetBinCenter(j),
 364                               v,
 365                               multhist->GetBinContent(i,j), multhist->GetBinError(i,j));
 366             //fDataObject->SetBinContent(i,j,v,multhist->GetBinContent(i,j));
 367             //fDataObject->SetBinError(i,j,v,multhist->GetBinError(i,j));
 368           }
 369         }
 370         */
 371         //if(nEvents)
 372         //  multhist->Scale(1/nEvents);
 373         //if(ringChar == 'O')
 374         //  multhist->RebinY(2);
 375         //removing edges
 376
 377         Int_t nNonZero = 0, nNonZeroInData = 0;
 378
 379         for(Int_t i =1;i<=multhistproj->GetNbinsX();i++) {
 380           if(multhistproj->GetBinContent(i) !=0)
 381             nNonZero++;
 382         }
 383         Int_t nBinsOld = fNbinsToCut;
 384         //if(det == 1 && ringChar =='I') {
 385         //  fNbinsToCut = 0;
 386         //      }
 387         TH1F* hRingMult = (TH1F*)fMultList[what]->FindObject(Form("hRingMult_FMD%d%c_%s",det,ringChar,fAnalysisNames[what].Data()));
 388
 389         for(Int_t i=1; i<=hRingMult->GetNbinsX(); i++) {
 390           if(multhistproj->GetBinContent(i)!=0) {
 391             nNonZeroInData++;
 392             if(nNonZeroInData <=fNbinsToCut || nNonZeroInData > (nNonZero - fNbinsToCut)) {
 393               continue;
 394             }
 395             Float_t oldweight = 0;
 396             Float_t oldwav    = 0;
 397             if(hRingMult->GetBinError(i)>0) {
 398               oldweight = 1/TMath::Power(hRingMult->GetBinError(i),2);
 399               oldwav    = oldweight*hRingMult->GetBinContent(i);
 400             }
 401             Float_t  weight = 1/TMath::Power(multhistproj->GetBinError(i),2);
 402             Float_t  wav    = oldwav + weight*multhistproj->GetBinContent(i);
 403             Float_t  sumofw = oldweight + weight;
 404             if(sumofw) {
 405               Float_t error = 1/TMath::Sqrt(sumofw);
 406
 407               hRingMult->SetBinContent(i,wav/sumofw);
 408               hRingMult->SetBinError(i,error);
 409
 410             }
 411           }
 412         }
 413         nNonZeroInData = 0;
 414
 415         TH1F* sumMultHist = (TH1F*)fMultList[what]->FindObject(Form("hMult_vtxbin%d_%s",v,fAnalysisNames[what].Data()));
 416
 417         for(Int_t i =1;i<=sumMultHist->GetNbinsX();i++) {
 418           if(multhistproj->GetBinContent(i) != 0 ) {
 419             nNonZeroInData++;
 420
 421             if(nNonZeroInData <=fNbinsToCut || nNonZeroInData > (nNonZero - fNbinsToCut)) {
 422               continue;
 423             }
 424             if(sumMultHist->GetBinContent(i) < SMALLNUMBER && TMath::Abs(multhistproj->GetBinContent(i)) > 0){
 425               sumMultHist->SetBinContent(i,multhistproj->GetBinContent(i));
 426               sumMultHist->SetBinError(i,multhistproj->GetBinError(i));
 427             }
 428             else {
 429
 430
 431               Float_t sumofw = (1/TMath::Power(multhistproj->GetBinError(i),2)) + (1/TMath::Power(sumMultHist->GetBinError(i),2));
 432               Float_t wav =  (multhistproj->GetBinContent(i)*(1/TMath::Power(multhistproj->GetBinError(i),2)) + sumMultHist->GetBinContent(i)*(1/TMath::Power(sumMultHist->GetBinError(i),2)))/sumofw;
 433               Float_t error = 1/TMath::Sqrt(sumofw);
 434               sumMultHist->SetBinContent(i, wav);
 435               sumMultHist->SetBinError(i, error);
 436             }
 437
 438
 439           }
 440         }
 441         fNbinsToCut = nBinsOld;
 442       }
 443
 444     }
 445   }
 446
 447   TH1F* sumMult  = new TH1F(Form("dNdeta_%s",fAnalysisNames[what].Data()),Form("dNdeta_%s",fAnalysisNames[what].Data()),nEtaBins,etaMin,etaMax);
 448   sumMult->Sumw2();
 449   fMultList[what]->Add(sumMult);
 450   TH1F* primMult  = new TH1F(Form("primary_%s",fAnalysisNames[what].Data()),Form("primary_%s",fAnalysisNames[what].Data()),nEtaBins,etaMin,etaMax);
 451   primMult->Sumw2();
 452   fMultList[what]->Add(primMult);
 453
 454   Float_t wav  = 0, sumofw=0, weight = 0, primwav  = 0, primsumofw=0, primweight = 0;//, etaofbin = 0;
 455   for(Int_t i =1; i<=sumMult->GetNbinsX();i++) {
 456
 457     for(Int_t v = 0; v<nVertexBins;v++) {
 458       if(what == kHits || what == kHitsTrVtx)
 459         hPrimVtx = (TH1F*)fMultList[what]->FindObject(Form("hMCHits_vtxbin%d_%s",v,fAnalysisNames[what].Data()));
 460       else
 461         hPrimVtx = (TH1F*)fList->FindObject(GetPrimName(what,0,0,v));
 462       TH1F* sumMultHist = (TH1F*)fMultList[what]->FindObject(Form("hMult_vtxbin%d_%s",v,fAnalysisNames[what].Data()));
 463       //etaofbin += sumMultHist->GetBinCenter(i);
 464       // if( hPrimVtx->GetBinContent(i)!=0 && hPrimVtx->GetBinError(i)>0.001 && sumMultHist->GetBinContent(i)!=0) {
 465       if( TMath::Abs(hPrimVtx->GetBinContent(i)) > 0) {
 466         primweight = 1/TMath::Power(hPrimVtx->GetBinError(i),2);
 467         primwav    = primwav + primweight*hPrimVtx->GetBinContent(i);
 468         primsumofw = primsumofw + primweight;
 469       }
 470       if(TMath::Abs(sumMultHist->GetBinContent(i)) > 0) {
 471
 472         weight = 1/TMath::Power(sumMultHist->GetBinError(i),2);
 473         wav    = wav + weight*sumMultHist->GetBinContent(i);
 474         sumofw = sumofw + weight;
 475       }
 476
 477     }
 478     if( primsumofw !=0 ) {// && sumofw !=0) {
 479       Float_t primerror = 1/TMath::Sqrt(primsumofw);
 480
 481       primMult->SetBinContent(i,primwav/primsumofw);
 482       primMult->SetBinError(i,primerror);
 483     }
 484     else {
 485       primMult->SetBinContent(i,0);
 486       primMult->SetBinError(i,0);
 487     }
 488
 489     if(sumofw) {
 490
 491       Float_t error = 1/TMath::Sqrt(sumofw);
 492
 493       sumMult->SetBinContent(i,wav/sumofw);
 494       sumMult->SetBinError(i,error);
 495     }
 496
 497
 498     wav = 0;
 499     sumofw = 0;
 500     weight = 0;
 501     primwav = 0;
 502     primsumofw = 0;
 503     primweight = 0;
 504
 505   }
 506
 507 }
 508 //_____________________________________________________________________
 509 void AliFMDDndeta::DrawDndeta(Analysis what, Int_t rebin, Bool_t realdata, TString pythiafile) {
 510   //Draw everything
 511   if(!fIsInit) {
 512     AliWarning("Not initialised - call Init to remedy");
 513     return;
 514   }
 515
 516   if(!fIsGenerated[what]) {
 517     AliWarning("Not generated - generate first then draw!");
 518     return;
 519   }
 520   AliFMDAnaParameters* pars =  AliFMDAnaParameters::Instance();
 521   SetNames(what);
 522   //UA5 data NSD
 523   Float_t x[19] = {0.125,0.375,0.625,0.875,1.125,1.375,1.625,1.875,2.125,2.375,
 524                    2.625,2.875,3.125,3.375,3.625,3.875,4.125,4.375,4.625};
 525   Float_t x2[19] = {-0.125,-0.375,-0.625,-0.875,-1.125,-1.375,-1.625,-1.875,
 526                     -2.125,-2.375,-2.625,-2.875,-3.125,-3.375,-3.625,-3.875,
 527                     -4.125,-4.375,-4.625};
 528
 529
 530   Float_t y[19] = {3.48, 3.38, 3.52, 3.68, 3.71, 3.86, 3.76, 3.66, 3.72 ,3.69,
 531                    3.56, 3.41, 3.15, 3.09, 2.74, 2.73, 2.32, 1.99, 1.69};
 532
 533   Float_t ey[19] = {0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,
 534                     0.07,0.07,0.08,0.08,0.09,0.09,0.1,0.13};
 535
 536   TGraphErrors* graph = new TGraphErrors(19,x,y,NULL,ey);
 537   TGraphErrors* graph2 = new TGraphErrors(19,x2,y,NULL,ey);
 538   //UA5 data INEL
 539   Float_t xinel[19] = {0.125, 0.375, 0.625, 0.875, 1.125, 1.375, 1.625, 1.875, 2.125, 2.375, 2.625, 2.875, 3.125, 3.375, 3.625, 3.875, 4.125, 4.375, 4.625};
 540   Float_t xinel2[19] = {-0.125, -0.375, -0.625, -0.875, -1.125, -1.375, -1.625, -1.875, -2.125, -2.375, -2.625, -2.875, -3.125, -3.375, -3.625, -3.875, -4.125, -4.375, -4.625};
 541   Float_t yinel[19] = {3.14, 3.04, 3.17, 3.33, 3.33, 3.53, 3.46, 3.41, 3.45, 3.39, 3.07, 3.07, 2.93, 2.93, 2.55, 2.48, 2.18, 1.91, 1.52};
 542   Float_t eyinel[19] = {0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.08, 0.08, 0.09, 0.09, 0.1, 0.13};
 543  TGraphErrors* graphinel = new TGraphErrors(19,xinel,yinel,NULL,eyinel);
 544  TGraphErrors* graphinel2 = new TGraphErrors(19,xinel2,yinel,NULL,eyinel);
 545
 546   graph->SetMarkerStyle(22);
 547   graph->SetMarkerColor(kBlue);
 548   graphinel->SetMarkerStyle(20);
 549   graphinel->SetMarkerColor(kGreen);
 550   graphinel2->SetMarkerStyle(24);
 551   graphinel2->SetMarkerColor(kGreen);
 552   //graph->Draw("P");
 553   graph2->SetMarkerStyle(26);
 554   graph2->SetMarkerColor(kBlue);
 555   graphinel->GetHistogram()->SetStats(kFALSE);
 556   graphinel2->GetHistogram()->SetStats(kFALSE);
 557
 558
 559   //End UA5 data
 560
 561   //Published ALICE data
 562   // Plot: p7742_d1x1y1
 563
 564   // INEL
 565
 566   TGraphAsymmErrors* p7742D1x1y1 = 0;
 567   double p7742D1x1y1xval[] = { -1.3, -1.1, -0.9, -0.7, -0.5, -0.3, -0.1, 0.1, 0.3,
 568     0.5, 0.7, 0.9, 1.1, 1.3 };
 569   double p7742D1x1y1xerrminus[] = { 0.09999999999999987, 0.09999999999999987, 0.09999999999999998, 0.10000000000000009, 0.09999999999999998, 0.10000000000000003, 0.1, 0.1, 0.09999999999999998,
 570     0.09999999999999998, 0.09999999999999998, 0.09999999999999998, 0.10000000000000009, 0.10000000000000009 };
 571   double p7742D1x1y1xerrplus[] = { 0.10000000000000009, 0.10000000000000009, 0.09999999999999998, 0.09999999999999998, 0.09999999999999998, 0.09999999999999998, 0.1, 0.1, 0.10000000000000003,
 572     0.09999999999999998, 0.10000000000000009, 0.09999999999999998, 0.09999999999999987, 0.09999999999999987 };
 573   double p7742D1x1y1yval[] = { 3.28, 3.28, 3.22, 3.12, 3.06, 3.02, 2.98, 3.02, 3.02,
 574     3.05, 3.15, 3.21, 3.26, 3.33 };
 575   double p7742D1x1y1yerrminus[] = { 0.06324555320336758, 0.06324555320336758, 0.06324555320336758, 0.06324555320336758, 0.06324555320336758, 0.05385164807134505, 0.05385164807134505, 0.05385164807134505, 0.05385164807134505,
 576     0.06324555320336758, 0.06324555320336758, 0.06324555320336758, 0.06324555320336758, 0.06324555320336758 };
 577   double p7742D1x1y1yerrplus[] = { 0.08246211251235322, 0.08246211251235322, 0.08246211251235322, 0.08246211251235322, 0.08246211251235322, 0.08246211251235322, 0.07280109889280519, 0.08246211251235322, 0.08246211251235322,
 578     0.08246211251235322, 0.08246211251235322, 0.08246211251235322, 0.08246211251235322, 0.08246211251235322 };
 579   int p7742D1x1y1numpoints = 14;
 580   p7742D1x1y1 = new TGraphAsymmErrors(p7742D1x1y1numpoints, p7742D1x1y1xval, p7742D1x1y1yval, p7742D1x1y1xerrminus, p7742D1x1y1xerrplus, p7742D1x1y1yerrminus, p7742D1x1y1yerrplus);
 581   p7742D1x1y1->SetName("/HepData/7742/d1x1y1");
 582   p7742D1x1y1->SetTitle("/HepData/7742/d1x1y1");
 583   p7742D1x1y1->SetMarkerStyle(21);
 584   p7742D1x1y1->SetMarkerColor(kRed);
 585   // p7742D1x1y1->Draw("Psame");
 586
 587
 588   //NSD
 589   TGraphAsymmErrors*  p7742D2x1y1 = 0;
 590  double p7742D2x1y1xval[] = { -1.3, -1.1, -0.9, -0.7, -0.5, -0.3, -0.1, 0.1, 0.3,
 591     0.5, 0.7, 0.9, 1.1, 1.3 };
 592   double p7742D2x1y1xerrminus[] = { 0.09999999999999987, 0.09999999999999987, 0.09999999999999998, 0.10000000000000009, 0.09999999999999998, 0.10000000000000003, 0.1, 0.1, 0.09999999999999998,
 593     0.09999999999999998, 0.09999999999999998, 0.09999999999999998, 0.10000000000000009, 0.10000000000000009 };
 594   double p7742D2x1y1xerrplus[] = { 0.10000000000000009, 0.10000000000000009, 0.09999999999999998, 0.09999999999999998, 0.09999999999999998, 0.09999999999999998, 0.1, 0.1, 0.10000000000000003,
 595     0.09999999999999998, 0.10000000000000009, 0.09999999999999998, 0.09999999999999987, 0.09999999999999987 };
 596   double p7742D2x1y1yval[] = { 3.9, 3.89, 3.81, 3.7, 3.64, 3.59, 3.53, 3.58, 3.59,
 597     3.61, 3.74, 3.8, 3.87, 3.95 };
 598   double p7742D2x1y1yerrminus[] = { 0.13341664064126335, 0.13152946437965907, 0.13152946437965907, 0.1216552506059644, 0.1216552506059644, 0.1216552506059644, 0.1216552506059644, 0.1216552506059644, 0.1216552506059644,
 599     0.1216552506059644, 0.1216552506059644, 0.13152946437965907, 0.13152946437965907, 0.13341664064126335 };
 600   double p7742D2x1y1yerrplus[] = { 0.13341664064126335, 0.13152946437965907, 0.13152946437965907, 0.1216552506059644, 0.1216552506059644, 0.1216552506059644, 0.1216552506059644, 0.1216552506059644, 0.1216552506059644,
 601     0.1216552506059644, 0.1216552506059644, 0.13152946437965907, 0.13152946437965907, 0.13341664064126335 };
 602   int p7742D2x1y1numpoints = 14;
 603
 604   p7742D2x1y1 = new TGraphAsymmErrors(p7742D2x1y1numpoints, p7742D2x1y1xval, p7742D2x1y1yval, p7742D2x1y1xerrminus, p7742D2x1y1xerrplus, p7742D2x1y1yerrminus, p7742D2x1y1yerrplus);
 605   p7742D2x1y1->SetName("/HepData/7742/d2x1y1");
 606   p7742D2x1y1->SetTitle("/HepData/7742/d2x1y1");
 607   p7742D2x1y1->SetMarkerStyle(21);
 608   p7742D2x1y1->SetMarkerColor(kRed);
 609   //p7742D2x1y1.Draw("AP");
 610
 611
 612
 613
 614
 615
 616   //End official data
 617
 618   // CMS published NSD data
 619
 620   TGraphAsymmErrors* p7743D8x1y1 = 0;
 621   double p7743D8x1y1xval[] = { -2.25, -1.75, -1.25, -0.75, -0.25, 0.25, 0.75, 1.25, 1.75,
 622     2.25 };
 623   double p7743D8x1y1xerrminus[] = { 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25 };
 624   double p7743D8x1y1xerrplus[] = { 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25 };
 625   double p7743D8x1y1yval[] = { 3.6, 3.73, 3.62, 3.54, 3.48, 3.48, 3.54, 3.62, 3.73,  3.6 };
 626   double p7743D8x1y1yerrminus[] = { 0.13, 0.14, 0.13, 0.13, 0.13, 0.13, 0.13, 0.13, 0.14,0.13 };
 627   double p7743D8x1y1yerrplus[] = { 0.13, 0.14, 0.13, 0.13, 0.13, 0.13, 0.13, 0.13, 0.14, 0.13 };
 628   int p7743D8x1y1numpoints = 10;
 629   p7743D8x1y1 = new TGraphAsymmErrors(p7743D8x1y1numpoints, p7743D8x1y1xval, p7743D8x1y1yval, p7743D8x1y1xerrminus, p7743D8x1y1xerrplus, p7743D8x1y1yerrminus, p7743D8x1y1yerrplus);
 630
 631   p7743D8x1y1->SetMarkerStyle(20);
 632   p7743D8x1y1->SetMarkerColor(kBlack);
 633
 634   //End CMS
 635
 636
 637   TH1I* hEvents   = (TH1I*)fList->FindObject(fEvents.Data());
 638
 639   TH1I* hMCEvents = (TH1I*)fList->FindObject(fPrimEvents.Data());
 640
 641    //SPD part
 642   TH1D* hSPDana = (TH1D*)fList->FindObject(Form("dNdeta_SPD_vtxbin%d_proj",5));
 643
 644   for(Int_t nn = 0; nn < pars->GetNvtxBins() ; nn++) {
 645     TH1D* hSPDanalysis = (TH1D*)fList->FindObject(Form("dNdeta_SPD_vtxbin%d_proj",nn));
 646     TH1D* hSPDanalysisTrVtx = (TH1D*)fList->FindObject(Form("dNdetaTrVtx_SPD_vtxbin%d_proj",nn));
 647     TH1D* hSPDanalysisNSD = (TH1D*)fList->FindObject(Form("dNdetaNSD_SPD_vtxbin%d_proj",nn));
 648
 649     Float_t nEventsSPD = (Float_t)hEvents->GetBinContent(nn+1);
 650     if(!nEventsSPD) continue;
 651     hSPDanalysis->Scale(1/nEventsSPD);
 652     hSPDanalysisTrVtx->Scale(1/nEventsSPD);
 653     hSPDanalysisNSD->Scale(1/nEventsSPD);
 654   }
 655   Float_t lowlimits[10] = {-1,-1.25,-1.5,-1.75,-1.9,-1.9,-1.9,-1.9,-1.9,-1.9};
 656   Float_t highlimits[10] = {1.9,1.9,1.9,1.9,1.9,1.9,1.75,1.5,1.25,1};
 657
 658   TH1F* hSPDcombi = new TH1F("SPDcombi","SPDcombi",hSPDana->GetNbinsX(),hSPDana->GetXaxis()->GetXmin(),hSPDana->GetXaxis()->GetXmax());
 659   TH1D* hSPDanalysis = 0;
 660   for(Int_t kk = 1; kk <=hSPDana->GetNbinsX(); kk++) {
 661
 662     Float_t weight = 0, wav=0,sumofw = 0;
 663
 664     // if(hSPDana->GetBinCenter(kk) < lowlimits[nn])
 665     //  continue;
 666     //if(hSPDana->GetBinCenter(kk) > highlimits[nn])
 667     //  continue;
 668
 669     for(Int_t nn =0; nn < pars->GetNvtxBins() ; nn++) {
 670       Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
 671       Float_t vtxZ1   = (delta*nn) - pars->GetVtxCutZ();
 672       Float_t vtxZ2   = (delta*(nn+1)) - pars->GetVtxCutZ();
 673
 674       if(vtxZ1<fVtxCut1 || vtxZ2 >fVtxCut2)
 675         continue;
 676
 677       //std::cout<<vtxZ1<<"   "<<vtxZ2<<std::endl;
 678
 679       if(what == kMult)
 680         hSPDanalysis = (TH1D*)fList->FindObject(Form("dNdeta_SPD_vtxbin%d_proj",nn));
 681       else if(what == kMultNSD)
 682         hSPDanalysis = (TH1D*)fList->FindObject(Form("dNdetaNSD_SPD_vtxbin%d_proj",nn));
 683       else if(what == kMultTrVtx)
 684         hSPDanalysis = (TH1D*)fList->FindObject(Form("dNdetaTrVtx_SPD_vtxbin%d_proj",nn));
 685       else continue;
 686       if(hSPDanalysis->GetBinCenter(kk) < lowlimits[nn])
 687         continue;
 688       if(hSPDanalysis->GetBinCenter(kk) > highlimits[nn])
 689         continue;
 690       //std::cout<<hSPDanalysis->GetBinCenter(kk)<<"  "<<lowlimits[nn]<<std::endl;
 691
 692       Float_t mult = hSPDanalysis->GetBinContent(kk);
 693       Float_t error = hSPDanalysis->GetBinError(kk);
 694       /*
 695       if(mult > 0 && hSPDanalysis->GetBinContent(kk-1) < SMALLNUMBER)
 696         continue;
 697       if(mult > 0 && hSPDanalysis->GetBinContent(kk-2) < SMALLNUMBER)
 698         continue;
 699
 700       if(mult > 0 && hSPDanalysis->GetBinContent(kk+1) < SMALLNUMBER)
 701         continue;
 702       if(mult > 0 && hSPDanalysis->GetBinContent(kk+2) < SMALLNUMBER)
 703         continue;
 704       */
 705       if(mult > 0) {
 706         weight = 1/TMath::Power(error,2);
 707         wav    = wav + weight*mult;
 708         sumofw = sumofw + weight;
 709       }
 710
 711     }
 712
 713     if(sumofw && wav) {
 714       Float_t errorTotal = 1/TMath::Sqrt(sumofw);
 715       hSPDcombi->SetBinContent(kk,wav/sumofw);
 716       hSPDcombi->SetBinError(kk,errorTotal);
 717     }
 718   }
 719
 720
 721   Float_t xb1 = hSPDcombi->GetNbinsX();
 722   Float_t xr1 = hSPDcombi->GetXaxis()->GetXmax() - hSPDcombi->GetXaxis()->GetXmin();
 723   hSPDcombi->Scale(xb1 / xr1);
 724   //hSPDcombi->Rebin(rebin);
 725   //hSPDcombi->Scale(1/(Float_t)rebin);
 726   //RebinHistogram(hSPDcombi,rebin);
 727   hSPDcombi->SetMarkerStyle(29);
 728   hSPDcombi->SetMarkerColor(kBlue);
 729   //hSPDcombi->DrawCopy("same");
 730
 731
 732   TCanvas* c1 = new TCanvas("cvertexbins","Overlaps",1400,800);
 733   c1->Divide(5,2);
 734   TCanvas* cCorrections = 0;
 735   TCanvas* cCorrectionsPhi = 0;
 736   if(fDrawAll) {
 737     cCorrections = new TCanvas("corrections","Correction vs Eta",1400,800);
 738     cCorrections->Divide(5,2);
 739
 740     cCorrectionsPhi = new TCanvas("correctionsPhi","Correction vs Phi",1400,800);
 741     cCorrectionsPhi->Divide(5,2);
 742   }
 743   //TCanvas* cphi = new TCanvas("dNdphi","dNdphi",1280,1024);
 744   // cphi->Divide(3,2);
 745
 746   Int_t nVertexBins  =  pars->GetNvtxBins();
 747
 748
 749   //Int_t npadphi = 0;
 750   TH1F*  hPrimVtx = 0;
 751   for(Int_t det = 1; det<=3; det++) {
 752     Int_t maxRing = (det == 1 ? 0 : 1);
 753     for(Int_t ring = 0;ring<=maxRing;ring++) {
 754       //  npadphi++;
 755
 756       for(Int_t v=0; v< nVertexBins; v++) {
 757         Char_t ringChar = (ring == 0 ? 'I' : 'O');
 758         Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
 759         Float_t vtxZ1   = (delta*v) - pars->GetVtxCutZ();
 760         Float_t vtxZ2   = (delta*(v+1)) - pars->GetVtxCutZ();
 761
 762         if(vtxZ1<fVtxCut1 || vtxZ2 >fVtxCut2)
 763           continue;
 764         Float_t nEvents = (Float_t)hEvents->GetBinContent(v+1);
 765         if(fDrawAll) {
 766         TH2F* hBgCor  = pars->GetBackgroundCorrection(det,ringChar,v);
 767         if(hBgCor)  {
 768           TH1D* hBgProj = hBgCor->ProjectionX(Form("hBgProj_FMD%d%c_vtxbin%d",det,ringChar,v));
 769           TH1D* hBgProjPhi = hBgCor->ProjectionY(Form("hBgProjPhi_FMD%d%c_vtxbin%d",det,ringChar,v));
 770
 771           hBgProjPhi->SetName(Form("FMD%d%c",det,ringChar));
 772           hBgProjPhi->SetTitle("");//Form("FMD%d",det));
 773           hBgProj->SetTitle("");
 774           //Float_t scalefactor = (hBgProj->GetXaxis()->GetXmax() - hBgProj->GetXaxis()->GetXmin()) / (Float_t)hBgProj->GetNbinsX();
 775           Float_t scalefactor = 1/(Float_t)hBgCor->GetNbinsY();
 776
 777           Float_t nFilledEtaBins = 0;
 778
 779           for(Int_t jj = 1; jj<=hBgProj->GetNbinsX(); jj++) {
 780             if(hBgProj->GetBinContent(jj) > 0.01)
 781               nFilledEtaBins++;
 782           }
 783
 784           Float_t scalefactorPhi = 1/nFilledEtaBins;
 785
 786
 787           hBgProj->Scale(scalefactor);
 788           cCorrections->cd(v+1);
 789
 790           hBgProj->SetMarkerColor(det + 2*ring);
 791           hBgProj->SetMarkerStyle(19 + det + 2*ring );
 792
 793           if((det + 2*ring) == 5)  hBgProj->SetMarkerColor(det + 2*ring+1);
 794           if((19 + det + 2*ring ) == 24)  hBgProj->SetMarkerStyle(29 );
 795
 796           if(det == 1) {
 797             hBgProj->GetYaxis()->SetRangeUser(0,4);
 798             hBgProj->SetStats(kFALSE);
 799             hBgProj->SetXTitle("#eta");
 800             hBgProj->DrawCopy("PE");
 801             TLatex* l = new TLatex(0.14,0.92,Form("Vtx range [%.1f, %.1f]",vtxZ1,vtxZ2));
 802             l->SetNDC(kTRUE);
 803             l->Draw();
 804           }
 805           else
 806             hBgProj->DrawCopy("PEsame");
 807
 808           hBgProjPhi->Scale(scalefactorPhi);
 809           cCorrectionsPhi->cd(v+1);
 810           hBgProjPhi->SetMarkerColor(det + 2*ring);
 811           hBgProjPhi->SetMarkerStyle(19 + det + 2*ring );
 812           if((det + 2*ring) == 5)  hBgProjPhi->SetMarkerColor(det + 2*ring+1);
 813           if((19 + det + 2*ring ) == 24)  hBgProjPhi->SetMarkerStyle(29 );
 814           if(det == 1) {
 815             hBgProjPhi->GetYaxis()->SetRangeUser(1,5);
 816             hBgProjPhi->SetStats(kFALSE);
 817
 818             hBgProjPhi->SetXTitle("#Phi");
 819             hBgProjPhi->DrawCopy("PE");
 820
 821
 822           }
 823           else
 824             hBgProjPhi->DrawCopy("PEsame");
 825
 826         }
 827
 828         }
 829
 830         TH1F* multhistproj = (TH1F*)fList->FindObject(GetAnalysisName(what,det,ringChar,v));
 831
 832         c1->cd(v+1);
 833
 834         if(det==1)  {
 835           multhistproj->SetMarkerStyle(20);
 836           multhistproj->SetMarkerColor(1);
 837         }
 838         if(det==2 && ringChar=='I') {
 839           multhistproj->SetMarkerStyle(21);
 840           multhistproj->SetMarkerColor(2);
 841         }
 842         if(det==2 && ringChar=='O') {
 843           multhistproj->SetMarkerStyle(3);
 844           multhistproj->SetMarkerColor(3);
 845         }
 846         if(det==3 && ringChar=='I') {
 847           multhistproj->SetMarkerStyle(22);
 848           multhistproj->SetMarkerColor(4);
 849         }
 850         if(det==3 && ringChar=='O') {
 851           multhistproj->SetMarkerStyle(23);
 852           multhistproj->SetMarkerColor(6);
 853         }
 854
 855
 856
 857         if(det == 1) {
 858           if(what == kHits || what == kHitsTrVtx)
 859             hPrimVtx = (TH1F*)fMultList[what]->FindObject(Form("hMCHits_vtxbin%d_%s",v,fAnalysisNames[what].Data()));
 860           else
 861             hPrimVtx = (TH1F*)fList->FindObject(GetPrimName(what,det,ringChar,v));
 862           // std::cout<<hPrimVtx<<"   "<<kHits<<"   "<<kHitsTrVtx<<"   "<<what<<std::endl;
 863           //std::cout<<Form("hMCHits_vtxbin%d_%s",v,fAnalysisNames[what].Data())<<std::endl;
 864           hPrimVtx->SetTitle("");
 865           TLatex* l = new TLatex(0.14,0.92,Form("Vtx range [%.1f, %.1f], %d events",vtxZ1,vtxZ2,(Int_t)nEvents));
 866           l->SetNDC(kTRUE);
 867           hPrimVtx->SetFillColor(kGray);
 868           hPrimVtx->SetLabelFont(132,"xyz");
 869           hPrimVtx->SetStats(0000);
 870           hPrimVtx->GetXaxis()->SetTitle("#eta");
 871           hPrimVtx->GetYaxis()->SetRangeUser(0,6);
 872           hPrimVtx->DrawCopy("E3");
 873           l->Draw();
 874           multhistproj->DrawCopy("same");
 875           TH1D* hSPDanavtxbin = 0;
 876           if(what == kMult || what == kMultTrVtx || what == kMultNSD)  {
 877
 878           if(what == kMult)
 879             hSPDanavtxbin = (TH1D*)fList->FindObject(Form("dNdeta_SPD_vtxbin%d_proj",v));
 880           if(what == kMultTrVtx)
 881             hSPDanavtxbin = (TH1D*)fList->FindObject(Form("dNdetaTrVtx_SPD_vtxbin%d_proj",v));
 882
 883           if(what == kMultNSD)
 884             hSPDanavtxbin = (TH1D*)fList->FindObject(Form("dNdetaNSD_SPD_vtxbin%d_proj",v));
 885
 886           hSPDanavtxbin->SetMarkerColor(kBlue);
 887           hSPDanavtxbin->SetMarkerStyle(kStar);
 888           hSPDanavtxbin->Scale(xb1 / xr1);
 889           hSPDanavtxbin->DrawCopy("same");
 890
 891           if(what != kMultNSD) {
 892             graphinel->Draw("sameP");
 893             graphinel2->Draw("sameP");
 894           }
 895           else
 896            {
 897              graph->Draw("sameP");
 898              graph2->Draw("sameP");
 899            }
 900           }
 901         }
 902         else
 903           multhistproj->DrawCopy("same");
 904
 905
 906       }
 907     }
 908   }
 909
 910   //Legends for corrections
 911   if(fDrawAll) {
 912     for(Int_t v=0; v< nVertexBins; v++) {
 913       TPad* pad= (TPad*)cCorrectionsPhi->cd(v+1);
 914       pad->BuildLegend(0.15,0.45,0.45,0.9);
 915       Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
 916       Float_t vtxZ1   = (delta*v) - pars->GetVtxCutZ();
 917       Float_t vtxZ2   = (delta*(v+1)) - pars->GetVtxCutZ();
 918       TLatex* l = new TLatex(0.14,0.92,Form("Vtx range [%.1f, %.1f]",vtxZ1,vtxZ2));
 919       l->SetNDC(kTRUE);
 920       l->Draw();
 921     }
 922   }
 923   for(Int_t v=0; v< nVertexBins; v++) {
 924     TH1F* sumMultHist = (TH1F*)fMultList[what]->FindObject(Form("hMult_vtxbin%d_%s",v,fAnalysisNames[what].Data()));
 925     c1->cd(v+1);
 926     sumMultHist->SetMarkerStyle(25);
 927     sumMultHist->DrawCopy("same");
 928
 929   }
 930   TH1F* primMult = (TH1F*)fMultList[what]->FindObject(Form("primary_%s",fAnalysisNames[what].Data()));
 931   TH1F* sumMult  = (TH1F*)fMultList[what]->FindObject(Form("dNdeta_%s",fAnalysisNames[what].Data()));
 932   sumMult->SetMarkerStyle(23);
 933   sumMult->SetMarkerColor(kRed);
 934   for(Int_t det = 1; det<=3;det++) {
 935     Int_t maxRing = (det == 1 ? 0 : 1);
 936     for(Int_t iring = 0; iring<=maxRing; iring++) {
 937       Char_t ringChar = (iring == 0 ? 'I' : 'O');
 938       TH1F* hRingMult= (TH1F*)fMultList[what]->FindObject(Form("hRingMult_FMD%d%c_%s",det,ringChar,fAnalysisNames[what].Data()));
 939       hRingMult->Add(primMult,-1);
 940       hRingMult->Divide(primMult);
 941     }
 942   }
 943
 944
 945   // End of SPD
 946
 947
 948   //TH1F*  hPrim = (TH1F*)fList->FindObject(fPrimEvents.Data());
 949   //  TFile testhit("/home/canute/ALICE/Simulations/TestOfAnalysis2/hitsdist.root");
 950   // TH1F* hHitCor = (TH1F*)testhit.Get("selectedHits");
 951   // hHitCor->Sumw2();
 952   // sumMult->Divide(hHitCor);
 953   TH1F*  hPrim = (TH1F*)fList->FindObject(fPrimdNdeta.Data());
 954   Float_t nPrimevents = hMCEvents->GetEntries();
 955   //std::cout<<hMCEvents->GetEntries()<<std::endl;
 956   Float_t xb = hPrim->GetNbinsX();
 957   Float_t xr = hPrim->GetXaxis()->GetXmax() - hPrim->GetXaxis()->GetXmin();
 958   //std::cout<<xb/xr<<std::endl;
 959   hPrim->Scale(xb / xr );
 960   if(nPrimevents > 0)
 961     hPrim->Scale(1/nPrimevents);
 962
 963
 964   //  primMult->Rebin(rebin);
 965   // primMult->Scale(1/rebin);
 966   // hPrim->Rebin(rebin);
 967   // hPrim->Scale(1/rebin);
 968   if(what != kHits && what != kHitsTrVtx)
 969     primMult = hPrim;
 970
 971   //  hReldif->Add(hPrim,-1);
 972   // hReldif->Divide(hPrim);
 973
 974
 975
 976   /////////////*******************New thing!!!
 977
 978   //TH3D* hist3d = (TH3D*)fDataObject->ProjectionXYZ("hist3d");
 979   // fDataObject->Sumw2();
 980   //TH2F* projeta = (TH2F*)fDataObject->Project3D("yx");
 981
 982   //  TProfile2D* projeta = fDataObject->Project3DProfile("yx");
 983   // projeta->SetSumw2();
 984   //TProfile* etaprofile = projeta->ProfileX("dNdeta_profile");
 985   //TProfile* etaprofile = projeta->ProfileX("dNdeta_profile");
 986   // TH1* etaprofile = projeta->ProjectionX("dNdeta_profile");
 987   //TProfile* etaprofile = fDataObject->ProfileX();
 988   /*
 989   TCanvas* ctest = new TCanvas("test","test",1200,600);
 990   ctest->Divide(3);
 991   ctest->cd(1);
 992   fDataObject->DrawCopy("box");
 993   ctest->cd(2);
 994   projeta->DrawCopy("colz");
 995   ctest->cd(3);
 996   etaprofile->DrawCopy();
 997   */
 998   //sumMult->Rebin(rebin);
 999   TH1F* hReldif = (TH1F*)sumMult->Clone("hReldif");
1000   if(rebin != 1) {
1001     RebinHistogram(sumMult,rebin);
1002     if(!realdata) {
1003       RebinHistogram(primMult,rebin);
1004       RebinHistogram(hReldif,rebin);
1005     }
1006   }
1007   hReldif->Add(primMult,-1);
1008   hReldif->Divide(primMult);
1009   TCanvas* c2 = new TCanvas("dN/deta","dN/deta",640,960);
1010   c2->Divide(1,2);//,0,0);
1011   c2->cd(2);
1012   gPad->Divide(1,2);
1013   gPad->cd(1);
1014   gPad->SetGridy();
1015   gPad->SetGridx();
1016   hReldif->SetTitle("");
1017   hReldif->GetYaxis()->SetTitle("Relative difference");
1018   hReldif->GetYaxis()->SetRangeUser(-0.2,0.2);
1019   hReldif->SetStats(0000);
1020   hReldif->GetXaxis()->SetTitle("#eta");
1021
1022   hReldif->DrawCopy();
1023   //SPD Rel dif
1024   TH1F* hReldifSPD = (TH1F*)hSPDcombi->Clone("hReldifSPD");
1025   if(rebin != 1) {
1026     RebinHistogram(hSPDcombi,rebin);
1027     if(!realdata) {
1028       RebinHistogram(hReldifSPD,rebin);
1029     }
1030   }
1031   hReldifSPD->Add(primMult,-1);
1032   hReldifSPD->Divide(primMult);
1033   hReldifSPD->DrawCopy("same");
1034
1035   TLine* zeroLine = new TLine(-4,0,6,0);
1036   zeroLine->SetLineWidth(2);
1037   zeroLine->SetLineColor(kBlack);
1038   zeroLine->Draw();
1039   hPrim->SetTitle("");
1040   hPrim->SetLabelFont(132,"xyz");
1041   hPrim->SetFillColor(kGray);
1042   primMult->SetFillColor(kBlue);
1043
1044
1045   c2->cd(1);
1046   gPad->SetGridy();
1047   gPad->SetGridx();
1048   hPrim->GetYaxis()->SetRangeUser(2,10);
1049   //hPrim->GetYaxis()->SetRangeUser(0,20);
1050   hPrim->SetStats(0000);
1051   hPrim->GetXaxis()->SetTitle("#eta");
1052   sumMult->SetTitle("");
1053   sumMult->SetStats(kFALSE);
1054   sumMult->SetXTitle("#eta");
1055   sumMult->SetYTitle("#frac{dN}{d#eta_{ch}}");
1056   // sumMult->GetYaxis()->SetRangeUser
1057   //sumMult->Scale(1/(Float_t)rebin);
1058   sumMult->DrawCopy("PE");
1059   //Syst errors
1060   TH1F* sumMultSystPos = (TH1F*)sumMult->Clone("systerrorsPos");
1061   TH1F* sumMultSystNeg = (TH1F*)sumMult->Clone("systerrorsNeg");
1062   for(Int_t jj = 1;jj <=sumMultSystPos->GetNbinsX();jj++) {
1063     if(sumMultSystPos->GetBinCenter(jj) < 0) {
1064       sumMultSystPos->SetBinError(jj,0);
1065       sumMultSystPos->SetBinContent(jj,0);
1066       continue;
1067     }
1068
1069     if(sumMultSystPos->GetBinContent(jj) > 0)
1070       sumMultSystPos->SetBinError(jj,TMath::Sqrt(TMath::Power(sumMultSystPos->GetBinError(jj),2)+TMath::Power(0.1*sumMultSystPos->GetBinContent(jj),2)));
1071   }
1072   for(Int_t jj = 1;jj <=sumMultSystNeg->GetNbinsX();jj++) {
1073     if(sumMultSystNeg->GetBinCenter(jj) > 0) {
1074       sumMultSystNeg->SetBinError(jj,0);
1075       sumMultSystNeg->SetBinContent(jj,0);
1076       continue;
1077     }
1078
1079     if(sumMultSystNeg->GetBinContent(jj) > 0)
1080       sumMultSystNeg->SetBinError(jj,TMath::Sqrt(TMath::Power(sumMultSystNeg->GetBinError(jj),2)+TMath::Power(0.1*sumMultSystNeg->GetBinContent(jj),2)));
1081   }
1082
1083
1084
1085   sumMultSystPos->SetFillColor(18);
1086   sumMultSystNeg->SetFillColor(18);
1087   //sumMultSystPos->DrawCopy("sameE5");
1088   //sumMultSystNeg->DrawCopy("sameE5");
1089   //End syst errors
1090
1091
1092
1093   if(what != kHits && what != kHitsTrVtx && hPrim->GetEntries())
1094     hPrim->DrawCopy("E3same");
1095   if(what == kHits || what == kHitsTrVtx)
1096     primMult->DrawCopy("E3same");
1097   hPrim->GetXaxis()->SetTitle("#eta");
1098
1099   TH1F* hMirror = new TH1F("hMirror","mirrored",sumMult->GetNbinsX(),sumMult->GetXaxis()->GetXmin(),sumMult->GetXaxis()->GetXmax());
1100
1101   for(Int_t i= (Int_t)0.5*sumMult->GetNbinsX(); i<=sumMult->GetNbinsX(); i++) {
1102     Float_t eta   = sumMult->GetBinCenter(i);
1103     Float_t value = sumMult->GetBinContent(i);
1104     Float_t error = sumMult->GetBinError(i);
1105     if(value > 0) {
1106       hMirror->SetBinContent(hMirror->FindBin(-1*eta),value);
1107       hMirror->SetBinError(hMirror->FindBin(-1*eta),error);
1108     }
1109   }
1110   TH1F* hReldifLR = (TH1F*)sumMult->Clone("hReldifLeftRight");
1111   hReldifLR->Add(hMirror,-1);
1112   hReldifLR->Divide(hMirror);
1113   c2->cd(2);
1114   gPad->cd(1);
1115   hReldifLR->GetYaxis()->SetRangeUser(-0.2,0.2);
1116
1117   hReldifLR->SetYTitle("Relative difference left-right");
1118   hReldifLR->SetLabelSize(2*hReldifLR->GetLabelSize("X"),"X");
1119   hReldifLR->SetLabelSize(2*hReldifLR->GetLabelSize("Y"),"Y");
1120   hReldifLR->SetTitleSize(2*hReldifLR->GetTitleSize("X"),"X");
1121   hReldifLR->SetTitleSize(2*hReldifLR->GetTitleSize("Y"),"Y");
1122   hReldifLR->SetTitleOffset(0.7*hReldifLR->GetTitleOffset("X"),"X");
1123   hReldifLR->SetTitleOffset(0.5*hReldifLR->GetTitleOffset("Y"),"Y");
1124
1125
1126   if(realdata)
1127     hReldifLR->DrawCopy();
1128   zeroLine->Draw();
1129   c2->cd(1);
1130   hMirror->SetMarkerStyle(25);
1131
1132   hPrim->SetMarkerStyle(8);
1133
1134   //graph2->Draw("P");
1135   TH1F* hPythiaMB = 0;
1136
1137   TFile* fpyt = 0;
1138
1139   if(!pythiafile.Contains("none"))
1140     fpyt = TFile::Open(pythiafile);
1141
1142   if(realdata ) {
1143     if(fpyt && !pythiafile.Contains("none")) {
1144       hPythiaMB = (TH1F*)fpyt->Get("hPythiaMB");
1145     }
1146     else
1147       std::cout<<"no pythia for this energy"<<std::endl;
1148   }
1149   if(what != kHits && what != kHitsTrVtx) {
1150     if(hPrim->GetEntries() != 0 && !realdata)
1151       hPrim->DrawCopy("PE3same");
1152     if(realdata) {
1153       //graph->Draw("PEsame");
1154       //graph2->Draw("PEsame");
1155       if(what != kMultNSD) {
1156         graphinel->Draw("PEsame");
1157         graphinel2->Draw("PEsame");
1158         p7742D1x1y1->Draw("PEsame");
1159       }
1160       else{
1161         graph->Draw("PEsame");
1162         graph2->Draw("PEsame");
1163         p7742D2x1y1->Draw("PEsame");
1164         p7743D8x1y1->Draw("PEsame");
1165
1166       }
1167
1168     }
1169
1170   }
1171
1172
1173   sumMult->DrawCopy("PEsame");
1174
1175
1176
1177   hSPDcombi->DrawCopy("same");
1178   if(realdata)
1179     hMirror->DrawCopy("PEsame");
1180
1181   TLegend* leg = new TLegend(0.3,0.2,0.7,0.45,"");
1182   if(!realdata) {
1183     if(what != kHits && what != kHitsTrVtx)
1184       leg->AddEntry(hPrim,"Primary","pf");
1185     else
1186       leg->AddEntry(primMult,"Hits","pf");
1187   }
1188
1189   if(what == kMult)
1190     leg->AddEntry(sumMult,"FMD INEL","p");
1191   else if(what == kMultTrVtx)
1192     leg->AddEntry(sumMult,"FMD TrVtx","p");
1193   else if(what == kMultNSD)
1194     leg->AddEntry(sumMult,"FMD NSD","p");
1195
1196
1197   if(realdata) {
1198   if(what == kMult)
1199     leg->AddEntry(hMirror,"Mirror FMD INEL","p");
1200   else if(what == kMultTrVtx)
1201     leg->AddEntry(hMirror,"Mirror FMD TrVtx","p");
1202   else if(what == kMultNSD)
1203     leg->AddEntry(hMirror,"Mirror FMD NSD","p");
1204
1205   if(what == kMultNSD) {
1206     leg->AddEntry(graph,"UA5 NSD","p");
1207     leg->AddEntry(graph2,"Mirror UA5 NSD","p"); }
1208   else if(what == kMult) {
1209     leg->AddEntry(graphinel,"UA5 INEL","p");
1210     leg->AddEntry(graphinel2,"Mirror UA5 INEL","p");
1211   }
1212     //leg->AddEntry(hPythiaMB,"Pythia MB","l");
1213   leg->AddEntry(hSPDcombi,"HHD SPD clusters","p");
1214   if(what == kMult)
1215     leg->AddEntry(p7742D1x1y1,"ALICE INEL published","p");
1216   if(what == kMultNSD) {
1217     leg->AddEntry(p7742D2x1y1,"ALICE NSD published","p");
1218     leg->AddEntry(p7743D8x1y1, "CMS NSD published","p");
1219   }
1220   }
1221   leg->Draw();
1222
1223   c2->cd(2);
1224   gPad->cd(2);
1225   TH1F* hRatioMultPythia = 0;
1226   TH1F* hRatioMultUA5 = 0;
1227   TH1F* hRatioMultUA5_rebin5 = new TH1F("hRatioMultUA5_rebin5","hRatioMultUA5_rebin5",sumMult->GetNbinsX(),sumMult->GetXaxis()->GetXmin(),sumMult->GetXaxis()->GetXmax());
1228
1229   Double_t xval=0, yval=0;
1230   Int_t ipos=0;
1231   for(Int_t bb =hRatioMultUA5_rebin5->FindBin(0); bb<=hRatioMultUA5_rebin5->GetNbinsX();bb++) {
1232
1233
1234
1235     xval=yval=0;
1236
1237     if(hRatioMultUA5_rebin5->GetBinCenter(bb) >= 0) {
1238       if(what == kMultNSD)
1239         graph->GetPoint(ipos,xval,yval);
1240       else
1241         graphinel->GetPoint(ipos,xval,yval);
1242
1243       if(yval>0) {
1244         hRatioMultUA5_rebin5->SetBinContent(bb,yval);
1245         hRatioMultUA5_rebin5->SetBinError(bb,graphinel->GetErrorY(ipos));
1246         if(hRatioMultUA5_rebin5->GetBinCenter(bb) < 4) {
1247           hRatioMultUA5_rebin5->SetBinContent(hRatioMultUA5_rebin5->FindBin(-1*hRatioMultUA5_rebin5->GetBinCenter(bb)),yval);
1248           hRatioMultUA5_rebin5->SetBinError(hRatioMultUA5_rebin5->FindBin(-1*hRatioMultUA5_rebin5->GetBinCenter(bb)),(what == kMultNSD ? graph->GetErrorY(ipos) : graph->GetErrorY(ipos)));
1249
1250         }
1251         ipos++;
1252       }
1253     }
1254   }
1255
1256
1257
1258   if(realdata) {
1259
1260     Float_t ratio = 1;
1261     if(hPythiaMB) {
1262       if(hPythiaMB->GetNbinsX() !=  sumMult->GetNbinsX())
1263         ratio = (Float_t)sumMult->GetNbinsX() / (Float_t)hPythiaMB->GetNbinsX();
1264     }
1265
1266     hRatioMultPythia = (TH1F*)sumMult->Clone("hRatioMultPythia");
1267     hRatioMultUA5    = (TH1F*)sumMult->Clone("hRatioMultUA5");
1268     if(ratio > 1) {
1269       hRatioMultPythia->Rebin((Int_t)ratio);
1270       hRatioMultPythia->Scale(1/ratio);
1271     }
1272     if(ratio < 1 && hPythiaMB) {
1273       hPythiaMB->Rebin((Int_t)(1/ratio));
1274       hPythiaMB->Scale(ratio);
1275     }
1276
1277     if(rebin !=5) {
1278       TGraphErrors* tmp1;
1279       if(what == kMultNSD)
1280       tmp1 = (TGraphErrors*)graph->Clone("UA5tmp");
1281     else
1282       tmp1 = (TGraphErrors*)graphinel->Clone("UA5tmp");
1283
1284     tmp1->Fit("pol8","Q0","Q0",1.5,6);
1285       TF1* hFit = tmp1->GetFunction("pol8");
1286       for(Int_t ii = hRatioMultUA5->GetNbinsX() / 2; ii<hRatioMultUA5->GetNbinsX();ii++) {
1287         if(hRatioMultUA5->GetBinContent(ii) > 0) {
1288           Float_t errorratio = hRatioMultUA5->GetBinError(ii) / hRatioMultUA5->GetBinContent(ii);
1289           hRatioMultUA5->SetBinContent(ii,
1290                                        hRatioMultUA5->GetBinContent(ii) /
1291                                        ((hRatioMultUA5->GetBinCenter(ii) > 4.8) ? hFit->Eval(hRatioMultUA5->GetBinCenter(ii-1)) : hFit->Eval(hRatioMultUA5->GetBinCenter(ii))));
1292           hRatioMultUA5->SetBinError(ii,hRatioMultUA5->GetBinContent(ii) * TMath::Sqrt(0.02*0.02 + TMath::Power(errorratio,2)));
1293
1294
1295         }
1296
1297       }
1298
1299       TGraphErrors* tmp2;
1300       if(what == kMultNSD)
1301         tmp2 = (TGraphErrors*)graph2->Clone("UA5tmp2");
1302       else
1303         tmp2 = (TGraphErrors*)graphinel2->Clone("UA5tmp2");
1304
1305       //tmp2->Fit("pol8","Q0","Q0",-3.7,-1.5);
1306       tmp2->Fit("pol7","Q0","Q0",-3.7,0);
1307       hFit = tmp2->GetFunction("pol7");
1308       for(Int_t ii = 0; ii<hRatioMultUA5->GetNbinsX()/2;ii++) {
1309         if(hRatioMultUA5->GetBinContent(ii) > 0) {
1310           Float_t errorratio = hRatioMultUA5->GetBinError(ii) / hRatioMultUA5->GetBinContent(ii);
1311           hRatioMultUA5->SetBinContent(ii,hRatioMultUA5->GetBinContent(ii) / hFit->Eval(hRatioMultUA5->GetBinCenter(ii))  );
1312           hRatioMultUA5->SetBinError(ii,hRatioMultUA5->GetBinContent(ii) * TMath::Sqrt(0.02*0.02 + TMath::Power(errorratio,2)));
1313
1314         }
1315
1316
1317         graphinel->GetHistogram()->SetStats(kFALSE);
1318         graphinel2->GetHistogram()->SetStats(kFALSE);
1319
1320       }
1321     }
1322     else hRatioMultUA5->Divide(hRatioMultUA5_rebin5);
1323       //}
1324
1325
1326
1327       gPad->SetGridx();
1328       gPad->SetGridy();
1329     TLine* oneLine = new TLine(-4,1,6,1);
1330     oneLine->SetLineWidth(2);
1331
1332     hRatioMultUA5->SetMarkerColor(kGreen);
1333     hRatioMultUA5->SetMarkerStyle(20);
1334     hRatioMultUA5->GetYaxis()->SetRangeUser(0.75,1.25);
1335
1336     hRatioMultUA5->SetYTitle("Ratio FMD/UA5");
1337     hRatioMultUA5->SetLabelSize(2*hRatioMultUA5->GetLabelSize("X"),"X");
1338     hRatioMultUA5->SetLabelSize(2*hRatioMultUA5->GetLabelSize("Y"),"Y");
1339     hRatioMultUA5->SetTitleSize(2*hRatioMultUA5->GetTitleSize("X"),"X");
1340     hRatioMultUA5->SetTitleSize(2*hRatioMultUA5->GetTitleSize("Y"),"Y");
1341     hRatioMultUA5->SetTitleOffset(0.7*hRatioMultUA5->GetTitleOffset("X"),"X");
1342     hRatioMultUA5->SetTitleOffset(0.5*hRatioMultUA5->GetTitleOffset("Y"),"Y");
1343
1344
1345
1346     hRatioMultUA5->DrawCopy();
1347     //hRatioMultPythia->DrawCopy("same");
1348     oneLine->Draw("same");
1349   }
1350
1351   TCanvas* cPaper = new TCanvas("FMD_SPD_UA5","FMD_SPD_UA5",800,600);
1352   cPaper->cd();
1353
1354   sumMult->DrawCopy();
1355   //sumMultSystPos->DrawCopy("sameE5");
1356   // sumMultSystNeg->DrawCopy("sameE5");
1357   //hTestdNdeta->DrawCopy("same");
1358   //hSPDcombi->DrawCopy("same");
1359   if(what == kMult)
1360     p7742D1x1y1->Draw("PEsame");
1361   if(what == kMultNSD) {
1362     p7742D2x1y1->Draw("PEsame");
1363     p7743D8x1y1->Draw("PEsame");
1364   }
1365
1366   TLegend* leg2 = new TLegend(0.3,0.2,0.7,0.45,"");
1367   if(what == kMult)
1368     leg2->AddEntry(sumMult,"FMD INEL","p");
1369   else if(what == kMultTrVtx)
1370     leg2->AddEntry(sumMult,"FMD TrVtx","p");
1371   else if(what == kMultNSD)
1372     leg2->AddEntry(sumMult,"FMD NSD","p");
1373
1374   if(realdata) {
1375
1376     /* if(what == kMult)
1377     leg2->AddEntry(hMirror,"Mirror FMD INEL","p");
1378   else if(what == kMultTrVtx)
1379     leg2->AddEntry(hMirror,"FMD TrVtx","p");
1380   else if(what == kMultNSD)
1381   leg2->AddEntry(hMirror,"FMD NSD","p");*/
1382
1383   if(what == kMultNSD) {
1384     leg2->AddEntry(graph,"UA5 NSD","p");
1385     leg2->AddEntry(graph2,"Mirror UA5 NSD","p"); }
1386   else if(what == kMult) {
1387     leg2->AddEntry(graphinel,"UA5 INEL","p");
1388     leg2->AddEntry(graphinel2,"Mirror UA5 INEL","p");
1389   }
1390     //leg2->AddEntry(hPythiaMB,"Pythia MB","l");
1391     //leg2->AddEntry(hSPDcombi,"HHD SPD clusters","p");
1392   if(what == kMult)
1393     leg2->AddEntry(p7742D1x1y1,"ALICE INEL published","p");
1394   if(what == kMultNSD) {
1395     leg2->AddEntry(p7742D2x1y1,"ALICE NSD published","p");
1396     leg2->AddEntry(p7743D8x1y1,"CMS NSD published","p");
1397   }
1398   }
1399   leg2->Draw();
1400
1401
1402
1403   if(what != kMultNSD) {
1404     graphinel->Draw("PEsame");
1405     graphinel2->Draw("PEsame");
1406   }
1407   else {
1408     graph->Draw("PEsame");
1409     graph2->Draw("PEsame");
1410   }
1411
1412   sumMult->DrawCopy("PEsame");
1413   c2->cd(1);
1414
1415   c2->Print("fmdana.png");
1416   TString species;
1417
1418   switch(what) {
1419   case kMult:
1420     species = "mult";
1421     break;
1422   case kMultTrVtx:
1423     species = "mult_TrVtx";
1424     break;
1425   case kHits:
1426     species = "hits";
1427     break;
1428   case kHitsTrVtx:
1429     species = "hits_TrVtx";
1430     break;
1431   case kMultNSD:
1432     species = "mult_NSD";
1433     break;
1434
1435   default:
1436     AliWarning("no valid Analysis entry!");
1437     break;
1438   }
1439   TFile fout(Form("fmd_dNdeta_%s.root",species.Data()),"recreate");
1440   if(hRatioMultPythia)
1441     hRatioMultPythia->Write();
1442   hPrim->Write();
1443   graph->Write("UA5nsd");
1444   graph2->Write("UA5mirrornsd");
1445   graphinel->Write("UA5inel");
1446   graphinel2->Write("UA5mirrorinel");
1447   sumMult->Write();
1448   hSPDcombi->Write();
1449   hReldif->Write();
1450   fMultList[what]->Write();
1451   c2->Write();
1452   //fDataObject->Write();
1453   fout.Close();
1454
1455 }
1456 //_____________________________________________________________________
1457 void AliFMDDndeta::CreateSharingEfficiency(const Char_t* filename, Bool_t store) {
1458   //Generate sharing efficiency
1459   Init(filename);
1460   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
1461   //pars->Init(kTRUE,AliFMDAnaParameters::kBackgroundCorrection);
1462   Int_t nVertexBins = pars->GetNvtxBins();
1463
1464   SetNames(kHits);
1465   // "nEvents";
1466   TH1I* hEvents          = (TH1I*)fList->FindObject(fEvents.Data());
1467   // "nMCEventsNoCuts"
1468   TH1I* hPrimEvents      = (TH1I*)fList->FindObject(fPrimEvents.Data());
1469
1470   SetNames(kHitsTrVtx);
1471   // "nEvents";
1472   TH1I* hEventsTrVtx     = (TH1I*)fList->FindObject(fEvents.Data());
1473   // "nMCEvents"
1474   TH1I* hPrimEventsTrVtx = (TH1I*)fList->FindObject(fPrimEvents.Data());
1475
1476   AliFMDAnaCalibSharingEfficiency* sharEff = new AliFMDAnaCalibSharingEfficiency();
1477
1478   for(Int_t det = 1; det<=3; det++) {
1479     Int_t maxRing = (det == 1 ? 0 : 1);
1480     for(Int_t ring = 0;ring<=maxRing;ring++) {
1481       Char_t ringChar = (ring == 0 ? 'I' : 'O');
1482       for(Int_t v=0; v< nVertexBins; v++) {
1483         // Get histograms like  hits_NoCuts_FMD%d%c_vtxbin%d_proj
1484         // - from AliFMDAnalysisTaskBackgroundCorrection.cxx
1485         TH1F* hHits = (TH1F*)fList->FindObject(GetAnalysisName(kHits,det, ringChar, v));
1486         // Get histograms like hits_FMD%d%c_vtxbin%d_proj
1487         // - from AliFMDAnalysisTaskBackgroundCorrection.cxx
1488         TH1F* hHitsTrVtx   = (TH1F*)fList->FindObject(GetAnalysisName(kHitsTrVtx,det, ringChar, v));
1489         // Get histograms like "hMCHits_nocuts_FMD%d%c_vtxbin%d"
1490         // - from AliFMDAnalysisTaskSharing.cxx
1491         TH1F* hMCHits      = (TH1F*)fList->FindObject(GetPrimName(kHits,det, ringChar, v));
1492         // Get histograms like "hMCHits_FMD%d%c_vtxbin%d"
1493         // - from AliFMDAnalysisTaskSharing.cxx
1494         TH1F* hMCHitsTrVtx = (TH1F*)fList->FindObject(GetPrimName(kHitsTrVtx,det, ringChar, v));
1495
1496         TH1F* hCorrection  = (TH1F*)hHits->Clone(Form("hCorrection_FMD%d%c_vtx%d",det, ringChar, v));
1497         TH1F* hCorrectionTrVtx  = (TH1F*)hHitsTrVtx->Clone(Form("hCorrection_FMD%d%c_vtx%d",det, ringChar, v));
1498         if(hEvents->GetBinContent(v+1))
1499           hCorrection->Scale(1./(Float_t)hEvents->GetBinContent(v+1));
1500         if(hPrimEvents->GetBinContent(v+1))
1501           hMCHits->Scale(1./(Float_t)hPrimEvents->GetBinContent(v+1));
1502         if(hEventsTrVtx->GetBinContent(v+1))
1503           hCorrectionTrVtx->Scale(1./(Float_t)hEventsTrVtx->GetBinContent(v+1));
1504         if(hPrimEventsTrVtx->GetBinContent(v+1))
1505           hMCHitsTrVtx->Scale(1./(Float_t)hPrimEventsTrVtx->GetBinContent(v+1));
1506
1507         hCorrection->Divide(hMCHits);
1508         hCorrectionTrVtx->Divide(hMCHitsTrVtx);
1509
1510         //sharEff->SetSharingEff(det,ringChar,v,hCorrection);
1511         sharEff->SetSharingEff(det,ringChar,v,hCorrection);
1512         sharEff->SetSharingEffTrVtx(det,ringChar,v,hCorrectionTrVtx);
1513
1514         //      std::cout<<hHits<<"  "<<hHitsTrVtx<<"   "<<hPrim<<"    "<<hPrimTrVtx<<std::endl;
1515
1516         }
1517       }
1518     }
1519
1520   if(store) {
1521     TFile fsharing(pars->GetPath(pars->GetSharingEffID()),"RECREATE");
1522     sharEff->Write(AliFMDAnaParameters::GetSharingEffID());
1523     fsharing.Close();
1524   }
1525
1526 }
1527 //_____________________________________________________________________
1528 void AliFMDDndeta::RebinHistogram(TH1F* hist, Int_t rebin) {
1529   //Clever rebinner
1530
1531   if(hist->GetNbinsX()%rebin) {
1532     std::cout<<"cannot rebin "<<hist->GetNbinsX()%rebin<< "is not zero"<<std::endl;
1533     return;
1534
1535   }
1536   TH1F* cloneHist = (TH1F*)hist->Clone();
1537   cloneHist->Rebin(rebin);
1538
1539   Int_t nBinsNew = hist->GetNbinsX() / rebin;
1540
1541   for(Int_t i =1;i<=nBinsNew; i++) {
1542     Float_t bincontent = 0;
1543     Float_t sumofw = 0;
1544     Float_t sumformean = 0;
1545     Int_t   nBins = 0;
1546     for(Int_t j=1; j<=rebin;j++) {
1547       if(hist->GetBinContent((i-1)*rebin + j) > 0) {
1548         bincontent = bincontent + hist->GetBinContent((i-1)*rebin + j);
1549         nBins++;
1550         Float_t weight = 1 / TMath::Power(hist->GetBinError((i-1)*rebin + j),2);
1551         sumofw = sumofw + weight;
1552         sumformean = sumformean + weight*hist->GetBinContent((i-1)*rebin + j);
1553
1554       }
1555     }
1556
1557     if(bincontent > 0 ) {
1558       cloneHist->SetBinContent(i,sumformean / sumofw);
1559       cloneHist->SetBinError(i,TMath::Sqrt(sumformean));//TMath::Sqrt(1/sumofw));
1560     }
1561   }
1562   hist->Rebin(rebin);
1563   for(Int_t i =1;i<=nBinsNew; i++) {
1564     hist->SetBinContent(i,cloneHist->GetBinContent(i));
1565   }
1566
1567
1568 }
1569 //_____________________________________________________________________
1570 //
1571 // EOF
1572 //
1573 // EOF