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  (0,7);
 598   //sumMult->Rebin(rebin);
 599   TH1F* hReldif = (TH1F*)sumMult->Clone("hReldif");
 600   if(rebin != 1) {
 601     RebinHistogram(sumMult,rebin);
 602     RebinHistogram(primMult,rebin);
 603     RebinHistogram(hReldif,rebin);
 604   }
 605   hReldif->Add(primMult,-1);
 606   hReldif->Divide(primMult);
 607   TCanvas* c2 = new TCanvas("dN/deta","dN/deta",640,960);
 608   c2->Divide(1,2);//,0,0);
 609   c2->cd(2);
 610   gPad->Divide(1,2);
 611   gPad->cd(1);
 612   gPad->SetGridy();
 613   gPad->SetGridx();
 614   hReldif->SetTitle("");
 615   hReldif->GetYaxis()->SetTitle("Relative difference");
 616   hReldif->GetYaxis()->SetRangeUser(-0.2,0.2);
 617   hReldif->SetStats(0000);
 618   hReldif->GetXaxis()->SetTitle("#eta");
 619   hReldif->DrawCopy();
 620   TLine* zeroLine = new TLine(-6,0,6,0);
 621   zeroLine->SetLineWidth(2);
 622   zeroLine->SetLineColor(kBlack);
 623   zeroLine->Draw();
 624   hPrim->SetTitle("");
 625   hPrim->SetLabelFont(132,"xyz");
 626   hPrim->SetFillColor(kGray);
 627   primMult->SetFillColor(kBlue);
 628
 629
 630   c2->cd(1);
 631   gPad->SetGridy();
 632   gPad->SetGridx();
 633   hPrim->GetYaxis()->SetRangeUser(2,10);
 634   //hPrim->GetYaxis()->SetRangeUser(0,20);
 635   hPrim->SetStats(0000);
 636   hPrim->GetXaxis()->SetTitle("#eta");
 637   sumMult->SetTitle("");
 638   sumMult->SetStats(kFALSE);
 639   // sumMult->GetYaxis()->SetRangeUser
 640   //sumMult->Scale(1/(Float_t)rebin);
 641   sumMult->DrawCopy("PE");
 642   if(what != kHits && what != kHitsTrVtx && hPrim->GetEntries())
 643     hPrim->DrawCopy("E3same");
 644   if(what == kHits || what == kHitsTrVtx)
 645     primMult->DrawCopy("E3same");
 646   hPrim->GetXaxis()->SetTitle("#eta");
 647
 648   TH1F* hMirror = new TH1F("hMirror","mirrored",sumMult->GetNbinsX(),sumMult->GetXaxis()->GetXmin(),sumMult->GetXaxis()->GetXmax());
 649
 650   for(Int_t i=0.5*sumMult->GetNbinsX(); i<=sumMult->GetNbinsX(); i++) {
 651     Float_t eta   = sumMult->GetBinCenter(i);
 652     Float_t value = sumMult->GetBinContent(i);
 653     Float_t error = sumMult->GetBinError(i);
 654     if(value > 0) {
 655       hMirror->SetBinContent(hMirror->FindBin(-1*eta),value);
 656       hMirror->SetBinError(hMirror->FindBin(-1*eta),error);
 657     }
 658   }
 659   TH1F* hReldifLR = (TH1F*)sumMult->Clone("hReldifLeftRight");
 660   hReldifLR->Add(hMirror,-1);
 661   hReldifLR->Divide(hMirror);
 662   c2->cd(2);
 663   gPad->cd(1);
 664   if(realdata)
 665     hReldifLR->DrawCopy("same");
 666   c2->cd(1);
 667   hMirror->SetMarkerStyle(25);
 668
 669   hPrim->SetMarkerStyle(8);
 670   //UA5 data NSD
 671   Float_t x[19] = {0.125,0.375,0.625,0.875,1.125,1.375,1.625,1.875,2.125,2.375,
 672                    2.625,2.875,3.125,3.375,3.625,3.875,4.125,4.375,4.625};
 673   Float_t x2[19] = {-0.125,-0.375,-0.625,-0.875,-1.125,-1.375,-1.625,-1.875,
 674                     -2.125,-2.375,-2.625,-2.875,-3.125,-3.375,-3.625,-3.875,
 675                     -4.125,-4.375,-4.625};
 676
 677
 678   Float_t y[19] = {3.48, 3.38, 3.52, 3.68, 3.71, 3.86, 3.76, 3.66, 3.72 ,3.69,
 679                    3.56, 3.41, 3.15, 3.09, 2.74, 2.73, 2.32, 1.99, 1.69};
 680
 681   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,
 682                     0.07,0.07,0.08,0.08,0.09,0.09,0.1,0.13};
 683
 684   TGraphErrors* graph = new TGraphErrors(19,x,y,NULL,ey);
 685   TGraphErrors* graph2 = new TGraphErrors(19,x2,y,NULL,ey);
 686   //UA5 data INEL
 687   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};
 688   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};
 689   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};
 690   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};
 691  TGraphErrors* graphinel = new TGraphErrors(19,xinel,yinel,NULL,eyinel);
 692  TGraphErrors* graphinel2 = new TGraphErrors(19,xinel2,yinel,NULL,eyinel);
 693
 694   graph->SetMarkerStyle(22);
 695   graph->SetMarkerColor(kBlue);
 696   graphinel->SetMarkerStyle(20);
 697   graphinel->SetMarkerColor(kGreen);
 698   graphinel2->SetMarkerStyle(24);
 699   graphinel2->SetMarkerColor(kGreen);
 700   //graph->Draw("P");
 701   graph2->SetMarkerStyle(26);
 702   graph2->SetMarkerColor(kBlue);
 703   //graph2->Draw("P");
 704   TH1F* hPythiaMB = 0;
 705   TFile fpyt("/home/canute/ALICE/FMDanalysis/FirstAnalysis/pythia_study/pythiahists.root","READ");
 706   if(realdata) {
 707
 708     hPythiaMB = (TH1F*)fpyt.Get("hPythiaMB");
 709     hPythiaMB->DrawCopy("same");
 710
 711     std::cout<<hPythiaMB<<std::endl;
 712   }
 713   if(what != kHits && what != kHitsTrVtx) {
 714     if(hPrim->GetEntries() != 0 && !realdata)
 715       hPrim->DrawCopy("PE3same");
 716     if(realdata) {
 717       //graph->Draw("PEsame");
 718       //graph2->Draw("PEsame");
 719       graphinel->Draw("PEsame");
 720       graphinel2->Draw("PEsame");
 721     }
 722
 723   }
 724   sumMult->DrawCopy("PEsame");
 725   if(realdata)
 726     hMirror->DrawCopy("PEsame");
 727   //  std::cout<<"FMD 1 "<<sumMult->Integral(sumMult->FindBin(3),sumMult->FindBin(5),"width")<<std::endl;
 728
 729   //c2->SaveAs("dNdeta_fmd.png");
 730   /*TFile* itsfile = TFile::Open(itsfilename);
 731   if(itsfile) {
 732     //itsfile->cd();
 733
 734     //TList* itslist = (TList*)itsfile->Get("ITS_analysis");
 735     //TH1F* itshist0 = (TH1F*)itslist->FindObject("dndeta_check_0");
 736     //TH1F* itshist1 = (TH1F*)itslist->FindObject("dndeta_check_1");
 737     //TH1F* itshist2 = (TH1F*)itslist->FindObject("dndeta_check_2");
 738     TH1F* itshist0 = (TH1F*)itsfile->Get("dndeta/dNdEta");
 739     TH1F* itshist1= (TH1F*)itsfile->Get("dndeta/dNdEta_1");
 740     TH1F* itshist2 = (TH1F*)itsfile->Get("dndeta/dNdEta_2");
 741     itshist0->SetMarkerStyle(27);
 742     itshist1->SetMarkerStyle(28);
 743     itshist2->SetMarkerStyle(30);
 744     itshist0->DrawCopy("same");
 745     itshist1->DrawCopy("same");
 746     itshist2->DrawCopy("same");
 747
 748   }
 749   */
 750
 751
 752   TLegend* leg = new TLegend(0.35,0.2,0.65,0.45,"");
 753   if(!realdata) {
 754     if(what != kHits && what != kHitsTrVtx)
 755       leg->AddEntry(hPrim,"Primary","pf");
 756     else
 757       leg->AddEntry(primMult,"Hits","pf");
 758   }
 759   //leg->AddEntry(primMult,"Analysed prim","f");
 760
 761
 762   leg->AddEntry(sumMult,"Analysis","p");
 763   if(realdata) {
 764     leg->AddEntry(hMirror,"Mirror Analysis","p");
 765     //leg->AddEntry(graph,"UA5 NSD","p");
 766     //leg->AddEntry(graph2,"Mirror UA5 NSD","p");
 767     leg->AddEntry(graphinel,"UA5 INEL","p");
 768     leg->AddEntry(graphinel2,"Mirror UA5 INEL","p");
 769     leg->AddEntry(hPythiaMB,"Pythia MB","l");
 770
 771   }
 772   leg->Draw();
 773
 774   c2->cd(2);
 775   gPad->cd(2);
 776   TH1F* hRatioMultPythia = 0;
 777   TH1F* hRatioMultUA5 = 0;
 778   if(realdata) {
 779
 780     Float_t ratio = 1;
 781     if(hPythiaMB->GetNbinsX() !=  sumMult->GetNbinsX())
 782       ratio = (Float_t)sumMult->GetNbinsX() / (Float_t)hPythiaMB->GetNbinsX();
 783     hRatioMultPythia = (TH1F*)sumMult->Clone("hRatioMultPythia");
 784     hRatioMultUA5    = (TH1F*)sumMult->Clone("hRatioMultUA5");
 785     if(ratio > 1) {
 786       hRatioMultPythia->Rebin(ratio);
 787       hRatioMultPythia->Scale(1/ratio);
 788     }
 789     if(ratio < 1) {
 790       hPythiaMB->Rebin(1/ratio);
 791       hPythiaMB->Scale(ratio);
 792     }
 793
 794     hRatioMultPythia->Divide(hPythiaMB);
 795
 796     for(Int_t j=1;j<=hRatioMultUA5->GetNbinsX(); j++) {
 797       Float_t data = hRatioMultUA5->GetBinContent(j);
 798       Float_t errordata = hRatioMultUA5->GetBinError(j);
 799       Float_t ua5  = 0;
 800       Float_t ua5error = 0;
 801
 802
 803       if(hRatioMultUA5->GetBinCenter(j) > 0) {
 804         Double_t* xv  = graphinel->GetX();
 805         Double_t* yv  = graphinel->GetY();
 806         Double_t* eyv = graphinel->GetEY();
 807
 808         for(Int_t kk =0; kk<graphinel->GetN();kk++) {
 809           if(xv[kk] < hRatioMultUA5->GetBinCenter(j) &&  xv[kk+1] > hRatioMultUA5->GetBinCenter(j)) {
 810             ua5 = yv[kk];
 811             ua5error = eyv[kk];
 812             }
 813         }
 814
 815           }
 816       else {
 817         Double_t* xv = graphinel2->GetX();
 818         Double_t* yv = graphinel2->GetY();
 819         Double_t* eyv = graphinel->GetEY();
 820
 821         for(Int_t kk =0; kk<graphinel2->GetN();kk++) {
 822           if(xv[kk+1] < hRatioMultUA5->GetBinCenter(j) &&  xv[kk] > hRatioMultUA5->GetBinCenter(j)) {
 823             ua5 = yv[kk];
 824             ua5error = eyv[kk];
 825           }
 826         }
 827
 828       }
 829
 830       Float_t ratio2 = 1;
 831       if(ua5> 0) {
 832         ratio2 =  data / ua5;
 833         Float_t errorratio = 0;
 834         if(ua5 != 0 && data !=0)
 835           errorratio = ratio2*TMath::Sqrt(TMath::Power(ua5error/ua5,2)+TMath::Power(errordata/data,2));
 836         hRatioMultUA5->SetBinContent(j,ratio2);
 837         hRatioMultUA5->SetBinError(j,errorratio);
 838       }
 839     }
 840
 841
 842
 843     gPad->SetGridx();
 844     gPad->SetGridy();
 845     TLine* oneLine = new TLine(-6,1,6,1);
 846     hRatioMultPythia->DrawCopy();
 847     hRatioMultUA5->SetMarkerColor(kGreen);
 848     hRatioMultUA5->SetMarkerStyle(20);
 849     hRatioMultUA5->DrawCopy("same");
 850     oneLine->Draw("same");
 851   }
 852   c2->cd(1);
 853   TString species;
 854
 855   switch(what) {
 856   case kMult:
 857     species.Form("mult");
 858     break;
 859   case kMultTrVtx:
 860     species.Form("mult_TrVtx");
 861     break;
 862   case kHits:
 863     species.Form("hits");
 864     break;
 865   case kHitsTrVtx:
 866     species.Form("hits");
 867     break;
 868   default:
 869     AliWarning("no valid Analysis entry!");
 870     break;
 871   }
 872   TFile fout(Form("fmd_dNdeta_%s.root",species.Data()),"recreate");
 873   if(hRatioMultPythia)
 874     hRatioMultPythia->Write();
 875   hPrim->Write();
 876   sumMult->Write();
 877   hReldif->Write();
 878   fMultList[what]->Write();
 879   c2->Write();
 880   fout.Close();
 881
 882 }
 883 //_____________________________________________________________________
 884 void AliFMDDndeta::CreateSharingEfficiency(const Char_t* filename, Bool_t store) {
 885   //Generate sharing efficiency
 886   Init(filename);
 887   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 888   //pars->Init(kTRUE,AliFMDAnaParameters::kBackgroundCorrection);
 889   Int_t nVertexBins = pars->GetNvtxBins();
 890
 891   SetNames(kHits);
 892   TH1I* hEvents          = (TH1I*)fList->FindObject(fEvents.Data());
 893   TH1I* hPrimEvents      = (TH1I*)fList->FindObject(fPrimEvents.Data());
 894
 895   SetNames(kHitsTrVtx);
 896   TH1I* hEventsTrVtx     = (TH1I*)fList->FindObject(fEvents.Data());
 897   TH1I* hPrimEventsTrVtx = (TH1I*)fList->FindObject(fPrimEvents.Data());
 898
 899   AliFMDAnaCalibSharingEfficiency* sharEff = new AliFMDAnaCalibSharingEfficiency();
 900
 901   for(Int_t det = 1; det<=3; det++) {
 902     Int_t maxRing = (det == 1 ? 0 : 1);
 903     for(Int_t ring = 0;ring<=maxRing;ring++) {
 904       Char_t ringChar = (ring == 0 ? 'I' : 'O');
 905       for(Int_t v=0; v< nVertexBins; v++) {
 906         TH1F* hHits = (TH1F*)fList->FindObject(GetAnalysisName(kHits,det, ringChar, v));
 907         TH1F* hHitsTrVtx   = (TH1F*)fList->FindObject(GetAnalysisName(kHitsTrVtx,det, ringChar, v));
 908         TH1F* hMCHits      = (TH1F*)fList->FindObject(GetPrimName(kHits,det, ringChar, v));
 909         TH1F* hMCHitsTrVtx = (TH1F*)fList->FindObject(GetPrimName(kHitsTrVtx,det, ringChar, v));
 910
 911         TH1F* hCorrection  = (TH1F*)hHits->Clone(Form("hCorrection_FMD%d%c_vtx%d"));
 912         TH1F* hCorrectionTrVtx  = (TH1F*)hHitsTrVtx->Clone(Form("hCorrection_FMD%d%c_vtx%d"));
 913
 914         hCorrection->Scale(1./(Float_t)hEvents->GetBinContent(v+1));
 915         hMCHits->Scale(1./(Float_t)hPrimEvents->GetBinContent(v+1));
 916         hCorrectionTrVtx->Scale(1./(Float_t)hEventsTrVtx->GetBinContent(v+1));
 917         hMCHitsTrVtx->Scale(1./(Float_t)hPrimEventsTrVtx->GetBinContent(v+1));
 918         hCorrection->Divide(hMCHits);
 919         hCorrectionTrVtx->Divide(hMCHitsTrVtx);
 920
 921         sharEff->SetSharingEff(det,ringChar,v,hCorrection);
 922         sharEff->SetSharingEffTrVtx(det,ringChar,v,hCorrectionTrVtx);
 923         //      std::cout<<hHits<<"  "<<hHitsTrVtx<<"   "<<hPrim<<"    "<<hPrimTrVtx<<std::endl;
 924
 925         }
 926       }
 927     }
 928
 929   if(store) {
 930     TFile fsharing(pars->GetPath(pars->GetSharingEffID()),"RECREATE");
 931     sharEff->Write(AliFMDAnaParameters::GetSharingEffID());
 932     fsharing.Close();
 933   }
 934
 935 }
 936 //_____________________________________________________________________
 937 void AliFMDDndeta::RebinHistogram(TH1F* hist, Int_t rebin) {
 938   //Clever rebinner
 939
 940   if(hist->GetNbinsX()%rebin) {
 941     std::cout<<"cannot rebin "<<hist->GetNbinsX()%rebin<< "is not zero"<<std::endl;
 942     return;
 943
 944   }
 945   TH1F* cloneHist = (TH1F*)hist->Clone();
 946   cloneHist->Rebin(rebin);
 947
 948   Int_t nBinsNew = hist->GetNbinsX() / rebin;
 949
 950   for(Int_t i =1;i<=nBinsNew; i++) {
 951     Float_t bincontent = 0;
 952     Float_t sumofw = 0;
 953     Float_t sumformean = 0;
 954     Int_t   nBins = 0;
 955     for(Int_t j=1; j<=rebin;j++) {
 956       if(hist->GetBinContent((i-1)*rebin + j) > 0) {
 957         bincontent = bincontent + hist->GetBinContent((i-1)*rebin + j);
 958         nBins++;
 959         Float_t weight = 1 / TMath::Power(hist->GetBinError((i-1)*rebin + j),2);
 960         sumofw = sumofw + weight;
 961         sumformean = sumformean + weight*hist->GetBinContent((i-1)*rebin + j);
 962
 963       }
 964     }
 965
 966     if(bincontent > 0 ) {
 967       cloneHist->SetBinContent(i,sumformean / sumofw);
 968       cloneHist->SetBinError(i,TMath::Sqrt(sumformean));//TMath::Sqrt(1/sumofw));
 969     }
 970   }
 971   hist->Rebin(rebin);
 972   for(Int_t i =1;i<=nBinsNew; i++) {
 973     hist->SetBinContent(i,cloneHist->GetBinContent(i));
 974   }
 975
 976
 977 }
 978 //_____________________________________________________________________
 979 //
 980 // EOF
 981 //