]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGLF/SPECTRA/PiKaPr/TestAOD/AnalysisBoth.C
TENDER becomes Tender
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / AnalysisBoth.C
index 0c352e1d934db05862ba539d8b39bb21645a0c06..d1a198e07235b8a59fd62c116a1c4c1bf4fc2380 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;
@@ -24,7 +26,11 @@ 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,
@@ -43,13 +49,19 @@ enum {
  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                                               
+ 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,Bool_t mcasdata=false);
 void AnalysisBoth (UInt_t options=0xF,TString outdate, TString outnamedata, TString outnamemc="",TString configfile="" )
 {
-       gStyle->SetOptStat(0);  
+gStyle->SetOptStat(0); 
        TH1::AddDirectory(kFALSE);
        gSystem->Load("libCore.so");  
        gSystem->Load("libPhysics.so");
@@ -61,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");
@@ -124,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);
        
@@ -168,8 +182,25 @@ void AnalysisBoth (UInt_t options=0xF,TString outdate, TString outnamedata, TStr
        }
        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);
@@ -229,8 +260,9 @@ void AnalysisBoth (UInt_t options=0xF,TString outdate, TString outnamedata, TStr
                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();
@@ -254,21 +286,30 @@ void AnalysisBoth (UInt_t options=0xF,TString outdate, TString outnamedata, TStr
                corrLeonardo[i]->Divide(corrLeonardo[i],rawspectramc[i],1,1,"B");
                
                
-               
-               //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]);
+               if((options&kusePIDcontaminatiofromfile)==kusePIDcontaminatiofromfile)
+               {
+                       CopyCorrectionFromFile(filenames[1],"conPID",contPID);
+                       CopyCorrectionFromFile(filenames[1],"conPIDprimary",contPIDpri);
 
-               //RecomputeErrors(contPIDpri[i]);
-               contPIDpri[i]->Divide(contPIDpri[i],rawspectramc[i],1,1,"B");
-       
-               eff[i]->Divide(eff[i],MCTruth[i],1,1,"B");
+               }
+               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();
@@ -282,7 +323,7 @@ void AnalysisBoth (UInt_t options=0xF,TString outdate, TString outnamedata, TStr
                contPID[i]->Sumw2();
                contPID[i]->Divide(contPID[i],rawspectramc[i],1,1,"B");
        
-               if(options&kuseprimaryPIDcont)
+               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"));
@@ -293,7 +334,21 @@ void AnalysisBoth (UInt_t options=0xF,TString outdate, TString outnamedata, TStr
                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");
@@ -323,6 +378,11 @@ void AnalysisBoth (UInt_t options=0xF,TString outdate, TString outnamedata, TStr
                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]);
+               }
        }
        outdate.ReplaceAll("/","_");
        configfile.ReplaceAll(".","_");
@@ -331,11 +391,23 @@ void AnalysisBoth (UInt_t options=0xF,TString outdate, TString outnamedata, TStr
                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)
+       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))
                        {
                                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]);
@@ -357,16 +429,33 @@ void AnalysisBoth (UInt_t options=0xF,TString outdate, TString outnamedata, TStr
                        {
                                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();
@@ -380,7 +469,6 @@ void AnalysisBoth (UInt_t options=0xF,TString outdate, TString outnamedata, TStr
        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++)
@@ -395,7 +483,7 @@ void AnalysisBoth (UInt_t options=0xF,TString outdate, TString outnamedata, TStr
                }
        }       
        allch->Scale(vertexcorrection);
-       spectraall->Scale(vertexcorrection/1.6);
+       spectraall->Scale(vertexcorrection/etacut);
 
        //spectraall->Scale(1.0/1.6);
        lout->Write("output",TObject::kSingleKey);      
@@ -467,7 +555,7 @@ Bool_t   OpenFile(TString dirname,TString outputname, Bool_t mcflag, Bool_t mcas
                if(!managerdata||!ecutsdata||!tcutsdata)
                        return false;
                if(managerdata->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries()!=ecutsdata->GetHistoCuts()->GetBinContent(3))
-                       cout<<"Please check DATA file possible problem with merging"<<" "<<anagerdata->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries()<<" "<<ecutsdata->GetHistoCuts()->GetBinContent(3)<<endl;
+                       cout<<"Please check DATA file possible problem with merging"<<" "<<managerdata->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries()<<" "<<ecutsdata->GetHistoCuts()->GetBinContent(3)<<endl;
 
        }
        return true;
@@ -502,7 +590,7 @@ 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));
+//                                     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++)
@@ -510,10 +598,10 @@ TH1F* GetOneHistFromPtDCAhisto(TString name,TString hnameout,AliSpectraBothHisto
                                                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<<"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;
+//                                     cout<<testyield<<" "<<TMath::Sqrt(testerror)<<" error2 "<<inyield<<" "<<TMath::Sqrt(inyield)<<endl;
                                                
                                        //histo->SetBinContent(ibin,inyield);
                                        //histo->SetBinError(ibin,TMath::Sqrt(inyield));
@@ -531,8 +619,8 @@ 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++)
@@ -544,7 +632,7 @@ void GetPtHistFromPtDCAhisto(TString hnamein, TString hnameout, AliSpectraBothHi
                        
                        
                        histo[index]=GetOneHistFromPtDCAhisto(nameinfinal,nameoutfinal,hman,dcacutxy);
-                       CleanHisto(histo[index],min[ipart],max[ipart]);
+                       CleanHisto(histo[index],minRanges[ipart],maxRanges[ipart]);
                }
        } 
 }
@@ -559,7 +647,7 @@ void CleanHisto(TH1F* h, Float_t minV, Float_t maxV,TH1* contpid=0x0)
                }       
                if(contpid)
                {
-                       if(contpid->GetBinContent(i)>0.201)
+                       if(contpid->GetBinContent(i)>fMaxContaminationPIDMC)
                        {
                                h->SetBinContent(i,0);
                                h->SetBinError(i,0);
@@ -604,6 +692,7 @@ void DCACorrectionMarek(AliSpectraBothHistoManager* hman_data, AliSpectraBothHis
                                                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);
@@ -618,9 +707,11 @@ 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;
+                                       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;
                                        hToFit->Sumw2();
@@ -678,11 +769,14 @@ void DCACorrectionMarek(AliSpectraBothHistoManager* hman_data, AliSpectraBothHis
                                        gPad->SetLogy();
         
                                        TObjArray *mc=0x0;
+                                       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);
                                        if(fitsettings&0x1)
@@ -691,6 +785,18 @@ void DCACorrectionMarek(AliSpectraBothHistoManager* hman_data, AliSpectraBothHis
                                                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
+                                       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)
@@ -699,7 +805,8 @@ void DCACorrectionMarek(AliSpectraBothHistoManager* hman_data, AliSpectraBothHis
                                        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(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;
@@ -816,17 +923,21 @@ void DCACorrectionMarek(AliSpectraBothHistoManager* hman_data, AliSpectraBothHis
                                                PrimMCPred->Scale(v1/PrimMCPred->Integral(PrimMCPred->GetXaxis()->FindBin(FitRange[0]),PrimMCPred->GetXaxis()->FindBin(FitRange[1])));
 
                                                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->SetLineColor(kBlack);
-                                               Leg1->AddEntry(result,"result","lp");
-                                               result->DrawClone("lhistsame");
+                                               Leg1->AddEntry(result,"fit result","l");
+                                               result->DrawClone("histsame");
                                        
                                                PrimMCPred->SetLineColor(kGreen+2);
                                                PrimMCPred->SetLineStyle(2);
                                                 PrimMCPred->SetLineWidth(3.0);
-                                               Leg1->AddEntry(PrimMCPred,"Prmi.","l");
-                                               PrimMCPred->DrawClone("lhistsame");
+                                               Leg1->AddEntry(PrimMCPred,"primaries","l");
+                                               PrimMCPred->DrawClone("histsame");
                                                if(fitsettings&0x1)
                                                {
 
@@ -835,8 +946,8 @@ void DCACorrectionMarek(AliSpectraBothHistoManager* hman_data, AliSpectraBothHis
                                                        secStMCPred->SetLineWidth(3.0);
 
                                                        secStMCPred->SetLineStyle(3);
-                                                       Leg1->AddEntry(secStMCPred,"Sec.WD","l");
-                                                       secStMCPred->DrawClone("lhistsame");
+                                                       Leg1->AddEntry(secStMCPred,"weak decays","l");
+                                                       secStMCPred->DrawClone("histsame");
 
                                                }
                                                if(fitsettings&0x2)
@@ -847,8 +958,8 @@ void DCACorrectionMarek(AliSpectraBothHistoManager* hman_data, AliSpectraBothHis
                                                        secMCPred->SetLineWidth(3.0);
 
                                                        secMCPred->SetLineStyle(4);     
-                                                       Leg1->AddEntry(secMCPred,"Sec.Mat","l");
-                                                       secMCPred->DrawClone("lhistsame");
+                                                       Leg1->AddEntry(secMCPred,"material","l");
+                                                       secMCPred->DrawClone("histsame");
            
                                                }   
                                        }                               
@@ -863,7 +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->Draw();
+                                       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();
@@ -938,11 +1053,11 @@ 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());
+               TFile *fGeanFlukaK= TFile::Open(fnameGeanFlukaK.Data());
                if (!fGeanFlukaK)
                {
-                       fnameGeanFlukaK="$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/TestAOD/GFCorrection/correctionForCrossSection.321.root";
-                       fGeanFlukaK= new TFile(fnameGeanFlukaK.Data());
+                       fnameGeanFlukaK="$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/TestAOD/correctionForCrossSection.321.root";
+                       fGeanFlukaK= TFile::Open(fnameGeanFlukaK.Data());
                        if (!fGeanFlukaK)
                                return;
                }
@@ -957,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)
@@ -1000,11 +1118,11 @@ void GFCorrection(TH1F **Spectra,Float_t tofpt,UInt_t options)
          Int_t kPos=0;
          Int_t kNeg=1;
           TString fnameGFProtons= "GFCorrection/correctionForCrossSection.root";
-         TFile* fGFProtons = new TFile (fnameGFProtons.Data());
+         TFile* fGFProtons = TFile::Open(fnameGFProtons.Data());
          if (!fGFProtons)
          { 
-               fnameGFProtons=="$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/TestAOD/GFCorrection/correctionForCrossSection.root";
-               TFile* fGFProtons = new TFile (fnameGFProtons.Data());
+               fnameGFProtons="$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/TestAOD/correctionForCrossSection.root";
+               fGFProtons = TFile::Open(fnameGFProtons.Data());
                if (!fGFProtons)
                        return;
          }
@@ -1026,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)
@@ -1222,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();
@@ -1238,9 +1358,9 @@ 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;
+                               conver=conver/etacut;
                                //cout<<maxy<<" "<<conver<<" "<<masstmp<<""<<spectra[i]->GetName()<<endl;
                                Float_t bincontent=allch->GetBinContent(j);
                                Float_t binerror=allch->GetBinError(j);
@@ -1331,7 +1451,7 @@ Float_t Normaliztionwithbin0integrals(UInt_t options)
        cout<<histodata->GetBinContent(2)<<" "<<histodata->GetBinContent(4)<<endl;
        if ((options&knormalizationwithbin0integralsdata)==knormalizationwithbin0integralsdata)
                return  dataevents+bin0mcRec->Integral(binmin,binmax);
-       else if ((options&kknormalizationwithbin0integralsMC)==knormalizationwithbin0integralsdata)
+       else if ((options&kknormalizationwithbin0integralsMC)==knormalizationwithbin0integralsMC)   
                return dataevents+bin0mcMC->Integral(binmin,binmax) ;
        else
                return 1;               
@@ -1343,7 +1463,7 @@ Bool_t ReadConfigFile(TString configfile)
        ifstream infile(configfile.Data());
        if(infile.is_open()==false)
                return false;
-       TString namesofSetting[28]={"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"};   
+       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)
@@ -1409,9 +1529,42 @@ Bool_t ReadConfigFile(TString configfile)
                        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
+               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;
@@ -1439,7 +1592,9 @@ Bool_t ReadConfigFile(TString configfile)
                        return false;
                }
        }
-       
+       for(int i=0;i<3;i++)
+               if(minRanges[i]>maxRanges[i])
+                       return false;
                
 
        return true;
@@ -1457,4 +1612,93 @@ void SubHistWithFullCorr(TH1F* h1, TH1F* h2, Float_t factor1=1.0, Float_t factor
                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();
+}