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 "TMath.h"
  16 #include "TLegend.h"
  17 #include "TLatex.h"
  18
  19 #define SMALLNUMBER 0.0001
  20
  21 ClassImp(AliFMDDndeta)
  22 //_____________________________________________________________________
  23
  24 AliFMDDndeta::AliFMDDndeta()
  25 : TObject(),
  26   fList(0),
  27   fMultList(),
  28   fNbinsToCut(2),
  29   fVtxCut1(-10),
  30   fVtxCut2(10),
  31   fIsInit(kFALSE),
  32   fIsGenerated(),
  33   fPrimEvents(),
  34   fEvents(),
  35   fPrimdNdeta()
  36 {
  37   fAnalysisNames[0] = "Hits";
  38   fAnalysisNames[1] = "HitsTrVtx";
  39   fAnalysisNames[2] = "dNdeta";
  40   fAnalysisNames[3] = "dNdetaTrVtx";
  41
  42   for(Int_t i=0; i<4;i++)
  43     fMultList[i] = new TList();
  44 }
  45 //_____________________________________________________________________
  46 void AliFMDDndeta::SetNames(Analysis what) {
  47   // Set names of histograms from analysis
  48
  49   switch(what) {
  50   case kHits :
  51     fPrimEvents.Form("nMCEventsNoCuts"); //was nMCEvents
  52     fEvents.Form("nEvents");
  53     fPrimdNdeta.Form("hMultvsEtaNoCuts");
  54     break;
  55   case kHitsTrVtx :
  56     fPrimEvents.Form("nMCEvents");
  57     fEvents.Form("nEvents");
  58     fPrimdNdeta.Form("hMultvsEtaNoCuts");
  59     break;
  60   case kMult :
  61     fPrimEvents.Form("nMCEventsNoCuts");
  62     fEvents.Form("nEvents");
  63     fPrimdNdeta.Form("hMultvsEtaNoCuts");
  64     break;
  65   case kMultTrVtx :
  66     fPrimEvents.Form("nMCEvents");
  67     fEvents.Form("nEvents");
  68     fPrimdNdeta.Form("hMultvsEta");
  69     break;
  70   default:
  71     break;
  72   }
  73 }
  74 //_____________________________________________________________________
  75 const char* AliFMDDndeta::GetAnalysisName(Analysis what, UShort_t det, Char_t ring, Int_t vtxbin) {
  76   //Get names of histograms
  77   const char* name = "";
  78
  79   switch(what) {
  80   case kHits :
  81     name = Form("hits_NoCuts_FMD%d%c_vtxbin%d_proj",det,ring,vtxbin); //NoCuts added
  82     break;
  83   case kHitsTrVtx :
  84     name = Form("hits_FMD%d%c_vtxbin%d_proj",det,ring,vtxbin);
  85     break;
  86   case kMult :
  87     name = Form("dNdeta_FMD%d%c_vtxbin%d_proj",det,ring,vtxbin);
  88     break;
  89   case kMultTrVtx :
  90     name = Form("dNdeta_FMD%d%c_TrVtx_vtxbin%d_proj",det,ring,vtxbin);
  91     break;
  92   default:
  93     break;
  94   }
  95   //std::cout<<name.Data()<<std::endl;
  96   return name;
  97 }
  98 //_____________________________________________________________________
  99 const char* AliFMDDndeta::GetPrimName(Analysis what, UShort_t det, Char_t ring, Int_t vtxbin) {
 100   //Get names of primaries
 101   const char* name = "";
 102
 103   switch(what) {
 104   case kHits :
 105     name = Form("hMCHits_nocuts_FMD%d%c_vtxbin%d",det,ring,vtxbin); //nocuts added
 106     break;
 107   case kHitsTrVtx :
 108     name = Form("hMCHits_FMD%d%c_vtxbin%d",det,ring,vtxbin);
 109     break;
 110   case kMult :
 111     name = Form("primmult_vtxbin%d",vtxbin);
 112     break;
 113   case kMultTrVtx :
 114     name = Form("primmult_NoCuts_vtxbin%d",vtxbin);
 115     break;
 116   default:
 117     break;
 118   }
 119   return name;
 120 }
 121 //_____________________________________________________________________
 122 void AliFMDDndeta::GenerateHits(Analysis what) {
 123   //Generate the hits distributions
 124   AliFMDAnaParameters* pars =  AliFMDAnaParameters::Instance();
 125   TH2F* hTmp         =  pars->GetBackgroundCorrection(1,'I',0);
 126   Int_t nVertexBins  =  pars->GetNvtxBins();
 127   Int_t nEtaBins     =  hTmp->GetNbinsX();
 128   Float_t etaMin     =  hTmp->GetXaxis()->GetXmin();
 129   Float_t etaMax     =  hTmp->GetXaxis()->GetXmax();
 130
 131   for(Int_t i = 0; i<nVertexBins;i++) {
 132     TH1F* hHits = new TH1F(Form("hMCHits_vtxbin%d_%s",i,fAnalysisNames[what]),Form("hHits_vtxbin%d_%s",i,fAnalysisNames[what]),nEtaBins,etaMin,etaMax);
 133     hHits->Sumw2();
 134     fMultList[what]->Add(hHits);
 135   }
 136
 137   for(Int_t det = 1; det<=3; det++) {
 138     Int_t maxRing = (det == 1 ? 0 : 1);
 139     for(Int_t ring = 0;ring<=maxRing;ring++) {
 140
 141       Char_t ringChar = (ring == 0 ? 'I' : 'O');
 142       for(Int_t v=0; v< nVertexBins; v++) {
 143         TH1F* hits = (TH1F*)fList->FindObject(GetPrimName(what,det,ringChar,v));
 144
 145         Int_t nNonZero = 0, nNonZeroInData = 0;
 146
 147         //removing edges
 148         for(Int_t i =1;i<=hits->GetNbinsX();i++) {
 149           if(hits->GetBinContent(i) !=0)
 150             nNonZero++;
 151         }
 152
 153         TH1F* sumMultHist = (TH1F*)fMultList[what]->FindObject(Form("hMCHits_vtxbin%d_%s",v,fAnalysisNames[what]));
 154         for(Int_t i =1;i<=sumMultHist->GetNbinsX();i++) {
 155           if(hits->GetBinContent(i) == 0 ) continue;
 156
 157           nNonZeroInData++;
 158
 159           if(nNonZeroInData <=fNbinsToCut || nNonZeroInData > (nNonZero - fNbinsToCut)) {
 160             continue;
 161           }
 162           if(sumMultHist->GetBinContent(i) < SMALLNUMBER && TMath::Abs(hits->GetBinContent(i)) > 0){
 163               sumMultHist->SetBinContent(i,hits->GetBinContent(i));
 164               sumMultHist->SetBinError(i,hits->GetBinError(i));
 165             }
 166           else {
 167
 168
 169             Float_t sumofw = (1/TMath::Power(hits->GetBinError(i),2)) + (1/TMath::Power(sumMultHist->GetBinError(i),2));
 170             Float_t wav =  (hits->GetBinContent(i)*(1/TMath::Power(hits->GetBinError(i),2)) + sumMultHist->GetBinContent(i)*(1/TMath::Power(sumMultHist->GetBinError(i),2)))/sumofw;
 171             Float_t error = 1/TMath::Sqrt(sumofw);
 172             sumMultHist->SetBinContent(i, wav);
 173             sumMultHist->SetBinError(i, error);
 174           }
 175         }
 176       }
 177     }
 178   }
 179 }
 180
 181 //_____________________________________________________________________
 182 void AliFMDDndeta::Init(const Char_t* filename) {
 183   //Initialize everything from file
 184   TFile* fin = TFile::Open(filename);
 185
 186   if(!fin) {
 187     AliWarning("No file - exiting !");
 188     return;
 189   }
 190   AliFMDAnaParameters* pars =  AliFMDAnaParameters::Instance();
 191   //pars->Init();
 192
 193   TList* list = (TList*)fin->Get(Form("%s/BackgroundCorrected",pars->GetDndetaAnalysisName()));
 194
 195   if(!list) //an old file ? Perhaps...
 196     list = (TList*)fin->Get("BackgroundCorrected");
 197
 198
 199
 200   Init(list);
 201
 202 }
 203 //_____________________________________________________________________
 204 void AliFMDDndeta::Init(TList* list) {
 205   //Initialize everything from TList
 206
 207   if(!list) {
 208     AliWarning("No list - exiting !");
 209     return;
 210   }
 211
 212   fList = (TList*)list->Clone("inputList");
 213
 214   fIsGenerated[kHits]      = kFALSE;
 215   fIsGenerated[kMult]      = kFALSE;
 216   fIsGenerated[kMultTrVtx] = kFALSE;
 217   fIsGenerated[kHitsTrVtx] = kFALSE;
 218
 219   fIsInit = kTRUE;
 220 }
 221 //_____________________________________________________________________
 222 void AliFMDDndeta::GenerateMult(Analysis what) {
 223   //Generate dNdeta
 224
 225   if(!fIsInit) {
 226     AliWarning("Not initialised - call Init to remedy");
 227     return;
 228   }
 229
 230   if(fIsGenerated[what]) {
 231     AliWarning("Already generated - have a look at the results!");
 232     return;
 233   }
 234   else fIsGenerated[what] = kTRUE;
 235
 236   SetNames(what);
 237
 238   if(what == kHits || what == kHitsTrVtx)
 239     GenerateHits(what);
 240
 241   TH1I* hMCEvents = (TH1I*)fList->FindObject(fPrimEvents.Data());
 242
 243   AliFMDAnaParameters* pars =  AliFMDAnaParameters::Instance();
 244
 245   TH2F* hTmp         =  pars->GetBackgroundCorrection(1,'I',0);
 246   Int_t nVertexBins  =  pars->GetNvtxBins();
 247   Int_t nEtaBins     =  hTmp->GetNbinsX();
 248   Float_t etaMin     =  hTmp->GetXaxis()->GetXmin();
 249   Float_t etaMax     =  hTmp->GetXaxis()->GetXmax();
 250
 251   for(Int_t i = 0; i<nVertexBins;i++) {
 252     TH1F* hMult = new TH1F(Form("hMult_vtxbin%d_%s",i,fAnalysisNames[what]),Form("hMult_vtxbin%d_%s",i,fAnalysisNames[what]),nEtaBins,etaMin,etaMax);
 253     hMult->Sumw2();
 254     fMultList[what]->Add(hMult);
 255   }
 256
 257   for(Int_t det = 1; det<=3;det++) {
 258     Int_t maxRing = (det == 1 ? 0 : 1);
 259     for(Int_t iring = 0; iring<=maxRing; iring++) {
 260       Char_t ringChar = (iring == 0 ? 'I' : 'O');
 261       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);
 262       fMultList[what]->Add(hRingMult);
 263     }
 264   }
 265   TH1I* hEvents = (TH1I*)fList->FindObject(fEvents.Data());
 266   TH1F*  hPrimVtx = 0;
 267
 268
 269
 270   for(Int_t det = 1; det<=3; det++) {
 271     Int_t maxRing = (det == 1 ? 0 : 1);
 272     for(Int_t ring = 0;ring<=maxRing;ring++) {
 273
 274       Char_t ringChar = (ring == 0 ? 'I' : 'O');
 275       for(Int_t v=0; v< nVertexBins; v++) {
 276         if(det == 1) {
 277           if(what == kHits || what == kHitsTrVtx)
 278             hPrimVtx = (TH1F*)fMultList[what]->FindObject(Form("hMCHits_vtxbin%d_%s",v,fAnalysisNames[what]));
 279           else
 280             hPrimVtx = (TH1F*)fList->FindObject(GetPrimName(what,det,ringChar,v));
 281
 282
 283           Float_t nPrimevents = hMCEvents->GetBinContent(v+1);
 284           Float_t xb = hPrimVtx->GetNbinsX();
 285           Float_t xr = hPrimVtx->GetXaxis()->GetXmax() - hPrimVtx->GetXaxis()->GetXmin();
 286           hPrimVtx->Scale(xb / xr );
 287           if(nPrimevents > 0)
 288             hPrimVtx->Scale(1/nPrimevents);
 289
 290         }
 291
 292         Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
 293         Float_t vtxZ1   = (delta*v) - pars->GetVtxCutZ();
 294         Float_t vtxZ2   = (delta*(v+1)) - pars->GetVtxCutZ();
 295
 296         if(vtxZ1< fVtxCut1 || vtxZ2 >fVtxCut2)
 297           continue;
 298         Float_t nEvents = (Float_t)hEvents->GetBinContent(v+1);
 299
 300         //TH1F* multhistproj = (TH1F*)fList->FindObject(Form("dNdeta_FMD%d%c_vtxbin%d_proj",det,ringChar,v));
 301         TH1F* multhistproj = (TH1F*)fList->FindObject(GetAnalysisName(what,det,ringChar,v));
 302         if(nEvents)
 303           multhistproj->Scale(1/nEvents);
 304         Float_t xnBins = multhistproj->GetNbinsX();
 305         Float_t xrange = multhistproj->GetXaxis()->GetXmax() - multhistproj->GetXaxis()->GetXmin();
 306         multhistproj->Scale(xnBins / xrange);
 307
 308
 309
 310         Int_t nNonZero = 0, nNonZeroInData = 0;
 311
 312         //removing edges
 313
 314         for(Int_t i =1;i<=multhistproj->GetNbinsX();i++) {
 315           if(multhistproj->GetBinContent(i) !=0)
 316             nNonZero++;
 317         }
 318         Int_t nBinsOld = fNbinsToCut;
 319         //      if(det == 2 && ringChar =='I') {
 320         //  fNbinsToCut = 1;
 321         //      }
 322         TH1F* hRingMult = (TH1F*)fMultList[what]->FindObject(Form("hRingMult_FMD%d%c_%s",det,ringChar,fAnalysisNames[what]));
 323
 324         for(Int_t i=1; i<=hRingMult->GetNbinsX(); i++) {
 325           if(multhistproj->GetBinContent(i)!=0) {
 326             nNonZeroInData++;
 327             if(nNonZeroInData <=fNbinsToCut || nNonZeroInData > (nNonZero - fNbinsToCut)) {
 328               continue;
 329             }
 330             Float_t oldweight = 0;
 331             Float_t oldwav    = 0;
 332             if(hRingMult->GetBinError(i)>0) {
 333               oldweight = 1/TMath::Power(hRingMult->GetBinError(i),2);
 334               oldwav    = oldweight*hRingMult->GetBinContent(i);
 335             }
 336             Float_t  weight = 1/TMath::Power(multhistproj->GetBinError(i),2);
 337             Float_t  wav    = oldwav + weight*multhistproj->GetBinContent(i);
 338             Float_t  sumofw = oldweight + weight;
 339             if(sumofw) {
 340               Float_t error = 1/TMath::Sqrt(sumofw);
 341
 342               hRingMult->SetBinContent(i,wav/sumofw);
 343               hRingMult->SetBinError(i,error);
 344
 345             }
 346           }
 347         }
 348         nNonZeroInData = 0;
 349
 350         TH1F* sumMultHist = (TH1F*)fMultList[what]->FindObject(Form("hMult_vtxbin%d_%s",v,fAnalysisNames[what]));
 351
 352         for(Int_t i =1;i<=sumMultHist->GetNbinsX();i++) {
 353           if(multhistproj->GetBinContent(i) != 0 ) {
 354             nNonZeroInData++;
 355
 356             if(nNonZeroInData <=fNbinsToCut || nNonZeroInData > (nNonZero - fNbinsToCut)) {
 357               continue;
 358             }
 359             if(sumMultHist->GetBinContent(i) < SMALLNUMBER && TMath::Abs(multhistproj->GetBinContent(i)) > 0){
 360               sumMultHist->SetBinContent(i,multhistproj->GetBinContent(i));
 361               sumMultHist->SetBinError(i,multhistproj->GetBinError(i));
 362             }
 363             else {
 364
 365
 366               Float_t sumofw = (1/TMath::Power(multhistproj->GetBinError(i),2)) + (1/TMath::Power(sumMultHist->GetBinError(i),2));
 367               Float_t wav =  (multhistproj->GetBinContent(i)*(1/TMath::Power(multhistproj->GetBinError(i),2)) + sumMultHist->GetBinContent(i)*(1/TMath::Power(sumMultHist->GetBinError(i),2)))/sumofw;
 368               Float_t error = 1/TMath::Sqrt(sumofw);
 369               sumMultHist->SetBinContent(i, wav);
 370               sumMultHist->SetBinError(i, error);
 371             }
 372
 373
 374           }
 375         }
 376         fNbinsToCut = nBinsOld;
 377       }
 378
 379     }
 380   }
 381
 382   TH1F* sumMult  = new TH1F(Form("dNdeta_%s",fAnalysisNames[what]),Form("dNdeta_%s",fAnalysisNames[what]),nEtaBins,etaMin,etaMax);
 383   sumMult->Sumw2();
 384   fMultList[what]->Add(sumMult);
 385   TH1F* primMult  = new TH1F(Form("primary_%s",fAnalysisNames[what]),Form("primary_%s",fAnalysisNames[what]),nEtaBins,etaMin,etaMax);
 386   primMult->Sumw2();
 387   fMultList[what]->Add(primMult);
 388
 389   Float_t wav  = 0, sumofw=0, weight = 0, primwav  = 0, primsumofw=0, primweight = 0;//, etaofbin = 0;
 390   for(Int_t i =1; i<=sumMult->GetNbinsX();i++) {
 391
 392     for(Int_t v = 0; v<nVertexBins;v++) {
 393       if(what == kHits || what == kHitsTrVtx)
 394         hPrimVtx = (TH1F*)fMultList[what]->FindObject(Form("hMCHits_vtxbin%d_%s",v,fAnalysisNames[what]));
 395       else
 396         hPrimVtx = (TH1F*)fList->FindObject(GetPrimName(what,0,0,v));
 397       TH1F* sumMultHist = (TH1F*)fMultList[what]->FindObject(Form("hMult_vtxbin%d_%s",v,fAnalysisNames[what]));
 398       //etaofbin += sumMultHist->GetBinCenter(i);
 399       // if( hPrimVtx->GetBinContent(i)!=0 && hPrimVtx->GetBinError(i)>0.001 && sumMultHist->GetBinContent(i)!=0) {
 400       if( TMath::Abs(hPrimVtx->GetBinContent(i)) > 0) {
 401         primweight = 1/TMath::Power(hPrimVtx->GetBinError(i),2);
 402         primwav    = primwav + primweight*hPrimVtx->GetBinContent(i);
 403         primsumofw = primsumofw + primweight;
 404       }
 405       if(TMath::Abs(sumMultHist->GetBinContent(i)) > 0) {
 406
 407         weight = 1/TMath::Power(sumMultHist->GetBinError(i),2);
 408         wav    = wav + weight*sumMultHist->GetBinContent(i);
 409         sumofw = sumofw + weight;
 410       }
 411
 412     }
 413     if( primsumofw !=0 ) {// && sumofw !=0) {
 414       Float_t primerror = 1/TMath::Sqrt(primsumofw);
 415
 416       primMult->SetBinContent(i,primwav/primsumofw);
 417       primMult->SetBinError(i,primerror);
 418     }
 419     else {
 420       primMult->SetBinContent(i,0);
 421       primMult->SetBinError(i,0);
 422     }
 423
 424     if(sumofw) {
 425
 426       Float_t error = 1/TMath::Sqrt(sumofw);
 427
 428       sumMult->SetBinContent(i,wav/sumofw);
 429       sumMult->SetBinError(i,error);
 430     }
 431
 432
 433     wav = 0;
 434     sumofw = 0;
 435     weight = 0;
 436     primwav = 0;
 437     primsumofw = 0;
 438     primweight = 0;
 439
 440   }
 441
 442 }
 443 //_____________________________________________________________________
 444 void AliFMDDndeta::DrawDndeta(Analysis what, Int_t rebin, Bool_t realdata) {
 445   //Draw everything
 446   if(!fIsInit) {
 447     AliWarning("Not initialised - call Init to remedy");
 448     return;
 449   }
 450
 451   if(!fIsGenerated[what]) {
 452     AliWarning("Not generated - generate first then draw!");
 453     return;
 454   }
 455
 456   SetNames(what);
 457
 458   AliFMDAnaParameters* pars =  AliFMDAnaParameters::Instance();
 459   TCanvas* c1 = new TCanvas("cvertexbins","Overlaps",1400,800);
 460   c1->Divide(5,2);
 461   Int_t nVertexBins  =  pars->GetNvtxBins();
 462
 463   TH1I* hEvents   = (TH1I*)fList->FindObject(fEvents.Data());
 464   TH1I* hMCEvents = (TH1I*)fList->FindObject(fPrimEvents.Data());
 465
 466   TH1F*  hPrimVtx = 0;
 467   for(Int_t det = 1; det<=3; det++) {
 468     Int_t maxRing = (det == 1 ? 0 : 1);
 469     for(Int_t ring = 0;ring<=maxRing;ring++) {
 470       for(Int_t v=0; v< nVertexBins; v++) {
 471         Char_t ringChar = (ring == 0 ? 'I' : 'O');
 472         Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
 473         Float_t vtxZ1   = (delta*v) - pars->GetVtxCutZ();
 474         Float_t vtxZ2   = (delta*(v+1)) - pars->GetVtxCutZ();
 475
 476         if(vtxZ1<fVtxCut1 || vtxZ2 >fVtxCut2)
 477           continue;
 478         Float_t nEvents = (Float_t)hEvents->GetBinContent(v+1);
 479
 480         TH1F* multhistproj = (TH1F*)fList->FindObject(GetAnalysisName(what,det,ringChar,v));
 481
 482         c1->cd(v+1);
 483
 484         if(det==1)  {
 485           multhistproj->SetMarkerStyle(20);
 486           multhistproj->SetMarkerColor(1);
 487         }
 488         if(det==2 && ringChar=='I') {
 489           multhistproj->SetMarkerStyle(21);
 490           multhistproj->SetMarkerColor(2);
 491         }
 492         if(det==2 && ringChar=='O') {
 493           multhistproj->SetMarkerStyle(3);
 494           multhistproj->SetMarkerColor(3);
 495         }
 496         if(det==3 && ringChar=='I') {
 497           multhistproj->SetMarkerStyle(22);
 498           multhistproj->SetMarkerColor(4);
 499         }
 500         if(det==3 && ringChar=='O') {
 501           multhistproj->SetMarkerStyle(23);
 502           multhistproj->SetMarkerColor(6);
 503         }
 504
 505
 506
 507         if(det == 1) {
 508           if(what == kHits || what == kHitsTrVtx)
 509             hPrimVtx = (TH1F*)fMultList[what]->FindObject(Form("hMCHits_vtxbin%d_%s",v,Form("primary_%s",fAnalysisNames[what])));
 510           else
 511             hPrimVtx = (TH1F*)fList->FindObject(GetPrimName(what,det,ringChar,v));
 512
 513           hPrimVtx->SetTitle("");
 514           TLatex* l = new TLatex(0.14,0.92,Form("Vtx range [%.1f, %.1f], %d events",vtxZ1,vtxZ2,(Int_t)nEvents));
 515           l->SetNDC(kTRUE);
 516           hPrimVtx->SetFillColor(kGray);
 517           hPrimVtx->SetLabelFont(132,"xyz");
 518           hPrimVtx->SetStats(0000);
 519           hPrimVtx->GetXaxis()->SetTitle("#eta");
 520           hPrimVtx->GetYaxis()->SetRangeUser(0,15);
 521           hPrimVtx->DrawCopy("E3");
 522           l->Draw();
 523           multhistproj->DrawCopy("same");
 524
 525         }
 526         else
 527           multhistproj->DrawCopy("same");
 528
 529
 530       }
 531     }
 532   }
 533
 534   for(Int_t v=0; v< nVertexBins; v++) {
 535     TH1F* sumMultHist = (TH1F*)fMultList[what]->FindObject(Form("hMult_vtxbin%d_%s",v,fAnalysisNames[what]));
 536     c1->cd(v+1);
 537     sumMultHist->SetMarkerStyle(25);
 538     sumMultHist->DrawCopy("same");
 539
 540   }
 541   TH1F* primMult = (TH1F*)fMultList[what]->FindObject(Form("primary_%s",fAnalysisNames[what]));
 542   TH1F* sumMult  = (TH1F*)fMultList[what]->FindObject(Form("dNdeta_%s",fAnalysisNames[what]));
 543   sumMult->SetMarkerStyle(23);
 544   sumMult->SetMarkerColor(kRed);
 545   for(Int_t det = 1; det<=3;det++) {
 546     Int_t maxRing = (det == 1 ? 0 : 1);
 547     for(Int_t iring = 0; iring<=maxRing; iring++) {
 548       Char_t ringChar = (iring == 0 ? 'I' : 'O');
 549       TH1F* hRingMult= (TH1F*)fMultList[what]->FindObject(Form("hRingMult_FMD%d%c_%s",det,ringChar,fAnalysisNames[what]));
 550       hRingMult->Add(primMult,-1);
 551       hRingMult->Divide(primMult);
 552     }
 553   }
 554
 555
 556   //TH1F*  hPrim = (TH1F*)fList->FindObject(fPrimEvents.Data());
 557   //  TFile testhit("/home/canute/ALICE/Simulations/TestOfAnalysis2/hitsdist.root");
 558   // TH1F* hHitCor = (TH1F*)testhit.Get("selectedHits");
 559   // hHitCor->Sumw2();
 560   // sumMult->Divide(hHitCor);
 561   TH1F*  hPrim = (TH1F*)fList->FindObject(fPrimdNdeta.Data());
 562   Float_t nPrimevents = hMCEvents->GetEntries();
 563   //std::cout<<hMCEvents->GetEntries()<<std::endl;
 564   Float_t xb = hPrim->GetNbinsX();
 565   Float_t xr = hPrim->GetXaxis()->GetXmax() - hPrim->GetXaxis()->GetXmin();
 566   //std::cout<<xb/xr<<std::endl;
 567   hPrim->Scale(xb / xr );
 568   if(nPrimevents > 0)
 569     hPrim->Scale(1/nPrimevents);
 570
 571   /*
 572   TFile testin("/home/canute/ALICE/FMDanalysis/GridAnalysis/fmd_dNdeta_hits.root","READ");
 573     // TFile testin("/home/canute/ALICE/FMDanalysis/productionData/fmd_dNdeta_hits.root","READ");
 574   TH1F* hcorrect = (TH1F*)testin.Get("hReldif");
 575   hcorrect->SetName("djarnis");
 576   std::cout<<hcorrect<<std::endl;
 577   for(Int_t bb = 1;bb<=hcorrect->GetNbinsX();bb++) {
 578     hcorrect->SetBinContent(bb,hcorrect->GetBinContent(bb)+1);
 579   }
 580
 581   sumMult->Divide(hcorrect);
 582   //delete hcorrect;
 583   */
 584
 585
 586   //  primMult->Rebin(rebin);
 587   // primMult->Scale(1/rebin);
 588   // hPrim->Rebin(rebin);
 589   // hPrim->Scale(1/rebin);
 590   if(what != kHits && what != kHitsTrVtx)
 591     primMult = hPrim;
 592
 593   //  hReldif->Add(hPrim,-1);
 594   // hReldif->Divide(hPrim);
 595
 596
 597
 598   //sumMult->Rebin(rebin);
 599   TH1F* hReldif = (TH1F*)sumMult->Clone("hReldif");
 600   if(rebin != 1) {
 601     RebinHistogram(sumMult,rebin);
 602     if(!realdata) {
 603       RebinHistogram(primMult,rebin);
 604       RebinHistogram(hReldif,rebin);
 605     }
 606   }
 607   hReldif->Add(primMult,-1);
 608   hReldif->Divide(primMult);
 609   TCanvas* c2 = new TCanvas("dN/deta","dN/deta",640,960);
 610   c2->Divide(1,2);//,0,0);
 611   c2->cd(2);
 612   gPad->Divide(1,2);
 613   gPad->cd(1);
 614   gPad->SetGridy();
 615   gPad->SetGridx();
 616   hReldif->SetTitle("");
 617   hReldif->GetYaxis()->SetTitle("Relative difference");
 618   hReldif->GetYaxis()->SetRangeUser(-0.2,0.2);
 619   hReldif->SetStats(0000);
 620   hReldif->GetXaxis()->SetTitle("#eta");
 621
 622   hReldif->DrawCopy();
 623   TLine* zeroLine = new TLine(-6,0,6,0);
 624   zeroLine->SetLineWidth(2);
 625   zeroLine->SetLineColor(kBlack);
 626   zeroLine->Draw();
 627   hPrim->SetTitle("");
 628   hPrim->SetLabelFont(132,"xyz");
 629   hPrim->SetFillColor(kGray);
 630   primMult->SetFillColor(kBlue);
 631
 632
 633   c2->cd(1);
 634   gPad->SetGridy();
 635   gPad->SetGridx();
 636   hPrim->GetYaxis()->SetRangeUser(2,10);
 637   //hPrim->GetYaxis()->SetRangeUser(0,20);
 638   hPrim->SetStats(0000);
 639   hPrim->GetXaxis()->SetTitle("#eta");
 640   sumMult->SetTitle("");
 641   sumMult->SetStats(kFALSE);
 642   // sumMult->GetYaxis()->SetRangeUser
 643   //sumMult->Scale(1/(Float_t)rebin);
 644   sumMult->DrawCopy("PE");
 645   if(what != kHits && what != kHitsTrVtx && hPrim->GetEntries())
 646     hPrim->DrawCopy("E3same");
 647   if(what == kHits || what == kHitsTrVtx)
 648     primMult->DrawCopy("E3same");
 649   hPrim->GetXaxis()->SetTitle("#eta");
 650
 651   TH1F* hMirror = new TH1F("hMirror","mirrored",sumMult->GetNbinsX(),sumMult->GetXaxis()->GetXmin(),sumMult->GetXaxis()->GetXmax());
 652
 653   for(Int_t i=0.5*sumMult->GetNbinsX(); i<=sumMult->GetNbinsX(); i++) {
 654     Float_t eta   = sumMult->GetBinCenter(i);
 655     Float_t value = sumMult->GetBinContent(i);
 656     Float_t error = sumMult->GetBinError(i);
 657     if(value > 0) {
 658       hMirror->SetBinContent(hMirror->FindBin(-1*eta),value);
 659       hMirror->SetBinError(hMirror->FindBin(-1*eta),error);
 660     }
 661   }
 662   TH1F* hReldifLR = (TH1F*)sumMult->Clone("hReldifLeftRight");
 663   hReldifLR->Add(hMirror,-1);
 664   hReldifLR->Divide(hMirror);
 665   c2->cd(2);
 666   gPad->cd(1);
 667   if(realdata)
 668     hReldifLR->DrawCopy("same");
 669   c2->cd(1);
 670   hMirror->SetMarkerStyle(25);
 671
 672   hPrim->SetMarkerStyle(8);
 673   //UA5 data NSD
 674   Float_t x[19] = {0.125,0.375,0.625,0.875,1.125,1.375,1.625,1.875,2.125,2.375,
 675                    2.625,2.875,3.125,3.375,3.625,3.875,4.125,4.375,4.625};
 676   Float_t x2[19] = {-0.125,-0.375,-0.625,-0.875,-1.125,-1.375,-1.625,-1.875,
 677                     -2.125,-2.375,-2.625,-2.875,-3.125,-3.375,-3.625,-3.875,
 678                     -4.125,-4.375,-4.625};
 679
 680
 681   Float_t y[19] = {3.48, 3.38, 3.52, 3.68, 3.71, 3.86, 3.76, 3.66, 3.72 ,3.69,
 682                    3.56, 3.41, 3.15, 3.09, 2.74, 2.73, 2.32, 1.99, 1.69};
 683
 684   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,
 685                     0.07,0.07,0.08,0.08,0.09,0.09,0.1,0.13};
 686
 687   TGraphErrors* graph = new TGraphErrors(19,x,y,NULL,ey);
 688   TGraphErrors* graph2 = new TGraphErrors(19,x2,y,NULL,ey);
 689   //UA5 data INEL
 690   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};
 691   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};
 692   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};
 693   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};
 694  TGraphErrors* graphinel = new TGraphErrors(19,xinel,yinel,NULL,eyinel);
 695  TGraphErrors* graphinel2 = new TGraphErrors(19,xinel2,yinel,NULL,eyinel);
 696
 697   graph->SetMarkerStyle(22);
 698   graph->SetMarkerColor(kBlue);
 699   graphinel->SetMarkerStyle(20);
 700   graphinel->SetMarkerColor(kGreen);
 701   graphinel2->SetMarkerStyle(24);
 702   graphinel2->SetMarkerColor(kGreen);
 703   //graph->Draw("P");
 704   graph2->SetMarkerStyle(26);
 705   graph2->SetMarkerColor(kBlue);
 706   //graph2->Draw("P");
 707   TH1F* hPythiaMB = 0;
 708   TFile fpyt("/home/canute/ALICE/FMDanalysis/FirstAnalysis/pythia_study/pythiahists.root","READ");
 709   if(realdata) {
 710
 711     hPythiaMB = (TH1F*)fpyt.Get("hPythiaMB");
 712     hPythiaMB->DrawCopy("same");
 713
 714     std::cout<<hPythiaMB<<std::endl;
 715   }
 716   if(what != kHits && what != kHitsTrVtx) {
 717     if(hPrim->GetEntries() != 0 && !realdata)
 718       hPrim->DrawCopy("PE3same");
 719     if(realdata) {
 720       //graph->Draw("PEsame");
 721       //graph2->Draw("PEsame");
 722       graphinel->Draw("PEsame");
 723       graphinel2->Draw("PEsame");
 724     }
 725
 726   }
 727   sumMult->DrawCopy("PEsame");
 728   if(realdata)
 729     hMirror->DrawCopy("PEsame");
 730   //  std::cout<<"FMD 1 "<<sumMult->Integral(sumMult->FindBin(3),sumMult->FindBin(5),"width")<<std::endl;
 731
 732   //c2->SaveAs("dNdeta_fmd.png");
 733   /*TFile* itsfile = TFile::Open(itsfilename);
 734   if(itsfile) {
 735     //itsfile->cd();
 736
 737     //TList* itslist = (TList*)itsfile->Get("ITS_analysis");
 738     //TH1F* itshist0 = (TH1F*)itslist->FindObject("dndeta_check_0");
 739     //TH1F* itshist1 = (TH1F*)itslist->FindObject("dndeta_check_1");
 740     //TH1F* itshist2 = (TH1F*)itslist->FindObject("dndeta_check_2");
 741     TH1F* itshist0 = (TH1F*)itsfile->Get("dndeta/dNdEta");
 742     TH1F* itshist1= (TH1F*)itsfile->Get("dndeta/dNdEta_1");
 743     TH1F* itshist2 = (TH1F*)itsfile->Get("dndeta/dNdEta_2");
 744     itshist0->SetMarkerStyle(27);
 745     itshist1->SetMarkerStyle(28);
 746     itshist2->SetMarkerStyle(30);
 747     itshist0->DrawCopy("same");
 748     itshist1->DrawCopy("same");
 749     itshist2->DrawCopy("same");
 750
 751   }
 752   */
 753
 754
 755   TLegend* leg = new TLegend(0.35,0.2,0.65,0.45,"");
 756   if(!realdata) {
 757     if(what != kHits && what != kHitsTrVtx)
 758       leg->AddEntry(hPrim,"Primary","pf");
 759     else
 760       leg->AddEntry(primMult,"Hits","pf");
 761   }
 762   //leg->AddEntry(primMult,"Analysed prim","f");
 763
 764
 765   leg->AddEntry(sumMult,"Analysis","p");
 766   if(realdata) {
 767     leg->AddEntry(hMirror,"Mirror Analysis","p");
 768     //leg->AddEntry(graph,"UA5 NSD","p");
 769     //leg->AddEntry(graph2,"Mirror UA5 NSD","p");
 770     leg->AddEntry(graphinel,"UA5 INEL","p");
 771     leg->AddEntry(graphinel2,"Mirror UA5 INEL","p");
 772     leg->AddEntry(hPythiaMB,"Pythia MB","l");
 773
 774   }
 775   leg->Draw();
 776
 777   c2->cd(2);
 778   gPad->cd(2);
 779   TH1F* hRatioMultPythia = 0;
 780   TH1F* hRatioMultUA5 = 0;
 781   if(realdata) {
 782
 783     Float_t ratio = 1;
 784     if(hPythiaMB->GetNbinsX() !=  sumMult->GetNbinsX())
 785       ratio = (Float_t)sumMult->GetNbinsX() / (Float_t)hPythiaMB->GetNbinsX();
 786     hRatioMultPythia = (TH1F*)sumMult->Clone("hRatioMultPythia");
 787     hRatioMultUA5    = (TH1F*)sumMult->Clone("hRatioMultUA5");
 788     if(ratio > 1) {
 789       hRatioMultPythia->Rebin(ratio);
 790       hRatioMultPythia->Scale(1/ratio);
 791     }
 792     if(ratio < 1) {
 793       hPythiaMB->Rebin(1/ratio);
 794       hPythiaMB->Scale(ratio);
 795     }
 796
 797     hRatioMultPythia->Divide(hPythiaMB);
 798
 799     for(Int_t j=1;j<=hRatioMultUA5->GetNbinsX(); j++) {
 800       Float_t data = hRatioMultUA5->GetBinContent(j);
 801       Float_t errordata = hRatioMultUA5->GetBinError(j);
 802       Float_t ua5  = 0;
 803       Float_t ua5error = 0;
 804
 805
 806       if(hRatioMultUA5->GetBinCenter(j) > 0) {
 807         Double_t* xv  = graphinel->GetX();
 808         Double_t* yv  = graphinel->GetY();
 809         Double_t* eyv = graphinel->GetEY();
 810
 811         for(Int_t kk =0; kk<graphinel->GetN();kk++) {
 812           if(xv[kk] < hRatioMultUA5->GetBinCenter(j) &&  xv[kk+1] > hRatioMultUA5->GetBinCenter(j)) {
 813             ua5 = yv[kk];
 814             ua5error = eyv[kk];
 815             }
 816         }
 817
 818           }
 819       else {
 820         Double_t* xv = graphinel2->GetX();
 821         Double_t* yv = graphinel2->GetY();
 822         Double_t* eyv = graphinel->GetEY();
 823
 824         for(Int_t kk =0; kk<graphinel2->GetN();kk++) {
 825           if(xv[kk+1] < hRatioMultUA5->GetBinCenter(j) &&  xv[kk] > hRatioMultUA5->GetBinCenter(j)) {
 826             ua5 = yv[kk];
 827             ua5error = eyv[kk];
 828           }
 829         }
 830
 831       }
 832
 833       Float_t ratio2 = 1;
 834       if(ua5> 0) {
 835         ratio2 =  data / ua5;
 836         Float_t errorratio = 0;
 837         if(ua5 != 0 && data !=0)
 838           errorratio = ratio2*TMath::Sqrt(TMath::Power(ua5error/ua5,2)+TMath::Power(errordata/data,2));
 839         hRatioMultUA5->SetBinContent(j,ratio2);
 840         hRatioMultUA5->SetBinError(j,errorratio);
 841       }
 842     }
 843
 844
 845
 846     gPad->SetGridx();
 847     gPad->SetGridy();
 848     TLine* oneLine = new TLine(-6,1,6,1);
 849     hRatioMultPythia->DrawCopy();
 850     hRatioMultUA5->SetMarkerColor(kGreen);
 851     hRatioMultUA5->SetMarkerStyle(20);
 852     hRatioMultUA5->DrawCopy("same");
 853     oneLine->Draw("same");
 854   }
 855   c2->cd(1);
 856   TString species;
 857
 858   switch(what) {
 859   case kMult:
 860     species.Form("mult");
 861     break;
 862   case kMultTrVtx:
 863     species.Form("mult_TrVtx");
 864     break;
 865   case kHits:
 866     species.Form("hits");
 867     break;
 868   case kHitsTrVtx:
 869     species.Form("hits");
 870     break;
 871   default:
 872     AliWarning("no valid Analysis entry!");
 873     break;
 874   }
 875   TFile fout(Form("fmd_dNdeta_%s.root",species.Data()),"recreate");
 876   if(hRatioMultPythia)
 877     hRatioMultPythia->Write();
 878   hPrim->Write();
 879   sumMult->Write();
 880   hReldif->Write();
 881   fMultList[what]->Write();
 882   c2->Write();
 883   fout.Close();
 884
 885 }
 886 //_____________________________________________________________________
 887 void AliFMDDndeta::CreateSharingEfficiency(const Char_t* filename, Bool_t store) {
 888   //Generate sharing efficiency
 889   Init(filename);
 890   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 891   //pars->Init(kTRUE,AliFMDAnaParameters::kBackgroundCorrection);
 892   Int_t nVertexBins = pars->GetNvtxBins();
 893
 894   SetNames(kHits);
 895   TH1I* hEvents          = (TH1I*)fList->FindObject(fEvents.Data());
 896   TH1I* hPrimEvents      = (TH1I*)fList->FindObject(fPrimEvents.Data());
 897
 898   SetNames(kHitsTrVtx);
 899   TH1I* hEventsTrVtx     = (TH1I*)fList->FindObject(fEvents.Data());
 900   TH1I* hPrimEventsTrVtx = (TH1I*)fList->FindObject(fPrimEvents.Data());
 901
 902   AliFMDAnaCalibSharingEfficiency* sharEff = new AliFMDAnaCalibSharingEfficiency();
 903
 904   for(Int_t det = 1; det<=3; det++) {
 905     Int_t maxRing = (det == 1 ? 0 : 1);
 906     for(Int_t ring = 0;ring<=maxRing;ring++) {
 907       Char_t ringChar = (ring == 0 ? 'I' : 'O');
 908       for(Int_t v=0; v< nVertexBins; v++) {
 909         TH1F* hHits = (TH1F*)fList->FindObject(GetAnalysisName(kHits,det, ringChar, v));
 910         TH1F* hHitsTrVtx   = (TH1F*)fList->FindObject(GetAnalysisName(kHitsTrVtx,det, ringChar, v));
 911         TH1F* hMCHits      = (TH1F*)fList->FindObject(GetPrimName(kHits,det, ringChar, v));
 912         TH1F* hMCHitsTrVtx = (TH1F*)fList->FindObject(GetPrimName(kHitsTrVtx,det, ringChar, v));
 913
 914         TH1F* hCorrection  = (TH1F*)hHits->Clone(Form("hCorrection_FMD%d%c_vtx%d"));
 915         TH1F* hCorrectionTrVtx  = (TH1F*)hHitsTrVtx->Clone(Form("hCorrection_FMD%d%c_vtx%d"));
 916
 917         hCorrection->Scale(1./(Float_t)hEvents->GetBinContent(v+1));
 918         hMCHits->Scale(1./(Float_t)hPrimEvents->GetBinContent(v+1));
 919         hCorrectionTrVtx->Scale(1./(Float_t)hEventsTrVtx->GetBinContent(v+1));
 920         hMCHitsTrVtx->Scale(1./(Float_t)hPrimEventsTrVtx->GetBinContent(v+1));
 921         hCorrection->Divide(hMCHits);
 922         hCorrectionTrVtx->Divide(hMCHitsTrVtx);
 923
 924         sharEff->SetSharingEff(det,ringChar,v,hCorrection);
 925         sharEff->SetSharingEffTrVtx(det,ringChar,v,hCorrectionTrVtx);
 926         //      std::cout<<hHits<<"  "<<hHitsTrVtx<<"   "<<hPrim<<"    "<<hPrimTrVtx<<std::endl;
 927
 928         }
 929       }
 930     }
 931
 932   if(store) {
 933     TFile fsharing(pars->GetPath(pars->GetSharingEffID()),"RECREATE");
 934     sharEff->Write(AliFMDAnaParameters::GetSharingEffID());
 935     fsharing.Close();
 936   }
 937
 938 }
 939 //_____________________________________________________________________
 940 void AliFMDDndeta::RebinHistogram(TH1F* hist, Int_t rebin) {
 941   //Clever rebinner
 942
 943   if(hist->GetNbinsX()%rebin) {
 944     std::cout<<"cannot rebin "<<hist->GetNbinsX()%rebin<< "is not zero"<<std::endl;
 945     return;
 946
 947   }
 948   TH1F* cloneHist = (TH1F*)hist->Clone();
 949   cloneHist->Rebin(rebin);
 950
 951   Int_t nBinsNew = hist->GetNbinsX() / rebin;
 952
 953   for(Int_t i =1;i<=nBinsNew; i++) {
 954     Float_t bincontent = 0;
 955     Float_t sumofw = 0;
 956     Float_t sumformean = 0;
 957     Int_t   nBins = 0;
 958     for(Int_t j=1; j<=rebin;j++) {
 959       if(hist->GetBinContent((i-1)*rebin + j) > 0) {
 960         bincontent = bincontent + hist->GetBinContent((i-1)*rebin + j);
 961         nBins++;
 962         Float_t weight = 1 / TMath::Power(hist->GetBinError((i-1)*rebin + j),2);
 963         sumofw = sumofw + weight;
 964         sumformean = sumformean + weight*hist->GetBinContent((i-1)*rebin + j);
 965
 966       }
 967     }
 968
 969     if(bincontent > 0 ) {
 970       cloneHist->SetBinContent(i,sumformean / sumofw);
 971       cloneHist->SetBinError(i,TMath::Sqrt(sumformean));//TMath::Sqrt(1/sumofw));
 972     }
 973   }
 974   hist->Rebin(rebin);
 975   for(Int_t i =1;i<=nBinsNew; i++) {
 976     hist->SetBinContent(i,cloneHist->GetBinContent(i));
 977   }
 978
 979
 980 }
 981 //_____________________________________________________________________
 982 //
 983 // EOF
 984 //