PWG2/FORWARD/analysis/AliFMDDndeta.cxx

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