//////////////////////////////////////////////////
//
//  This macro was written by Domenico Colella (domenico.colella@cern.ch).
//   -- first version    [12 November 2013]
//   -- modified version [3 November 2013]: added inv. mass distr. without PID info
//
//   ------------------------
//   ------ Arguments -------
//   ------------------------
//   --  icasType          =  0) Xi- 1) Xi+ 2) Omega- 3) Omega+
//   --  collidingsystem   =  0) PbPb  1) pp  2) pPb
//   --  isMC              =  kTRUE if running on MC production 
//   --  fileDir           =  "Input file directory"
//   --  filein            =  "Input file name"
//
//
//   -------------------------------------
//   ------ QATask output content --------
//   -------------------------------------
//   The output produced by the QATask is a CFContainer with 4 steps and 21 variables.
//   The meaning of each variable within the container are listed here:
//   --  0   = Max DCA Cascade Daughters                 pp: 2.0     PbPb: 0.3     pPb: 2.0
//   --  1   = Min DCA Bach To PV                        pp: 0.01    PbPb: 0.03    pPb: 0.03
//   --  2   = Min Cascade Cosine Of PA                  pp: 0.98    PbPb: 0.999   pPb: 0.95
//   --  3   = Min Cascade Radius Fid. Vol.              pp: 0.2     PbPb: 0.9     pPb: 0.4
//   --  4   = Window Invariant Mass Lambda              pp: 0.008   PbPb: 0.0008  pPb: 0.010
//   --  5   = Max DCA V0 Daughters                      pp: 1.5     PbPb: 1.0     pPb: 2.0
//   --  6   = Min V0 Cosine Of PA To PV                 pp: pT dep. PbPb: 0.98    pPb: 0.95
//   --  7   = Min V0 Radius Fid. Vol.                   pp: 0.2     PbPb: 0.9     pPb: 1.0
//   --  8   = Min DCA V0 To PV                          pp: 0.01    PbPb: 0.05    pPb: 0.05
//   --  9   = Min DCA Pos To PV                         pp: 0.05    PbPb: 0.1     pPb: 0.02
//   --  10  = Min DCA Neg To PV                         pp: 0.05    PbPb: 0.1     pPb: 0.02
//   --  11  = Invariant Mass distribution for Xi
//   --  12  = Invariant Mass distribution for Omega
//   --  13  = Transverse Momentum distribution
//   --  14  = Rapidity distribution for Xi
//   --  15  = Rapidity distribution for Omega
//   --  16  = Proper length distribution for the cascade
//   --  17  = Proper length distribution for the V0
//   --  18  = Min V0 Cosine Of PA To Xi Vertex         pp: pT dep. PbPb: pT dep.
//   --  19  = Centrality
//   --  20  = ESD track multiplicity
//   The last two variables are empty in the case proton-proton collisions.
//   In case of MC production one more CFContainer is produced, containing infos on the 
//   generated particles. As the previous container, this one is composed by 4 steps, one
//   for each cascade and 7 variables:
//    -- 0   = Total momentum
//    -- 1   = Transverse momentum
//    -- 2   = Rapidity
//    -- 3   = Pseudo-rapidity
//    -- 4   = Theta angle
//    -- 5   = Phi angle
//    -- 6   = Centrality
//    The previous container is still produced with the informations from the reconstructed
//    particles.
//
//
//   -----------------------------------
//   ------ Present Macro Checks -------
//   -----------------------------------
//   Using this macro many checks on the cascade topological reconstruction procedure
//   can be performed. In particular, the shape and the limit for the topological 
//   variable distributions as well as other kinematical variable distributions. The
//   reconstruction of the cascades are performed using two classes AliCascadeVertexer.cxx 
//   and AliV0vertexer.cxx contained in /STEER/ESD/ folder in Aliroot.
//   In the following are listed the contents of each page of the produced pdf:
//   
//   -- [Page 1] Distributions for the variables: 
//                DCA cascade daughters,  Bachelor IP to PV, 
//                Cascade cosine of PA,   Cascade radius of fiducial volume, 
//                Invariant mass Lambda,  DCA V0 daughters.
//   -- [Page 2] Distributions for the variables:
//                V0 cosine of PA to PV,  Min V0 Radius fiducial volume, 
//                Min DCA V0 To PV,       Min DCA positive To PV, 
//                Min DCA negative To PV, V0 cosine of PA to XiV
//   -- [Page 3] Distributions for the variables;
//                InvMass,                Transverse momentum,
//                Rapidity,               Cascade proper length, 
//                V0 proper length.
//   -- [Page 4] Check on the invariant mass distribution fit.
//   -- [Page 5] Check on the invariant mass distribution, without the PID info, fit. 
//   -- [Page 6] Only in case of PbPb or pPb collisions, the event centrality
//               distribution.
//   -- [Page 7] Only in case of MC production, distributions for the MC generated
//               particles, of the variables:
//                Total momentum,         Transverse momentum,
//                Rapidity,               Pseudo-rapidity,
//                Theta angle,            Phi angle,
//
//
//////////////////////////////////////////////////////




class AliCFContainer;

void PostProcessQAMultistrange(Int_t   icasType        = 0,                             // 0) Xi- 1) Xi+ 2) Omega- 3) Omega+
                               Int_t   collidingsystem = 0,                             // 0) PbPb  1) pp 2) pPb
                               Bool_t  isMC            = kFALSE,                         // kTRUE-->MC and kFALSE-->Exp.
                               Char_t *fileDir         = ".",                           // Input file directory
                               Char_t *filein          = "AnalysisResults.root"         // Input file name
                              ) {


     //___________________
     //DEFINE DRAW OPTIONS
     gStyle->SetOptStat(1110);
     gStyle->SetOptStat(kFALSE);
     gStyle->SetOptTitle(kFALSE);
     gStyle->SetFrameLineWidth(2.5);
     gStyle->SetCanvasColor(0);
     gStyle->SetPadColor(0);
     gStyle->SetHistLineWidth(2.5);
     gStyle->SetLabelSize(0.05, "x");
     gStyle->SetLabelSize(0.05, "y");
     gStyle->SetTitleSize(0.05, "x");
     gStyle->SetTitleSize(0.05, "y");
     gStyle->SetTitleOffset(1.1, "x");
     gStyle->SetPadBottomMargin(0.14);

     //_______________________
     //SOURCE USEFUL LIBRARIES
     gSystem->Load("libANALYSIS.so");
     gSystem->Load("libANALYSISalice.so");
     gSystem->Load("libCORRFW.so");

     //_________________________________
     //SOURCE THE FILE AND THE CONTAINER
     TFile *f1 = new TFile(Form("%s/%s",fileDir,filein));
     AliCFContainer *cf = (AliCFContainer*) (f1->Get("PWGLFStrangeness.outputCheckCascade/fCFContCascadeCuts"));  
     TList *hlist = (TList*) f1->Get("PWGLFStrangeness.outputCheckCascade/fListHistMultistrangeQA");
 
     //____________
     //DEEFINE TEXT
     TLatex* t1 = new TLatex(0.6,0.55,"#color[3]{OK!!}");
     t1->SetTextSize(0.1);
     t1->SetNDC();
     TLatex* t2 = new TLatex(0.6,0.55,"#color[2]{NOT OK!!}");
     t2->SetTextSize(0.1);
     t2->SetNDC();
     t2->SetTextColor(2);
     TLatex* tcasc;
     if      (icasType == 0) tcasc = new TLatex(0.8,0.7,"#color[1]{#Xi^{-}}");
     else if (icasType == 1) tcasc = new TLatex(0.8,0.7,"#color[1]{#Xi^{+}}");
     else if (icasType == 2) tcasc = new TLatex(0.8,0.7,"#color[1]{#Omega^{-}}");
     else if (icasType == 3) tcasc = new TLatex(0.8,0.7,"#color[1]{#Omega^{+}}");
     tcasc->SetTextSize(0.15);
     tcasc->SetNDC();
     tcasc->SetTextColor(2);
     TLatex* tpdgmass;
     if      (icasType == 0) tpdgmass = new TLatex(0.55,0.7,"#color[1]{PDG mass: 1.321 GeV/c^{2}}");
     else if (icasType == 1) tpdgmass = new TLatex(0.55,0.7,"#color[1]{PDG mass: 1.321 GeV/c^{2}}");
     else if (icasType == 2) tpdgmass = new TLatex(0.55,0.7,"#color[1]{PDG mass: 1.672 GeV/c^{2}}");
     else if (icasType == 3) tpdgmass = new TLatex(0.55,0.7,"#color[1]{PDG mass: 1.672 GeV/c^{2}}");
     tpdgmass->SetTextSize(0.07);
     tpdgmass->SetNDC();
     tpdgmass->SetTextColor(2);

     
     //________________________________ 
     //DEFINE 1st CANVAS AND DRAW PLOTS
     TCanvas *c1 = new TCanvas("c1","",1200,800);
     c1->Divide(2,3); 
       //Pad 1: DCA cascade daughters
       c1->cd(1);
       gPad->SetLogy();
       TH1D *hvar0 = cf->ShowProjection(0,icasType);
       hvar0->Draw("histo");
       Double_t x0;
       if      (collidingsystem == 0) x0 = 0.3;
       else if (collidingsystem == 1) x0 = 2.0;
       else if (collidingsystem == 2) x0 = 2.0;
       TLine *line0 = new TLine(x0,0.,x0,hvar0->GetBinContent(hvar0->GetMaximumBin()));
       line0->SetLineColor(kRed);
       line0->SetLineStyle(9);
       line0->SetLineWidth(2.0);
       line0->Draw("same");
          Bool_t check_0 = checkOverTheLimit(hvar0, x0);
          if (check_0) { cout<<"The cut is OK!!"<<endl; t1->Draw(); }
          else         { cout<<"The cut is NOT OK!!"<<endl; t2->Draw(); }
       tcasc->Draw();
       //Pad 2: Bachelor IP to PV
       c1->cd(2);
       gPad->SetLogy();
       TH1D *hvar1 = cf->ShowProjection(1,icasType);
       hvar1->GetXaxis()->SetRangeUser(0.,0.24);
       hvar1->Draw("histo");
       Double_t x1;
       if      (collidingsystem == 0) x1 = 0.03;
       else if (collidingsystem == 1) x1 = 0.01;
       else if (collidingsystem == 2) x1 = 0.03;
       TLine *line1 = new TLine(x1,0.,x1,hvar1->GetBinContent(hvar1->GetMaximumBin()));
       line1->SetLineColor(kRed);
       line1->SetLineStyle(9);
       line1->SetLineWidth(2.0);
       line1->Draw("same");
          Bool_t check_1 = checkUnderTheLimit(hvar1, x1);
          if (check_1) { cout<<"The cut is OK!!"<<endl; t1->Draw(); }
          else         { cout<<"The cut is NOT OK!!"<<endl; t2->Draw(); }
       //Pad 3: Cascade cosine of Pointing Angle
       c1->cd(3);
       gPad->SetLogy();
       TH1D *hvar2 = cf->ShowProjection(2,icasType);
       Double_t max2 = hvar2->GetBinContent(hvar2->GetMaximumBin());
       hvar2->GetYaxis()->SetRangeUser(0.01,max2*1.5);
       hvar2->Draw("histo");
       Double_t x2;
       if      (collidingsystem == 0) x2 = 0.999;
       else if (collidingsystem == 1) x2 = 0.98;
       else if (collidingsystem == 2) x2 = 0.95;
       TLine *line2 = new TLine(x2,0.,x2,hvar2->GetBinContent(hvar2->GetMaximumBin()));
       line2->SetLineColor(kRed);
       line2->SetLineStyle(9);
       line2->SetLineWidth(2.0);
       line2->Draw("same");
       line1->Draw("same");
          Bool_t check_2 = checkUnderTheLimit(hvar2, x2);
          if (check_2) { cout<<"The cut is OK!!"<<endl; t1->Draw(); }
          else         { cout<<"The cut is NOT OK!!"<<endl; t2->Draw(); }
       //Pad 4: Cascade radius of fiducial volume
       c1->cd(4);
       gPad->SetLogy();
       TH1D *hvar3 = cf->ShowProjection(3,icasType);
       hvar3->GetXaxis()->SetRangeUser(0.,3.8);
       hvar3->Draw("histo");
       Double_t x3;
       if      (collidingsystem == 0) x3 = 0.9;
       else if (collidingsystem == 1) x3 = 0.2;
       else if (collidingsystem == 2) x3 = 0.4;
       TLine *line3 = new TLine(x3,0.,x3,hvar3->GetBinContent(hvar3->GetMaximumBin()));
       line3->SetLineColor(kRed);
       line3->SetLineStyle(9);
       line3->SetLineWidth(2.0);
       line3->Draw("same");
          Bool_t check_3 = checkUnderTheLimit(hvar3, x3);
          if (check_3) { cout<<"The cut is OK!!"<<endl; t1->Draw(); }
          else         { cout<<"The cut is NOT OK!!"<<endl; t2->Draw(); }
       //Pad 5: Invariant mass Lambda
       c1->cd(5);
       TH1D *hvar4 = cf->ShowProjection(4,icasType);
       hvar4->Draw("histo");
       Double_t x41;
       if      (collidingsystem < 2)  x41 = 1.116 + 0.008;
       else if (collidingsystem == 2) x41 = 1.116 + 0.010;
       TLine *line41 = new TLine(x41,0.,x41,hvar4->GetBinContent(hvar4->GetMaximumBin()));
       line41->SetLineColor(kRed);
       line41->SetLineStyle(9);
       line41->SetLineWidth(2.0);
       line41->Draw("same");
       Double_t x42;
       if      (collidingsystem < 2)  x42 = 1.115 - 0.008;
       else if (collidingsystem == 2) x42 = 1.115 - 0.010;
       TLine *line42 = new TLine(x42,0.,x42,hvar4->GetBinContent(hvar4->GetMaximumBin()));
       line42->SetLineColor(kRed);
       line42->SetLineStyle(9);
       line42->SetLineWidth(2.0);
       line42->Draw("same");
          Bool_t check_4_1 = checkUnderTheLimit(hvar4, x42);
          Bool_t check_4_2 = checkOverTheLimit(hvar4, x41);
          if (check_4_1 && check_4_2) { cout<<"The cut is OK!!"<<endl; t1->Draw(); }
          else                        { cout<<"The cut is NOT OK!!"<<endl; t2->Draw(); }
       //Pad 6: DCA V0 daughters
       c1->cd(6);
       gPad->SetLogy();
       TH1D *hvar5 = cf->ShowProjection(5,icasType);
       hvar5->Draw("histo");
       Double_t x5;
       if      (collidingsystem == 0) x5 = 1.0;
       else if (collidingsystem == 1) x5 = 1.5;
       else if (collidingsystem == 2) x5 = 2.0;
       TLine *line5 = new TLine(x5,0.,x5,hvar5->GetBinContent(hvar5->GetMaximumBin()));
       line5->SetLineColor(kRed);
       line5->SetLineStyle(9);
       line5->SetLineWidth(2.0);
       line5->Draw("same");
          Bool_t check_5 = checkOverTheLimit(hvar5, x5);
          if (check_5) { cout<<"The cut is OK!!"<<endl; t1->Draw(); }
          else         { cout<<"The cut is NOT OK!!"<<endl; t2->Draw(); }
     c1->SaveAs("fig_lf_Multistrange.pdf(");
    
     //________________________________
     //DEFINE 2nd CANVAS AND DRAW PLOTS
     TCanvas *c2 = new TCanvas("c2","",1200,800);
     c2->Divide(2,3);
       //Pad 1: V0 cosine of Pointing Angle to PV
       c2->cd(1);
       gPad->SetLogy();
       TH1D *hvar6 = cf->ShowProjection(6,icasType);
       Double_t max6 = hvar6->GetBinContent(hvar6->GetMaximumBin());
       hvar6->GetYaxis()->SetRangeUser(0.01,max6*1.5);
       hvar6->Draw("histo");
       tcasc->Draw();
       //Pad 2: Min V0 Radius Fid. Vol.  
       c2->cd(2);
       gPad->SetLogy();
       TH1D *hvar7 = cf->ShowProjection(7,icasType);
       hvar7->GetXaxis()->SetRangeUser(0.,3.0);
       hvar7->Draw("histo");
       Double_t x7;
       if      (collidingsystem == 0) x7 = 0.9;
       else if (collidingsystem == 1) x7 = 0.2;
       else if (collidingsystem == 2) x7 = 0.4;
       TLine *line7 = new TLine(x7,0.,x7,hvar7->GetBinContent(hvar7->GetMaximumBin()));
       line7->SetLineColor(kRed);
       line7->SetLineStyle(9);
       line7->SetLineWidth(2.0);
       line7->Draw("same");
          Bool_t check_7 = checkUnderTheLimit(hvar7, x7);
          if (check_7) { cout<<"The cut is OK!!"<<endl; t1->Draw(); }
          else         { cout<<"The cut is NOT OK!!"<<endl; t2->Draw(); }
       //Pad3: Min DCA V0 To PV
       c2->cd(3);
       gPad->SetLogy();
       TH1D *hvar8 = cf->ShowProjection(8,icasType);
       hvar8->GetXaxis()->SetRangeUser(0.,0.3);
       hvar8->Draw("histo");
       Double_t x8;
       if      (collidingsystem == 0) x8 = 0.05;
       else if (collidingsystem == 1) x8 = 0.01;
       else if (collidingsystem == 2) x8 = 0.05;
       TLine *line8 = new TLine(x8,0.,x8,hvar8->GetBinContent(hvar8->GetMaximumBin()));
       line8->SetLineColor(kRed);
       line8->SetLineStyle(9);
       line8->SetLineWidth(2.0);
       line8->Draw("same");
          Bool_t check_8 = checkUnderTheLimit(hvar8, x8);
          if (check_8) { cout<<"The cut is OK!!"<<endl; t1->Draw(); }
          else         { cout<<"The cut is NOT OK!!"<<endl; t2->Draw(); }
       //Pad 4: Min DCA Pos To PV
       c2->cd(4);
       gPad->SetLogy();
       TH1D *hvar9 = cf->ShowProjection(9,icasType);
       hvar9->GetXaxis()->SetRangeUser(0.,0.2);
       hvar9->Draw("histo");
       Double_t x9;
       if      (collidingsystem == 0) x9 = 0.1;
       else if (collidingsystem == 1) x9 = 0.05;
       else if (collidingsystem == 2) x9 = 0.02;
       TLine *line9 = new TLine(x9,0.,x9,hvar9->GetBinContent(hvar9->GetMaximumBin()));
       line9->SetLineColor(kRed);
       line9->SetLineStyle(9);
       line9->SetLineWidth(2.0);
       line9->Draw("same");
          Bool_t check_9 = checkUnderTheLimit(hvar9, x9);
          if (check_9) { cout<<"The cut is OK!!"<<endl; t1->Draw(); }
          else         { cout<<"The cut is NOT OK!!"<<endl; t2->Draw(); }
       //Pad 5: Min DCA Neg To PV
       c2->cd(5);
       gPad->SetLogy();
       TH1D *hvar10 = cf->ShowProjection(10,icasType);
       hvar10->GetXaxis()->SetRangeUser(0.,0.2);
       hvar10->Draw("histo");
       Double_t x10;
       if      (collidingsystem == 0) x10 = 0.1;
       else if (collidingsystem == 1) x10 = 0.05;
       else if (collidingsystem == 2) x10 = 0.02;
       TLine *line10 = new TLine(x10,0.,x10,hvar10->GetBinContent(hvar10->GetMaximumBin()));
       line10->SetLineColor(kRed);
       line10->SetLineStyle(9);
       line10->SetLineWidth(2.0);
       line10->Draw("same");
          Bool_t check_10 = checkUnderTheLimit(hvar10, x10);
          if (check_10) { cout<<"The cut is OK!!"<<endl; t1->Draw(); }
          else         { cout<<"The cut is NOT OK!!"<<endl; t2->Draw(); }
       //Pad 6: V0 cosine of Pointing Angle to Xi vtx
       c2->cd(6);
       gPad->SetLogy();
       TH1D *hvar20 = cf->ShowProjection(18,icasType);
       Double_t max20 = hvar20->GetBinContent(hvar20->GetMaximumBin());
       hvar20->GetYaxis()->SetRangeUser(0.01,max20*1.5);
       hvar20->Draw("histo");
     c2->SaveAs("fig_lf_Multistrange.pdf");

     //________________________________
     //DEFINE 3rd CANVAS AND DRAW PLOTS
     TCanvas *c3 = new TCanvas("c3","",1200,800);
     c3->Divide(2,2);
       //Pad 1: Transverse momentum
       c3->cd(1);
       TH1D *hvar13 = cf->ShowProjection(13,icasType);
       hvar13->Draw("histo");
       tcasc->Draw();
       //Pad 2: Y
       c3->cd(2);
       TH1D *hvar14 = cf->ShowProjection(14+icasType/2,icasType);
       hvar14->Draw("histo");
       //Pad 3: Cascade proper length
       c3->cd(3);
       TH1D *hvar18;
       hvar18 = cf->ShowProjection(16,icasType);
       hvar18->GetXaxis()->SetRangeUser(0.,90.);
       hvar18->Draw("histo");
       //Pad 4: V0 proper length 
       c3->cd(4);
       TH1D *hvar19;
       hvar19 = cf->ShowProjection(17,icasType);
       hvar19->GetXaxis()->SetRangeUser(0.,90.);
       hvar19->Draw("histo");
       //Pad 6
       // empty 
     c3->SaveAs("fig_lf_Multistrange.pdf");

     //________________________________ 
     //DEFINE 4th CANVAS AND DRAW PLOTS
     TCanvas *c4 = new TCanvas("c4","",1200,800);
     c4->Divide(2,2);
       //Pad1: invariant mass fit
       c4->cd(1);
       TH1D *hvar23 = cf->ShowProjection(11+icasType/2,icasType);
       hvar23->Draw("histo");
       tcasc->Draw();
        // - SOME PARAMETER VALUE
        Bool_t kfitgauss = kFALSE;
        Bool_t kfitleft  = kFALSE;
        Bool_t kfitright = kFALSE;
        // - SOME DEFINITIONS
        Float_t lowlimmass;
        Float_t uplimmass;
        Float_t lowgausslim;
        Float_t upgausslim;
        if (icasType==0||icasType==1) {
            lowlimmass=1.30;
            uplimmass=1.34;
            lowgausslim=1.312;
            upgausslim=1.332;
        } else {
            lowlimmass=1.645;
            uplimmass=1.70;
            lowgausslim=1.668;
            upgausslim=1.678;
        }
        TF1*  fitinvmass = new TF1("fitinvmass","gaus(0)+pol2(3)",lowlimmass,uplimmass);
        fitinvmass->SetParName(0, "cnstntG");
        fitinvmass->SetParName(1, "meanG");
        fitinvmass->SetParName(2, "sigmaG");
        fitinvmass->SetParLimits(0,0.,500000.);
        if (icasType==0||icasType==1) {
            fitinvmass->SetParameter(1, 1.32171);
            fitinvmass->SetParLimits(1, 1.31,1.33);
            fitinvmass->SetParLimits(2,0.001,0.005);
        } else {
            fitinvmass->SetParameter(1, 1.67245);
            fitinvmass->SetParLimits(1, 1.664,1.68);
            fitinvmass->SetParLimits(2,0.0008,0.006);
        }
        hvar23->Fit("fitinvmass","rimeN");
        fitinvmass->SetLineColor(kRed);
        fitinvmass->Draw("same");
        Float_t meanGauss   = fitinvmass->GetParameter(1);
        Float_t sigmaGauss  = fitinvmass->GetParameter(2);
       cout<<"Mean: "<<meanGauss<<endl;
       cout<<"Sigma: "<<sigmaGauss<<endl;
       //Pad2: Text
       c4->cd(2);
       Float_t refwidth = 0.002;
       if (icasType > 1) refwidth = 0.0025;
       TPaveText *pave1 = new TPaveText(0.05,0.3,0.95,0.7);
       pave1->SetFillColor(0);
       pave1->SetTextSize(0.04);
       pave1->SetTextAlign(12);
       pave1->AddText("Invariant mass distribution #color[4]{WITH} PID on dauther tracks");
       if (icasType < 2) pave1->AddText("PDG mass: 1.32171 GeV/c^{2}");
       else              pave1->AddText("PDG mass: 1.67245 GeV/c^{2}");
       pave1->AddText(Form("#color[1]{Mass from Fit: %.5f #pm %.5f GeV/c^{2}}",meanGauss,sigmaGauss));
       if (sigmaGauss > refwidth - 0.0003 && sigmaGauss < refwidth + 0.0003) pave1->AddText("#color[3]{OK!! The width is compatible with standard.}");
       else                                                                  pave1->AddText("#color[2]{NOT OK!! Problem.}");
       pave1->Draw();
       cout<<"   "<<refwidth - 0.0003<<"<"<<sigmaGauss<<"<"<<refwidth + 0.0003<<endl;
       //Pad3: invariant mass fit
       c4->cd(3);
       TH1D *hvar24;
       if      (icasType == 0) hvar24 = (TH1D*) hlist->FindObject(Form("fHistMassXiMinus"));
       else if (icasType == 1) hvar24 = (TH1D*) hlist->FindObject(Form("fHistMassXiPlus"));
       else if (icasType == 2) hvar24 = (TH1D*) hlist->FindObject(Form("fHistMassOmegaMinus"));
       else if (icasType == 3) hvar24 = (TH1D*) hlist->FindObject(Form("fHistMassOmegaPlus"));
       hvar24->Draw("histo");
        TF1*  fitinvmassnoPID = new TF1("fitinvmassnoPID","gaus(0)+pol2(3)",lowlimmass,uplimmass);
        fitinvmassnoPID->SetParName(0, "cnstntG");
        fitinvmassnoPID->SetParName(1, "meanG");
        fitinvmassnoPID->SetParName(2, "sigmaG");
        fitinvmassnoPID->SetParLimits(0,0.,500000.);
        if (icasType==0||icasType==1) {
            fitinvmassnoPID->SetParameter(1, 1.32171);
            fitinvmassnoPID->SetParLimits(1, 1.31,1.33);
            fitinvmassnoPID->SetParLimits(2,0.001,0.005);
        } else {
            fitinvmassnoPID->SetParameter(1, 1.67245);
            fitinvmassnoPID->SetParLimits(1, 1.664,1.68);
            fitinvmassnoPID->SetParLimits(2,0.0008,0.006);
        }
        hvar24->Fit("fitinvmassnoPID","rimeN");
        fitinvmassnoPID->SetLineColor(kRed);
        fitinvmassnoPID->Draw("same");
        Float_t meanGaussnoPID   = fitinvmassnoPID->GetParameter(1);
        Float_t sigmaGaussnoPID  = fitinvmassnoPID->GetParameter(2);
       cout<<"Mean: "<<meanGaussnoPID<<endl;
       cout<<"Sigma: "<<sigmaGaussnoPID<<endl;
       //Pad4: Text
       c4->cd(4);
       Float_t refwidth = 0.002;
       if (icasType > 1) refwidth = 0.0025;
       TPaveText *pave2 = new TPaveText(0.05,0.3,0.95,0.7);
       pave2->SetFillColor(0);
       pave2->SetTextSize(0.04);
       pave2->SetTextAlign(12);
       pave2->AddText("Invariant mass distribution #color[4]{WITHOUT} PID on dauther tracks");
       if (icasType < 2) pave2->AddText("PDG mass: 1.32171 GeV/c^{2}");
       else              pave2->AddText("PDG mass: 1.67245 GeV/c^{2}");
       pave2->AddText(Form("#color[1]{Mass from Fit: %.5f #pm %.5f GeV/c^{2}}",meanGaussnoPID,sigmaGaussnoPID));
       if (sigmaGaussnoPID > refwidth - 0.0003 && sigmaGaussnoPID < refwidth + 0.0003) pave2->AddText("#color[3]{OK!! The width is compatible with standard.}");
       else                                                                            pave2->AddText("#color[2]{NOT OK!! Problem.}");
       pave2->Draw();
       cout<<"   "<<refwidth - 0.0003<<"<"<<sigmaGaussnoPID<<"<"<<refwidth + 0.0003<<endl;
     c4->SaveAs("fig_lf_Multistrange.pdf");


     //________________________________ 
     //DEFINE 6th CANVAS AND DRAW PLOTS
     TCanvas *c6 = new TCanvas("c6","",1200,400);//1200,270);
     c6->Divide(2,1);
       //Pad 1: Event Selection
       c6->cd(1);
       TH1D* hvar26 = hlist->FindObject(Form("fHistEventSel"));
       hvar26->Draw("histo");
       //Pad 2: Centrality/Multiplicity distribution
       c6->cd(2);
       TH1D *hvar16 = cf->ShowProjection(19,icasType);
       hvar16->Draw("histo");
     if      (!isMC) c6->SaveAs("fig_lf_Multistrange.pdf)");
     else if (isMC) c6->SaveAs("fig_lf_Multistrange.pdf");


     //_______________________________
     //CHECK ON MONTE CARLO PRODUCTION
     if (isMC) { 
           
            AliCFContainer *cfMC = (AliCFContainer*) (f1->Get("PWGLFStrangeness.outputCheckCascade/fCFContCascadeMCgen"));
            //DEFINE 6th CANVAS AND DRAW PLOTS
            TCanvas *c7 = new TCanvas("c7","",1200,800);
            c7->Divide(2,3);
            //Pad 1: Total Momentum
            c7->cd(1);
            TH1D *hvar17 = cfMC->ShowProjection(0,icasType);
            hvar17->Draw("histo");
            tcasc->Draw();
            //Pad 2: Transverse Momentum
            c7->cd(2);
            TH1D *hvar18 = cfMC->ShowProjection(1,icasType);
            hvar18->Draw("histo");
            //Pad 3: Rapidity (y)
            c7->cd(3);
            TH1D *hvar19 = cfMC->ShowProjection(2,icasType);
            hvar19->Draw("histo");
            //Pad 4: Pseudo-rapidity (eta)
            c7->cd(4);
            TH1D *hvar20 = cfMC->ShowProjection(3,icasType);
            hvar20->Draw("histo");
            //Pad 5: Theta
            c7->cd(5);
            TH1D *hvar21 = cfMC->ShowProjection(4,icasType);
            hvar21->Draw("histo");
            //Pad 6: Phi
            c7->cd(6);
            TH1D *hvar22 = cfMC->ShowProjection(5,icasType);
            hvar22->Draw("histo");
            c7->SaveAs("fig_lf_Multistrange.pdf)");
     }


}



//______________________
Bool_t checkUnderTheLimit(TH1D *lHist, Double_t limit) {

      Int_t binlimit = lHist->FindBin(limit);
      Bool_t checkOk = kTRUE;
      for (Int_t i = 1; i < binlimit; i++) {
           Int_t content = 0;
           content = lHist->GetBinContent(i);
           if (content != 0) checkOk = kFALSE;
           //cout<<"Content bin "<<i<<": "<<content<<endl;
      }
      return checkOk;

}

//______________________
Bool_t checkOverTheLimit(TH1D *lHist, Double_t limit) {

      Int_t binlimit = lHist->FindBin(limit);
      Int_t lastbin  = lHist->GetNbinsX();
      Bool_t checkOk = kTRUE;
      for (Int_t i = binlimit; i < lastbin+1; i++) {
           Int_t content = 0;
           content = lHist->GetBinContent(i);
           if (content != 0) checkOk = kFALSE;
           //cout<<"Content bin "<<i<<": "<<content<<endl;
      }
      return checkOk;

}