FMD/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 "TPad.h"
  13 #include "iostream"
  14 #include "TMath.h"
  15 #include "TLegend.h"
  16 #include "TLatex.h"
  17 ClassImp(AliFMDDndeta)
  18 //_____________________________________________________________________
  19
  20 AliFMDDndeta::AliFMDDndeta()
  21 : TObject(),
  22   fList(0),
  23   fMultList(),
  24   fNbinsToCut(2),
  25   fVtxCut(10),
  26   fIsInit(kFALSE),
  27   fIsGenerated(),
  28   fPrimEvents()
  29 {
  30
  31 }
  32 //_____________________________________________________________________
  33 void AliFMDDndeta::SetNames(Analysis what) {
  34
  35   switch(what) {
  36   case kHits :
  37     fPrimEvents.Form("nMCEventsNoCuts"); //was nMCEvents
  38     fEvents.Form("nEvents");
  39     fPrimdNdeta.Form("hMultvsEtaNoCuts");
  40     break;
  41   case kHitsTrVtx :
  42     fPrimEvents.Form("nMCEvents");
  43     fEvents.Form("nEvents");
  44     fPrimdNdeta.Form("hMultvsEtaNoCuts");
  45     break;
  46   case kMult :
  47     fPrimEvents.Form("nMCEventsNoCuts");
  48     fEvents.Form("nEvents");
  49     fPrimdNdeta.Form("hMultvsEtaNoCuts");
  50     break;
  51   case kMultTrVtx :
  52     fPrimEvents.Form("nMCEvents");
  53     fEvents.Form("nEvents");
  54     fPrimdNdeta.Form("hMultvsEta");
  55     break;
  56   default:
  57     break;
  58   }
  59 }
  60 //_____________________________________________________________________
  61 const char* AliFMDDndeta::GetAnalysisName(Analysis what, UShort_t det, Char_t ring, Int_t vtxbin) {
  62
  63   const char* name;
  64
  65   switch(what) {
  66   case kHits :
  67     name = Form("hits_NoCuts_FMD%d%c_vtxbin%d_proj",det,ring,vtxbin); //NoCuts added
  68     break;
  69   case kHitsTrVtx :
  70     name = Form("hits_FMD%d%c_vtxbin%d_proj",det,ring,vtxbin);
  71     break;
  72   case kMult :
  73     name = Form("dNdeta_FMD%d%c_vtxbin%d_proj",det,ring,vtxbin);
  74     break;
  75   case kMultTrVtx :
  76     name = Form("dNdeta_FMD%d%c_TrVtx_vtxbin%d_proj",det,ring,vtxbin);
  77     break;
  78   default:
  79     break;
  80   }
  81   //std::cout<<name.Data()<<std::endl;
  82   return name;
  83 }
  84 //_____________________________________________________________________
  85 const char* AliFMDDndeta::GetPrimName(Analysis what, UShort_t det, Char_t ring, Int_t vtxbin) {
  86
  87   const char* name;
  88
  89   switch(what) {
  90   case kHits :
  91     name = Form("hMCHits_nocuts_FMD%d%c_vtxbin%d",det,ring,vtxbin); //nocuts added
  92     break;
  93   case kHitsTrVtx :
  94     name = Form("hMCHits_FMD%d%c_vtxbin%d",det,ring,vtxbin);
  95     break;
  96   case kMult :
  97     name = Form("primmult_vtxbin%d",vtxbin);
  98     break;
  99   case kMultTrVtx :
 100     name = Form("primmult_NoCuts_vtxbin%d",vtxbin);
 101     break;
 102   default:
 103     break;
 104   }
 105   return name;
 106 }
 107 //_____________________________________________________________________
 108 void AliFMDDndeta::GenerateHits() {
 109
 110   AliFMDAnaParameters* pars =  AliFMDAnaParameters::Instance();
 111   TH2F* hTmp         =  pars->GetBackgroundCorrection(1,'I',0);
 112   Int_t nVertexBins  =  pars->GetNvtxBins();
 113   Int_t nEtaBins     =  hTmp->GetNbinsX();
 114   Float_t etaMin     =  hTmp->GetXaxis()->GetXmin();
 115   Float_t etaMax     =  hTmp->GetXaxis()->GetXmax();
 116
 117   for(Int_t i = 0; i<nVertexBins;i++) {
 118     TH1F* hHits = new TH1F(Form("hMCHits_vtxbin%d",i),Form("hHits_vtxbin%d",i),nEtaBins,etaMin,etaMax);
 119     hHits->Sumw2();
 120     fMultList.Add(hHits);
 121   }
 122
 123   for(Int_t det = 1; det<=3; det++) {
 124     Int_t maxRing = (det == 1 ? 0 : 1);
 125     for(Int_t ring = 0;ring<=maxRing;ring++) {
 126
 127       Char_t ringChar = (ring == 0 ? 'I' : 'O');
 128       for(Int_t v=0; v< nVertexBins; v++) {
 129         TH1F* hits = (TH1F*)fList->FindObject(GetPrimName(kHits,det,ringChar,v));
 130
 131         Int_t nNonZero = 0, nNonZeroInData = 0;
 132
 133         //removing edges
 134         for(Int_t i =1;i<=hits->GetNbinsX();i++) {
 135           if(hits->GetBinContent(i) !=0)
 136             nNonZero++;
 137         }
 138
 139         TH1F* sumMultHist = (TH1F*)fMultList.FindObject(Form("hMCHits_vtxbin%d",v));
 140         for(Int_t i =1;i<=sumMultHist->GetNbinsX();i++) {
 141           if(hits->GetBinContent(i) == 0 ) continue;
 142
 143           nNonZeroInData++;
 144
 145           if(nNonZeroInData <=fNbinsToCut || nNonZeroInData > (nNonZero - fNbinsToCut)) {
 146             continue;
 147           }
 148           if(sumMultHist->GetBinContent(i) == 0 && hits->GetBinContent(i) != 0){
 149               sumMultHist->SetBinContent(i,hits->GetBinContent(i));
 150               sumMultHist->SetBinError(i,hits->GetBinError(i));
 151             }
 152           else {
 153
 154
 155             Float_t sumofw = (1/TMath::Power(hits->GetBinError(i),2)) + (1/TMath::Power(sumMultHist->GetBinError(i),2));
 156             Float_t wav =  (hits->GetBinContent(i)*(1/TMath::Power(hits->GetBinError(i),2)) + sumMultHist->GetBinContent(i)*(1/TMath::Power(sumMultHist->GetBinError(i),2)))/sumofw;
 157             Float_t error = 1/TMath::Sqrt(sumofw);
 158             sumMultHist->SetBinContent(i, wav);
 159             sumMultHist->SetBinError(i, error);
 160           }
 161         }
 162       }
 163     }
 164   }
 165 }
 166
 167 //_____________________________________________________________________
 168 void AliFMDDndeta::Init(const Char_t* filename) {
 169
 170   TFile* fin = TFile::Open(filename);
 171
 172   if(!fin) {
 173     AliWarning("No file - aborting !");
 174     return;
 175   }
 176   fList = (TList*)fin->Get("BackgroundCorrected");
 177   AliFMDAnaParameters* pars =  AliFMDAnaParameters::Instance();
 178   pars->Init();
 179
 180   fIsGenerated[kHits]      = kFALSE;
 181   fIsGenerated[kMult]      = kFALSE;
 182   fIsGenerated[kMultTrVtx] = kFALSE;
 183
 184   fIsInit = kTRUE;
 185 }
 186 //_____________________________________________________________________
 187 void AliFMDDndeta::GenerateMult(Analysis what) {
 188
 189
 190   if(!fIsInit) {
 191     AliWarning("Not initialised - call Init to remedy");
 192     return;
 193   }
 194
 195   if(fIsGenerated[what]) {
 196     AliWarning("Already generated - have a look at the results!");
 197     return;
 198   }
 199   else fIsGenerated[what] = kTRUE;
 200
 201   SetNames(what);
 202
 203   if(what == kHits)
 204     GenerateHits();
 205
 206   TH1I* hMCEvents = (TH1I*)fList->FindObject(fPrimEvents.Data());
 207
 208   AliFMDAnaParameters* pars =  AliFMDAnaParameters::Instance();
 209
 210   TH2F* hTmp         =  pars->GetBackgroundCorrection(1,'I',0);
 211   Int_t nVertexBins  =  pars->GetNvtxBins();
 212   Int_t nEtaBins     =  hTmp->GetNbinsX();
 213   Float_t etaMin     =  hTmp->GetXaxis()->GetXmin();
 214   Float_t etaMax     =  hTmp->GetXaxis()->GetXmax();
 215
 216   for(Int_t i = 0; i<nVertexBins;i++) {
 217     TH1F* hMult = new TH1F(Form("hMult_vtxbin%d",i),Form("hMult_vtxbin%d",i),nEtaBins,etaMin,etaMax);
 218     hMult->Sumw2();
 219     fMultList.Add(hMult);
 220   }
 221
 222   for(Int_t det = 1; det<=3;det++) {
 223     Int_t maxRing = (det == 1 ? 0 : 1);
 224     for(Int_t iring = 0; iring<=maxRing; iring++) {
 225       Char_t ringChar = (iring == 0 ? 'I' : 'O');
 226       TH1F* hRingMult= new TH1F(Form("hRingMult_FMD%d%c",det,ringChar),Form("hRingMult_FMD%d%c",det,ringChar),nEtaBins,etaMin,etaMax);
 227       fMultList.Add(hRingMult);
 228     }
 229   }
 230   TH1I* hEvents = (TH1I*)fList->FindObject(fEvents.Data());
 231   TH1F*  hPrimVtx = 0;
 232
 233
 234
 235   for(Int_t det = 1; det<=3; det++) {
 236     Int_t maxRing = (det == 1 ? 0 : 1);
 237     for(Int_t ring = 0;ring<=maxRing;ring++) {
 238
 239       Char_t ringChar = (ring == 0 ? 'I' : 'O');
 240       for(Int_t v=0; v< nVertexBins; v++) {
 241         if(det == 1) {
 242           if(what == kHits)
 243             hPrimVtx = (TH1F*)fMultList.FindObject(Form("hMCHits_vtxbin%d",v));
 244           else
 245             hPrimVtx = (TH1F*)fList->FindObject(GetPrimName(what,det,ringChar,v));
 246
 247
 248           Float_t nPrimevents = hMCEvents->GetBinContent(v+1);
 249           Float_t xb = hPrimVtx->GetNbinsX();
 250           Float_t xr = hPrimVtx->GetXaxis()->GetXmax() - hPrimVtx->GetXaxis()->GetXmin();
 251           hPrimVtx->Scale(xb / xr );
 252           if(nPrimevents > 0)
 253             hPrimVtx->Scale(1/nPrimevents);
 254
 255         }
 256
 257         Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
 258         Float_t vtxZ1   = (delta*v) - pars->GetVtxCutZ();
 259         Float_t vtxZ2   = (delta*(v+1)) - pars->GetVtxCutZ();
 260
 261         if(vtxZ1<(-1*fVtxCut) || vtxZ2 >fVtxCut)
 262           continue;
 263         Float_t nEvents = (Float_t)hEvents->GetBinContent(v+1);
 264
 265         //TH1F* multhistproj = (TH1F*)fList->FindObject(Form("dNdeta_FMD%d%c_vtxbin%d_proj",det,ringChar,v));
 266         TH1F* multhistproj = (TH1F*)fList->FindObject(GetAnalysisName(what,det,ringChar,v));
 267         if(nEvents)
 268           multhistproj->Scale(1/nEvents);
 269         Float_t xnBins = multhistproj->GetNbinsX();
 270         Float_t xrange = multhistproj->GetXaxis()->GetXmax() - multhistproj->GetXaxis()->GetXmin();
 271         multhistproj->Scale(xnBins / xrange);
 272
 273
 274
 275         Int_t nNonZero = 0, nNonZeroInData = 0;
 276
 277         //removing edges
 278
 279         for(Int_t i =1;i<=multhistproj->GetNbinsX();i++) {
 280           if(multhistproj->GetBinContent(i) !=0)
 281             nNonZero++;
 282         }
 283         Int_t nBinsOld = fNbinsToCut;
 284         if(det == 2 && ringChar =='I') {
 285           fNbinsToCut = 0;
 286         }
 287         TH1F* hRingMult = (TH1F*)fMultList.FindObject(Form("hRingMult_FMD%d%c",det,ringChar));
 288
 289         for(Int_t i=1; i<=hRingMult->GetNbinsX(); i++) {
 290           if(multhistproj->GetBinContent(i)!=0) {
 291             nNonZeroInData++;
 292             if(nNonZeroInData <=fNbinsToCut || nNonZeroInData > (nNonZero - fNbinsToCut)) {
 293               continue;
 294             }
 295             Float_t oldweight = 0;
 296             Float_t oldwav    = 0;
 297             if(hRingMult->GetBinError(i)>0) {
 298               oldweight = 1/TMath::Power(hRingMult->GetBinError(i),2);
 299               oldwav    = oldweight*hRingMult->GetBinContent(i);
 300             }
 301             Float_t  weight = 1/TMath::Power(multhistproj->GetBinError(i),2);
 302             Float_t  wav    = oldwav + weight*multhistproj->GetBinContent(i);
 303             Float_t  sumofw = oldweight + weight;
 304             if(sumofw) {
 305               Float_t error = 1/TMath::Sqrt(sumofw);
 306
 307               hRingMult->SetBinContent(i,wav/sumofw);
 308               hRingMult->SetBinError(i,error);
 309
 310             }
 311           }
 312         }
 313         nNonZeroInData = 0;
 314
 315         TH1F* sumMultHist = (TH1F*)fMultList.FindObject(Form("hMult_vtxbin%d",v));
 316         for(Int_t i =1;i<=sumMultHist->GetNbinsX();i++) {
 317           if(multhistproj->GetBinContent(i) != 0 ) {
 318             nNonZeroInData++;
 319
 320             if(nNonZeroInData <=fNbinsToCut || nNonZeroInData > (nNonZero - fNbinsToCut)) {
 321               continue;
 322             }
 323             if(sumMultHist->GetBinContent(i) == 0 && multhistproj->GetBinContent(i) != 0){
 324               sumMultHist->SetBinContent(i,multhistproj->GetBinContent(i));
 325               sumMultHist->SetBinError(i,multhistproj->GetBinError(i));
 326             }
 327             else {
 328
 329
 330               Float_t sumofw = (1/TMath::Power(multhistproj->GetBinError(i),2)) + (1/TMath::Power(sumMultHist->GetBinError(i),2));
 331               Float_t wav =  (multhistproj->GetBinContent(i)*(1/TMath::Power(multhistproj->GetBinError(i),2)) + sumMultHist->GetBinContent(i)*(1/TMath::Power(sumMultHist->GetBinError(i),2)))/sumofw;
 332               Float_t error = 1/TMath::Sqrt(sumofw);
 333               sumMultHist->SetBinContent(i, wav);
 334               sumMultHist->SetBinError(i, error);
 335             }
 336
 337
 338           }
 339         }
 340         fNbinsToCut = nBinsOld;
 341       }
 342
 343     }
 344   }
 345
 346   TH1F* sumMult  = new TH1F("hSumMult","hSumMult",nEtaBins,etaMin,etaMax);
 347   sumMult->Sumw2();
 348   fMultList.Add(sumMult);
 349   TH1F* primMult  = new TH1F("hPrimMult","hPrimMult",nEtaBins,etaMin,etaMax);
 350   primMult->Sumw2();
 351   fMultList.Add(primMult);
 352
 353   Float_t wav  = 0, sumofw=0, weight = 0, primwav  = 0, primsumofw=0, primweight = 0;//, etaofbin = 0;
 354   for(Int_t i =1; i<=sumMult->GetNbinsX();i++) {
 355
 356     for(Int_t v = 0; v<nVertexBins;v++) {
 357       if(what == kHits)
 358         hPrimVtx = (TH1F*)fMultList.FindObject(Form("hMCHits_vtxbin%d",v));
 359       else
 360         hPrimVtx = (TH1F*)fList->FindObject(GetPrimName(what,0,0,v));
 361       TH1F* sumMultHist = (TH1F*)fMultList.FindObject(Form("hMult_vtxbin%d",v));
 362       //etaofbin += sumMultHist->GetBinCenter(i);
 363       // if( hPrimVtx->GetBinContent(i)!=0 && hPrimVtx->GetBinError(i)>0.001 && sumMultHist->GetBinContent(i)!=0) {
 364       if( hPrimVtx->GetBinContent(i)!=0) {
 365         primweight = 1/TMath::Power(hPrimVtx->GetBinError(i),2);
 366         primwav    = primwav + primweight*hPrimVtx->GetBinContent(i);
 367         primsumofw = primsumofw + primweight;
 368       }
 369       if(sumMultHist->GetBinContent(i)!=0) {
 370
 371         weight = 1/TMath::Power(sumMultHist->GetBinError(i),2);
 372         wav    = wav + weight*sumMultHist->GetBinContent(i);
 373         sumofw = sumofw + weight;
 374       }
 375
 376     }
 377     if( primsumofw !=0 ) {// && sumofw !=0) {
 378       Float_t primerror = 1/TMath::Sqrt(primsumofw);
 379
 380       primMult->SetBinContent(i,primwav/primsumofw);
 381       primMult->SetBinError(i,primerror);
 382     }
 383     else {
 384       primMult->SetBinContent(i,0);
 385       primMult->SetBinError(i,0);
 386     }
 387
 388     if(sumofw) {
 389
 390       Float_t error = 1/TMath::Sqrt(sumofw);
 391
 392       sumMult->SetBinContent(i,wav/sumofw);
 393       sumMult->SetBinError(i,error);
 394     }
 395
 396
 397     wav = 0;
 398     sumofw = 0;
 399     weight = 0;
 400     primwav = 0;
 401     primsumofw = 0;
 402     primweight = 0;
 403
 404   }
 405
 406 }
 407 //_____________________________________________________________________
 408 void AliFMDDndeta::DrawDndeta(Analysis what, Int_t rebin) {
 409
 410   if(!fIsInit) {
 411     AliWarning("Not initialised - call Init to remedy");
 412     return;
 413   }
 414
 415   if(!fIsGenerated[what]) {
 416     AliWarning("Not generated - generate first then draw!");
 417     return;
 418   }
 419
 420   SetNames(what);
 421
 422   AliFMDAnaParameters* pars =  AliFMDAnaParameters::Instance();
 423   TCanvas* c1 = new TCanvas("cvertexbins","Overlaps",1400,800);
 424   c1->Divide(5,2);
 425   Int_t nVertexBins  =  pars->GetNvtxBins();
 426
 427   TH1I* hEvents   = (TH1I*)fList->FindObject(fEvents.Data());
 428   TH1I* hMCEvents = (TH1I*)fList->FindObject(fPrimEvents.Data());
 429
 430   TH1F*  hPrimVtx = 0;
 431   for(Int_t det = 1; det<=3; det++) {
 432     Int_t maxRing = (det == 1 ? 0 : 1);
 433     for(Int_t ring = 0;ring<=maxRing;ring++) {
 434       for(Int_t v=0; v< nVertexBins; v++) {
 435         Char_t ringChar = (ring == 0 ? 'I' : 'O');
 436         Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
 437         Float_t vtxZ1   = (delta*v) - pars->GetVtxCutZ();
 438         Float_t vtxZ2   = (delta*(v+1)) - pars->GetVtxCutZ();
 439
 440         if(vtxZ1<(-1*fVtxCut) || vtxZ2 >fVtxCut)
 441           continue;
 442         Float_t nEvents = (Float_t)hEvents->GetBinContent(v+1);
 443
 444         TH1F* multhistproj = (TH1F*)fList->FindObject(GetAnalysisName(what,det,ringChar,v));
 445
 446         c1->cd(v+1);
 447
 448         if(det==1)  {
 449           multhistproj->SetMarkerStyle(20);
 450           multhistproj->SetMarkerColor(1);
 451         }
 452         if(det==2 && ringChar=='I') {
 453           multhistproj->SetMarkerStyle(21);
 454           multhistproj->SetMarkerColor(2);
 455         }
 456         if(det==2 && ringChar=='O') {
 457           multhistproj->SetMarkerStyle(3);
 458           multhistproj->SetMarkerColor(3);
 459         }
 460         if(det==3 && ringChar=='I') {
 461           multhistproj->SetMarkerStyle(22);
 462           multhistproj->SetMarkerColor(4);
 463         }
 464         if(det==3 && ringChar=='O') {
 465           multhistproj->SetMarkerStyle(23);
 466           multhistproj->SetMarkerColor(6);
 467         }
 468
 469
 470
 471         if(det == 1) {
 472           if(what == kHits)
 473             hPrimVtx = (TH1F*)fMultList.FindObject(Form("hMCHits_vtxbin%d",v));
 474           else
 475             hPrimVtx = (TH1F*)fList->FindObject(GetPrimName(what,det,ringChar,v));
 476
 477           hPrimVtx->SetTitle("");
 478           TLatex* l = new TLatex(0.14,0.92,Form("Vtx range [%.1f, %.1f], %d events",vtxZ1,vtxZ2,(Int_t)nEvents));
 479           l->SetNDC(kTRUE);
 480           hPrimVtx->SetFillColor(kGray);
 481           hPrimVtx->SetLabelFont(132,"xyz");
 482           hPrimVtx->SetStats(0000);
 483           hPrimVtx->GetXaxis()->SetTitle("#eta");
 484           hPrimVtx->GetYaxis()->SetRangeUser(0,15);
 485           hPrimVtx->DrawCopy("E3");
 486           l->Draw();
 487           multhistproj->DrawCopy("same");
 488
 489         }
 490         else
 491           multhistproj->DrawCopy("same");
 492
 493
 494       }
 495     }
 496   }
 497
 498   for(Int_t v=0; v< nVertexBins; v++) {
 499     TH1F* sumMultHist = (TH1F*)fMultList.FindObject(Form("hMult_vtxbin%d",v));
 500     c1->cd(v+1);
 501     sumMultHist->SetMarkerStyle(25);
 502     sumMultHist->DrawCopy("same");
 503
 504   }
 505   TH1F* primMult = (TH1F*)fMultList.FindObject("hPrimMult");
 506   TH1F* sumMult  = (TH1F*)fMultList.FindObject("hSumMult");
 507   sumMult->SetMarkerStyle(23);
 508   sumMult->SetMarkerColor(kRed);
 509   for(Int_t det = 1; det<=3;det++) {
 510     Int_t maxRing = (det == 1 ? 0 : 1);
 511     for(Int_t iring = 0; iring<=maxRing; iring++) {
 512       Char_t ringChar = (iring == 0 ? 'I' : 'O');
 513       TH1F* hRingMult= (TH1F*)fMultList.FindObject(Form("hRingMult_FMD%d%c",det,ringChar));
 514       hRingMult->Add(primMult,-1);
 515       hRingMult->Divide(primMult);
 516     }
 517   }
 518
 519
 520   //TH1F*  hPrim = (TH1F*)fList->FindObject(fPrimEvents.Data());
 521   //  TFile testhit("/home/canute/ALICE/Simulations/TestOfAnalysis2/hitsdist.root");
 522   // TH1F* hHitCor = (TH1F*)testhit.Get("selectedHits");
 523   // hHitCor->Sumw2();
 524   // sumMult->Divide(hHitCor);
 525   TH1F*  hPrim = (TH1F*)fList->FindObject(fPrimdNdeta.Data());
 526   Float_t nPrimevents = hMCEvents->GetEntries();
 527   //std::cout<<hMCEvents->GetEntries()<<std::endl;
 528   Float_t xb = hPrim->GetNbinsX();
 529   Float_t xr = hPrim->GetXaxis()->GetXmax() - hPrim->GetXaxis()->GetXmin();
 530   //std::cout<<xb/xr<<std::endl;
 531   hPrim->Scale(xb / xr );
 532   if(nPrimevents > 0)
 533     hPrim->Scale(1/nPrimevents);
 534
 535     TFile testin("/home/canute/ALICE/FMDanalysis/GridAnalysis/fmd_dNdeta_hits.root","READ");
 536     // TFile testin("/home/canute/ALICE/FMDanalysis/productionData/fmd_dNdeta_hits.root","READ");
 537   TH1F* hcorrect = (TH1F*)testin.Get("hReldif");
 538   hcorrect->SetName("djarnis");
 539   std::cout<<hcorrect<<std::endl;
 540   for(Int_t bb = 1;bb<=hcorrect->GetNbinsX();bb++) {
 541     hcorrect->SetBinContent(bb,hcorrect->GetBinContent(bb)+1);
 542   }
 543
 544   sumMult->Divide(hcorrect);
 545   //delete hcorrect;
 546
 547
 548   sumMult->Rebin(rebin);
 549   sumMult->Scale(1/rebin);
 550   primMult->Rebin(rebin);
 551   primMult->Scale(1/rebin);
 552   hPrim->Rebin(rebin);
 553   hPrim->Scale(1/rebin);
 554   if(what != kHits)
 555     primMult = hPrim;
 556   TH1F* hReldif = (TH1F*)sumMult->Clone("hReldif");
 557   hReldif->Add(primMult,-1);
 558   hReldif->Divide(primMult);
 559   //  hReldif->Add(hPrim,-1);
 560   // hReldif->Divide(hPrim);
 561
 562   TCanvas* c2 = new TCanvas("dN/deta","dN/deta",640,960);
 563   c2->Divide(1,2);//,0,0);
 564   c2->cd(2);
 565   gPad->SetGridy();
 566   gPad->SetGridx();
 567   hReldif->SetTitle("");
 568   hReldif->GetYaxis()->SetTitle("Relative difference");
 569   hReldif->GetYaxis()->SetRangeUser(-0.2,0.2);
 570   hReldif->SetStats(0000);
 571   hReldif->GetXaxis()->SetTitle("#eta");
 572   hReldif->DrawCopy();
 573   TLine* zeroLine = new TLine(-6,0,6,0);
 574   zeroLine->SetLineWidth(2);
 575   zeroLine->SetLineColor(kBlack);
 576   zeroLine->Draw();
 577   hPrim->SetTitle("");
 578   hPrim->SetLabelFont(132,"xyz");
 579   hPrim->SetFillColor(kGray);
 580   primMult->SetFillColor(kBlue);
 581   c2->cd(1);
 582   gPad->SetGridy();
 583   gPad->SetGridx();
 584   hPrim->GetYaxis()->SetRangeUser(2,10);
 585   //hPrim->GetYaxis()->SetRangeUser(0,20);
 586   hPrim->SetStats(0000);
 587   hPrim->GetXaxis()->SetTitle("#eta");
 588
 589   sumMult->DrawCopy("PE");
 590   hPrim->DrawCopy("E3same");
 591   //primMult->DrawCopy("E3same");
 592   hPrim->GetXaxis()->SetTitle("#eta");
 593
 594   hPrim->SetMarkerStyle(8);
 595
 596   hPrim->DrawCopy("PE3same");
 597   sumMult->DrawCopy("PEsame");
 598   std::cout<<"FMD 1 "<<sumMult->Integral(sumMult->FindBin(3),sumMult->FindBin(5),"width")<<std::endl;
 599   /*TFile* itsfile = TFile::Open(itsfilename);
 600   if(itsfile) {
 601     //itsfile->cd();
 602
 603     //TList* itslist = (TList*)itsfile->Get("ITS_analysis");
 604     //TH1F* itshist0 = (TH1F*)itslist->FindObject("dndeta_check_0");
 605     //TH1F* itshist1 = (TH1F*)itslist->FindObject("dndeta_check_1");
 606     //TH1F* itshist2 = (TH1F*)itslist->FindObject("dndeta_check_2");
 607     TH1F* itshist0 = (TH1F*)itsfile->Get("dndeta/dNdEta");
 608     TH1F* itshist1= (TH1F*)itsfile->Get("dndeta/dNdEta_1");
 609     TH1F* itshist2 = (TH1F*)itsfile->Get("dndeta/dNdEta_2");
 610     itshist0->SetMarkerStyle(27);
 611     itshist1->SetMarkerStyle(28);
 612     itshist2->SetMarkerStyle(30);
 613     itshist0->DrawCopy("same");
 614     itshist1->DrawCopy("same");
 615     itshist2->DrawCopy("same");
 616
 617   }
 618   */
 619
 620
 621   TLegend* leg = new TLegend(0.35,0.2,0.65,0.35,"");
 622   leg->AddEntry(hPrim,"Primary","pf");
 623   //leg->AddEntry(primMult,"Analysed prim","f");
 624   leg->AddEntry(sumMult,"Analysis","p");
 625   leg->Draw();
 626   TString species;
 627
 628   switch(what) {
 629   case kMult:
 630     species.Form("mult");
 631     break;
 632   case kMultTrVtx:
 633     species.Form("mult_TrVtx");
 634     break;
 635   case kHits:
 636     species.Form("hits");
 637     break;
 638   default:
 639     AliWarning("no valid Analysis entry!");
 640     break;
 641   }
 642   TFile fout(Form("fmd_dNdeta_%s.root",species.Data()),"recreate");
 643   sumMult->Write();
 644   hReldif->Write();
 645   fMultList.Write();
 646   fout.Close();
 647
 648 }
 649 //_____________________________________________________________________
 650 void AliFMDDndeta::CreateSharingEfficiency(const Char_t* filename, Bool_t store) {
 651
 652   Init(filename);
 653   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 654   //pars->Init(kTRUE,AliFMDAnaParameters::kBackgroundCorrection);
 655   Int_t nVertexBins = pars->GetNvtxBins();
 656
 657   SetNames(kHits);
 658   TH1I* hEvents          = (TH1I*)fList->FindObject(fEvents.Data());
 659   TH1I* hPrimEvents      = (TH1I*)fList->FindObject(fPrimEvents.Data());
 660
 661   SetNames(kHitsTrVtx);
 662   TH1I* hEventsTrVtx     = (TH1I*)fList->FindObject(fEvents.Data());
 663   TH1I* hPrimEventsTrVtx = (TH1I*)fList->FindObject(fPrimEvents.Data());
 664
 665   AliFMDAnaCalibSharingEfficiency* sharEff = new AliFMDAnaCalibSharingEfficiency();
 666
 667   for(Int_t det = 1; det<=3; det++) {
 668     Int_t maxRing = (det == 1 ? 0 : 1);
 669     for(Int_t ring = 0;ring<=maxRing;ring++) {
 670       Char_t ringChar = (ring == 0 ? 'I' : 'O');
 671       for(Int_t v=0; v< nVertexBins; v++) {
 672         TH1F* hHits = (TH1F*)fList->FindObject(GetAnalysisName(kHits,det, ringChar, v));
 673         TH1F* hHitsTrVtx   = (TH1F*)fList->FindObject(GetAnalysisName(kHitsTrVtx,det, ringChar, v));
 674         TH1F* hMCHits      = (TH1F*)fList->FindObject(GetPrimName(kHits,det, ringChar, v));
 675         TH1F* hMCHitsTrVtx = (TH1F*)fList->FindObject(GetPrimName(kHitsTrVtx,det, ringChar, v));
 676
 677         TH1F* hCorrection  = (TH1F*)hHits->Clone(Form("hCorrection_FMD%d%c_vtx%d"));
 678         TH1F* hCorrectionTrVtx  = (TH1F*)hHitsTrVtx->Clone(Form("hCorrection_FMD%d%c_vtx%d"));
 679
 680         hCorrection->Scale(1./(Float_t)hEvents->GetBinContent(v+1));
 681         hMCHits->Scale(1./(Float_t)hPrimEvents->GetBinContent(v+1));
 682         hCorrectionTrVtx->Scale(1./(Float_t)hEventsTrVtx->GetBinContent(v+1));
 683         hMCHitsTrVtx->Scale(1./(Float_t)hPrimEventsTrVtx->GetBinContent(v+1));
 684         hCorrection->Divide(hMCHits);
 685         hCorrectionTrVtx->Divide(hMCHitsTrVtx);
 686
 687         sharEff->SetSharingEff(det,ringChar,v,hCorrection);
 688         sharEff->SetSharingEffTrVtx(det,ringChar,v,hCorrectionTrVtx);
 689         //      std::cout<<hHits<<"  "<<hHitsTrVtx<<"   "<<hPrim<<"    "<<hPrimTrVtx<<std::endl;
 690
 691         }
 692       }
 693     }
 694
 695   if(store) {
 696     TFile fsharing(pars->GetPath(pars->GetSharingEffID()),"RECREATE");
 697     sharEff->Write(AliFMDAnaParameters::GetSharingEffID());
 698     fsharing.Close();
 699   }
 700
 701 }
 702
 703 //_____________________________________________________________________
 704 //
 705 // EOF
 706 //