PWGLF/SPECTRA/PiKaPr/TestAOD/AnalysisBoth.C

   1 class AliSpectraBothHistoManager;
   2 class AliSpectraBothEventCuts;
   3 class AliSpectraBothTrackCuts;
   4 TString Charge[]={"Pos","Neg"};
   5 TString Sign[]={"Plus","Minus"};
   6 TString Particle[]={"Pion","Kaon","Proton"};
   7 AliSpectraBothHistoManager* managerdata=0x0;
   8 AliSpectraBothEventCuts* ecutsdata=0x0;
   9 AliSpectraBothTrackCuts* tcutsdata=0x0;
  10
  11 AliSpectraBothHistoManager* managermc=0x0;
  12 AliSpectraBothEventCuts* ecutsmc=0x0;
  13 AliSpectraBothTrackCuts* tcutsmc=0x0;
  14
  15 Float_t TOFMatchingScalling[2]={-1,-1};
  16 Int_t Color[3]={1,2,4};
  17 Int_t Marker[6]={20,21,22,24,25,26};
  18 Double_t Range[3]={0.3,0.3,0.5}; // LowPt range for pi k p
  19
  20
  21 Double_t FitRange[2]={-2.0,2.0};
  22 Double_t CutRange[2]={-3.0,3.0};
  23 Double_t minptformaterial[6]={0.0,0.2,0.0,0.0,0.2,0.0};
  24 Double_t maxptformaterial[6]={0.0,0.6,1.3,0.0,0.6,0.0};
  25 Double_t minptforWD[6]={0.2,100.0,0.3,0.2,100.0,0.3};
  26 Double_t maxptforWD[6]={1.5,-100.0,2.0,1.5,-100.0,2.0};
  27
  28
  29 enum ECharge_t {
  30   kPositive,
  31   kNegative,
  32   kNCharges
  33 };
  34 enum {
  35  kdodca=0x1, //dca fits are made
  36  kgeantflukaKaon=0x2,// geant fluka correction is used for kaons
  37  kgeantflukaProton=0x4, // geant fluka correction is used for protons and antiprotons
  38  knormalizationtoeventspassingPhySel=0x8,// spectra are divided by number of events passing physic selection
  39  kveretxcorrectionandbadchunkscorr=0x10, // correction for difference in z vertex distribution in data and MC and correction for bad chunks is applied
  40  kmcisusedasdata=0x20, // the result of the looping over MC is used as data input
  41  kdonotusedcacuts=0x40, // allows to use the constant dca cut for all pt bins not the pt dependet defined in stardrad track cuts 2011
  42  kuseprimaryPIDcont=0x80, //pid contamination is calculated using only primiary particle in this case K should use dca fits
  43  knormalizationwithbin0integralsdata=0x100, // the normalization factor is calcualte using integral over z vertex distributions (in this case reconstructed vertex disitrbution uses z vertex for data)
  44  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
  45  kuserangeonfigfile=0x400, // use of config file for dca fit settings
  46  kskipconcutonspectra=0x800 //do not use conPID<02 cut  useful for syst. studies
  47 };
  48
  49 Bool_t OpenFile(TString dirname, TString outputname, Bool_t mcflag,Bool_t mcasdata=false);
  50 void AnalysisBoth (UInt_t options=0xF,TString outdate, TString outnamedata, TString outnamemc="",TString configfile="" )
  51 {
  52         gStyle->SetOptStat(0);
  53         TH1::AddDirectory(kFALSE);
  54         gSystem->Load("libCore.so");
  55         gSystem->Load("libPhysics.so");
  56         gSystem->Load("libTree");
  57         gSystem->Load("libMatrix");
  58         gSystem->Load("libSTEERBase");
  59         gSystem->Load("libESD");
  60         gSystem->Load("libAOD");
  61         gSystem->Load("libANALYSIS");
  62         gSystem->Load("libOADB");
  63         gSystem->Load("libANALYSISalice");
  64         gSystem->Load("libTENDER");
  65         gSystem->Load("libCORRFW");
  66         gSystem->Load("libPWGTools");
  67         gSystem->Load("libPWGLFspectra");
  68
  69         gROOT->LoadMacro("$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/TestAOD/QAPlotsBoth.C");
  70         Double_t mass[3];
  71         mass[0]   = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
  72         mass[1]   = TDatabasePDG::Instance()->GetParticle("K+")->Mass();
  73         mass[2] = TDatabasePDG::Instance()->GetParticle("proton")->Mass();
  74
  75         TFormula* dcacutxy=0x0;
  76         if(!(options&kdonotusedcacuts))
  77         {
  78
  79                 AliESDtrackCuts* esdtrackcuts= AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE);
  80                 TString formulastring(esdtrackcuts->GetMaxDCAToVertexXYPtDep());
  81                 formulastring.ReplaceAll("pt","x");
  82                 dcacutxy=new TFormula("dcacutxy",formulastring.Data());
  83         }
  84         if(options&kuserangeonfigfile)
  85                 if(!ReadConfigFile(configfile))
  86                         return;
  87         TList* lout=new TList();
  88
  89
  90         TString indirname=Form("/output/train%s",outdate.Data());
  91         //TString indirname("/output/train24072012");
  92         if(outnamemc.Length()==0)
  93         outnamemc=outnamedata;
  94         cout<<indirname.Data()<<" "<<outnamemc.Data()<<endl;
  95         // Do the job
  96
  97
  98         OpenFile(indirname,outnamemc,true);
  99         OpenFile(indirname,outnamedata,false,((Bool_t)(options&kmcisusedasdata)));
 100         if(!managermc||!managerdata)
 101         {
 102                 cout<<managermc<<" "<<managerdata<<endl;
 103                 return;
 104         }
 105         TH1F* rawspectradata[6];
 106         TH1F* rawspectramc[6];
 107         TH1F* MCTruth[6];
 108         TH1F* eff[6];
 109         TH1F* contallMC[6];
 110         TH1F* contPID[6];
 111         TH1F* contWD[6];
 112         TH1F* contMat[6];
 113         TH1F* confinal[6];
 114         TH1F* contPIDpri[6];
 115         TH1F* contSecMC[6];
 116
 117         TH1F* contfit[12];
 118         TH1F* contWDfit[12];
 119         TH1F* contMatfit[12];
 120         TH1F* primaryfit[12];
 121
 122
 123
 124         TH1F* spectra[6];
 125         TH1F* spectraLeonardo[6];
 126
 127         TH1F* corrLeonardo[6];
 128         //GetSpectra(managerdata,rawspectradata,true);
 129         //GetSpectra(managermc,rawspectramc,true,true);
 130
 131         GetPtHistFromPtDCAhisto("hHistPtRecSigma","SpectraMC",managermc,rawspectramc,dcacutxy);
 132         GetPtHistFromPtDCAhisto("hHistPtRecSigma","SpectraDATA",managerdata,rawspectradata,dcacutxy);
 133         GetPtHistFromPtDCAhisto("hHistPtRecTruePrimary","eff",managermc,eff,dcacutxy);
 134         GetPtHistFromPtDCAhisto("hHistPtRecTrue","conPID",managermc,contPID,dcacutxy);
 135         GetPtHistFromPtDCAhisto("hHistPtRecSigmaSecondaryWeakDecay","conWD",managermc,contWD,dcacutxy);
 136         GetPtHistFromPtDCAhisto("hHistPtRecSigmaSecondaryMaterial","conMat",managermc,contMat,dcacutxy);
 137         GetPtHistFromPtDCAhisto("hHistPtRecSigmaPrimary","conPIDprimary",managermc,contPIDpri,dcacutxy);
 138
 139
 140         Double_t neventsmcall = 1 ;  //if loop over MC is done after or befor events cuts this will be changed
 141         Double_t neventsdata =  1;
 142         Double_t neventsmc =  1;
 143
 144         //Normaliztion of MCtruth depends if the loop was done after of before ESD event cuts.
 145         //In currect code this cannot be check on the level of macro.
 146         //If the loop was done before MC should be done to all processed events (NumberOfProcessedEvents())
 147         //If loop was done after MC should be normalized to all accepted events (NumberOfEvents())
 148         // The option one will be alaways use.
 149
 150         neventsmcall= ecutsmc->NumberOfProcessedEvents();
 151         if(options&knormalizationtoeventspassingPhySel)
 152         {
 153                 //neventsmcall= ecutsmc->NumberOfProcessedEvents();
 154                  neventsdata=ecutsdata->NumberOfPhysSelEvents();
 155                  neventsmc=ecutsmc->NumberOfPhysSelEvents();
 156         }
 157         else if ((options&knormalizationwithbin0integralsdata)||(options&knormalizationwithbin0integralsMC))
 158         {
 159                 neventsdata=Normaliztionwithbin0integrals(options);
 160                 neventsmc=ecutsmc->NumberOfPhysSelEvents();
 161         }
 162         else
 163         {
 164                 neventsdata=ecutsdata->NumberOfEvents(); //number of accepted events
 165                  neventsmc=ecutsmc->NumberOfEvents();
 166                 neventsmcall= ecutsmc->NumberOfEvents();
 167
 168         }
 169         GetMCTruth(MCTruth);
 170         cout<<neventsdata<<" Events"<<endl;
 171
 172
 173         TH1F* allgen=((TH1F*)managermc->GetPtHistogram1D("hHistPtGen",1,1))->Clone();
 174         allgen->SetName("AllGen");
 175         TH1F* allrecMC=GetOneHistFromPtDCAhisto("hHistPtRec","rawallMC",managermc,dcacutxy);
 176         TH1F* alleff=GetOneHistFromPtDCAhisto("hHistPtRecPrimary","effall",managermc,dcacutxy);
 177         TH1F* allrecdata=GetOneHistFromPtDCAhisto("hHistPtRec","rawalldata",managerdata,dcacutxy);
 178
 179
 180
 181
 182         TH1F* spectraall=(TH1F*)allrecdata->Clone("recNch");
 183         spectraall->SetTitle("recNch");
 184         TH1F* contall=(TH1F*)allrecMC->Clone("contall");
 185         contall->SetTitle("contall");
 186         //contall->Add(alleff,-1);
 187         SubHistWithFullCorr(contall,alleff);
 188         alleff->Divide(alleff,allgen,1,1,"B");
 189         contall->Divide(contall,allrecMC,1,1,"B");
 190
 191         GetCorrectedSpectra(spectraall,allrecdata,alleff,contall);
 192         Divideby2pipt(spectraall);
 193
 194         allrecdata->Scale(1./neventsdata,"width");
 195         allgen->Scale(1./neventsmcall,"width");
 196         allrecMC->Scale(1./neventsmc,"width");
 197         spectraall->Scale(1./neventsdata,"width");
 198
 199
 200         lout->Add(allgen);
 201         lout->Add(allrecMC);
 202         lout->Add(alleff);
 203         lout->Add(allrecdata);
 204         lout->Add(spectraall);
 205         lout->Add(contall);
 206
 207         for (int i=0;i<6;i++)
 208         {
 209
 210
 211                 TString tmpname(rawspectramc[i]->GetTitle());
 212                 tmpname.ReplaceAll("SpectraMC","%s");
 213                 contallMC[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contallMC"));
 214                 contallMC[i]->SetTitle(Form(tmpname.Data(),"contallMC"));
 215                 contfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contfit"));
 216                 contfit[i]->SetTitle(Form(tmpname.Data(),"contfit"));
 217                 contWDfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contWDfit"));
 218                 contWDfit[i]->SetTitle(Form(tmpname.Data(),"contWDfit"));
 219                 contMatfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contMatfit"));
 220                 contMatfit[i]->SetTitle(Form(tmpname.Data(),"contMatfit"));
 221                 primaryfit[i]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"primaryfit"));
 222                 primaryfit[i]->SetTitle(Form(tmpname.Data(),"primaryfit"));
 223                 contfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contfitonMC"));
 224                 contfit[i+6]->SetTitle(Form(tmpname.Data(),"contfitonMC"));
 225                 contWDfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contWDfitonMC"));
 226                 contWDfit[i+6]->SetTitle(Form(tmpname.Data(),"contWDfitonMC"));
 227                 contMatfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"contMatfitonMC"));
 228                 contMatfit[i+6]->SetTitle(Form(tmpname.Data(),"contMatfitonMC"));
 229                 primaryfit[i+6]=(TH1F*)rawspectramc[i]->Clone(Form(tmpname.Data(),"primaryfitMC"));
 230                 primaryfit[i+6]->SetTitle(Form(tmpname.Data(),"primaryfitMC"));
 231
 232
 233
 234                 contfit[i]->Reset();
 235                 contWDfit[i]->Reset();
 236                 contMatfit[i]->Reset();
 237                 primaryfit[i]->Reset();
 238
 239
 240                 contfit[i+6]->Reset();
 241                 contWDfit[i+6]->Reset();
 242                 contMatfit[i+6]->Reset();
 243                 primaryfit[i+6]->Reset();
 244
 245                 SetBintoOne(primaryfit[i]);
 246                 SetBintoOne(primaryfit[i+6]);
 247                 spectra[i]=(TH1F*)rawspectradata[i]->Clone(Form(tmpname.Data(),"SpectraFinal"));
 248                 spectra[i]->SetTitle(Form(tmpname.Data(),"SpectraFinal"));
 249                 spectraLeonardo[i]=(TH1F*)rawspectradata[i]->Clone(Form(tmpname.Data(),"SpectraFinalLeonardo"));
 250                 spectraLeonardo[i]->SetTitle(Form(tmpname.Data(),"SpectraFinalLeonardo"));
 251
 252                 corrLeonardo[i]=(TH1F*)MCTruth[i]->Clone(Form(tmpname.Data(),"CorrFactLeonardo"));
 253                 corrLeonardo[i]->SetTitle(Form(tmpname.Data(),"CorrFactLeonardo"));
 254                 corrLeonardo[i]->Divide(corrLeonardo[i],rawspectramc[i],1,1,"B");
 255
 256
 257
 258                 //contallMC[i]->Add(eff[i],-1.0);
 259                 SubHistWithFullCorr(contallMC[i],eff[i]);
 260                 //RecomputeErrors(contallMC[i]);
 261                 contallMC[i]->Sumw2();
 262                 contallMC[i]->Divide(contallMC[i],rawspectramc[i],1,1,"B");
 263
 264                 // contamintaion from PID but only primaries
 265                 //contPIDpri[i]->Add(eff[i],-1.0);
 266                 SubHistWithFullCorr(contPIDpri[i],eff[i]);
 267
 268                 //RecomputeErrors(contPIDpri[i]);
 269                 contPIDpri[i]->Divide(contPIDpri[i],rawspectramc[i],1,1,"B");
 270
 271                 eff[i]->Divide(eff[i],MCTruth[i],1,1,"B");
 272
 273
 274                 contPID[i]->Sumw2();
 275                 rawspectramc[i]->Sumw2();
 276                 //contPID[i]->Add(contPID[i],rawspectramc[i],-1,1);
 277                 SubHistWithFullCorr(contPID[i],rawspectramc[i]);
 278                 contPID[i]->Scale(-1.0);
 279
 280                 //RecomputeErrors(contPID[i]);
 281                 contPID[i]->ResetStats();
 282                 contPID[i]->Sumw2();
 283                 contPID[i]->Divide(contPID[i],rawspectramc[i],1,1,"B");
 284
 285                 if(options&kuseprimaryPIDcont)
 286                         confinal[i]=(TH1F*)contPIDpri[i]->Clone(Form(tmpname.Data(),"confinal"));
 287                 else
 288                         confinal[i]=(TH1F*)contPID[i]->Clone(Form(tmpname.Data(),"confinal"));
 289                 confinal[i]->SetTitle(Form(tmpname.Data(),"confinal"));
 290
 291                 contSecMC[i]=(TH1F*)contWD[i]->Clone(Form(tmpname.Data(),"conSecMC"));
 292                 contSecMC[i]->Add(contMat[i]);
 293                 contWD[i]->Divide(contWD[i],rawspectramc[i],1,1,"B");
 294                 contMat[i]->Divide(contMat[i],rawspectramc[i],1,1,"B");
 295                 contSecMC[i]->Divide(contSecMC[i],rawspectramc[i],1,1,"B");
 296
 297
 298                 rawspectradata[i]->Scale(1./neventsdata,"width");
 299                 rawspectramc[i]->Scale(1./neventsmc,"width");
 300                 MCTruth[i]->Scale(1./neventsmcall,"width");
 301                 spectraLeonardo[i]->Scale(1./neventsdata,"width");
 302
 303
 304
 305                 lout->Add(rawspectradata[i]);
 306                 lout->Add(rawspectramc[i]);
 307                 lout->Add(MCTruth[i]);
 308                 lout->Add(eff[i]);
 309                 lout->Add(contallMC[i]);
 310                 lout->Add(contPID[i]);
 311                 lout->Add(contWD[i]);
 312                 lout->Add(contMat[i]);
 313                 lout->Add(contfit[i]);
 314                 lout->Add(contWDfit[i]);
 315                 lout->Add(contMatfit[i]);
 316                 lout->Add(primaryfit[i]);
 317                 lout->Add(contfit[i+6]);
 318                 lout->Add(contWDfit[i+6]);
 319                 lout->Add(contMatfit[i+6]);
 320                 lout->Add(primaryfit[i+6]);
 321                 lout->Add(spectra[i]);
 322                 lout->Add(spectraLeonardo[i]);
 323                 lout->Add(confinal[i]);
 324                 lout->Add(contPIDpri[i]);
 325                 lout->Add(contSecMC[i]);
 326         }
 327         outdate.ReplaceAll("/","_");
 328         configfile.ReplaceAll(".","_");
 329         TFile* fout=0x0;
 330         if(configfile.Length()>0&&(options&kuserangeonfigfile))
 331                 fout=new TFile(Form("./results/ResMY_%s_%s_%#X_%s.root",outnamemc.Data(),outdate.Data(),options,configfile.Data()),"RECREATE");
 332         else
 333                 fout=new TFile(Form("./results/ResMY_%s_%s_%#X.root",outnamemc.Data(),outdate.Data(),options),"RECREATE");
 334         if (options&kdodca)
 335                 DCACorrectionMarek(managerdata,managermc,dcacutxy,fout,contfit,contWDfit,contMatfit,primaryfit);
 336         for (int i=0;i<6;i++)
 337         {
 338                         if(options&kdodca)
 339                         {
 340                                 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
 341                                         confinal[i]->Add(contfit[i]);
 342                                 GetCorrectedSpectra(spectra[i],rawspectradata[i],eff[i],confinal[i]);
 343                         }
 344                         else
 345                         {
 346                                 GetCorrectedSpectra(spectra[i],rawspectradata[i],eff[i],contallMC[i]);
 347                         }
 348                         GetCorrectedSpectraLeonardo(spectraLeonardo[i],corrLeonardo[i],primaryfit[i],primaryfit[i+6]);
 349                         if(options&kskipconcutonspectra)
 350                                 continue;
 351                         if(options&kuseprimaryPIDcont)
 352                         {
 353                                 CleanHisto(spectra[i],-1,100,contPIDpri[i]);
 354                                 CleanHisto(spectraLeonardo[i],-1,100,contPIDpri[i]);
 355                         }
 356                         else
 357                         {
 358                                 CleanHisto(spectra[i],-1,100,contPID[i]);
 359                                 CleanHisto(spectraLeonardo[i],-1,100,contPID[i]);
 360                         }
 361         }
 362
 363         GFCorrection(spectra,tcutsdata->GetPtTOFMatching(),options);
 364         GFCorrection(spectraLeonardo,tcutsdata->GetPtTOFMatching(),options);
 365          MatchingTOFEff(spectra,lout);
 366           MatchingTOFEff(spectraLeonardo);
 367         TH1F* allch=GetSumAllCh(spectra,mass);
 368         lout->Add(allch);
 369
 370 //      lout->ls();
 371         fout->cd();
 372         TList* listqa=new TList();
 373         TList* canvaslist=new TList();
 374         Float_t vertexcorrection=1.0;
 375         Float_t corrbadchunksvtx=1.0;
 376         if ((options&knormalizationwithbin0integralsdata)||(options&knormalizationwithbin0integralsMC))
 377                 corrbadchunksvtx=QAPlotsBoth(managerdata,managermc,ecutsdata,ecutsmc,tcutsdata,tcutsmc,listqa,canvaslist,0);
 378         else
 379                 corrbadchunksvtx=QAPlotsBoth(managerdata,managermc,ecutsdata,ecutsmc,tcutsdata,tcutsmc,listqa,canvaslist,1);
 380         if (options&kveretxcorrectionandbadchunkscorr)
 381                 vertexcorrection=corrbadchunksvtx;
 382         cout<<" VTX corr="<<vertexcorrection<<endl;
 383         Double_t ycut=tcutsdata->GetY();
 384         if(TMath::Abs(ycut)>0.0)
 385                 vertexcorrection=vertexcorrection/(2.0*ycut);
 386         for (int i=0;i<6;i++)
 387         {
 388                 spectra[i]->Scale(vertexcorrection);
 389                 spectraLeonardo[i]->Scale(vertexcorrection);
 390                 if(TMath::Abs(ycut)>0.0)
 391                 {
 392                         rawspectradata[i]->Scale(1.0/(2.0*ycut));
 393                         rawspectramc[i]->Scale(1.0/(2.0*ycut));
 394                         MCTruth[i]->Scale(1.0/(2.0*ycut));
 395                 }
 396         }
 397         allch->Scale(vertexcorrection);
 398         spectraall->Scale(vertexcorrection/1.6);
 399
 400         //spectraall->Scale(1.0/1.6);
 401         lout->Write("output",TObject::kSingleKey);
 402         listqa->Write("outputQA",TObject::kSingleKey);
 403         canvaslist->Write("outputcanvas",TObject::kSingleKey);
 404
 405         fout->Close();
 406         //Normaliztionwithbin0integrals();
 407
 408 }
 409
 410 Bool_t   OpenFile(TString dirname,TString outputname, Bool_t mcflag, Bool_t mcasdata)
 411 {
 412
 413
 414         TString nameFile = Form("./%s/AnalysisResults%s.root",dirname.Data(),(mcflag?"MC":"DATA"));
 415         TFile *file = TFile::Open(nameFile.Data());
 416         if(!file)
 417         {
 418                 cout<<"no file"<<endl;
 419                 return false;
 420         }
 421         TString sname=Form("OutputBothSpectraTask_%s_%s",(mcflag?"MC":"Data"),outputname.Data());
 422         if(mcasdata)
 423         {
 424                 cout<<"using MC as data "<<endl;
 425                 sname=Form("OutputBothSpectraTask_%s_%s","MC",outputname.Data());
 426         }
 427         file->ls();
 428         TDirectoryFile *dir=(TDirectoryFile*)file->Get(sname.Data());
 429         if(!dir)
 430         {
 431         //      cout<<"no dir "<<sname.Data()<<endl;
 432                 if(mcasdata)
 433                 {
 434                         cout<<"using MC as data "<<endl;
 435                         sname=Form("OutputAODSpectraTask_%s_%s","MC",outputname.Data());
 436                 }
 437                 else
 438                         sname=Form("OutputAODSpectraTask_%s_%s",(mcflag?"MC":"Data"),outputname.Data());
 439         //      cout<<"trying "<<sname.Data()<<endl;
 440                 dir=(TDirectoryFile*)file->Get(sname.Data());
 441                 if(!dir)
 442                 {
 443                         cout<<"no dir "<<sname.Data()<<endl;
 444                         return false;
 445                 }
 446         }
 447         cout << " -- Info about " <<(mcflag?"MC":"DATA") <<" -- "<< endl;
 448         if(mcflag)
 449         {
 450                 managermc= (AliSpectraBothHistoManager*) dir->Get("SpectraHistos");
 451                 ecutsmc = (AliSpectraBothEventCuts*) dir->Get("Event Cuts");
 452                 tcutsmc = (AliSpectraBothTrackCuts*) dir->Get("Track Cuts");
 453                 ecutsmc->PrintCuts();
 454                 tcutsmc->PrintCuts();
 455                 if(!managermc||!ecutsmc||!tcutsmc)
 456                         return false;
 457                 if(managermc->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries()!=ecutsmc->GetHistoCuts()->GetBinContent(3))
 458                         cout<<"Please check MC file possible problem with merging"<<" "<<managermc->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries()<<" "<<ecutsmc->GetHistoCuts()->GetBinContent(3)<<endl;
 459         }
 460         else
 461         {
 462                 managerdata= (AliSpectraBothHistoManager*) dir->Get("SpectraHistos");
 463                 ecutsdata = (AliSpectraBothEventCuts*) dir->Get("Event Cuts");
 464                 tcutsdata = (AliSpectraBothTrackCuts*) dir->Get("Track Cuts");
 465                 ecutsdata->PrintCuts();
 466                 tcutsdata->PrintCuts();
 467                 if(!managerdata||!ecutsdata||!tcutsdata)
 468                         return false;
 469                 if(managerdata->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries()!=ecutsdata->GetHistoCuts()->GetBinContent(3))
 470                         cout<<"Please check DATA file possible problem with merging"<<" "<<anagerdata->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries()<<" "<<ecutsdata->GetHistoCuts()->GetBinContent(3)<<endl;
 471
 472         }
 473         return true;
 474 }
 475
 476  void GetMCTruth(TH1F** MCTruth)
 477  {
 478         for(Int_t icharge=0;icharge<2;icharge++)
 479         {
 480                 for(Int_t ipart=0;ipart<3;ipart++)
 481                 {
 482                         Int_t index=ipart+3*icharge;
 483                         TString hname=Form("hHistPtGenTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
 484                         MCTruth[index]=(TH1F*)((TH1F*)managermc->GetPtHistogram1D(hname.Data(),1,1))->Clone();
 485                         MCTruth[index]->SetName(Form("MCTruth_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
 486                         MCTruth[index]->SetTitle(Form("MCTruth_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
 487                         MCTruth[index]->Sumw2();
 488                 }
 489         }
 490 }
 491
 492 TH1F* GetOneHistFromPtDCAhisto(TString name,TString hnameout,AliSpectraBothHistoManager* hman,TFormula* dcacutxy)
 493 {
 494                         histo =(TH1F*)((TH1F*) hman->GetPtHistogram1D(name.Data(),-1,-1))->Clone();
 495                         histo->SetName(hnameout.Data());
 496                         histo->SetTitle(hnameout.Data());
 497
 498                         if(dcacutxy)
 499                         {
 500                                 for(int ibin=1;ibin<histo->GetNbinsX();ibin++)
 501                                 {
 502                                         Double_t lowedge=histo->GetBinLowEdge(ibin);
 503                                         Float_t cut=dcacutxy->Eval(lowedge);
 504                                         TH1F* dcahist=(TH1F*)hman->GetDCAHistogram1D(name.Data(),lowedge,lowedge);
 505                                         //Float_t inyield=dcahist->Integral(dcahist->GetXaxis()->FindBin(-1.0*cut),dcahist->GetXaxis()->FindBin(cut));
 506                                         Float_t testyield=0.0;
 507                                         Float_t testerror=0.0;
 508                                         for (int itest=dcahist->GetXaxis()->FindBin(-1.0*cut);itest<=dcahist->GetXaxis()->FindBin(cut);itest++)
 509                                         {
 510                                                 testyield+=dcahist->GetBinContent(itest);
 511                                                 testerror+=dcahist->GetBinError(itest)*dcahist->GetBinError(itest);
 512                                         }
 513                                         //cout<<"corr data "<<histo->GetBinContent(ibin)<<" "<<inyield<<" "<<dcahist->Integral()<<" "<<hnameout.Data()<<endl;
 514                                         //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;
 515
 516                                         //cout<<testyield<<" "<<TMath::Sqrt(testerror)<<" error2 "<<inyield<<" "<<TMath::Sqrt(inyield)<<endl;
 517
 518                                         //histo->SetBinContent(ibin,inyield);
 519                                         //histo->SetBinError(ibin,TMath::Sqrt(inyield));
 520                                         histo->SetBinContent(ibin,testyield);
 521                                         histo->SetBinError(ibin,TMath::Sqrt(testerror));
 522
 523                                 }
 524                         }
 525                         histo->Sumw2();
 526                         return histo;
 527 }
 528
 529
 530
 531
 532 void GetPtHistFromPtDCAhisto(TString hnamein, TString hnameout, AliSpectraBothHistoManager* hman,TH1F** histo,TFormula* dcacutxy)
 533 {
 534         Float_t min[3]={0.3,0.3,0.4};
 535         Float_t max[3]={1.5,1.2,2.2};
 536         for(Int_t icharge=0;icharge<2;icharge++)
 537         {
 538                 for(Int_t ipart=0;ipart<3;ipart++)
 539                 {
 540                         Int_t index=ipart+3*icharge;
 541                         Printf("Getting %s",hnamein.Data());
 542                         TString nameinfinal=Form("%s%s%s",hnamein.Data(),Particle[ipart].Data(),Sign[icharge].Data());
 543                         TString nameoutfinal=Form("%s%s%s",hnameout.Data(),Particle[ipart].Data(),Sign[icharge].Data());
 544
 545
 546                         histo[index]=GetOneHistFromPtDCAhisto(nameinfinal,nameoutfinal,hman,dcacutxy);
 547                         CleanHisto(histo[index],min[ipart],max[ipart]);
 548                 }
 549         }
 550 }
 551 void CleanHisto(TH1F* h, Float_t minV, Float_t maxV,TH1* contpid=0x0)
 552 {
 553         for (int i=0;i<=h->GetNbinsX();i++)
 554         {
 555                 if(h->GetXaxis()->GetBinCenter(i)<minV||h->GetXaxis()->GetBinCenter(i)>maxV)
 556                 {
 557                         h->SetBinContent(i,0);
 558                         h->SetBinError(i,0);
 559                 }
 560                 if(contpid)
 561                 {
 562                         if(contpid->GetBinContent(i)>0.201)
 563                         {
 564                                 h->SetBinContent(i,0);
 565                                 h->SetBinError(i,0);
 566                         }
 567                 }
 568         }
 569 }
 570
 571
 572 void DCACorrectionMarek(AliSpectraBothHistoManager* hman_data, AliSpectraBothHistoManager* hman_mc,TFormula* fun,TFile *fout,TH1F** hcon,TH1F** hconWD,TH1F** hconMat,TH1F** hprimary)
 573 {
 574   printf("\n\n-> DCA Correction");
 575
 576
 577   Printf("\DCACorr");
 578   TString sample[2]={"data","mc"};
 579   ofstream debug("debugDCA.txt");
 580   TList* listofdcafits=new TList();
 581   for(Int_t icharge=0;icharge<2;icharge++)
 582   {
 583                 for(Int_t ipart=0;ipart<3;ipart++)
 584                 {
 585                         Int_t index=ipart+3*icharge;
 586                         for(Int_t isample=0;isample<2;isample++)
 587                         {
 588
 589                                 for(Int_t ibin_data=7;ibin_data<40;ibin_data++)
 590                                 {
 591
 592                                         Double_t lowedge=hcon[index]->GetBinLowEdge(ibin_data);
 593                                         Double_t binwidth=hcon[index]->GetBinWidth(ibin_data);
 594                                         debug<<"NEW "<<Particle[ipart].Data()<<" "<<Sign[icharge].Data()<<" "<<lowedge<<endl;
 595                                         if(fun)
 596                                         {
 597                                                 CutRange[1]=fun->Eval(lowedge);
 598                                                 CutRange[0]=-1.0*CutRange[1];
 599                                         }
 600                                         debug<<"cut  "<<CutRange[1]<<" "<<CutRange[0]<<endl;
 601                                         Short_t fitsettings=DCAfitsettings(lowedge+0.5*binwidth,index);
 602                                         debug<<"settings "<< fitsettings<<endl;
 603                                         if(fitsettings==0)
 604                                                 continue;
 605
 606                                         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);
 607                                         TLegend* Leg1 = new TLegend(0.6,0.6,0.85,0.85,"","NDC");
 608                                         Leg1->SetFillStyle(kFALSE);
 609                                         Leg1->SetLineColor(kWhite);
 610                                         Leg1->SetBorderSize(0);
 611
 612
 613                                         if(isample==0)
 614                                                 TH1F *hToFit =(TH1F*) ((TH1F*)hman_data->GetDCAHistogram1D(Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone();
 615                                         if(isample==1)
 616                                                 TH1F *hToFit =(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone();
 617                                         debug<<Particle[ipart].Data()<<" "<<Sign[icharge].Data()<<" "<<lowedge<<endl;
 618                                         TH1F *hmc1=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigmaPrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone();
 619                                         TH1F *hmc2=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigmaSecondaryWeakDecay%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone();
 620                                         TH1F *hmc3=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigmaSecondaryMaterial%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge))->Clone();
 621                                         Double_t minentries=1;
 622                                         debug<<hToFit->GetEntries()<<" "<<hmc1->GetEntries()<<" "<<hmc2->GetEntries()<<" "<<hmc3->GetEntries()<<endl;
 623                                         debug<<((fitsettings&0x1)&&hmc2->GetEntries()<=minentries)<<" "<<((fitsettings&0x2)&&hmc3->GetEntries()<=minentries)<<endl;
 624                                         if(hToFit->GetEntries()<=minentries || hmc1->GetEntries()<=minentries || ((fitsettings&0x1)&&hmc2->GetEntries()<=minentries) || ((fitsettings&0x2)&&hmc3->GetEntries()<=minentries))
 625                                                 continue;
 626                                         hToFit->Sumw2();
 627                                         hmc1->Sumw2();
 628                                         hmc2->Sumw2();
 629                                         hmc3->Sumw2();
 630
 631                                         Float_t corrforrebinning[4]={1.0,1.0,1.0,1.0};
 632                                         Int_t binCutRange[]={hmc1->GetXaxis()->FindBin(CutRange[0]),hmc1->GetXaxis()->FindBin(CutRange[1])};
 633
 634
 635                                         if(hmc3->GetNbinsX()>300)
 636                                         {
 637
 638                                                 corrforrebinning[0]=hToFit->Integral(binCutRange[0],binCutRange[1]);
 639                                                 corrforrebinning[1]=hmc1->Integral(binCutRange[0],binCutRange[1]);
 640                                                 corrforrebinning[2]=hmc2->Integral(binCutRange[0],binCutRange[1]);
 641                                                 corrforrebinning[3]=hmc3->Integral(binCutRange[0],binCutRange[1]);
 642
 643                                                 hToFit->Rebin(30);
 644                                                 hmc1->Rebin(30);
 645                                                 hmc2->Rebin(30);
 646                                                 hmc3->Rebin(30);
 647
 648                                                 binCutRange[0]=hmc1->GetXaxis()->FindBin(CutRange[0]);
 649                                                 binCutRange[1]=hmc1->GetXaxis()->FindBin(CutRange[1]);
 650
 651
 652                                                 //after rebbing we lose resolution of the dca this correction also us to do obtain inside used dca
 653
 654                                                 if(hToFit->Integral(binCutRange[0],binCutRange[1])>0.0)
 655                                                         corrforrebinning[0]=corrforrebinning[0]/hToFit->Integral(binCutRange[0],binCutRange[1]);
 656                                                 else
 657                                                         corrforrebinning[0]=1.0;
 658                                                 if(hmc1->Integral(binCutRange[0],binCutRange[1])>0.0)
 659                                                         corrforrebinning[1]=corrforrebinning[1]/hmc1->Integral(binCutRange[0],binCutRange[1]);
 660                                                 else
 661                                                         corrforrebinning[1]=1.0;
 662                                                 if(hmc2->Integral(binCutRange[0],binCutRange[1])>0.0)
 663                                                         corrforrebinning[2]=corrforrebinning[2]/hmc2->Integral(binCutRange[0],binCutRange[1]);
 664                                                 else
 665                                                         corrforrebinning[2]=1.0;
 666                                                 if(hmc3->Integral(binCutRange[0],binCutRange[1])>0.0)
 667                                                         corrforrebinning[3]=corrforrebinning[3]/hmc3->Integral(binCutRange[0],binCutRange[1]);
 668                                                 else
 669                                                         corrforrebinning[3]=1.0;
 670
 671
 672
 673                                         }
 674
 675                                         cDCA->cd();
 676                                         gPad->SetGridy();
 677                                         gPad->SetGridx();
 678                                         gPad->SetLogy();
 679
 680                                         TObjArray *mc=0x0;
 681                                         if(fitsettings==3)
 682                                                 mc = new TObjArray(3);        // MC histograms are put in this array
 683                                         else
 684                                                 mc = new TObjArray(2);
 685
 686
 687                                         mc->Add(hmc1);
 688                                         if(fitsettings&0x1)
 689                                                 mc->Add(hmc2);
 690                                         if(fitsettings&0x2)
 691                                                 mc->Add(hmc3);
 692                                         TFractionFitter* fit = new TFractionFitter(hToFit,mc); // initialise
 693                                         fit->Constrain(0,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
 694                                         if(fitsettings&0x1)
 695                                                 fit->Constrain(1,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
 696                                         if(fitsettings&0x2)
 697                                                 fit->Constrain(1+(fitsettings&0x1),0.0,1.0);               // constrain fraction 1 to be between 0 and 1
 698
 699                                         Int_t binFitRange[]={hmc1->GetXaxis()->FindBin(FitRange[0]),hmc1->GetXaxis()->FindBin(FitRange[1])};
 700                                         fit->SetRangeX(binFitRange[0],binFitRange[1]);
 701                                         hToFit->GetXaxis()->SetRange(binFitRange[0],binFitRange[1]);
 702                                         hToFit->SetTitle(Form("DCA distr - %s %s %s %lf",sample[isample].Data(),Particle[ipart].Data(),Sign[icharge].Data(),lowedge));
 703                                         Int_t status = fit->Fit();               // perform the fit
 704                                         cout << "fit status: " << status << endl;
 705                                         debug<<"fit status: " << status << endl;
 706
 707                                         if (status == 0)
 708                                         {
 709                                                 Double_t v1=0.0,v2=0.0,v3=0.0;
 710                                                 Double_t ev1=0.0,ev2=0.0,ev3=0.0;
 711                                                 Double_t cov=0.0;
 712
 713                                                 // check on fit status
 714                                                 TH1F* result = (TH1F*) fit->GetPlot();
 715                                                 TH1F* PrimMCPred=(TH1F*)fit->GetMCPrediction(0);
 716
 717                                                 TH1F* secStMCPred=0X0;
 718
 719                                                 TH1F* secMCPred=0x0;
 720
 721                                                 //first method, use directly the fit result
 722                                                 fit->GetResult(0,v1,ev1);
 723
 724                                                 if(fitsettings&0x1)
 725                                                 {
 726                                                         fit->GetResult(1,v2,ev2);
 727                                                         secStMCPred=(TH1F*)fit->GetMCPrediction(1);
 728
 729                                                 }
 730                                                 if(fitsettings&0x2)
 731                                                 {
 732                                                         fit->GetResult(1+(fitsettings&0x1),v3,ev3);
 733                                                         secMCPred=(TH1F*)fit->GetMCPrediction(1+(fitsettings&0x1));
 734                                                         if(fitsettings&0x1)
 735                                                                 cov=fit->GetFitter()->GetCovarianceMatrixElement(1,2);
 736
 737
 738                                                 }
 739
 740
 741                                                 debug<<v1<<" "<<ev1<<" "<<v2<<" "<<ev2<<" "<<v3<<" "<<ev3<<" "<<" "<<cov<<endl;
 742
 743                                                 // becuase dca cut range is not a fit range the results from TFractionFitter should be rescale
 744                                                 // This can be done in two ways or use input histograms or output histograms
 745                                                 // The difference between those two methods should be on the level of statistical error
 746                                                 // I use output histograms
 747
 748                                                 // Method 1 input histo
 749
 750                                                 Float_t normalizationdata=hToFit->Integral(hToFit->GetXaxis()->FindBin(CutRange[0]),hToFit->GetXaxis()->FindBin(CutRange[1]))/hToFit->Integral(binFitRange[0],binFitRange[1]);
 751                                                 normalizationdata*=corrforrebinning[0];
 752
 753                                                 Float_t normalizationmc1=(hmc1->Integral(binCutRange[0],binCutRange[1])/hmc1->Integral(binFitRange[0],binFitRange[1]))/normalizationdata;
 754                                                 Float_t normalizationmc2=0.0;
 755                                                 if(fitsettings&0x1)
 756                                                         normalizationmc2=(hmc2->Integral(binCutRange[0],binCutRange[1])/hmc2->Integral(binFitRange[0],binFitRange[1]))/normalizationdata;
 757                                                 Float_t normalizationmc3=0.0;
 758                                                 if(fitsettings&0x2)
 759                                                         normalizationmc3=(hmc3->Integral(binCutRange[0],binCutRange[1])/hmc3->Integral(binFitRange[0],binFitRange[1]))/normalizationdata;
 760
 761                                                 normalizationmc1*=corrforrebinning[1];
 762                                                 normalizationmc2*=corrforrebinning[2];
 763                                                 normalizationmc3*=corrforrebinning[3];
 764
 765                                                 debug<<"After Nor"<<endl;
 766                                                 debug<<v1*normalizationmc1<<" "<<ev1*normalizationmc1<<" "<<v2*normalizationmc2<<" "<<ev2*normalizationmc2<<" "<<v3*normalizationmc3<<" "<<ev3*normalizationmc3<<" "<<endl;
 767                                                 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;
 768                                                 debug<<"addtional info"<<endl;
 769
 770
 771                                                //Method 2 output histo
 772
 773                                                 Float_t normalizationdata1=result->Integral(binCutRange[0],binCutRange[1])/result->Integral(binFitRange[0],binFitRange[1]);
 774
 775
 776                                                 // if the cut range is bigger the fit range we should calculate the normalization factor for data using the data histogram
 777                                                 // because result histogram has entries only in fits range
 778                                                 if(FitRange[0]>CutRange[0]||FitRange[1]<CutRange[1])
 779                                                         normalizationdata1=normalizationdata;
 780
 781                                                 normalizationdata1*=corrforrebinning[0];
 782
 783
 784                                                 Float_t normalizationmc11=(PrimMCPred->Integral(binCutRange[0],binCutRange[1])/PrimMCPred->Integral(binFitRange[0],binFitRange[1]))/normalizationdata1;
 785                                                 Float_t normalizationmc21=0.0;
 786                                                 if(fitsettings&0x1)
 787                                                         normalizationmc21=(secStMCPred->Integral(binCutRange[0],binCutRange[1])/secStMCPred->Integral(binFitRange[0],binFitRange[1]))/normalizationdata1;
 788                                                 Float_t normalizationmc31=0.0;
 789                                                 if(fitsettings&0x2)
 790                                                         normalizationmc31=(secMCPred->Integral(binCutRange[0],binCutRange[1])/secMCPred->Integral(binFitRange[0],binFitRange[1]))/normalizationdata1;
 791
 792                                                 normalizationmc11*=corrforrebinning[1];
 793                                                 normalizationmc21*=corrforrebinning[2];
 794                                                 normalizationmc31*=corrforrebinning[3];
 795
 796                                                 debug<<"After Nor 2"<<endl;
 797                                                 debug<<v1*normalizationmc11<<" "<<ev1*normalizationmc11<<" "<<v2*normalizationmc21<<" "<<ev2*normalizationmc21<<" "<<v3*normalizationmc31<<" "<<ev3*normalizationmc31<<endl;
 798
 799                                                 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;
 800
 801
 802                                                 hconWD[index+6*isample]->SetBinContent(ibin_data,v2*normalizationmc21);
 803                                                 hconWD[index+6*isample]->SetBinError(ibin_data,ev2*normalizationmc21);
 804                                                 hconMat[index+6*isample]->SetBinContent(ibin_data,v3*normalizationmc31);
 805                                                 hconMat[index+6*isample]->SetBinError(ibin_data,ev3*normalizationmc31);
 806                                                 hprimary[index+6*isample]->SetBinContent(ibin_data,v1*normalizationmc11);
 807                                                 hprimary[index+6*isample]->SetBinError(ibin_data,ev1*normalizationmc11);
 808                                                 hcon[index+6*isample]->SetBinContent(ibin_data,v2*normalizationmc21+v3*normalizationmc31);
 809                                                 hcon[index+6*isample]->SetBinError(ibin_data,TMath::Sqrt(ev2*ev2*normalizationmc21*normalizationmc21+ev3*ev3*normalizationmc31*normalizationmc31+cov*normalizationmc31*normalizationmc21));
 810
 811
 812
 813                                                 //Drawing section
 814                                                 result->Scale(1.0/result->Integral(result->GetXaxis()->FindBin(FitRange[0]),result->GetXaxis()->FindBin(FitRange[1])));
 815                                                 hToFit->Scale(1.0/hToFit->Integral(hToFit->GetXaxis()->FindBin(FitRange[0]),hToFit->GetXaxis()->FindBin(FitRange[1])));
 816                                                 PrimMCPred->Scale(v1/PrimMCPred->Integral(PrimMCPred->GetXaxis()->FindBin(FitRange[0]),PrimMCPred->GetXaxis()->FindBin(FitRange[1])));
 817
 818                                                 hToFit->SetMinimum(0.0001);
 819                                                 hToFit->DrawClone("E1x0");
 820                                                 result->SetTitle("Fit result");
 821                                                 result->SetLineColor(kBlack);
 822                                                 Leg1->AddEntry(result,"result","lp");
 823                                                 result->DrawClone("lhistsame");
 824
 825                                                 PrimMCPred->SetLineColor(kGreen+2);
 826                                                 PrimMCPred->SetLineStyle(2);
 827                                                  PrimMCPred->SetLineWidth(3.0);
 828                                                 Leg1->AddEntry(PrimMCPred,"Prmi.","l");
 829                                                 PrimMCPred->DrawClone("lhistsame");
 830                                                 if(fitsettings&0x1)
 831                                                 {
 832
 833                                                         secStMCPred->Scale(v2/secStMCPred->Integral(secStMCPred->GetXaxis()->FindBin(FitRange[0]),secStMCPred->GetXaxis()->FindBin(FitRange[1])));
 834                                                         secStMCPred->SetLineColor(kRed);
 835                                                         secStMCPred->SetLineWidth(3.0);
 836
 837                                                         secStMCPred->SetLineStyle(3);
 838                                                         Leg1->AddEntry(secStMCPred,"Sec.WD","l");
 839                                                         secStMCPred->DrawClone("lhistsame");
 840
 841                                                 }
 842                                                 if(fitsettings&0x2)
 843                                                 {
 844
 845                                                         secMCPred->Scale(v3/secMCPred->Integral(secMCPred->GetXaxis()->FindBin(FitRange[0]),secMCPred->GetXaxis()->FindBin(FitRange[1])));
 846                                                         secMCPred->SetLineColor(kBlue);
 847                                                         secMCPred->SetLineWidth(3.0);
 848
 849                                                         secMCPred->SetLineStyle(4);
 850                                                         Leg1->AddEntry(secMCPred,"Sec.Mat","l");
 851                                                         secMCPred->DrawClone("lhistsame");
 852
 853                                                 }
 854                                         }
 855                                         else
 856                                         {
 857                                                 hconWD[index+6*isample]->SetBinContent(ibin_data,0.0);
 858                                                 hconWD[index+6*isample]->SetBinError(ibin_data,0.0);
 859                                                 hconMat[index+6*isample]->SetBinContent(ibin_data,0.0);
 860                                                 hconMat[index+6*isample]->SetBinError(ibin_data,0.0);
 861                                                 hcon[index+6*isample]->SetBinContent(ibin_data,0.0);
 862                                                 hcon[index+6*isample]->SetBinError(ibin_data,0.0);
 863                                                 hprimary[index+6*isample]->SetBinContent(ibin_data,1.0);
 864                                                 hprimary[index+6*isample]->SetBinError(ibin_data,0.0);
 865                                         }
 866                                         Leg1->Draw();
 867                                         listofdcafits->Add(cDCA);
 868
 869                                         //cDCA->Write();
 870                                         delete hToFit;
 871                                 }
 872
 873
 874                         }
 875
 876                 }
 877         }
 878         fout->cd();
 879         listofdcafits->Write("DCAfits",TObject::kSingleKey);
 880 }
 881
 882 void RecomputeErrors(TH1* h)
 883 {
 884         for (int i=0; i<=h->GetXaxis()->GetNbins(); i++)
 885         {
 886                 cout<<h->GetBinContent(i)<<" "<<h->GetBinError(i)<<" error "<<TMath::Sqrt(h->GetBinContent(i))<<endl;
 887                 h->SetBinError(i,TMath::Sqrt(h->GetBinContent(i)));
 888         }
 889         h->Sumw2();
 890 }
 891 void SetBintoOne(TH1* h)
 892 {
 893         for (int i=0;i<=h->GetXaxis()->GetNbins(); i++)
 894         {
 895                 h->SetBinContent(i,1);
 896                 h->SetBinError(i,0);
 897         }
 898 }
 899
 900
 901 void GetCorrectedSpectra(TH1F* corr,TH1F* raw,TH1F* eff, TH1F* con)
 902 {
 903         for (int i=0;i<=corr->GetXaxis()->GetNbins(); i++)
 904         {
 905                 corr->SetBinContent(i,1);
 906                 corr->SetBinError(i,0);
 907         }
 908         corr->Sumw2();
 909         corr->Add(con,-1);
 910         corr->Sumw2();
 911         corr->Divide(eff);
 912         corr->Sumw2();
 913         corr->Multiply(raw);
 914         corr->Sumw2();
 915 }
 916 void GetCorrectedSpectraLeonardo(TH1F* spectra,TH1F* correction, TH1F* hprimaryData,TH1F* hprimaryMC)
 917 {
 918         spectra->Sumw2();
 919         spectra->Multiply(correction);
 920         spectra->Sumw2();
 921         hprimaryData->Sumw2();
 922         spectra->Multiply(hprimaryData);
 923         hprimaryMC->Sumw2();
 924         spectra->Divide(hprimaryMC);
 925 }
 926
 927 void GFCorrection(TH1F **Spectra,Float_t tofpt,UInt_t options)
 928 {
 929         if (options&kgeantflukaKaon)
 930         {
 931                  //Geant/Fluka Correction
 932                 Printf("\nGF correction for Kaons");
 933                 //Getting GF For Kaons in TPC
 934                 TGraph *gGFCorrectionKaonPlus=new TGraph();
 935                 gGFCorrectionKaonPlus->SetName("gGFCorrectionKaonPlus");
 936                 gGFCorrectionKaonPlus->SetTitle("gGFCorrectionKaonPlus");
 937                 TGraph *gGFCorrectionKaonMinus=new TGraph();
 938                 gGFCorrectionKaonMinus->SetName("gGFCorrectionKaonMinus");
 939                 gGFCorrectionKaonMinus->SetTitle("gGFCorrectionKaonMinus");
 940                  TString fnameGeanFlukaK="GFCorrection/correctionForCrossSection.321.root";
 941                 TFile *fGeanFlukaK= new TFile(fnameGeanFlukaK.Data());
 942                 if (!fGeanFlukaK)
 943                 {
 944                         fnameGeanFlukaK="$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/TestAOD/GFCorrection/correctionForCrossSection.321.root";
 945                         fGeanFlukaK= new TFile(fnameGeanFlukaK.Data());
 946                         if (!fGeanFlukaK)
 947                                 return;
 948                 }
 949                 TH1F *hGeantFlukaKPos=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionParticles");
 950                 TH1F *hGeantFlukaKNeg=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionAntiParticles");
 951                 //getting GF func for Kaons with TOF
 952                 TF1 *fGFKPosTracking;
 953                 fGFKPosTracking = TrackingEff_geantflukaCorrection(3,kPositive);
 954                 TF1 *fGFKNegTracking;
 955                  fGFKNegTracking = TrackingEff_geantflukaCorrection(3,kNegative);
 956                  TF1 *fGFKPosMatching;
 957                  fGFKPosMatching = TOFmatchMC_geantflukaCorrection(3,kPositive);
 958                  TF1 *fGFKNegMatching;
 959                  fGFKNegMatching = TOFmatchMC_geantflukaCorrection(3,kNegative);
 960                  for(Int_t binK=0;binK<=Spectra[1]->GetNbinsX();binK++)
 961                 {
 962                         if(Spectra[1]->GetBinCenter(binK)<tofpt)
 963                         {//use TPC GeantFlukaCorrection
 964                                 Float_t FlukaCorrKPos=hGeantFlukaKPos->GetBinContent(hGeantFlukaKPos->FindBin(Spectra[1]->GetBinCenter(binK)));
 965                                 Float_t FlukaCorrKNeg=hGeantFlukaKNeg->GetBinContent(hGeantFlukaKNeg->FindBin(Spectra[4]->GetBinCenter(binK)));
 966                         //      Printf("TPC Geant/Fluka: pt:%f  Pos:%f  Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPos,FlukaCorrKNeg);
 967                                 Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPos);
 968                                 Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNeg);
 969                                 Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPos);
 970                                 Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNeg);
 971                                 gGFCorrectionKaonPlus->SetPoint(binK,Spectra[1]->GetBinCenter(binK),FlukaCorrKPos);
 972                                 gGFCorrectionKaonMinus->SetPoint(binK,Spectra[4]->GetBinCenter(binK),FlukaCorrKNeg);
 973                         }
 974                         else
 975                         {
 976                                 gGFCorrectionKaonPlus->SetPoint(binK,Spectra[1]->GetBinCenter(binK),0);
 977                                 gGFCorrectionKaonMinus->SetPoint(binK,Spectra[4]->GetBinCenter(binK),0);
 978                                 Float_t FlukaCorrKPosTracking=fGFKPosTracking->Eval(Spectra[1]->GetBinCenter(binK));
 979                                 Float_t FlukaCorrKNegTracking=fGFKNegTracking->Eval(Spectra[1]->GetBinCenter(binK));
 980                         //      Printf("TPC/TOF Geant/Fluka Tracking: pt:%f  Pos:%f  Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPosTracking,FlukaCorrKNegTracking);
 981                                 Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPosTracking);
 982                                 Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNegTracking);
 983                                 Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPosTracking);
 984                                 Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNegTracking);
 985                                 Float_t FlukaCorrKPosMatching=fGFKPosMatching->Eval(Spectra[1]->GetBinCenter(binK));
 986                                 Float_t FlukaCorrKNegMatching=fGFKNegMatching->Eval(Spectra[1]->GetBinCenter(binK));
 987                         //      Printf("TPC/TOF Geant/Fluka Matching: pt:%f  Pos:%f  Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPosMatching,FlukaCorrKNegMatching);
 988                                 Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPosMatching);
 989                                 Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNegMatching);
 990                                 Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPosMatching);
 991                                 Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNegMatching);
 992                         }
 993                 }
 994           }
 995           if(!(options&kgeantflukaProton))
 996                 return;
 997           //Geant Fluka for P in TPC
 998           Printf("\nGF correction for Protons");
 999           const Int_t kNCharge=2;
1000           Int_t kPos=0;
1001           Int_t kNeg=1;
1002           TString fnameGFProtons= "GFCorrection/correctionForCrossSection.root";
1003           TFile* fGFProtons = new TFile (fnameGFProtons.Data());
1004           if (!fGFProtons)
1005           {
1006                 fnameGFProtons=="$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/TestAOD/GFCorrection/correctionForCrossSection.root";
1007                 TFile* fGFProtons = new TFile (fnameGFProtons.Data());
1008                 if (!fGFProtons)
1009                         return;
1010           }
1011
1012
1013
1014           TH2D * hCorrFluka[kNCharge];
1015           TH2D * hCorrFluka[2];
1016           hCorrFluka[kPos] = (TH2D*)fGFProtons->Get("gHistCorrectionForCrossSectionProtons");
1017           hCorrFluka[kNeg] = (TH2D*)fGFProtons->Get("gHistCorrectionForCrossSectionAntiProtons");
1018           //getting GF func for Kaons with TPCTOF
1019           TF1 *fGFpPosTracking;
1020           fGFpPosTracking = TrackingEff_geantflukaCorrection(4,kPositive);
1021           TF1 *fGFpNegTracking;
1022           fGFpNegTracking = TrackingEff_geantflukaCorrection(4,kNegative);
1023           TF1 *fGFpPosMatching;
1024           fGFpPosMatching = TOFmatchMC_geantflukaCorrection(4,kPositive);
1025           TF1 *fGFpNegMatching;
1026           fGFpNegMatching = TOFmatchMC_geantflukaCorrection(4,kNegative);
1027
1028
1029                 Int_t nbins = Spectra[2]->GetNbinsX();
1030
1031                 for(Int_t ibin = 0; ibin < nbins; ibin++)
1032                 {
1033                         if(Spectra[2]->GetBinCenter(ibin)<tofpt)
1034                         {//use TPC GeantFlukaCorrection
1035                                 for(Int_t icharge = 0; icharge < kNCharge; icharge++)
1036                                 {
1037                                         Int_t nbinsy=hCorrFluka[icharge]->GetNbinsY();
1038                                         Float_t pt = Spectra[2]->GetBinCenter(ibin);
1039                                         Float_t minPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(1);
1040                                         Float_t maxPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(nbinsy+1);
1041                                         if (pt < minPtCorrection)
1042                                                 pt = minPtCorrection+0.0001;
1043                                         if (pt > maxPtCorrection)
1044                                                 pt = maxPtCorrection;
1045                                         Float_t correction = hCorrFluka[icharge]->GetBinContent(1,hCorrFluka[icharge]->GetYaxis()->FindBin(pt));
1046
1047                                         if (correction > 0.0)
1048                                         {// If the bin is empty this is a  0
1049                         //                      cout<<icharge<<" "<<ibin<<" "<<correction<<endl;
1050                                                 Spectra[icharge*3+2]->SetBinContent(ibin,Spectra[icharge*3+2]->GetBinContent(ibin)*correction);
1051                                                 Spectra[icharge*3+2]->SetBinError(ibin,Spectra[icharge*3+2]->GetBinError(ibin)*correction);
1052                                         }
1053                                         else if (Spectra[icharge*3+2]->GetBinContent(ibin) > 0.0)
1054                                         {
1055                                                 // If we are skipping a non-empty bin, we notify the user
1056                                                 cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for " <<"protons Pt:"<< pt<< endl;
1057                                                 cout << " Bin content: " << Spectra[icharge*3+2]->GetBinContent(ibin)  << endl;
1058                                         }
1059                                 }
1060                         }
1061                         else
1062                         {
1063                                 Float_t FlukaCorrpPosTracking=fGFpPosTracking->Eval(Spectra[2]->GetBinCenter(ibin));
1064                                 Float_t FlukaCorrpNegTracking=fGFpNegTracking->Eval(Spectra[2]->GetBinCenter(ibin));
1065                         //      Printf("TPC/TOF Geant/Fluka Tracking: pt:%f  Pos:%f  Neg:%f",Spectra[2]->GetBinCenter(ibin),FlukaCorrpPosTracking,FlukaCorrpNegTracking);
1066                                 Spectra[2]->SetBinContent(ibin,Spectra[2]->GetBinContent(ibin)*FlukaCorrpPosTracking);
1067                                 Spectra[5]->SetBinContent(ibin,Spectra[5]->GetBinContent(ibin)*FlukaCorrpNegTracking);
1068                                 Spectra[2]->SetBinError(ibin,Spectra[2]->GetBinError(ibin)*FlukaCorrpPosTracking);
1069                                 Spectra[5]->SetBinError(ibin,Spectra[5]->GetBinError(ibin)*FlukaCorrpNegTracking);
1070                                 Float_t FlukaCorrpPosMatching=fGFpPosMatching->Eval(Spectra[2]->GetBinCenter(ibin));
1071                                 Float_t FlukaCorrpNegMatching=fGFpNegMatching->Eval(Spectra[2]->GetBinCenter(ibin));
1072                         //      Printf("TPC/TOF Geant/Fluka Matching: pt:%f  Pos:%f  Neg:%f",Spectra[2]->GetBinCenter(ibin),FlukaCorrpPosMatching,FlukaCorrpNegMatching);
1073                                 Spectra[2]->SetBinContent(ibin,Spectra[2]->GetBinContent(ibin)*FlukaCorrpPosMatching);
1074                                 Spectra[5]->SetBinContent(ibin,Spectra[5]->GetBinContent(ibin)*FlukaCorrpNegMatching);
1075                                 Spectra[2]->SetBinError(ibin,Spectra[2]->GetBinError(ibin)*FlukaCorrpPosMatching);
1076                                 Spectra[5]->SetBinError(ibin,Spectra[5]->GetBinError(ibin)*FlukaCorrpNegMatching);
1077                         }
1078                 }
1079          fGeanFlukaK->Close();
1080          delete fGeanFlukaK;
1081 }
1082
1083
1084 ///////////
1085 TF1 *
1086 TrackingEff_geantflukaCorrection(Int_t ipart, Int_t icharge)
1087 {
1088
1089   if (ipart == 3 && icharge == kNegative) {
1090     TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionKaMinus(x)", 0., 5.);
1091     return f;
1092   }
1093   else if (ipart == 4 && icharge == kNegative) {
1094     TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionPrMinus(x)", 0., 5.);
1095   }
1096   else
1097     TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionNull(x)", 0., 5.);
1098
1099   return f;
1100 }
1101
1102 Double_t
1103 TrackingPtGeantFlukaCorrectionNull(Double_t pTmc)
1104 {
1105   return 1.;
1106 }
1107
1108 Double_t
1109 TrackingPtGeantFlukaCorrectionPrMinus(Double_t pTmc)
1110 {
1111   return (1 - 0.129758 *TMath::Exp(-pTmc*0.679612));
1112 }
1113
1114 Double_t
1115 TrackingPtGeantFlukaCorrectionKaMinus(Double_t pTmc)
1116 {
1117   return TMath::Min((0.972865 + 0.0117093*pTmc), 1.);
1118 }
1119 ///////////////////////////////////////////
1120 TF1 *
1121 TOFmatchMC_geantflukaCorrection(Int_t ipart, Int_t icharge)
1122 {
1123
1124   if (ipart == 3 && icharge == kNegative) {
1125     TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionKaMinus(x)", 0., 5.);
1126     return f;
1127   }
1128   else if (ipart == 4 && icharge == kNegative) {
1129     TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionPrMinus(x)", 0., 5.);
1130   }
1131   else
1132     TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionNull(x)", 0., 5.);
1133
1134   return f;
1135 }
1136
1137
1138 Double_t
1139 MatchingPtGeantFlukaCorrectionNull(Double_t pTmc)
1140 {
1141   return 1.;
1142 }
1143
1144 Double_t
1145 MatchingPtGeantFlukaCorrectionPrMinus(Double_t pTmc)
1146 {
1147   Float_t ptTPCoutP =pTmc*(1-6.81059e-01*TMath::Exp(-pTmc*4.20094));
1148   return (TMath::Power(1 - 0.129758*TMath::Exp(-ptTPCoutP*0.679612),0.07162/0.03471));
1149 }
1150
1151 Double_t
1152 MatchingPtGeantFlukaCorrectionKaMinus(Double_t pTmc)
1153 {
1154   Float_t ptTPCoutK=pTmc*(1- 3.37297e-03/pTmc/pTmc - 3.26544e-03/pTmc);
1155   return TMath::Min((TMath::Power(0.972865 + 0.0117093*ptTPCoutK,0.07162/0.03471)), 1.);
1156 }
1157 void MatchingTOFEff(TH1F** Spectra, TList* list=0x0)
1158 {
1159           if(TOFMatchingScalling[0]<0.0&&TOFMatchingScalling[1]<0.0)
1160           {
1161                 TH1F *hMatcEffPos_data=(TH1F*)tcutsdata->GetHistoNMatchedPos();
1162                   hMatcEffPos_data->Sumw2();
1163                   //hMatcEffPos_data->Divide((TH1F*)tcutsdata->GetHistoNSelectedPos());
1164                   hMatcEffPos_data->Divide(hMatcEffPos_data,(TH1F*)tcutsdata->GetHistoNSelectedPos(),1,1,"B");
1165                   hMatcEffPos_data->SetTitle("Matching Eff Pos - data");
1166                   TH1F *hMatcEffNeg_data=(TH1F*)tcutsdata->GetHistoNMatchedNeg();
1167                   hMatcEffNeg_data->Sumw2();
1168                   //hMatcEffNeg_data->Divide((TH1F*)tcutsdata->GetHistoNSelectedNeg());
1169                   hMatcEffNeg_data->Divide(hMatcEffNeg_data,(TH1F*)tcutsdata->GetHistoNSelectedNeg(),1,1,"B");
1170                   hMatcEffNeg_data->SetTitle("Matching Eff Neg - data");
1171                   TH1F *hMatcEffPos_mc=(TH1F*)tcutsmc->GetHistoNMatchedPos();
1172                   hMatcEffPos_mc->Sumw2();
1173                   //hMatcEffPos_mc->Divide((TH1F*)tcutsmc->GetHistoNSelectedPos());
1174                   hMatcEffPos_mc->Divide(hMatcEffPos_mc,(TH1F*)tcutsmc->GetHistoNSelectedPos(),1,1,"B");
1175                   hMatcEffPos_mc->SetTitle("Matching Eff Pos - mc");
1176                   TH1F *hMatcEffNeg_mc=(TH1F*)tcutsmc->GetHistoNMatchedNeg();
1177                   hMatcEffNeg_mc->Sumw2();
1178                   //hMatcEffNeg_mc->Divide((TH1F*)tcutsmc->GetHistoNSelectedNeg());
1179                   hMatcEffNeg_mc->Divide(hMatcEffNeg_mc,(TH1F*)tcutsmc->GetHistoNSelectedNeg(),1,1,"B");
1180                   hMatcEffNeg_mc->SetTitle("Matching Eff Neg - mc");
1181
1182
1183                   hMatcEffPos_data->Divide(hMatcEffPos_mc);
1184                   hMatcEffNeg_data->Divide(hMatcEffNeg_mc);
1185                   hMatcEffPos_data->SetName("MatchingTOFPos");
1186                   hMatcEffNeg_data->SetName("MatchingTOFNeg");
1187
1188
1189                   TF1 *pol0MatchPos_data=new TF1("pol0MatchPos_data","pol0",.6,5);
1190                   hMatcEffPos_data->Fit("pol0MatchPos_data","MNR");
1191                   TF1 *pol0MatchNeg_data=new TF1("pol0MatchNeg_data","pol0",.6,5);
1192                   hMatcEffNeg_data->Fit("pol0MatchNeg_data","MNR");
1193
1194                 list->Add(hMatcEffPos_data);
1195                   list->Add(hMatcEffNeg_data);
1196
1197
1198                   TOFMatchingScalling[0]=pol0MatchPos_data->GetParameter(0);
1199                   TOFMatchingScalling[1]=pol0MatchNeg_data->GetParameter(0);
1200
1201           }
1202           //Correction spectra for matching efficiency
1203           //For the moment I'm using the inclusive correction
1204           for(Int_t ipart=0;ipart<3;ipart++)
1205           {
1206
1207                 for(Int_t ibin=1;ibin<Spectra[ipart]->GetNbinsX();ibin++)
1208                 {
1209                         Float_t ptspectra=Spectra[ipart]->GetBinCenter(ibin);
1210                         if(ptspectra<tcutsdata->GetPtTOFMatching())
1211                                 continue;
1212                   //Spectra[ipart]->SetBinContent(ibin,( Spectra[ipart]->GetBinContent(ibin)/hMatcEffPos_data->GetBinContent(hMatcEffPos_data->FindBin(ptspectra))));
1213                   //Spectra[ipart+3]->SetBinContent(ibin,( Spectra[ipart+3]->GetBinContent(ibin)/hMatcEffNeg_data->GetBinContent(hMatcEffNeg_data->FindBin(ptspectra))));
1214                         Spectra[ipart]->SetBinContent(ibin,( Spectra[ipart]->GetBinContent(ibin)/TOFMatchingScalling[0]));
1215                         Spectra[ipart+3]->SetBinContent(ibin,( Spectra[ipart+3]->GetBinContent(ibin)/TOFMatchingScalling[1]));
1216                 }
1217           }
1218 }
1219 Double_t eta2y(Double_t pt, Double_t mass, Double_t eta)
1220 {
1221   Double_t mt = TMath::Sqrt(mass * mass + pt * pt);
1222   return TMath::ASinH(pt / mt * TMath::SinH(eta));
1223 }
1224
1225 TH1* GetSumAllCh(TH1F** spectra, Double_t* mass  )
1226 {
1227         TH1F* allch=(((TH1F*))spectra[0]->Clone("allCh"));
1228         allch->Reset();
1229         for (int i=0;i<6;i++)
1230         {
1231                 Double_t masstmp=mass[i%3];
1232                 for (int j=1;j<spectra[i]->GetXaxis()->GetNbins();j++)
1233                 {
1234                         Float_t value=spectra[i]->GetBinContent(j);
1235                         Float_t error=spectra[i]->GetBinError(j);
1236                         if(value>0.0)
1237                         {
1238                                 Float_t pt=spectra[i]->GetXaxis()->GetBinCenter(j);
1239                                 Float_t mt2=masstmp*masstmp+pt*pt;
1240                                 Float_t mt=TMath::Sqrt(mt2);
1241                                 Float_t maxy=eta2y(pt,masstmp,0.8);
1242                                 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))));
1243                                 conver=conver/1.6;
1244                                 //cout<<maxy<<" "<<conver<<" "<<masstmp<<""<<spectra[i]->GetName()<<endl;
1245                                 Float_t bincontent=allch->GetBinContent(j);
1246                                 Float_t binerror=allch->GetBinError(j);
1247                                 bincontent+=conver*value;
1248                                 binerror=TMath::Sqrt(binerror*binerror+conver*conver*error*error);
1249                                 allch->SetBinContent(j,bincontent);
1250                                 allch->SetBinError(j,binerror);
1251                         }
1252
1253                 }
1254         }
1255         Divideby2pipt(allch);
1256         allch->SetTitle("N_{ch};p_{T} (GeV/c);1/(2 #pi p_{T})dN/p_{T} |#eta|<0.8");
1257         return allch;
1258 }
1259
1260 void Divideby2pipt(TH1* hist)
1261 {
1262
1263         for (int i=1;i<hist->GetXaxis()->GetNbins();i++)
1264         {
1265                 Float_t value=hist->GetBinContent(i);
1266                 Float_t error=hist->GetBinError(i);
1267                 Float_t pt=hist->GetXaxis()->GetBinCenter(i);
1268                 hist->SetBinContent(i,value/(2.0*TMath::Pi()*pt));
1269                 hist->SetBinError(i,error/(2.0*TMath::Pi()*pt));
1270
1271         }
1272 }
1273
1274 Short_t DCAfitsettings (Float_t pt, Int_t type)
1275 {
1276         Short_t value=0x0;
1277         if (pt<maxptformaterial[type]&&pt>minptformaterial[type])
1278                 value=value+2;
1279         if (pt<maxptforWD[type]&&pt>minptforWD[type])
1280                 value=value+1;
1281         return value;
1282
1283 }
1284
1285 Float_t Normaliztionwithbin0integrals(UInt_t options)
1286 {
1287
1288         TH1F* bin0mcRec=(TH1F*)ecutsmc->GetHistoVtxGenerated()->Clone("Bin0_rec");
1289         TH1F* bin0mcMC=(TH1F*)ecutsmc->GetHistoVtxGenerated()->Clone("Bin0_MC");
1290
1291         TH1F* vertexmc=ecutsmc->GetHistoVtxAftSelwithoutZvertexCut();
1292         TH1F* vertexmcMCz=ecutsmc->GetHistoVtxAftSelwithoutZvertexCutusingMCz();
1293         TH1F* vertexdata=ecutsdata->GetHistoVtxAftSelwithoutZvertexCut();
1294
1295         TH1I* histodata=ecutsdata->GetHistoCuts();
1296         TH1I* histomc=ecutsmc->GetHistoCuts();
1297
1298         Float_t dataevents=(Float_t)histodata->GetBinContent(3);
1299         //cout<<histodata->GetBinContent(2)<<endl;
1300         Float_t databin0events=((Float_t)histodata->GetBinContent(2))-((Float_t)histodata->GetBinContent(4));
1301
1302         bin0mcRec->Sumw2();
1303         bin0mcMC->Sumw2();
1304
1305         bin0mcRec->Add(vertexmc,-1);
1306         bin0mcMC->Add(vertexmcMCz,-1);
1307
1308         bin0mcRec->Divide(vertexmc);
1309         bin0mcMC->Divide(vertexmcMCz);
1310
1311         bin0mcRec->Multiply(vertexdata);
1312         bin0mcMC->Multiply(vertexdata);
1313
1314         Float_t bin0mcRecN=0.0;
1315         Float_t bin0mcMCN=0.0;
1316
1317         for (int i=0;i<=bin0mcRec->GetXaxis()->GetNbins();i++)
1318         {
1319                 bin0mcRecN+=bin0mcRec->GetBinContent(i);
1320                 bin0mcMCN+=bin0mcMC->GetBinContent(i);
1321
1322         }
1323         bin0mcRec->Scale(databin0events/bin0mcRecN);
1324         bin0mcMC->Scale(databin0events/bin0mcMCN);
1325
1326         Int_t binmin=bin0mcRec->GetXaxis()->FindBin(-10);
1327         Int_t binmax=bin0mcRec->GetXaxis()->FindBin(10)-1;
1328         cout<<  bin0mcRec->GetXaxis()->GetBinLowEdge(binmin)<<" "<<bin0mcRec->GetXaxis()->GetBinUpEdge(binmax)<<endl;
1329         cout<<bin0mcRecN<<" "<<bin0mcMCN<<" "<<databin0events<<endl;
1330         cout<<dataevents<<" normalization "<<dataevents+bin0mcRec->Integral(binmin,binmax)<<" "<<dataevents+bin0mcMC->Integral(binmin,binmax)<<endl;
1331         cout<<histodata->GetBinContent(2)<<" "<<histodata->GetBinContent(4)<<endl;
1332         if ((options&knormalizationwithbin0integralsdata)==knormalizationwithbin0integralsdata)
1333                 return  dataevents+bin0mcRec->Integral(binmin,binmax);
1334         else if ((options&kknormalizationwithbin0integralsMC)==knormalizationwithbin0integralsdata)
1335                 return dataevents+bin0mcMC->Integral(binmin,binmax) ;
1336         else
1337                 return 1;
1338 }
1339
1340
1341 Bool_t ReadConfigFile(TString configfile)
1342 {
1343         ifstream infile(configfile.Data());
1344         if(infile.is_open()==false)
1345                 return false;
1346         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"};
1347
1348         char buffer[256];
1349         while (infile.eof()==false)
1350         {
1351                 buffer[0]='#';
1352                 while (buffer[0]=='#'&&infile.eof()==false)
1353                         infile.getline(buffer,256);
1354                 TString tmpstring(buffer);
1355                 cout<<buffer<<endl;
1356                 if(tmpstring.Contains(namesofSetting[0]))
1357                         CutRange[0]=(tmpstring.Remove(0,namesofSetting[0].Length()+1)).Atof();
1358                 else if (tmpstring.Contains(namesofSetting[1]))
1359                         CutRange[1]=(tmpstring.Remove(0,namesofSetting[1].Length()+1)).Atof();
1360                 else if (tmpstring.Contains(namesofSetting[2]))
1361                         FitRange[0]=(tmpstring.Remove(0,namesofSetting[2].Length()+1)).Atof();
1362                 else if (tmpstring.Contains(namesofSetting[3]))
1363                         FitRange[1]=(tmpstring.Remove(0,namesofSetting[3].Length()+1)).Atof();
1364                 else if (tmpstring.Contains(namesofSetting[4]))
1365                         minptformaterial[0]=(tmpstring.Remove(0,namesofSetting[4].Length()+1)).Atof();
1366                 else if (tmpstring.Contains(namesofSetting[5]))
1367                         maxptformaterial[0]=(tmpstring.Remove(0,namesofSetting[5].Length()+1)).Atof();
1368                 else if (tmpstring.Contains(namesofSetting[6]))
1369                         minptformaterial[1]=(tmpstring.Remove(0,namesofSetting[6].Length()+1)).Atof();
1370                 else if (tmpstring.Contains(namesofSetting[7]))
1371                         maxptformaterial[1]=(tmpstring.Remove(0,namesofSetting[7].Length()+1)).Atof();
1372                 else if (tmpstring.Contains(namesofSetting[8]))
1373                         minptformaterial[2]=(tmpstring.Remove(0,namesofSetting[8].Length()+1)).Atof();
1374                 else if (tmpstring.Contains(namesofSetting[9]))
1375                         maxptformaterial[2]=(tmpstring.Remove(0,namesofSetting[9].Length()+1)).Atof();
1376                 else if (tmpstring.Contains(namesofSetting[10]))
1377                         minptformaterial[3]=(tmpstring.Remove(0,namesofSetting[10].Length()+1)).Atof();
1378                 else if (tmpstring.Contains(namesofSetting[11]))
1379                         maxptformaterial[3]=(tmpstring.Remove(0,namesofSetting[11].Length()+1)).Atof();
1380                 else if (tmpstring.Contains(namesofSetting[12]))
1381                         minptformaterial[4]=(tmpstring.Remove(0,namesofSetting[12].Length()+1)).Atof();
1382                 else if (tmpstring.Contains(namesofSetting[13]))
1383                         maxptformaterial[4]=(tmpstring.Remove(0,namesofSetting[13].Length()+1)).Atof();
1384                 else if (tmpstring.Contains(namesofSetting[14]))
1385                         minptformaterial[5]=(tmpstring.Remove(0,namesofSetting[14].Length()+1)).Atof();
1386                 else if (tmpstring.Contains(namesofSetting[15]))
1387                         maxptformaterial[5]=(tmpstring.Remove(0,namesofSetting[15].Length()+1)).Atof();
1388                 else if (tmpstring.Contains(namesofSetting[16]))
1389                         minptforWD[0]=(tmpstring.Remove(0,namesofSetting[16].Length()+1)).Atof();
1390                 else if (tmpstring.Contains(namesofSetting[17]))
1391                         maxptforWD[0]=(tmpstring.Remove(0,namesofSetting[17].Length()+1)).Atof();
1392                 else if (tmpstring.Contains(namesofSetting[18]))
1393                         minptforWD[1]=(tmpstring.Remove(0,namesofSetting[18].Length()+1)).Atof();
1394                 else if (tmpstring.Contains(namesofSetting[19]))
1395                         maxptforWD[1]=(tmpstring.Remove(0,namesofSetting[19].Length()+1)).Atof();
1396                 else if (tmpstring.Contains(namesofSetting[20]))
1397                         minptforWD[2]=(tmpstring.Remove(0,namesofSetting[20].Length()+1)).Atof();
1398                 else if (tmpstring.Contains(namesofSetting[21]))
1399                         maxptforWD[2]=(tmpstring.Remove(0,namesofSetting[21].Length()+1)).Atof();
1400                 else if (tmpstring.Contains(namesofSetting[22]))
1401                         minptforWD[3]=(tmpstring.Remove(0,namesofSetting[22].Length()+1)).Atof();
1402                 else if (tmpstring.Contains(namesofSetting[23]))
1403                         maxptforWD[3]=(tmpstring.Remove(0,namesofSetting[23].Length()+1)).Atof();
1404                 else if (tmpstring.Contains(namesofSetting[24]))
1405                         minptforWD[4]=(tmpstring.Remove(0,namesofSetting[24].Length()+1)).Atof();
1406                 else if (tmpstring.Contains(namesofSetting[25]))
1407                         maxptforWD[4]=(tmpstring.Remove(0,namesofSetting[25].Length()+1)).Atof();
1408                 else if (tmpstring.Contains(namesofSetting[26]))
1409                         minptforWD[5]=(tmpstring.Remove(0,namesofSetting[26].Length()+1)).Atof();
1410                 else if (tmpstring.Contains(namesofSetting[27]))
1411                         maxptforWD[5]=(tmpstring.Remove(0,namesofSetting[27].Length()+1)).Atof();
1412                 else
1413                         continue;
1414
1415         }
1416         for(int i=0;i<6;i++)
1417                 cout<<minptformaterial[i]<<" "<<maxptformaterial[i]<<" "<<minptforWD[i]<<" "<<maxptforWD[i]<<endl;
1418         cout<<FitRange[0]<<" "<<FitRange[1]<<" "<<CutRange[0]<<CutRange[1]<<endl;
1419         if(FitRange[0]>=FitRange[1])
1420         {
1421                 cout<<"A"<<endl;
1422                 return false;
1423         }
1424         if(CutRange[0]>=CutRange[1])
1425         {
1426                 cout<<"B"<<endl;
1427                 return false;
1428         }
1429         for(int i=0;i<6;i++)
1430         {
1431                 if((minptformaterial[i]>maxptformaterial[i]&&minptformaterial[i]>0.0)||minptformaterial[i]<0.0||maxptformaterial[i]<0.0)
1432                 {
1433                         cout<<"C"<<endl;
1434                         return false;
1435                 }
1436                 if((minptforWD[i]>maxptforWD[i]&&minptforWD[i]>0.0)||minptforWD[i]<0.0||maxptforWD[i]<0.0)
1437                 {
1438                         cout<<"D"<<endl;
1439                         return false;
1440                 }
1441         }
1442
1443
1444
1445         return true;
1446 }
1447
1448 void SubHistWithFullCorr(TH1F* h1, TH1F* h2, Float_t factor1=1.0, Float_t factor2=1.0)
1449 {
1450         if(h1->GetNbinsX()!=h2->GetNbinsX())
1451                 return;
1452         for (int i=0;i<=h1->GetNbinsX();i++)
1453         {
1454                 Float_t tmpvalue=factor1*h1->GetBinContent(i)-factor2*h2->GetBinContent(i);
1455                 Float_t tmperror=TMath::Abs(factor1*factor1*h1->GetBinError(i)*h1->GetBinError(i)-factor2*factor2*h2->GetBinError(i)*h2->GetBinError(i));
1456                 h1->SetBinContent(i,tmpvalue);
1457                 h1->SetBinError(i,TMath::Sqrt(tmperror));
1458         }
1459
1460 }