TENDER becomes Tender
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / AnalysisBoth.C
index 30a0569..d1a198e 100644 (file)
@@ -30,7 +30,7 @@ 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,
@@ -51,14 +51,17 @@ enum {
  kuserangeonfigfile=0x400, // use of config file for dca fit settings
  kskipconcutonspectra=0x800, //do not use conPID<02 cut  useful for syst. studies
  kuseTOFmatchingcorrection=0x1000, // if set tof matching correction is applied.
- kuseTOFcorrforPIDsignalmatching=0x2000 // rescale the for spectra by the factor given in config files                 
-                                                       
+ 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");
@@ -70,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");
@@ -133,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);
        
@@ -177,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);
@@ -238,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();
@@ -263,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();
@@ -291,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"));
@@ -302,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");
@@ -332,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(".","_");
@@ -340,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]);
@@ -644,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();
@@ -988,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/correctionForCrossSection.321.root";
-                       fGeanFlukaK= new TFile(fnameGeanFlukaK.Data());
+                       fGeanFlukaK= TFile::Open(fnameGeanFlukaK.Data());
                        if (!fGeanFlukaK)
                                return;
                }
@@ -1007,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)
@@ -1050,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/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;
          }
@@ -1076,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)
@@ -1393,7 +1463,7 @@ Bool_t ReadConfigFile(TString configfile)
        ifstream infile(configfile.Data());
        if(infile.is_open()==false)
                return false;
-       TString namesofSetting[38]={"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"};   
+       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)
@@ -1479,7 +1549,13 @@ Bool_t ReadConfigFile(TString configfile)
                        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  
+               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;
 
 
@@ -1570,4 +1646,59 @@ void TOFPIDsignalmatchingApply(TH1* h, Float_t 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();
 }