TENDER becomes Tender
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / AnalysisBoth.C
index 595ccaf..d1a198e 100644 (file)
@@ -1,9 +1,11 @@
+
 class AliSpectraBothHistoManager;
 class AliSpectraBothEventCuts; 
 class AliSpectraBothTrackCuts;
 TString Charge[]={"Pos","Neg"};
 TString Sign[]={"Plus","Minus"};
 TString Particle[]={"Pion","Kaon","Proton"};
+TString symboles[]={"#pi^{+}","K^{+}","p","pi^{-}","K^{-}","#bar{p}"}; 
 AliSpectraBothHistoManager* managerdata=0x0;
 AliSpectraBothEventCuts* ecutsdata=0x0; 
 AliSpectraBothTrackCuts* tcutsdata=0x0;
@@ -17,22 +19,49 @@ Int_t Color[3]={1,2,4};
 Int_t Marker[6]={20,21,22,24,25,26};
 Double_t Range[3]={0.3,0.3,0.5}; // LowPt range for pi k p
 
+
+Double_t FitRange[2]={-2.0,2.0};
+Double_t CutRange[2]={-3.0,3.0};
+Double_t minptformaterial[6]={0.0,0.2,0.0,0.0,0.2,0.0};
+Double_t maxptformaterial[6]={0.0,0.6,1.3,0.0,0.6,0.0};
+Double_t minptforWD[6]={0.2,100.0,0.3,0.2,100.0,0.3};
+Double_t maxptforWD[6]={1.5,-100.0,2.0,1.5,-100.0,2.0};
+Double_t minRanges[3]={0.3,0.3,0.45};
+Double_t maxRanges[3]={1.5,1.2,2.2};
+Double_t TOFPIDsignalmatching[]={-1.0,-1.0,-1.0};
+Double_t fMaxContaminationPIDMC=0.2;
+TString filenames[]={"eff.root","pid.root","sec.root"};
+
 enum ECharge_t {
   kPositive,
   kNegative,
   kNCharges
 };
 enum {
- kdodca=0x1,
- kgeantflukaKaon=0x2,
- kgeantflukaProton=0x4,
- knormalizationtoeventspassingPhySel=0x8,
- kveretxcorrectionandbadchunkscorr=0x10
+ kdodca=0x1, //dca fits are made 
+ kgeantflukaKaon=0x2,// geant fluka correction is used for kaons 
+ kgeantflukaProton=0x4, // geant fluka correction is used for protons and antiprotons
+ knormalizationtoeventspassingPhySel=0x8,// spectra are divided by number of events passing physic selection   
+ kveretxcorrectionandbadchunkscorr=0x10, // correction for difference in z vertex distribution in data and MC and correction for bad chunks is applied
+ kmcisusedasdata=0x20, // the result of the looping over MC is used as data input 
+ kdonotusedcacuts=0x40, // allows to use the constant dca cut for all pt bins not the pt dependet defined in stardrad track cuts 2011
+ kuseprimaryPIDcont=0x80, //pid contamination is calculated using only primiary particle in this case K should use dca fits 
+ knormalizationwithbin0integralsdata=0x100, // the normalization factor is calcualte using integral over z vertex distributions (in this case reconstructed vertex disitrbution uses z vertex for data) 
+ knormalizationwithbin0integralsMC=0x200, //in this case reconstructed vertex disitrbution uses z vertex for data, those to options will be use only if knormalizationtoeventspassingPhySel is not set
+ kuserangeonfigfile=0x400, // use of config file for dca fit settings
+ kskipconcutonspectra=0x800, //do not use conPID<02 cut  useful for syst. studies
+ kuseTOFmatchingcorrection=0x1000, // if set tof matching correction is applied.
+ kuseTOFcorrforPIDsignalmatching=0x2000, // rescale the for spectra by the factor given in config files        
+ kuseeffcorrectionfromfile=0x4000, //use the efficiency from the file specfied in config file                  
+ kusePIDcontaminatiofromfile=0x8000, //use the PID contamination from the file specfied in config file
+ kuseseccontaminatiofromfile=0x10000 //use the secondary contamination from the file specfied in config file
+                                                       
 };     
 
-Bool_t OpenFile(TString dirname, TString outputname, Bool_t mcflag);
-void AnalysisBoth (UInt_t options=0xF,TString outdate, TString outnamedata, TString outnamemc="" )
+Bool_t OpenFile(TString dirname, TString outputname, Bool_t mcflag,Bool_t mcasdata=false);
+void AnalysisBoth (UInt_t options=0xF,TString outdate, TString outnamedata, TString outnamemc="",TString configfile="" )
 {
+gStyle->SetOptStat(0); 
        TH1::AddDirectory(kFALSE);
        gSystem->Load("libCore.so");  
        gSystem->Load("libPhysics.so");
@@ -44,7 +73,7 @@ void AnalysisBoth (UInt_t options=0xF,TString outdate, TString outnamedata, TStr
        gSystem->Load("libANALYSIS");
        gSystem->Load("libOADB");
        gSystem->Load("libANALYSISalice");
-       gSystem->Load("libTENDER");
+       gSystem->Load("libTender");
        gSystem->Load("libCORRFW");
        gSystem->Load("libPWGTools");
        gSystem->Load("libPWGLFspectra");
@@ -55,13 +84,20 @@ void AnalysisBoth (UInt_t options=0xF,TString outdate, TString outnamedata, TStr
        mass[1]   = TDatabasePDG::Instance()->GetParticle("K+")->Mass();
        mass[2] = TDatabasePDG::Instance()->GetParticle("proton")->Mass();
 
-       AliESDtrackCuts* esdtrackcuts= AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE);
-       TString formulastring(esdtrackcuts->GetMaxDCAToVertexXYPtDep());
-       formulastring.ReplaceAll("pt","x");
-       TFormula* dcacutxy=new TFormula("dcacutxy",formulastring.Data());
-
+       TFormula* dcacutxy=0x0;
+       if(!(options&kdonotusedcacuts))
+       {
+       
+               AliESDtrackCuts* esdtrackcuts= AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE);
+               TString formulastring(esdtrackcuts->GetMaxDCAToVertexXYPtDep());        
+               formulastring.ReplaceAll("pt","x");
+               dcacutxy=new TFormula("dcacutxy",formulastring.Data());
+       }
+       if(options&kuserangeonfigfile)
+               if(!ReadConfigFile(configfile))
+                       return;         
        TList* lout=new TList();
-
+       
 
        TString indirname=Form("/output/train%s",outdate.Data());
        //TString indirname("/output/train24072012");
@@ -72,7 +108,7 @@ void AnalysisBoth (UInt_t options=0xF,TString outdate, TString outnamedata, TStr
 
 
        OpenFile(indirname,outnamemc,true);
-       OpenFile(indirname,outnamedata,false);
+       OpenFile(indirname,outnamedata,false,((Bool_t)(options&kmcisusedasdata)));
        if(!managermc||!managerdata)
        {
                cout<<managermc<<" "<<managerdata<<endl;
@@ -87,6 +123,8 @@ void AnalysisBoth (UInt_t options=0xF,TString outdate, TString outnamedata, TStr
        TH1F* contWD[6];
        TH1F* contMat[6];
        TH1F* confinal[6];
+       TH1F* contPIDpri[6];
+       TH1F* contSecMC[6];
        
        TH1F* contfit[12];
        TH1F* contWDfit[12];
@@ -98,7 +136,9 @@ void AnalysisBoth (UInt_t options=0xF,TString outdate, TString outnamedata, TStr
        TH1F* spectra[6];
        TH1F* spectraLeonardo[6];
        
-       TH1F* corrLeonardo[6]; 
+       TH1F* corrLeonardo[6];
+        TH1F* doubleconuntsdata[3]; 
+        TH1F* doubleconuntsMC[3]; 
        //GetSpectra(managerdata,rawspectradata,true);
        //GetSpectra(managermc,rawspectramc,true,true);
        
@@ -108,31 +148,59 @@ void AnalysisBoth (UInt_t options=0xF,TString outdate, TString outnamedata, TStr
        GetPtHistFromPtDCAhisto("hHistPtRecTrue","conPID",managermc,contPID,dcacutxy);
        GetPtHistFromPtDCAhisto("hHistPtRecSigmaSecondaryWeakDecay","conWD",managermc,contWD,dcacutxy);
        GetPtHistFromPtDCAhisto("hHistPtRecSigmaSecondaryMaterial","conMat",managermc,contMat,dcacutxy);
+       GetPtHistFromPtDCAhisto("hHistPtRecSigmaPrimary","conPIDprimary",managermc,contPIDpri,dcacutxy);
+
        
-       
-//     Double_t neventsdata =  ecutsdata->NumberOfEvents();
        Double_t neventsmcall = 1 ;  //if loop over MC is done after or befor events cuts this will be changed 
        Double_t neventsdata =  1;
        Double_t neventsmc =  1;
 
+       //Normaliztion of MCtruth depends if the loop was done after of before ESD event cuts.
+       //In currect code this cannot be check on the level of macro.
+       //If the loop was done before MC should be done to all processed events (NumberOfProcessedEvents())
+       //If loop was done after MC should be normalized to all accepted events (NumberOfEvents()) 
+       // The option one will be alaways use.
+       
+       neventsmcall= ecutsmc->NumberOfProcessedEvents();
        if(options&knormalizationtoeventspassingPhySel)
        {
-               neventsmcall= ecutsmc->NumberOfProcessedEvents();
+               //neventsmcall= ecutsmc->NumberOfProcessedEvents();
                 neventsdata=ecutsdata->NumberOfPhysSelEvents();
                 neventsmc=ecutsmc->NumberOfPhysSelEvents();
        }
+       else if ((options&knormalizationwithbin0integralsdata)||(options&knormalizationwithbin0integralsMC))
+       {
+               neventsdata=Normaliztionwithbin0integrals(options);
+               neventsmc=ecutsmc->NumberOfPhysSelEvents();
+       }
        else
        {
-               neventsdata=ecutsdata->NumberOfEvents();
+               neventsdata=ecutsdata->NumberOfEvents(); //number of accepted events
                 neventsmc=ecutsmc->NumberOfEvents();
                neventsmcall= ecutsmc->NumberOfEvents();
 
-
        }
        GetMCTruth(MCTruth);
-       
-       
-       
+       cout<<neventsdata<<" Events"<<endl;
+        cout<< neventsmc<<" Events "<<endl;
+       if(neventsdata<1)
+       {
+               cout<<"No events DATA"<<endl;
+               return;
+       }
+       if(neventsmc<1&&(((options&kusePIDcontaminatiofromfile)!=kusePIDcontaminatiofromfile)||((options&kuseseccontaminatiofromfile)!=kuseseccontaminatiofromfile)||((options&kuseseccontaminatiofromfile)!=kuseseccontaminatiofromfile)))
+       {
+               cout<<"No events MC"<<endl;
+               return;
+       }
+       else if(neventsmc<1)
+       {
+               neventsmc=1; //dumpy normalization
+               neventsmcall=1; //dumpy normalization
+
+       }
+
+       cout<<"A "<<neventsmc<<endl;
        TH1F* allgen=((TH1F*)managermc->GetPtHistogram1D("hHistPtGen",1,1))->Clone();
        allgen->SetName("AllGen");
        TH1F* allrecMC=GetOneHistFromPtDCAhisto("hHistPtRec","rawallMC",managermc,dcacutxy);
@@ -143,8 +211,11 @@ void AnalysisBoth (UInt_t options=0xF,TString outdate, TString outnamedata, TStr
        
        
        TH1F* spectraall=(TH1F*)allrecdata->Clone("recNch");
+       spectraall->SetTitle("recNch");
        TH1F* contall=(TH1F*)allrecMC->Clone("contall");
-       contall->Add(alleff,-1);
+       contall->SetTitle("contall");
+       //contall->Add(alleff,-1);
+       SubHistWithFullCorr(contall,alleff);
        alleff->Divide(alleff,allgen,1,1,"B");
        contall->Divide(contall,allrecMC,1,1,"B");
        
@@ -162,24 +233,36 @@ void AnalysisBoth (UInt_t options=0xF,TString outdate, TString outnamedata, TStr
        lout->Add(alleff);
        lout->Add(allrecdata);
        lout->Add(spectraall);
-       lout->Add(contall);     
+       lout->Add(contall);
+       
        for (int i=0;i<6;i++)
        {
        
                
                TString tmpname(rawspectramc[i]->GetTitle());
                tmpname.ReplaceAll("SpectraMC","%s");
-               contallMC[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contallMC")); 
+               contallMC[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contallMC"));
+               contallMC[i]->SetTitle(Form(tmpname.Data(),"contallMC"));
                contfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contfit"));
+               contfit[i]->SetTitle(Form(tmpname.Data(),"contfit"));
                contWDfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contWDfit"));
+               contWDfit[i]->SetTitle(Form(tmpname.Data(),"contWDfit"));
                contMatfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contMatfit"));
+               contMatfit[i]->SetTitle(Form(tmpname.Data(),"contMatfit"));
                primaryfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"primaryfit"));
-               
+               primaryfit[i]->SetTitle(Form(tmpname.Data(),"primaryfit"));
                contfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contfitonMC"));
+               contfit[i+6]->SetTitle(Form(tmpname.Data(),"contfitonMC"));
                contWDfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contWDfitonMC"));
+               contWDfit[i+6]->SetTitle(Form(tmpname.Data(),"contWDfitonMC"));
                contMatfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contMatfitonMC"));
+               contMatfit[i+6]->SetTitle(Form(tmpname.Data(),"contMatfitonMC"));
                primaryfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"primaryfitMC"));
+               primaryfit[i+6]->SetTitle(Form(tmpname.Data(),"primaryfitMC"));
+               
+       
                
+                       
                contfit[i]->Reset();
                contWDfit[i]->Reset();
                contMatfit[i]->Reset();
@@ -194,37 +277,78 @@ void AnalysisBoth (UInt_t options=0xF,TString outdate, TString outnamedata, TStr
                SetBintoOne(primaryfit[i]);
                SetBintoOne(primaryfit[i+6]);           
                spectra[i]=(TH1F*)rawspectradata[i]->Clone(Form(tmpname.Data(),"SpectraFinal"));
-               
+               spectra[i]->SetTitle(Form(tmpname.Data(),"SpectraFinal"));
                spectraLeonardo[i]=(TH1F*)rawspectradata[i]->Clone(Form(tmpname.Data(),"SpectraFinalLeonardo"));
+               spectraLeonardo[i]->SetTitle(Form(tmpname.Data(),"SpectraFinalLeonardo"));
+
                corrLeonardo[i]=(TH1F*)MCTruth[i]->Clone(Form(tmpname.Data(),"CorrFactLeonardo"));
-               
+               corrLeonardo[i]->SetTitle(Form(tmpname.Data(),"CorrFactLeonardo"));                     
                corrLeonardo[i]->Divide(corrLeonardo[i],rawspectramc[i],1,1,"B");
                
                
-               
-               contallMC[i]->Add(eff[i],-1.0);
-               RecomputeErrors(contallMC[i]);
-               contallMC[i]->Sumw2(); 
-               contallMC[i]->Divide(contallMC[i],rawspectramc[i],1,1,"B");
-               
-               eff[i]->Divide(eff[i],MCTruth[i],1,1,"B");
+               if((options&kusePIDcontaminatiofromfile)==kusePIDcontaminatiofromfile)
+               {
+                       CopyCorrectionFromFile(filenames[1],"conPID",contPID);
+                       CopyCorrectionFromFile(filenames[1],"conPIDprimary",contPIDpri);
+
+               }
+               else
+               {
+                       //contallMC[i]->Add(eff[i],-1.0);
+                       SubHistWithFullCorr(contallMC[i],eff[i]);
+                       //RecomputeErrors(contallMC[i]);
+                       contallMC[i]->Sumw2(); 
+                       contallMC[i]->Divide(contallMC[i],rawspectramc[i],1,1,"B");
+                       // contamintaion from PID but only primaries
+                       //contPIDpri[i]->Add(eff[i],-1.0);
+                       SubHistWithFullCorr(contPIDpri[i],eff[i]);
+                       //RecomputeErrors(contPIDpri[i]);
+                       contPIDpri[i]->Divide(contPIDpri[i],rawspectramc[i],1,1,"B");
+               }
+
+               if((options&kuseeffcorrectionfromfile)==kuseeffcorrectionfromfile)
+                       CopyCorrectionFromFile(filenames[0],"eff",eff);
+               else
+                       eff[i]->Divide(eff[i],MCTruth[i],1,1,"B");
                
                
                contPID[i]->Sumw2();
                rawspectramc[i]->Sumw2();
-               contPID[i]->Add(contPID[i],rawspectramc[i],-1,1);
-               RecomputeErrors(contPID[i]);
+               //contPID[i]->Add(contPID[i],rawspectramc[i],-1,1);
+               SubHistWithFullCorr(contPID[i],rawspectramc[i]);
+               contPID[i]->Scale(-1.0);
+
+               //RecomputeErrors(contPID[i]);
                contPID[i]->ResetStats();
                contPID[i]->Sumw2();
                contPID[i]->Divide(contPID[i],rawspectramc[i],1,1,"B");
-               
-               confinal[i]=(TH1F*)contPID[i]->Clone(Form(tmpname.Data(),"confinal"));
-
-
+       
+               if((options&kuseprimaryPIDcont)==kuseprimaryPIDcont)
+                       confinal[i]=(TH1F*)contPIDpri[i]->Clone(Form(tmpname.Data(),"confinal"));
+               else    
+                       confinal[i]=(TH1F*)contPID[i]->Clone(Form(tmpname.Data(),"confinal"));
+               confinal[i]->SetTitle(Form(tmpname.Data(),"confinal"));
+
+               contSecMC[i]=(TH1F*)contWD[i]->Clone(Form(tmpname.Data(),"conSecMC"));
+               contSecMC[i]->Add(contMat[i]);
                contWD[i]->Divide(contWD[i],rawspectramc[i],1,1,"B");
                contMat[i]->Divide(contMat[i],rawspectramc[i],1,1,"B");
+               contSecMC[i]->Divide(contSecMC[i],rawspectramc[i],1,1,"B");
+               if(i>2)
+               {
+                       doubleconuntsdata[i-3]=(TH1F*)rawspectradata[i]->Clone(Form("DoublecountsDATA%s",Particle[i-3].Data()));
+                       doubleconuntsdata[i-3]->Reset();
+                       doubleconuntsMC[i-3]=(TH1F*)rawspectramc[i]->Clone(Form("DoublecountsMC%s",Particle[i-3].Data()));
+                       doubleconuntsMC[i-3]->Reset();
+
+                       CalculateDoubleCounts(doubleconuntsdata[i-3],rawspectradata,i-3,kTRUE);
+                       CalculateDoubleCounts(doubleconuntsMC[i-3],rawspectramc,i-3,kFALSE);
+
+                       
        
+               }
        
+
        
                rawspectradata[i]->Scale(1./neventsdata,"width");
                rawspectramc[i]->Scale(1./neventsmc,"width");
@@ -244,21 +368,49 @@ void AnalysisBoth (UInt_t options=0xF,TString outdate, TString outnamedata, TStr
                lout->Add(contfit[i]);
                lout->Add(contWDfit[i]);
                lout->Add(contMatfit[i]);
+               lout->Add(primaryfit[i]);       
                lout->Add(contfit[i+6]);
                lout->Add(contWDfit[i+6]);
                lout->Add(contMatfit[i+6]);
+               lout->Add(primaryfit[i+6]);     
                lout->Add(spectra[i]);
                lout->Add(spectraLeonardo[i]);
                lout->Add(confinal[i]);
+               lout->Add(contPIDpri[i]);
+               lout->Add(contSecMC[i]);
+               if(i>2)
+               {
+                       lout->Add(doubleconuntsdata[i-3]);
+                       lout->Add(doubleconuntsMC[i-3]);
+               }
        }
-       TFile* fout=new TFile(Form("./results/ResMY_%s_%s.root",outnamemc.Data(),outdate.Data()),"RECREATE");
-       if (options&kdodca)
+       outdate.ReplaceAll("/","_");
+       configfile.ReplaceAll(".","_");
+       TFile* fout=0x0;
+       if(configfile.Length()>0&&(options&kuserangeonfigfile))
+               fout=new TFile(Form("./results/ResMY_%s_%s_%#X_%s.root",outnamemc.Data(),outdate.Data(),options,configfile.Data()),"RECREATE");
+       else
+               fout=new TFile(Form("./results/ResMY_%s_%s_%#X.root",outnamemc.Data(),outdate.Data(),options),"RECREATE");
+       if (((options&kdodca)==kdodca)&&((options&kuseseccontaminatiofromfile)!=kuseseccontaminatiofromfile))
                DCACorrectionMarek(managerdata,managermc,dcacutxy,fout,contfit,contWDfit,contMatfit,primaryfit);
+       else if ((options&kuseseccontaminatiofromfile)==kuseseccontaminatiofromfile)
+       {
+               CopyCorrectionFromFile(filenames[2],"contfit",contfit);
+               CopyCorrectionFromFile(filenames[2],"contWDfit",contWDfit);
+               CopyCorrectionFromFile(filenames[2],"contMatfit",contMatfit);
+               CopyCorrectionFromFile(filenames[2],"primaryfit",primaryfit);
+
+       }
+       else 
+               cout<<"Secondary from MC"<<endl;
+               
+
        for (int i=0;i<6;i++)
        {
-                       if(options&kdodca)
+                       if(((options&kdodca)==kdodca)||((options&kuseseccontaminatiofromfile)==kuseseccontaminatiofromfile))
                        {
-                               confinal[i]->Add(contfit[i]);
+                               if((options&kuseprimaryPIDcont)||(i!=1&&i!=4)) //if we do not use cont PId only for primary  and this is a kaon that do not use fit
+                                       confinal[i]->Add(contfit[i]);
                                GetCorrectedSpectra(spectra[i],rawspectradata[i],eff[i],confinal[i]);
                        }
                        else
@@ -266,36 +418,72 @@ void AnalysisBoth (UInt_t options=0xF,TString outdate, TString outnamedata, TStr
                                GetCorrectedSpectra(spectra[i],rawspectradata[i],eff[i],contallMC[i]);  
                        }
                        GetCorrectedSpectraLeonardo(spectraLeonardo[i],corrLeonardo[i],primaryfit[i],primaryfit[i+6]);
-                       CleanHisto(spectra[i],-1,100,contPID[i]);
-                       CleanHisto(spectraLeonardo[i],-1,100,contPID[i]);                               
+                       if(options&kskipconcutonspectra)
+                               continue;
+                       if(options&kuseprimaryPIDcont)
+                       {
+                               CleanHisto(spectra[i],-1,100,contPIDpri[i]);
+                               CleanHisto(spectraLeonardo[i],-1,100,contPIDpri[i]);            
+                       }
+                       else
+                       {
+                               CleanHisto(spectra[i],-1,100,contPID[i]);
+                               CleanHisto(spectraLeonardo[i],-1,100,contPID[i]);
+                       }
+                       // Apply correction for wrongly simulated TOF signal in MC 
+                       if(options&kuseTOFcorrforPIDsignalmatching)
+                               TOFPIDsignalmatchingApply(spectra[i],TOFPIDsignalmatching[i%3]);
+                                                       
        }
        
        GFCorrection(spectra,tcutsdata->GetPtTOFMatching(),options);
        GFCorrection(spectraLeonardo,tcutsdata->GetPtTOFMatching(),options);
-        MatchingTOFEff(spectra,lout);
-         MatchingTOFEff(spectraLeonardo);
-       TH1F* allch=GetSumAllCh(spectra,mass);
-       lout->Add(allch);       
+       
+       Double_t ycut=tcutsdata->GetY();
+       Double_t etacut=tcutsdata->GetEtaMax()-tcutsdata->GetEtaMin();
+       if(etacut<0.00001)
+       {
+               cout<<"using eta window 1.6"<<endl; 
+               etacut=1.6;
+
+       }
 
+       TH1F* allch=GetSumAllCh(spectra,mass,etacut);
+       lout->Add(allch);       
+               if(options&kuseTOFmatchingcorrection)
+       {       
+               MatchingTOFEff(spectra,lout);
+               MatchingTOFEff(spectraLeonardo);
+               TOFMatchingForNch(spectraall);
+       }
 //     lout->ls();
        fout->cd();     
        TList* listqa=new TList();
        TList* canvaslist=new TList();
        Float_t vertexcorrection=1.0;
-       Float_t corrbadchunksvtx=QAPlotsBoth(managerdata,managermc,ecutsdata,ecutsmc,tcutsdata,tcutsmc,listqa,canvaslist);
+       Float_t corrbadchunksvtx=1.0;
+       if ((options&knormalizationwithbin0integralsdata)||(options&knormalizationwithbin0integralsMC))
+               corrbadchunksvtx=QAPlotsBoth(managerdata,managermc,ecutsdata,ecutsmc,tcutsdata,tcutsmc,listqa,canvaslist,0);
+       else
+               corrbadchunksvtx=QAPlotsBoth(managerdata,managermc,ecutsdata,ecutsmc,tcutsdata,tcutsmc,listqa,canvaslist,1);
        if (options&kveretxcorrectionandbadchunkscorr)
                vertexcorrection=corrbadchunksvtx;
        cout<<" VTX corr="<<vertexcorrection<<endl;
-       Double_t ycut=tcutsdata->GetY();
        if(TMath::Abs(ycut)>0.0)
                vertexcorrection=vertexcorrection/(2.0*ycut);
        for (int i=0;i<6;i++)
        {
                spectra[i]->Scale(vertexcorrection);
                spectraLeonardo[i]->Scale(vertexcorrection);
+               if(TMath::Abs(ycut)>0.0)
+               {
+                       rawspectradata[i]->Scale(1.0/(2.0*ycut));
+                       rawspectramc[i]->Scale(1.0/(2.0*ycut));
+                       MCTruth[i]->Scale(1.0/(2.0*ycut));
+               }
        }       
        allch->Scale(vertexcorrection);
-       spectraall->Scale(vertexcorrection/1.6);
+       spectraall->Scale(vertexcorrection/etacut);
 
        //spectraall->Scale(1.0/1.6);
        lout->Write("output",TObject::kSingleKey);      
@@ -303,11 +491,13 @@ void AnalysisBoth (UInt_t options=0xF,TString outdate, TString outnamedata, TStr
        canvaslist->Write("outputcanvas",TObject::kSingleKey);
 
        fout->Close();
+       //Normaliztionwithbin0integrals();
 
 }
 
-Bool_t   OpenFile(TString dirname,TString outputname, Bool_t mcflag)
+Bool_t   OpenFile(TString dirname,TString outputname, Bool_t mcflag, Bool_t mcasdata)
 {
+       
 
        TString nameFile = Form("./%s/AnalysisResults%s.root",dirname.Data(),(mcflag?"MC":"DATA"));
        TFile *file = TFile::Open(nameFile.Data());
@@ -317,12 +507,23 @@ Bool_t   OpenFile(TString dirname,TString outputname, Bool_t mcflag)
                return false;
        }       
        TString sname=Form("OutputBothSpectraTask_%s_%s",(mcflag?"MC":"Data"),outputname.Data());
+       if(mcasdata)
+       {
+               cout<<"using MC as data "<<endl;
+               sname=Form("OutputBothSpectraTask_%s_%s","MC",outputname.Data());
+       }
        file->ls();
        TDirectoryFile *dir=(TDirectoryFile*)file->Get(sname.Data());
        if(!dir)
        {
-       //      cout<<"no dir "<<sname.Data()<<endl;
-               sname=Form("OutputAODSpectraTask_%s_%s",(mcflag?"MC":"Data"),outputname.Data());
+       //      cout<<"no dir "<<sname.Data()<<endl;    
+               if(mcasdata)
+               {
+                       cout<<"using MC as data "<<endl;
+                       sname=Form("OutputAODSpectraTask_%s_%s","MC",outputname.Data());
+               }
+               else    
+                       sname=Form("OutputAODSpectraTask_%s_%s",(mcflag?"MC":"Data"),outputname.Data());
        //      cout<<"trying "<<sname.Data()<<endl;
                dir=(TDirectoryFile*)file->Get(sname.Data());
                if(!dir)
@@ -341,6 +542,8 @@ Bool_t   OpenFile(TString dirname,TString outputname, Bool_t mcflag)
                tcutsmc->PrintCuts();
                if(!managermc||!ecutsmc||!tcutsmc)
                        return false;
+               if(managermc->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries()!=ecutsmc->GetHistoCuts()->GetBinContent(3))
+                       cout<<"Please check MC file possible problem with merging"<<" "<<managermc->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries()<<" "<<ecutsmc->GetHistoCuts()->GetBinContent(3)<<endl;
        }
        else
        {
@@ -351,6 +554,9 @@ Bool_t   OpenFile(TString dirname,TString outputname, Bool_t mcflag)
                tcutsdata->PrintCuts();
                if(!managerdata||!ecutsdata||!tcutsdata)
                        return false;
+               if(managerdata->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries()!=ecutsdata->GetHistoCuts()->GetBinContent(3))
+                       cout<<"Please check DATA file possible problem with merging"<<" "<<managerdata->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries()<<" "<<ecutsdata->GetHistoCuts()->GetBinContent(3)<<endl;
+
        }
        return true;
 }
@@ -383,12 +589,25 @@ TH1F* GetOneHistFromPtDCAhisto(TString name,TString hnameout,AliSpectraBothHisto
                                {
                                        Double_t lowedge=histo->GetBinLowEdge(ibin);
                                        Float_t cut=dcacutxy->Eval(lowedge);
-                                       TH1F* dcahist=(TH1F*)hman->GetDCAHistogram1D(name.Data(),lowedge,lowedge));
-                                       Float_t inyield=dcahist->Integral(dcahist->GetXaxis()->FindBin(-1.0*cut),dcahist->GetXaxis()->FindBin(cut));
-                                       cout<<"corr data "<<histo->GetBinContent(ibin)<<" "<<inyield<<" "<<dcahist->Integral()<<" "<<hnameout.Data()<<endl;
-                                       cout<<"test dca "<<lowedge<<" "<<dcacutxy->Eval(lowedge)<<" "<<dcacutxy->Eval(histo->GetXaxis()->GetBinUpEdge(ibin))<<" "<<dcahist->GetBinLowEdge(dcahist->GetXaxis()->FindBin(-1.0*cut))<<" "<<dcahist->GetXaxis()->GetBinUpEdge(dcahist->GetXaxis()->FindBin(-1.0*cut))<<endl;
-                                       histo->SetBinContent(ibin,inyield);
-                                       histo->SetBinError(ibin,TMath::Sqrt(inyield));
+                                       TH1F* dcahist=(TH1F*)hman->GetDCAHistogram1D(name.Data(),lowedge,lowedge);
+//                                     Float_t inyield=dcahist->Integral(dcahist->GetXaxis()->FindBin(-1.0*cut),dcahist->GetXaxis()->FindBin(cut));
+                                       Float_t testyield=0.0;
+                                       Float_t testerror=0.0;  
+                                       for (int itest=dcahist->GetXaxis()->FindBin(-1.0*cut);itest<=dcahist->GetXaxis()->FindBin(cut);itest++)
+                                       {
+                                               testyield+=dcahist->GetBinContent(itest);
+                                               testerror+=dcahist->GetBinError(itest)*dcahist->GetBinError(itest);
+                                       }
+//                                     cout<<"corr data "<<histo->GetBinContent(ibin)<<" "<<inyield<<" "<<dcahist->Integral()<<" "<<hnameout.Data()<<endl;
+//                                     cout<<"test dca "<<lowedge<<" "<<dcacutxy->Eval(lowedge)<<" "<<dcacutxy->Eval(histo->GetXaxis()->GetBinUpEdge(ibin))<<" "<<dcahist->GetBinLowEdge(dcahist->GetXaxis()->FindBin(-1.0*cut))<<" "<<dcahist->GetXaxis()->GetBinUpEdge(dcahist->GetXaxis()->FindBin(-1.0*cut))<<endl;
+
+//                                     cout<<testyield<<" "<<TMath::Sqrt(testerror)<<" error2 "<<inyield<<" "<<TMath::Sqrt(inyield)<<endl;
+                                               
+                                       //histo->SetBinContent(ibin,inyield);
+                                       //histo->SetBinError(ibin,TMath::Sqrt(inyield));
+                                       histo->SetBinContent(ibin,testyield);
+                                       histo->SetBinError(ibin,TMath::Sqrt(testerror));
+
                                }
                        }
                        histo->Sumw2();
@@ -400,39 +619,20 @@ TH1F* GetOneHistFromPtDCAhisto(TString name,TString hnameout,AliSpectraBothHisto
        
 void GetPtHistFromPtDCAhisto(TString hnamein, TString hnameout, AliSpectraBothHistoManager* hman,TH1F** histo,TFormula* dcacutxy)
 {
-       Float_t min[3]={0.3,0.3,0.4};
-       Float_t max[3]={1.5,1.2,2.2};
+       //Float_t min[3]={0.3,0.3,0.45};
+       //Float_t max[3]={1.5,1.2,2.2};
        for(Int_t icharge=0;icharge<2;icharge++)
        {
                for(Int_t ipart=0;ipart<3;ipart++)
                {
                        Int_t index=ipart+3*icharge;
-                       //TString hname=Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data());
-                       printf("Getting %s",hnamein.Data());
+                       Printf("Getting %s",hnamein.Data());
                        TString nameinfinal=Form("%s%s%s",hnamein.Data(),Particle[ipart].Data(),Sign[icharge].Data());
                        TString nameoutfinal=Form("%s%s%s",hnameout.Data(),Particle[ipart].Data(),Sign[icharge].Data());
                        
                        
                        histo[index]=GetOneHistFromPtDCAhisto(nameinfinal,nameoutfinal,hman,dcacutxy);
-                       /*histo[index] =(TH1F*)((TH1F*) hman->GetPtHistogram1D(Form("%s%s%s",hnamein.Data(),Particle[ipart].Data(),Sign[icharge].Data()),-1,-1))->Clone();
-                       histo[index]->SetName(Form("%s%s%s",hnameout.Data(),Particle[ipart].Data(),Sign[icharge].Data()));
-                       histo[index]->SetTitle(Form("%s%s%s",hnameout.Data(),Particle[ipart].Data(),Sign[icharge].Data()));
-                 
-                       if(dcacutxy)
-                       {
-                               for(int ibin=1;ibin<histo[index]->GetNbinsX();ibin++)
-                               {
-                                       Double_t lowedge=histo[index]->GetBinLowEdge(ibin);
-                                       Float_t cut=dcacutxy->Eval(lowedge);
-                                       TH1F* dcahist=(TH1F*)hman->GetDCAHistogram1D(Form("%s%s%s",hnamein.Data(),Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge));
-                                       Float_t inyield=dcahist->Integral(dcahist->GetXaxis()->FindBin(-1.0*cut),dcahist->GetXaxis()->FindBin(cut));
-                                       cout<<"corr data "<<histo[index]->GetBinContent(ibin)<<" "<<inyield<<" "<<dcahist->Integral()<<endl;
-                                       histo[index]->SetBinContent(ibin,inyield);
-                                       histo[index]->SetBinError(ibin,TMath::Sqrt(inyield));
-                               }
-                       }*/
-                       CleanHisto(histo[index],min[ipart],max[ipart]);
-                       // histo[index]->Sumw2();
+                       CleanHisto(histo[index],minRanges[ipart],maxRanges[ipart]);
                }
        } 
 }
@@ -447,7 +647,7 @@ void CleanHisto(TH1F* h, Float_t minV, Float_t maxV,TH1* contpid=0x0)
                }       
                if(contpid)
                {
-                       if(contpid->GetBinContent(i)>0.2)
+                       if(contpid->GetBinContent(i)>fMaxContaminationPIDMC)
                        {
                                h->SetBinContent(i,0);
                                h->SetBinError(i,0);
@@ -460,18 +660,9 @@ void CleanHisto(TH1F* h, Float_t minV, Float_t maxV,TH1* contpid=0x0)
 void DCACorrectionMarek(AliSpectraBothHistoManager* hman_data, AliSpectraBothHistoManager* hman_mc,TFormula* fun,TFile *fout,TH1F** hcon,TH1F** hconWD,TH1F** hconMat,TH1F** hprimary)
 {
   printf("\n\n-> DCA Correction");  
-  Double_t FitRange[2]={-2.0,2.0};
-  Double_t CutRange[2]={-1.0,1.0};
-  Double_t minptformaterial[6]={0.0,100.0,0.0,0.0,100.0,0.0};
-  Double_t maxptformaterial[6]={0.0,-100.0,1.3,0.0,-100.0,0.0};
-  Bool_t usefit[3]={true,false,true};
-  Double_t maxbinforfit[6]={1.5,0,2.0,1.5,0,2.0};
+  
   Printf("\DCACorr");
- // TH1F *hcorrection[2];
- /* TCanvas *ccorrection=new TCanvas("DCAcorrection","DCAcorrection",700,500);
-  TCanvas *cRatiocorrection=new TCanvas("DCARatiocorrection","DCARatiocorrection",700,500);
-  cRatiocorrection->Divide(2,1);
-  ccorrection->Divide(2,1);*/
   TString sample[2]={"data","mc"};
   ofstream debug("debugDCA.txt");
   TList* listofdcafits=new TList();
@@ -483,25 +674,31 @@ void DCACorrectionMarek(AliSpectraBothHistoManager* hman_data, AliSpectraBothHis
                        for(Int_t isample=0;isample<2;isample++)
                        {
 
-                               /*hcorrection[isample]=(TH1F*)Spectra[index]->Clone();
-                               hcorrection[isample]->Reset("all");*/
                                for(Int_t ibin_data=7;ibin_data<40;ibin_data++)
                                {       
                                                
                                        Double_t lowedge=hcon[index]->GetBinLowEdge(ibin_data);
                                        Double_t binwidth=hcon[index]->GetBinWidth(ibin_data);
-                                       debug<<"NEW "<<Particle[ipart].Data()<<" "<<Sign[icharge].Data()<<" "<<lowedge<<endl;                                           
-                                       CutRange[1]=fun->Eval(lowedge);
-                                       CutRange[0]=-1.0*CutRange[1];
+                                       debug<<"NEW "<<Particle[ipart].Data()<<" "<<Sign[icharge].Data()<<" "<<lowedge<<endl;
+                                       if(fun)
+                                       {                                               
+                                               CutRange[1]=fun->Eval(lowedge);
+                                               CutRange[0]=-1.0*CutRange[1];
+                                       }       
                                        debug<<"cut  "<<CutRange[1]<<" "<<CutRange[0]<<endl;            
-                                       Bool_t useMaterial=kFALSE;
-                                       cout<<"try to fit "<<lowedge<<" "<<maxbinforfit[index]<<endl;
-                                       if(lowedge>maxbinforfit[index])
-                                               continue;
-                                       if(lowedge>minptformaterial[index]&&lowedge<maxptformaterial[index])
-                                               useMaterial=kTRUE;
-         
+                                       Short_t fitsettings=DCAfitsettings(lowedge+0.5*binwidth,index);
+                                       debug<<"settings "<< fitsettings<<endl;
+                                       if(fitsettings==0)
+                                               continue;       
+                                               
                                        TCanvas *cDCA=new TCanvas(Form("cDCA%d%s%s%sbin%d",index,sample[isample].Data(),Particle[ipart].Data(),Sign[icharge].Data(),ibin_data),Form("cDCA%d%s%s%sbin%d",index,sample[isample].Data(),Particle[ipart].Data(),Sign[icharge].Data(),ibin_data),1700,1500);
+                                       cDCA->SetMargin(0.1,0.02,0.1,0.02);
+                                       TLegend* Leg1 = new TLegend(0.6,0.6,0.85,0.85,"","NDC");
+                                       Leg1->SetFillStyle(kFALSE);
+                                       Leg1->SetLineColor(kWhite);
+                                       Leg1->SetBorderSize(0);
+
+                                       
                                        if(isample==0)
                                                TH1F *hToFit =(TH1F*) ((TH1F*)hman_data->GetDCAHistogram1D(Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone();
                                        if(isample==1)
@@ -510,51 +707,58 @@ void DCACorrectionMarek(AliSpectraBothHistoManager* hman_data, AliSpectraBothHis
                                        TH1F *hmc1=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigmaPrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone();
                                        TH1F *hmc2=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigmaSecondaryWeakDecay%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone();
                                        TH1F *hmc3=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigmaSecondaryMaterial%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone();
-                                       Double_t minentries=1;
-                                       debug<<" Entries "<<isample<<" "<<hToFit->GetEntries()<<" "<<hmc1->GetEntries()<<" "<<hmc2->GetEntries()<<" "<<hmc3->GetEntries()<<endl;
-                                       debug<<"2 Entries "<<isample<<" "<<hToFit->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))<<" "<<hmc1->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))<<" "<<hmc2->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))<<" "<<hmc3->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))<<endl;
-                                       debug<< CutRange[0]<<" "<<CutRange[1]<<" "<<lowedge<<endl;
-                                       if(hToFit->GetEntries()<=minentries || hmc1->GetEntries()<=minentries || hmc2->GetEntries()<=minentries || hmc3->GetEntries()<=minentries)
+                                       Double_t minentries=100;
+                                       debug<<hToFit->GetEntries()<<" "<<hmc1->GetEntries()<<" "<<hmc2->GetEntries()<<" "<<hmc3->GetEntries()<<endl;
+                                       debug<<((fitsettings&0x1)&&hmc2->GetEntries()<=minentries)<<" "<<((fitsettings&0x2)&&hmc3->GetEntries()<=minentries)<<endl;
+                                       cout<<hToFit->GetEntries()<<" "<<hmc1->GetEntries()<<" "<<hmc2->GetEntries()<<" "<<hmc3->GetEntries()<<endl;
+                                        cout<<((fitsettings&0x1)&&hmc2->GetEntries()<=minentries)<<" "<<((fitsettings&0x2)&&hmc3->GetEntries()<=minentries)<<endl;
+                                       if(hToFit->GetEntries()<=minentries || hmc1->GetEntries()<=minentries || ((fitsettings&0x1)&&hmc2->GetEntries()<=minentries) || ((fitsettings&0x2)&&hmc3->GetEntries()<=minentries))
                                                continue;
-                                       //Data and MC can have different stat
                                        hToFit->Sumw2();
                                        hmc1->Sumw2();
                                        hmc2->Sumw2();
                                        hmc3->Sumw2();
                                        
                                        Float_t corrforrebinning[4]={1.0,1.0,1.0,1.0};  
+                                       Int_t binCutRange[]={hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1])};                    
+
                                        
-       
                                        if(hmc3->GetNbinsX()>300)
                                        {
                                        
-                                               corrforrebinning[0]=hToFit->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]));
-                                               corrforrebinning[1]=hmc1->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]));
-                                               corrforrebinning[2]=hmc2->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]));
-                                               corrforrebinning[3]=hmc3->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]));
+                                               corrforrebinning[0]=hToFit->Integral(binCutRange[0],binCutRange[1]);
+                                               corrforrebinning[1]=hmc1->Integral(binCutRange[0],binCutRange[1]);
+                                               corrforrebinning[2]=hmc2->Integral(binCutRange[0],binCutRange[1]);
+                                               corrforrebinning[3]=hmc3->Integral(binCutRange[0],binCutRange[1]);
 
                                                hToFit->Rebin(30);
                                                hmc1->Rebin(30);
                                                hmc2->Rebin(30);
                                                hmc3->Rebin(30);
-                                               if(hToFit->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))>0.0)
-                                                       corrforrebinning[0]=corrforrebinning[0]/hToFit->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]));
+
+                                               binCutRange[0]=hmc1->GetXaxis()->FindBin(CutRange[0]);
+                                               binCutRange[1]=hmc1->GetXaxis()->FindBin(CutRange[1]);
+
+
+                                               //after rebbing we lose resolution of the dca this correction also us to do obtain inside used dca
+
+                                               if(hToFit->Integral(binCutRange[0],binCutRange[1])>0.0)
+                                                       corrforrebinning[0]=corrforrebinning[0]/hToFit->Integral(binCutRange[0],binCutRange[1]);
                                                else    
                                                        corrforrebinning[0]=1.0;
-                                               if(hmc1->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))>0.0)
-                                                       corrforrebinning[1]=corrforrebinning[1]/hmc1->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]));
+                                               if(hmc1->Integral(binCutRange[0],binCutRange[1])>0.0)
+                                                       corrforrebinning[1]=corrforrebinning[1]/hmc1->Integral(binCutRange[0],binCutRange[1]);
                                                else    
                                                        corrforrebinning[1]=1.0;
-                                               if(hmc2->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))>0.0)
-                                                       corrforrebinning[2]=corrforrebinning[2]/hmc2->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]));
+                                               if(hmc2->Integral(binCutRange[0],binCutRange[1])>0.0)
+                                                       corrforrebinning[2]=corrforrebinning[2]/hmc2->Integral(binCutRange[0],binCutRange[1]);
                                                else    
                                                        corrforrebinning[2]=1.0;
-                                               if(hmc3->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))>0.0)
-                                                       corrforrebinning[3]=corrforrebinning[3]/hmc3->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]));
+                                               if(hmc3->Integral(binCutRange[0],binCutRange[1])>0.0)
+                                                       corrforrebinning[3]=corrforrebinning[3]/hmc3->Integral(binCutRange[0],binCutRange[1]);
                                                else    
                                                        corrforrebinning[3]=1.0;
                                                        
-                                               debug<<" cor bin "<<corrforrebinning[0]<<" "<<corrforrebinning[1]<<" "<<corrforrebinning[2]<<" "<<corrforrebinning[3]<<endl;
 
 
                                        }
@@ -565,101 +769,151 @@ void DCACorrectionMarek(AliSpectraBothHistoManager* hman_data, AliSpectraBothHis
                                        gPad->SetLogy();
         
                                        TObjArray *mc=0x0;
-                                       if(useMaterial)
+                                       Int_t Npar=3;   
+                                       if(fitsettings==3)
                                                mc = new TObjArray(3);        // MC histograms are put in this array
                                        else
+                                       {
                                                mc = new TObjArray(2);
+                                               Npar=2;
+                                       }
+                                               
                                        mc->Add(hmc1);
-                                       mc->Add(hmc2);
-                                       if(useMaterial)
+                                       if(fitsettings&0x1)
+                                               mc->Add(hmc2);
+                                       if(fitsettings&0x2)
                                                mc->Add(hmc3);
                                        TFractionFitter* fit = new TFractionFitter(hToFit,mc); // initialise
                                        fit->Constrain(0,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
-                                       fit->Constrain(1,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
-                                       if(useMaterial)
-                                               fit->Constrain(2,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
-                                       fit->SetRangeX(hToFit->GetXaxis()->FindBin(FitRange[0]),hToFit->GetXaxis()->FindBin(FitRange[1]));
-                                       hToFit->GetXaxis()->SetRange(hToFit->GetXaxis()->FindBin(FitRange[0]),hToFit->GetXaxis()->FindBin(FitRange[1]));
-                                       hToFit->SetTitle(Form("DCA distr - %s %s %s %lf",sample[isample].Data(),Particle[ipart].Data(),Sign[icharge].Data(),lowedge));
+                                       Double_t defaultStep = 0.01;
+                                       for (int iparm = 0; iparm < Npar; ++iparm) 
+                                       {
+                                               TString name("frac"); name += iparm;
+                                               if(iparm==0)
+                                                       fit->GetFitter()->SetParameter(iparm,name.Data(),0.85,0.01,0.0,1.0);
+                                               else if (Npar==2)
+                                                       fit->GetFitter()->SetParameter(iparm,name.Data(),0.15,0.01,0.0,1.0);
+                                               else
+                                                       fit->GetFitter()->SetParameter(iparm,name.Data(),0.075,0.01,0.0,1.0);
+       
+                                       }
+                                       if(fitsettings&0x1)
+                                               fit->Constrain(1,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
+                                       if(fitsettings&0x2)
+                                               fit->Constrain(1+(fitsettings&0x1),0.0,1.0);               // constrain fraction 1 to be between 0 and 1
+                                       
+                                       Int_t binFitRange[]={hmc1->GetXaxis()->FindBin(FitRange[0]),hmc1->GetXaxis()->FindBin(FitRange[1])};                    
+                                       fit->SetRangeX(binFitRange[0],binFitRange[1]);
+                                       hToFit->GetXaxis()->SetRange(binFitRange[0],binFitRange[1]);
+                               //      hToFit->SetTitle(Form("DCA distr - %s %s %s %lf",sample[isample].Data(),Particle[ipart].Data(),Sign[icharge].Data(),lowedge));
+                                       hToFit->SetTitle("");
                                        Int_t status = fit->Fit();               // perform the fit
                                        cout << "fit status: " << status << endl;
                                        debug<<"fit status: " << status << endl;
                
-                                       if (status == 0 && usefit[ipart])
-                                       {                          // check on fit status
+                                       if (status == 0)
+                                       {       
+                                               Double_t v1=0.0,v2=0.0,v3=0.0;
+                                               Double_t ev1=0.0,ev2=0.0,ev3=0.0;
+                                               Double_t cov=0.0;
+                                                                  
+                                               // check on fit status
                                                TH1F* result = (TH1F*) fit->GetPlot();
                                                TH1F* PrimMCPred=(TH1F*)fit->GetMCPrediction(0);
-                                               TH1F* secStMCPred=(TH1F*)fit->GetMCPrediction(1);
+                                               
+                                               TH1F* secStMCPred=0X0;
                                        
                                                TH1F* secMCPred=0x0;
-                                               if(useMaterial)
-                                                       secMCPred=(TH1F*)fit->GetMCPrediction(2);
            
-                                               Double_t v1=0,v2=0,v3=0;
-                                               Double_t ev1=0,ev2=0,ev3=0;
-                                               Double_t cov=0.0;
                                                //first method, use directly the fit result
                                                fit->GetResult(0,v1,ev1);
-                                               fit->GetResult(1,v2,ev2);
-                                               if(useMaterial)
+
+                                               if(fitsettings&0x1)
                                                {
-                                                       fit->GetResult(2,v3,ev3);
-                                                       fit->GetFitter()->GetCovarianceMatrixElement(1,2);
+                                                       fit->GetResult(1,v2,ev2);
+                                                       secStMCPred=(TH1F*)fit->GetMCPrediction(1);
+
+                                               }
+                                               if(fitsettings&0x2)
+                                               {
+                                                       fit->GetResult(1+(fitsettings&0x1),v3,ev3);
+                                                       secMCPred=(TH1F*)fit->GetMCPrediction(1+(fitsettings&0x1));
+                                                       if(fitsettings&0x1)
+                                                               cov=fit->GetFitter()->GetCovarianceMatrixElement(1,2);
+
+
                                                }
-                                               debug<<v1<<" "<<ev1<<" "<<v2<<" "<<ev2<<" "<<v3<<" "<<ev3<<" "<<endl;
+                                               
                                        
-           
-          
-                                               Float_t normalizationdata=hToFit->Integral(hToFit->GetXaxis()->FindBin(CutRange[0]),hToFit->GetXaxis()->FindBin(CutRange[1]))/hToFit->Integral(hToFit->GetXaxis()->FindBin(FitRange[0]),hToFit->GetXaxis()->FindBin(FitRange[1]));
+                                               debug<<v1<<" "<<ev1<<" "<<v2<<" "<<ev2<<" "<<v3<<" "<<ev3<<" "<<" "<<cov<<endl;
+                                       
+                                               // becuase dca cut range is not a fit range the results from TFractionFitter should be rescale
+                                               // This can be done in two ways or use input histograms or output histograms
+                                               // The difference between those two methods should be on the level of statistical error         
+                                               // I use output histograms 
                                                
-                                               Float_t normalizationmc1=hmc1->Integral(hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1]))/hmc1->Integral(hmc1->GetXaxis()->FindBin(FitRange[0]),hmc1->GetXaxis()->FindBin(FitRange[1]))/normalizationdata;
-                                               Float_t normalizationmc2=hmc2->Integral(hmc2->GetXaxis()->FindBin(CutRange[0]),hmc2->GetXaxis()->FindBin(CutRange[1]))/hmc2->Integral(hmc1->GetXaxis()->FindBin(FitRange[0]),hmc2->GetXaxis()->FindBin(FitRange[1]))/normalizationdata;
-                                               Float_t normalizationmc3=hmc3->Integral(hmc3->GetXaxis()->FindBin(CutRange[0]),hmc3->GetXaxis()->FindBin(CutRange[1]))/hmc3->Integral(hmc1->GetXaxis()->FindBin(FitRange[0]),hmc3->GetXaxis()->FindBin(FitRange[1]))/normalizationdata;
+                                               // Method 1 input histo         
+
+                                               Float_t normalizationdata=hToFit->Integral(hToFit->GetXaxis()->FindBin(CutRange[0]),hToFit->GetXaxis()->FindBin(CutRange[1]))/hToFit->Integral(binFitRange[0],binFitRange[1]);
+                                               normalizationdata*=corrforrebinning[0];
+
+                                               Float_t normalizationmc1=(hmc1->Integral(binCutRange[0],binCutRange[1])/hmc1->Integral(binFitRange[0],binFitRange[1]))/normalizationdata;
+                                               Float_t normalizationmc2=0.0;
+                                               if(fitsettings&0x1)
+                                                       normalizationmc2=(hmc2->Integral(binCutRange[0],binCutRange[1])/hmc2->Integral(binFitRange[0],binFitRange[1]))/normalizationdata;
+                                               Float_t normalizationmc3=0.0;
+                                               if(fitsettings&0x2)
+                                                       normalizationmc3=(hmc3->Integral(binCutRange[0],binCutRange[1])/hmc3->Integral(binFitRange[0],binFitRange[1]))/normalizationdata;
+
+                                               normalizationmc1*=corrforrebinning[1];
+                                               normalizationmc2*=corrforrebinning[2];
+                                               normalizationmc3*=corrforrebinning[3];
+
                                                debug<<"After Nor"<<endl;
                                                debug<<v1*normalizationmc1<<" "<<ev1*normalizationmc1<<" "<<v2*normalizationmc2<<" "<<ev2*normalizationmc2<<" "<<v3*normalizationmc3<<" "<<ev3*normalizationmc3<<" "<<endl;
                                                debug<<1.0-v1*normalizationmc1<<" "<<ev1*normalizationmc1<<" "<<v2*normalizationmc2+v3*normalizationmc3<<" "<<TMath::Sqrt(ev2*ev2*normalizationmc2*normalizationmc2+ev3*ev3*normalizationmc3*normalizationmc3+cov*normalizationmc3*normalizationmc2)<<endl;
-                                               Float_t normalizationdata1=result->Integral(result->GetXaxis()->FindBin(CutRange[0]),result->GetXaxis()->FindBin(CutRange[1]))/result->Integral(result->GetXaxis()->FindBin(FitRange[0]),result->GetXaxis()->FindBin(FitRange[1]));
+                                               debug<<"addtional info"<<endl;
+
+
+                                              //Method 2 output histo          
+
+                                               Float_t normalizationdata1=result->Integral(binCutRange[0],binCutRange[1])/result->Integral(binFitRange[0],binFitRange[1]);
                                                
 
+                                               // if the cut range is bigger the fit range we should calculate the normalization factor for data using the data histogram 
+                                               // because result histogram has entries only in fits range       
+                                               if(FitRange[0]>CutRange[0]||FitRange[1]<CutRange[1])    
+                                                       normalizationdata1=normalizationdata;
+                                       
                                                normalizationdata1*=corrforrebinning[0];
 
 
-                                               Float_t normalizationmc11=PrimMCPred->Integral(PrimMCPred->GetXaxis()->FindBin(CutRange[0]),PrimMCPred->GetXaxis()->FindBin(CutRange[1]))/PrimMCPred->Integral(PrimMCPred->GetXaxis()->FindBin(FitRange[0]),PrimMCPred->GetXaxis()->FindBin(FitRange[1]))/normalizationdata1;
-                                               Float_t normalizationmc21=secStMCPred->Integral(secStMCPred->GetXaxis()->FindBin(CutRange[0]),secStMCPred->GetXaxis()->FindBin(CutRange[1]))/secStMCPred->Integral(secStMCPred->GetXaxis()->FindBin(FitRange[0]),secStMCPred->GetXaxis()->FindBin(FitRange[1]))/normalizationdata1;
-                                               Float_t normalizationmc31=0;
-                                               if(useMaterial)
-                                                       normalizationmc31=secMCPred->Integral(secMCPred->GetXaxis()->FindBin(CutRange[0]),secMCPred->GetXaxis()->FindBin(CutRange[1]))/secMCPred->Integral(secMCPred->GetXaxis()->FindBin(FitRange[0]),secMCPred->GetXaxis()->FindBin(FitRange[1]))/normalizationdata1;
+                                               Float_t normalizationmc11=(PrimMCPred->Integral(binCutRange[0],binCutRange[1])/PrimMCPred->Integral(binFitRange[0],binFitRange[1]))/normalizationdata1;
+                                               Float_t normalizationmc21=0.0;
+                                               if(fitsettings&0x1)
+                                                       normalizationmc21=(secStMCPred->Integral(binCutRange[0],binCutRange[1])/secStMCPred->Integral(binFitRange[0],binFitRange[1]))/normalizationdata1;
+                                               Float_t normalizationmc31=0.0;
+                                               if(fitsettings&0x2)
+                                                       normalizationmc31=(secMCPred->Integral(binCutRange[0],binCutRange[1])/secMCPred->Integral(binFitRange[0],binFitRange[1]))/normalizationdata1;
                                                
                                                normalizationmc11*=corrforrebinning[1];
                                                normalizationmc21*=corrforrebinning[2];
                                                normalizationmc31*=corrforrebinning[3];
 
-                               debug<<"After Nor 2"<<endl;
+                                               debug<<"After Nor 2"<<endl;
                                                debug<<v1*normalizationmc11<<" "<<ev1*normalizationmc11<<" "<<v2*normalizationmc21<<" "<<ev2*normalizationmc21<<" "<<v3*normalizationmc31<<" "<<ev3*normalizationmc31<<endl;
                                                
                                                debug<<1.0-v1*normalizationmc11<<" "<<ev1*normalizationmc11<<" "<<v2*normalizationmc21+v3*normalizationmc31<<" "<<TMath::Sqrt(ev2*ev2*normalizationmc21*normalizationmc21+ev3*ev3*normalizationmc31*normalizationmc31+cov*normalizationmc31*normalizationmc21)<<endl;
                                        
-                                               debug<<CutRange[0]<<" "<<CutRange[1]<<endl;             
-                                               debug<<" Entries "<<isample<<" "<<hToFit->GetEntries()<<" "<<hmc1->GetEntries()<<" "<<hmc2->GetEntries()<<" "<<hmc3->GetEntries()<<endl;
-                                               debug<<"2 Entries "<<isample<<" "<<hToFit->Integral(result->GetXaxis()->FindBin(CutRange[0]),result->GetXaxis()->FindBin(CutRange[1]))<<" "<<hmc1->Integral(result->GetXaxis()->FindBin(CutRange[0]),result->GetXaxis()->FindBin(CutRange[1]))<<" "<<hmc2->Integral(result->GetXaxis()->FindBin(CutRange[0]),result->GetXaxis()->FindBin(CutRange[1]))<<" "<<hmc3->Integral(result->GetXaxis()->FindBin(CutRange[0]),result->GetXaxis()->FindBin(CutRange[1]))<<endl;
-
                                                
                                                hconWD[index+6*isample]->SetBinContent(ibin_data,v2*normalizationmc21);
                                                hconWD[index+6*isample]->SetBinError(ibin_data,ev2*normalizationmc21);
                                                hconMat[index+6*isample]->SetBinContent(ibin_data,v3*normalizationmc31);
                                                hconMat[index+6*isample]->SetBinError(ibin_data,ev3*normalizationmc31);
                                                hprimary[index+6*isample]->SetBinContent(ibin_data,v1*normalizationmc11);
-                                               hprimary[index+6*isample]->SetBinContent(ibin_data,v1*normalizationmc11);
-                                               if(useMaterial)
-                                               {
-                                                       hcon[index+6*isample]->SetBinContent(ibin_data,v2*normalizationmc21+v3*normalizationmc31);
-                                                       hcon[index+6*isample]->SetBinError(ibin_data,TMath::Sqrt(ev2*ev2*normalizationmc21*normalizationmc21+ev3*ev3*normalizationmc31*normalizationmc31+cov*normalizationmc31*normalizationmc21));
-                                               }
-                                               else
-                                               {
-                                                       hcon[index+6*isample]->SetBinContent(ibin_data,v2*normalizationmc21);
-                                                       hcon[index+6*isample]->SetBinError(ibin_data,ev2*normalizationmc21);
-                                               }
+                                               hprimary[index+6*isample]->SetBinError(ibin_data,ev1*normalizationmc11);
+                                               hcon[index+6*isample]->SetBinContent(ibin_data,v2*normalizationmc21+v3*normalizationmc31);
+                                               hcon[index+6*isample]->SetBinError(ibin_data,TMath::Sqrt(ev2*ev2*normalizationmc21*normalizationmc21+ev3*ev3*normalizationmc31*normalizationmc31+cov*normalizationmc31*normalizationmc21));
                                                
                                                
                                                
@@ -667,26 +921,48 @@ void DCACorrectionMarek(AliSpectraBothHistoManager* hman_data, AliSpectraBothHis
                                                result->Scale(1.0/result->Integral(result->GetXaxis()->FindBin(FitRange[0]),result->GetXaxis()->FindBin(FitRange[1])));
                                                hToFit->Scale(1.0/hToFit->Integral(hToFit->GetXaxis()->FindBin(FitRange[0]),hToFit->GetXaxis()->FindBin(FitRange[1])));
                                                PrimMCPred->Scale(v1/PrimMCPred->Integral(PrimMCPred->GetXaxis()->FindBin(FitRange[0]),PrimMCPred->GetXaxis()->FindBin(FitRange[1])));
-                                               secStMCPred->Scale(v2/secStMCPred->Integral(secStMCPred->GetXaxis()->FindBin(FitRange[0]),secStMCPred->GetXaxis()->FindBin(FitRange[1])));
-                                               if(useMaterial)
-                                                       secMCPred->Scale(v3/secMCPred->Integral(secMCPred->GetXaxis()->FindBin(FitRange[0]),secMCPred->GetXaxis()->FindBin(FitRange[1])));          
-                                                  
-                                               result->SetLineColor(kBlack);
-                                               PrimMCPred->SetLineColor(kGreen+2);
-                                               secStMCPred->SetLineColor(kRed);
-                                               
+
                                                hToFit->SetMinimum(0.0001);
+                                               hToFit->GetXaxis()->SetTitle("DCA_{xy} (cm)");
+                                               hToFit->GetYaxis()->SetTitle("N_{counts}/N_{counts}(-3cm;3cm)");
+                                               hToFit->GetYaxis()->SetTitleOffset(1.3);
                                                hToFit->DrawClone("E1x0");
+                                               Leg1->AddEntry(hToFit,"data","p");
                                                result->SetTitle("Fit result");
-                                               result->DrawClone("lhistsame");
-                                               PrimMCPred->DrawClone("lhistsame");
-                                               secStMCPred->DrawClone("lhistsame");
-                                               if(useMaterial)
+                                               result->SetLineColor(kBlack);
+                                               Leg1->AddEntry(result,"fit result","l");
+                                               result->DrawClone("histsame");
+                                       
+                                               PrimMCPred->SetLineColor(kGreen+2);
+                                               PrimMCPred->SetLineStyle(2);
+                                                PrimMCPred->SetLineWidth(3.0);
+                                               Leg1->AddEntry(PrimMCPred,"primaries","l");
+                                               PrimMCPred->DrawClone("histsame");
+                                               if(fitsettings&0x1)
                                                {
-                                                       secMCPred->SetLineColor(kBlue); 
-                                                       secMCPred->DrawClone("lhistsame");
+
+                                                       secStMCPred->Scale(v2/secStMCPred->Integral(secStMCPred->GetXaxis()->FindBin(FitRange[0]),secStMCPred->GetXaxis()->FindBin(FitRange[1])));
+                                                       secStMCPred->SetLineColor(kRed);
+                                                       secStMCPred->SetLineWidth(3.0);
+
+                                                       secStMCPred->SetLineStyle(3);
+                                                       Leg1->AddEntry(secStMCPred,"weak decays","l");
+                                                       secStMCPred->DrawClone("histsame");
+
                                                }
-                                       }
+                                               if(fitsettings&0x2)
+                                               {
+                                                       
+                                                       secMCPred->Scale(v3/secMCPred->Integral(secMCPred->GetXaxis()->FindBin(FitRange[0]),secMCPred->GetXaxis()->FindBin(FitRange[1])));
+                                                       secMCPred->SetLineColor(kBlue);
+                                                       secMCPred->SetLineWidth(3.0);
+
+                                                       secMCPred->SetLineStyle(4);     
+                                                       Leg1->AddEntry(secMCPred,"material","l");
+                                                       secMCPred->DrawClone("histsame");
+           
+                                               }   
+                                       }                               
                                        else
                                        {
                                                hconWD[index+6*isample]->SetBinContent(ibin_data,0.0);
@@ -698,6 +974,11 @@ void DCACorrectionMarek(AliSpectraBothHistoManager* hman_data, AliSpectraBothHis
                                                hprimary[index+6*isample]->SetBinContent(ibin_data,1.0);
                                                hprimary[index+6*isample]->SetBinError(ibin_data,0.0);
                                        }
+                                       Leg1->DrawClone();
+                                       TLatex* texttitle=new TLatex();
+                                       texttitle->SetNDC();
+                                       texttitle->SetTextSize(0.04);
+                                       texttitle->DrawLatex(0.12,0.92,Form("%s %.2f<#it{p}_{T} < %.2f (GeV/#it{c})",symboles[index].Data(),lowedge,binwidth+lowedge));
                                        listofdcafits->Add(cDCA);
                                        
                                        //cDCA->Write();
@@ -716,7 +997,10 @@ void DCACorrectionMarek(AliSpectraBothHistoManager* hman_data, AliSpectraBothHis
 void RecomputeErrors(TH1* h)
 {
        for (int i=0; i<=h->GetXaxis()->GetNbins(); i++)
+       {
+               cout<<h->GetBinContent(i)<<" "<<h->GetBinError(i)<<" error "<<TMath::Sqrt(h->GetBinContent(i))<<endl;
                h->SetBinError(i,TMath::Sqrt(h->GetBinContent(i)));
+       }
        h->Sumw2();     
 }
 void SetBintoOne(TH1* h)
@@ -769,9 +1053,14 @@ void GFCorrection(TH1F **Spectra,Float_t tofpt,UInt_t options)
                gGFCorrectionKaonMinus->SetName("gGFCorrectionKaonMinus");
                gGFCorrectionKaonMinus->SetTitle("gGFCorrectionKaonMinus");
                 TString fnameGeanFlukaK="GFCorrection/correctionForCrossSection.321.root";
-               TFile *fGeanFlukaK= new TFile(fnameGeanFlukaK.Data());
-               if (!fGeanFlukaK)               
-                       return;
+               TFile *fGeanFlukaK= TFile::Open(fnameGeanFlukaK.Data());
+               if (!fGeanFlukaK)
+               {
+                       fnameGeanFlukaK="$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/TestAOD/correctionForCrossSection.321.root";
+                       fGeanFlukaK= TFile::Open(fnameGeanFlukaK.Data());
+                       if (!fGeanFlukaK)
+                               return;
+               }
                TH1F *hGeantFlukaKPos=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionParticles");
                TH1F *hGeantFlukaKNeg=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionAntiParticles");
                //getting GF func for Kaons with TOF
@@ -783,6 +1072,9 @@ void GFCorrection(TH1F **Spectra,Float_t tofpt,UInt_t options)
                 fGFKPosMatching = TOFmatchMC_geantflukaCorrection(3,kPositive);
                 TF1 *fGFKNegMatching;
                 fGFKNegMatching = TOFmatchMC_geantflukaCorrection(3,kNegative);
+                if(tcutsdata->GetUseTypeDependedTOFCut())
+                       tofpt=tcutsdata->GetPtTOFMatchingKaon();
+
                 for(Int_t binK=0;binK<=Spectra[1]->GetNbinsX();binK++)
                {
                        if(Spectra[1]->GetBinCenter(binK)<tofpt)
@@ -825,7 +1117,18 @@ void GFCorrection(TH1F **Spectra,Float_t tofpt,UInt_t options)
          const Int_t kNCharge=2;
          Int_t kPos=0;
          Int_t kNeg=1;
-         TFile* fGFProtons = new TFile ("GFCorrection/correctionForCrossSection.root");
+          TString fnameGFProtons= "GFCorrection/correctionForCrossSection.root";
+         TFile* fGFProtons = TFile::Open(fnameGFProtons.Data());
+         if (!fGFProtons)
+         { 
+               fnameGFProtons="$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/TestAOD/correctionForCrossSection.root";
+               fGFProtons = TFile::Open(fnameGFProtons.Data());
+               if (!fGFProtons)
+                       return;
+         }
+               
+
+
          TH2D * hCorrFluka[kNCharge];
          TH2D * hCorrFluka[2];
          hCorrFluka[kPos] = (TH2D*)fGFProtons->Get("gHistCorrectionForCrossSectionProtons");
@@ -841,8 +1144,10 @@ void GFCorrection(TH1F **Spectra,Float_t tofpt,UInt_t options)
          fGFpNegMatching = TOFmatchMC_geantflukaCorrection(4,kNegative);
          
         
-               Int_t nbins = Spectra[2]->GetNbinsX();
-               
+          Int_t nbins = Spectra[2]->GetNbinsX();
+          if(tcutsdata->GetUseTypeDependedTOFCut())
+               tofpt=tcutsdata->GetPtTOFMatchingProton();
+
                for(Int_t ibin = 0; ibin < nbins; ibin++)
                {
                        if(Spectra[2]->GetBinCenter(ibin)<tofpt)
@@ -973,17 +1278,25 @@ void MatchingTOFEff(TH1F** Spectra, TList* list=0x0)
 {
          if(TOFMatchingScalling[0]<0.0&&TOFMatchingScalling[1]<0.0)
          {
-                 TH1F *hMatcEffPos_data=(TH1F*)tcutsdata->GetHistoNMatchedPos();
-                 hMatcEffPos_data->Divide((TH1F*)tcutsdata->GetHistoNSelectedPos());
+               TH1F *hMatcEffPos_data=(TH1F*)tcutsdata->GetHistoNMatchedPos();
+                 hMatcEffPos_data->Sumw2();
+                 //hMatcEffPos_data->Divide((TH1F*)tcutsdata->GetHistoNSelectedPos());
+                 hMatcEffPos_data->Divide(hMatcEffPos_data,(TH1F*)tcutsdata->GetHistoNSelectedPos(),1,1,"B");
                  hMatcEffPos_data->SetTitle("Matching Eff Pos - data");
                  TH1F *hMatcEffNeg_data=(TH1F*)tcutsdata->GetHistoNMatchedNeg();
-                 hMatcEffNeg_data->Divide((TH1F*)tcutsdata->GetHistoNSelectedNeg());
+                 hMatcEffNeg_data->Sumw2();
+                 //hMatcEffNeg_data->Divide((TH1F*)tcutsdata->GetHistoNSelectedNeg());
+                 hMatcEffNeg_data->Divide(hMatcEffNeg_data,(TH1F*)tcutsdata->GetHistoNSelectedNeg(),1,1,"B");
                  hMatcEffNeg_data->SetTitle("Matching Eff Neg - data");
                  TH1F *hMatcEffPos_mc=(TH1F*)tcutsmc->GetHistoNMatchedPos();
-                 hMatcEffPos_mc->Divide((TH1F*)tcutsmc->GetHistoNSelectedPos());
+                 hMatcEffPos_mc->Sumw2();
+                 //hMatcEffPos_mc->Divide((TH1F*)tcutsmc->GetHistoNSelectedPos());
+                 hMatcEffPos_mc->Divide(hMatcEffPos_mc,(TH1F*)tcutsmc->GetHistoNSelectedPos(),1,1,"B");
                  hMatcEffPos_mc->SetTitle("Matching Eff Pos - mc");
                  TH1F *hMatcEffNeg_mc=(TH1F*)tcutsmc->GetHistoNMatchedNeg();
-                 hMatcEffNeg_mc->Divide((TH1F*)tcutsmc->GetHistoNSelectedNeg());
+                 hMatcEffNeg_mc->Sumw2();
+                 //hMatcEffNeg_mc->Divide((TH1F*)tcutsmc->GetHistoNSelectedNeg());
+                 hMatcEffNeg_mc->Divide(hMatcEffNeg_mc,(TH1F*)tcutsmc->GetHistoNSelectedNeg(),1,1,"B");
                  hMatcEffNeg_mc->SetTitle("Matching Eff Neg - mc");
 
 
@@ -1004,6 +1317,7 @@ void MatchingTOFEff(TH1F** Spectra, TList* list=0x0)
                        
                  TOFMatchingScalling[0]=pol0MatchPos_data->GetParameter(0);
                  TOFMatchingScalling[1]=pol0MatchNeg_data->GetParameter(0);
+
          }
          //Correction spectra for matching efficiency
          //For the moment I'm using the inclusive correction
@@ -1028,7 +1342,7 @@ Double_t eta2y(Double_t pt, Double_t mass, Double_t eta)
   return TMath::ASinH(pt / mt * TMath::SinH(eta));
 }
 
-TH1* GetSumAllCh(TH1F** spectra, Double_t* mass  )
+TH1* GetSumAllCh(TH1F** spectra, Double_t* mass,Double_t etacut)
 {
        TH1F* allch=(((TH1F*))spectra[0]->Clone("allCh"));
        allch->Reset();
@@ -1044,10 +1358,10 @@ TH1* GetSumAllCh(TH1F** spectra, Double_t* mass  )
                                Float_t pt=spectra[i]->GetXaxis()->GetBinCenter(j);
                                Float_t mt2=masstmp*masstmp+pt*pt;              
                                Float_t mt=TMath::Sqrt(mt2);
-                               Float_t maxy=eta2y(pt,masstmp,0.8);
+                               Float_t maxy=eta2y(pt,masstmp,etacut/2.0);
                                Float_t conver=maxy*(TMath::Sqrt(1-masstmp*masstmp/(mt2*TMath::CosH(maxy)*TMath::CosH(maxy)))+TMath::Sqrt(1-masstmp*masstmp/(mt2*TMath::CosH(0.0)*TMath::CosH(0.0))));
-                               conver=conver/1.6;
-                               cout<<maxy<<" "<<conver<<" "<<masstmp<<""<<spectra[i]->GetName()<<endl;
+                               conver=conver/etacut;
+                               //cout<<maxy<<" "<<conver<<" "<<masstmp<<""<<spectra[i]->GetName()<<endl;
                                Float_t bincontent=allch->GetBinContent(j);
                                Float_t binerror=allch->GetBinError(j);
                                bincontent+=conver*value;
@@ -1076,3 +1390,315 @@ void Divideby2pipt(TH1* hist)
 
        }
 }
+
+Short_t DCAfitsettings (Float_t pt, Int_t type)
+{
+       Short_t value=0x0;
+       if (pt<maxptformaterial[type]&&pt>minptformaterial[type])
+               value=value+2;
+       if (pt<maxptforWD[type]&&pt>minptforWD[type])
+               value=value+1;
+       return value;   
+
+} 
+
+Float_t Normaliztionwithbin0integrals(UInt_t options)
+{
+       
+       TH1F* bin0mcRec=(TH1F*)ecutsmc->GetHistoVtxGenerated()->Clone("Bin0_rec");
+       TH1F* bin0mcMC=(TH1F*)ecutsmc->GetHistoVtxGenerated()->Clone("Bin0_MC");
+
+       TH1F* vertexmc=ecutsmc->GetHistoVtxAftSelwithoutZvertexCut(); 
+       TH1F* vertexmcMCz=ecutsmc->GetHistoVtxAftSelwithoutZvertexCutusingMCz(); 
+       TH1F* vertexdata=ecutsdata->GetHistoVtxAftSelwithoutZvertexCut();
+
+       TH1I* histodata=ecutsdata->GetHistoCuts();
+       TH1I* histomc=ecutsmc->GetHistoCuts();
+
+       Float_t dataevents=(Float_t)histodata->GetBinContent(3);
+       //cout<<histodata->GetBinContent(2)<<endl;
+       Float_t databin0events=((Float_t)histodata->GetBinContent(2))-((Float_t)histodata->GetBinContent(4));   
+
+       bin0mcRec->Sumw2();
+       bin0mcMC->Sumw2();
+               
+       bin0mcRec->Add(vertexmc,-1);
+       bin0mcMC->Add(vertexmcMCz,-1);
+       
+       bin0mcRec->Divide(vertexmc);
+       bin0mcMC->Divide(vertexmcMCz);
+       
+       bin0mcRec->Multiply(vertexdata);
+       bin0mcMC->Multiply(vertexdata);
+       
+       Float_t bin0mcRecN=0.0;
+       Float_t bin0mcMCN=0.0;
+
+       for (int i=0;i<=bin0mcRec->GetXaxis()->GetNbins();i++)
+       {
+               bin0mcRecN+=bin0mcRec->GetBinContent(i);
+               bin0mcMCN+=bin0mcMC->GetBinContent(i);
+
+       }
+       bin0mcRec->Scale(databin0events/bin0mcRecN);
+       bin0mcMC->Scale(databin0events/bin0mcMCN);              
+       
+       Int_t binmin=bin0mcRec->GetXaxis()->FindBin(-10);
+       Int_t binmax=bin0mcRec->GetXaxis()->FindBin(10)-1;
+       cout<<  bin0mcRec->GetXaxis()->GetBinLowEdge(binmin)<<" "<<bin0mcRec->GetXaxis()->GetBinUpEdge(binmax)<<endl;
+       cout<<bin0mcRecN<<" "<<bin0mcMCN<<" "<<databin0events<<endl;    
+       cout<<dataevents<<" normalization "<<dataevents+bin0mcRec->Integral(binmin,binmax)<<" "<<dataevents+bin0mcMC->Integral(binmin,binmax)<<endl;
+       cout<<histodata->GetBinContent(2)<<" "<<histodata->GetBinContent(4)<<endl;
+       if ((options&knormalizationwithbin0integralsdata)==knormalizationwithbin0integralsdata)
+               return  dataevents+bin0mcRec->Integral(binmin,binmax);
+       else if ((options&kknormalizationwithbin0integralsMC)==knormalizationwithbin0integralsMC)   
+               return dataevents+bin0mcMC->Integral(binmin,binmax) ;
+       else
+               return 1;               
+}
+
+Bool_t ReadConfigFile(TString configfile)
+{
+       ifstream infile(configfile.Data());
+       if(infile.is_open()==false)
+               return false;
+       TString namesofSetting[]={"CutRangeMin","CutRangeMax","FitRangeMin","FitRangeMax","MinMatPionPlus","MaxMatPionPlus","MinMatKaonPlus","MaxMatKaonPlus","MinMatProtonPlus","MaxMatProtonPlus","MinMatPionMinus","MaxMatPionMinus","MinMatKaonMinus","MaxMatKaonMinus","MinMatProtonMinus","MaxMatProtonMinus","MinWDPionPlus","MaxWDPionPlus","MinWDKaonPlus","MaxWDKaonPlus","MinWDProtonPlus","MaxWDProtonPlus","MinWDPionMinus","MaxWDPionMinus","MinWDKaonMinus","MaxWDKaonMinus","MinWDProtonMinus","MaxWDProtonMinus","MaxContaminationPIDMC","MinPions","MaxPions","MinKaons","MaxKaons","MinProtons","MaxProtons","TOFPIDsignalmatchPion","TOFPIDsignalmatchKaon","TOFPIDsignalmatchProton","NamefileEff","NameFilePIDcon","NameFileSeccon"};     
+
+       char buffer[256];
+       while (infile.eof()==false)
+       {
+               buffer[0]='#'; 
+               while (buffer[0]=='#'&&infile.eof()==false)
+                       infile.getline(buffer,256);
+               TString tmpstring(buffer);
+               cout<<buffer<<endl;
+               if(tmpstring.Contains(namesofSetting[0]))
+                       CutRange[0]=(tmpstring.Remove(0,namesofSetting[0].Length()+1)).Atof();
+               else if (tmpstring.Contains(namesofSetting[1]))
+                       CutRange[1]=(tmpstring.Remove(0,namesofSetting[1].Length()+1)).Atof();  
+               else if (tmpstring.Contains(namesofSetting[2]))
+                       FitRange[0]=(tmpstring.Remove(0,namesofSetting[2].Length()+1)).Atof();
+               else if (tmpstring.Contains(namesofSetting[3]))
+                       FitRange[1]=(tmpstring.Remove(0,namesofSetting[3].Length()+1)).Atof();
+               else if (tmpstring.Contains(namesofSetting[4]))
+                       minptformaterial[0]=(tmpstring.Remove(0,namesofSetting[4].Length()+1)).Atof();          
+               else if (tmpstring.Contains(namesofSetting[5]))
+                       maxptformaterial[0]=(tmpstring.Remove(0,namesofSetting[5].Length()+1)).Atof();  
+               else if (tmpstring.Contains(namesofSetting[6]))
+                       minptformaterial[1]=(tmpstring.Remove(0,namesofSetting[6].Length()+1)).Atof();          
+               else if (tmpstring.Contains(namesofSetting[7]))
+                       maxptformaterial[1]=(tmpstring.Remove(0,namesofSetting[7].Length()+1)).Atof();  
+               else if (tmpstring.Contains(namesofSetting[8]))
+                       minptformaterial[2]=(tmpstring.Remove(0,namesofSetting[8].Length()+1)).Atof();          
+               else if (tmpstring.Contains(namesofSetting[9]))
+                       maxptformaterial[2]=(tmpstring.Remove(0,namesofSetting[9].Length()+1)).Atof();
+               else if (tmpstring.Contains(namesofSetting[10]))
+                       minptformaterial[3]=(tmpstring.Remove(0,namesofSetting[10].Length()+1)).Atof();         
+               else if (tmpstring.Contains(namesofSetting[11]))
+                       maxptformaterial[3]=(tmpstring.Remove(0,namesofSetting[11].Length()+1)).Atof();
+               else if (tmpstring.Contains(namesofSetting[12]))
+                       minptformaterial[4]=(tmpstring.Remove(0,namesofSetting[12].Length()+1)).Atof();         
+               else if (tmpstring.Contains(namesofSetting[13]))
+                       maxptformaterial[4]=(tmpstring.Remove(0,namesofSetting[13].Length()+1)).Atof();
+               else if (tmpstring.Contains(namesofSetting[14]))
+                       minptformaterial[5]=(tmpstring.Remove(0,namesofSetting[14].Length()+1)).Atof();         
+               else if (tmpstring.Contains(namesofSetting[15]))
+                       maxptformaterial[5]=(tmpstring.Remove(0,namesofSetting[15].Length()+1)).Atof();
+               else if (tmpstring.Contains(namesofSetting[16]))
+                       minptforWD[0]=(tmpstring.Remove(0,namesofSetting[16].Length()+1)).Atof();               
+               else if (tmpstring.Contains(namesofSetting[17]))
+                       maxptforWD[0]=(tmpstring.Remove(0,namesofSetting[17].Length()+1)).Atof();       
+               else if (tmpstring.Contains(namesofSetting[18]))
+                       minptforWD[1]=(tmpstring.Remove(0,namesofSetting[18].Length()+1)).Atof();               
+               else if (tmpstring.Contains(namesofSetting[19]))
+                       maxptforWD[1]=(tmpstring.Remove(0,namesofSetting[19].Length()+1)).Atof();       
+               else if (tmpstring.Contains(namesofSetting[20]))
+                       minptforWD[2]=(tmpstring.Remove(0,namesofSetting[20].Length()+1)).Atof();               
+               else if (tmpstring.Contains(namesofSetting[21]))
+                       maxptforWD[2]=(tmpstring.Remove(0,namesofSetting[21].Length()+1)).Atof();
+               else if (tmpstring.Contains(namesofSetting[22]))
+                       minptforWD[3]=(tmpstring.Remove(0,namesofSetting[22].Length()+1)).Atof();               
+               else if (tmpstring.Contains(namesofSetting[23]))
+                       maxptforWD[3]=(tmpstring.Remove(0,namesofSetting[23].Length()+1)).Atof();
+               else if (tmpstring.Contains(namesofSetting[24]))
+                       minptforWD[4]=(tmpstring.Remove(0,namesofSetting[24].Length()+1)).Atof();               
+               else if (tmpstring.Contains(namesofSetting[25]))
+                       maxptforWD[4]=(tmpstring.Remove(0,namesofSetting[25].Length()+1)).Atof();
+               else if (tmpstring.Contains(namesofSetting[26]))
+                       minptforWD[5]=(tmpstring.Remove(0,namesofSetting[26].Length()+1)).Atof();               
+               else if (tmpstring.Contains(namesofSetting[27]))
+                       maxptforWD[5]=(tmpstring.Remove(0,namesofSetting[27].Length()+1)).Atof();                               
+               else if (tmpstring.Contains(namesofSetting[28]))
+                       fMaxContaminationPIDMC=(tmpstring.Remove(0,namesofSetting[28].Length()+1)).Atof();              
+               else if (tmpstring.Contains(namesofSetting[29]))
+                       minRanges[0]=(tmpstring.Remove(0,namesofSetting[29].Length()+1)).Atof();
+               else if (tmpstring.Contains(namesofSetting[30]))
+                       maxRanges[0]=(tmpstring.Remove(0,namesofSetting[30].Length()+1)).Atof();                
+               else if (tmpstring.Contains(namesofSetting[31]))
+                       minRanges[1]=(tmpstring.Remove(0,namesofSetting[31].Length()+1)).Atof();
+               else if (tmpstring.Contains(namesofSetting[32]))
+                       maxRanges[1]=(tmpstring.Remove(0,namesofSetting[32].Length()+1)).Atof();                
+               else if (tmpstring.Contains(namesofSetting[33]))
+                       minRanges[2]=(tmpstring.Remove(0,namesofSetting[33].Length()+1)).Atof();
+               else if (tmpstring.Contains(namesofSetting[34]))
+                       maxRanges[2]=(tmpstring.Remove(0,namesofSetting[34].Length()+1)).Atof();                
+               else if (tmpstring.Contains(namesofSetting[35]))
+                       TOFPIDsignalmatching[0]=(tmpstring.Remove(0,namesofSetting[35].Length()+1)).Atof();             
+               else if (tmpstring.Contains(namesofSetting[36]))
+                       TOFPIDsignalmatching[1]=(tmpstring.Remove(0,namesofSetting[36].Length()+1)).Atof();
+               else if (tmpstring.Contains(namesofSetting[37]))
+                       TOFPIDsignalmatching[2]=(tmpstring.Remove(0,namesofSetting[37].Length()+1)).Atof();     
+               else if (tmpstring.Contains(namesofSetting[38]))
+                       filenames[0]=tmpstring.Remove(0,namesofSetting[38].Length()+1);
+               else if (tmpstring.Contains(namesofSetting[39]))
+                       filenames[1]=tmpstring.Remove(0,namesofSetting[39].Length()+1);
+               else if (tmpstring.Contains(namesofSetting[40]))
+                       filenames[2]=tmpstring.Remove(0,namesofSetting[40].Length()+1);
+               else 
+                       continue;
+
+
+//     Double_t minRanges[3]={0.3,0.3,0.45};
+//Double_t maxRanges[3]={1.5,1.2,2.2};
+//Double_t fMaxContaminationPIDMC=0.2;
+
+
+
+       }
+       for(int i=0;i<6;i++)
+               cout<<minptformaterial[i]<<" "<<maxptformaterial[i]<<" "<<minptforWD[i]<<" "<<maxptforWD[i]<<endl;
+       cout<<FitRange[0]<<" "<<FitRange[1]<<" "<<CutRange[0]<<CutRange[1]<<endl;
+       if(FitRange[0]>=FitRange[1])    
+       {
+               cout<<"A"<<endl;                                
+               return false;
+       }
+       if(CutRange[0]>=CutRange[1])
+       {       
+               cout<<"B"<<endl;                                
+               return false;
+       }
+       for(int i=0;i<6;i++)
+       {
+               if((minptformaterial[i]>maxptformaterial[i]&&minptformaterial[i]>0.0)||minptformaterial[i]<0.0||maxptformaterial[i]<0.0)
+               {
+                       cout<<"C"<<endl;
+                       return false;
+               }
+               if((minptforWD[i]>maxptforWD[i]&&minptforWD[i]>0.0)||minptforWD[i]<0.0||maxptforWD[i]<0.0)
+               {
+                       cout<<"D"<<endl;
+                       return false;
+               }
+       }
+       for(int i=0;i<3;i++)
+               if(minRanges[i]>maxRanges[i])
+                       return false;
+               
+
+       return true;
+}
+
+void SubHistWithFullCorr(TH1F* h1, TH1F* h2, Float_t factor1=1.0, Float_t factor2=1.0)
+{
+       if(h1->GetNbinsX()!=h2->GetNbinsX())
+               return;
+       for (int i=0;i<=h1->GetNbinsX();i++)
+       {
+               Float_t tmpvalue=factor1*h1->GetBinContent(i)-factor2*h2->GetBinContent(i);
+               Float_t tmperror=TMath::Abs(factor1*factor1*h1->GetBinError(i)*h1->GetBinError(i)-factor2*factor2*h2->GetBinError(i)*h2->GetBinError(i));
+               h1->SetBinContent(i,tmpvalue);
+               h1->SetBinError(i,TMath::Sqrt(tmperror));
+       }               
+       
+}
+
+void TOFMatchingForNch(TH1* h)
+{
+        if(TOFMatchingScalling[0]>0.0&&TOFMatchingScalling[1]>0.0)
+        {
+               Float_t factor=0.5*TOFMatchingScalling[0]+0.5*TOFMatchingScalling[1];
+               for(Int_t ibin=1;ibin<h->GetNbinsX();ibin++)
+               {
+                       Float_t ptspectra=h->GetBinCenter(ibin);
+                       if(ptspectra<tcutsdata->GetPtTOFMatching())
+                               continue;
+                       h->SetBinContent(ibin,(h->GetBinContent(ibin)/factor));
+               }
+
+        }
+         else
+               return;                 
+
+
+}
+void TOFPIDsignalmatchingApply(TH1* h, Float_t factor)
+{
+       if(factor<0.0)
+               return;
+       for(Int_t ibin=1;ibin<h->GetNbinsX();ibin++)
+       {
+               Float_t ptspectra=h->GetBinCenter(ibin);
+               if(ptspectra<tcutsdata->GetPtTOFMatching())
+                       continue;
+               h->SetBinContent(ibin,(h->GetBinContent(ibin)*factor));
+
+       }
+
+}
+void CalculateDoubleCounts(TH1* doubleconunts,TH1** rawspectra,Int_t ipar, Bool_t dataflag)
+{
+       TH2F* tmphist=0x0;      
+       if (dataflag)
+               tmphist=(TH2F*)managerdata->GetGenMulvsRawMulHistogram("hHistDoubleCounts");    
+       else
+               tmphist=(TH2F*)managermc->GetGenMulvsRawMulHistogram("hHistDoubleCounts");
+
+       if(!tmphist)    
+               return;
+       TH1F* tmphist1=(TH1F*)rawspectra[ipar]->Clone("test");
+       tmphist1->Add(rawspectra[ipar+3]);
+       doubleconunts->Add(tmphist->ProjectionX(doubleconunts->GetName(),1,1));
+       if(ipar!=2)
+               doubleconunts->Add(tmphist->ProjectionX("pi+k",2,2));                   
+       if(ipar!=1)
+               doubleconunts->Add(tmphist->ProjectionX("pi+p",3,3));
+       if(ipar!=0)
+               doubleconunts->Add(tmphist->ProjectionX("k+p",4,4));
+       doubleconunts->Divide(doubleconunts,tmphist1,1,1,"B");
+
+       delete tmphist1;
+       
+
+}
+
+void CopyCorrectionFromFile(TString filename,TString correctionname,TH1F** corrtab)
+{
+       TFile* ftmp=TFile::Open(filename.Data());
+       TString tmp("tmp");
+       TList* ltmp=(TList*)ftmp->Get("output");
+
+       for(Int_t icharge=0;icharge<2;icharge++)
+        {
+                for(Int_t ipart=0;ipart<3;ipart++)
+                {
+                        Int_t index=ipart+3*icharge;
+                        tmp.Form("%s%s%s",correctionname.Data(),Particle[ipart].Data(),Sign[icharge].Data());
+                       TH1* histtmp=(TH1*)ltmp->FindObject(tmp.Data());
+                       if(!histtmp)
+                               continue;
+                       for(Int_t ibin=1;ibin<=histtmp->GetNbinsX();ibin++)
+                       {
+                               Float_t content=histtmp->GetBinContent(ibin);
+                               Float_t error=histtmp->GetBinError(ibin);
+                               Int_t bin=corrtab[index]->FindBin(histtmp->GetBinCenter(ibin));
+                               corrtab[index]->SetBinContent(bin,content);
+                               corrtab[index]->SetBinError(bin,error);
+                               
+                       }  
+               }
+       }
+       delete ltmp;
+       ftmp->Close();
+}