PWGLF/SPECTRA/ChargedHadrons/dNdPt/macros/GenerateCorrMatr_pp.C

   1 Double_t GenerateCorrMatr_pp(const char* mcfile, const char* mctask, const char* idstring ,const char* outfile, const char* gifdir = 0)
   2 //Double_t GenerateCorrMatr_pp()
   3 {
   4 //tmp setting
   5 // const char* mcfile = "/d/alice09/mknichel/train/V007.MC_pp/2011-05-05_2347.7147/mergedPeriods/MC_pp/7TeV/LHC10f6a/mknichel_dNdPtpp_TPCITS.root";
   6 // const char* mctask = "mknichel_dNdPtpp_TPCITS";
   7 // const char* idstring = "TPCITS";
   8 // const char* outfile = "/u/mknichel/alice/dNdPt_pp/temp/corrMatr_LHC10f6a_TPCITS.root";
   9 // const char* gifdir = "/u/mknichel/alice/dNdPt_pp/temp";
  10
  11 TString fname (mcfile);
  12
  13 TString id (idstring);
  14 TString taskname = "mknichel_dNdPtpp_" + id;
  15 TString objname = "dNdPtAnalysis_" + id;
  16
  17     // settings vor zVertex cut (event and track level)
  18     Double_t zVert = 10.0;
  19
  20     // setting on eta cut (track level)
  21     Double_t eta = 0.8;
  22
  23     // strangeness scaling factor (for secondaries from strange decays)
  24     Double_t sscale = 2.0;  //ignored
  25
  26
  27     // Set pt binning
  28     const Int_t ptNbins = 31;
  29     Double_t bins[32] = {0.0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0,2.4,2.8, 3.0,  4.0, 50.0, 100.0 };
  30     Double_t* binsPt = new Double_t[32];
  31     for (int ibin=0; ibin<32; ibin++) {binsPt[ibin] = bins[ibin];}
  32
  33
  34
  35   Int_t ptNbinsTrackEventCorr = 36;
  36   Int_t etaNbins = 30;
  37   Int_t zvNbins = 12;
  38
  39   Double_t binsMultDefault[28] = {-0.5, 0.5 , 1.5 , 2.5 , 3.5 , 4.5 , 5.5 , 6.5 , 7.5 , 8.5,
  40                                      9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5,
  41                                      19.5,20.5, 21.5, 22.5, 23.5, 24.5, 29.5, 149.5};
  42
  43   Double_t binsPtTrackEventCorrDefault[37] = {0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.2,2.4,2.6,3.0,4.0,50.0};
  44
  45   Double_t binsPtDefault[69] = {0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.2,2.4,2.6,2.8,3.0,3.2,3.4,3.6,3.8,4.0,4.5,5.0,5.5,6.0,6.5,7.0,8.0,9.0,10.0,11.0,12.0,13.0,14.0,15.0,16.0,18.0,20.0,22.0,24.0,26.0,28.0,30.0,32.0,34.0,36.0,40.0,45.0,50.0};
  46
  47   Double_t binsEtaDefault[31] = {-1.5,-1.4,-1.3,-1.2,-1.1,-1.0,-0.9,-0.8,-0.7,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,0.,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5};
  48
  49   Double_t binsZvDefault[13] = {-30.,-25.,-20.,-15.,-10.,-5.,0.,5.,10.,15.,20.,25.,30.};
  50
  51     // binning for corrections (note difference for efficiency and secondaries)
  52     if (fname.Contains("2.76TeV") || fname.Contains("900GeV")) {
  53
  54     const Int_t ptNbinsSecCorr  = 28;
  55     const Int_t ptNbinsEffCorr  = 31;
  56     const Int_t etaNbinsCorr = 4;
  57     const Int_t zvNbinsCorr  = 4;
  58     const Int_t multNbinsCorr = 2;
  59     Double_t ptBinsSecCorr[29] = {0.0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0, 1.2, 1.4, 1.6, 2.0, 2.4, 3.0, 50.0, 100.0 };
  60     Double_t ptBinsEffCorr[32] = {0.0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0, 1.2, 1.4, 1.6, 2.0, 2.4, 2.8, 3.0, 3.4, 4.0, 50.0, 100.0 };
  61     Double_t etaBinsCorr[5] = {-0.8,-0.4,0.,0.4,0.8};
  62
  63     Double_t zvBinsCorr[5] = {-10.,-5.,0.,5.,10.};
  64     Double_t multBinsCorr[3] = {-0.5,0.5,500.};
  65
  66
  67     } else {
  68
  69     const Int_t ptNbinsSecCorr  = 38;
  70     const Int_t ptNbinsEffCorr  = 50;
  71     const Int_t etaNbinsCorr = 8;
  72     const Int_t zvNbinsCorr  = 4;
  73     const Int_t multNbinsCorr = 2;
  74     Double_t ptBinsSecCorr[39] = {0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.2,2.4,2.6,3.0,3.4,4.0,50.0,100.0};
  75     Double_t ptBinsEffCorr[51] = {0.0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 8.0, 10.0, 50.0, 100.0 };
  76     Double_t etaBinsCorr[9] = {-0.8,-0.6,-0.4,-0.2,0.,0.2,0.4,0.6,0.8};
  77
  78     Double_t zvBinsCorr[5] = {-10.,-5.,0.,5.,10.};
  79     Double_t multBinsCorr[3] = {-0.5,0.5,500.};
  80
  81     }
  82
  83
  84     Int_t binsTrackEffCorr[3]={zvNbinsCorr,ptNbinsEffCorr,etaNbinsCorr};
  85     Int_t binsTrackSecCorr[3]={zvNbinsCorr,ptNbinsSecCorr,etaNbinsCorr};
  86
  87     THnSparse* hSTrackEffCorr = new THnSparseF("hSTrackEffCorr","Zv:pT:eta",3,binsTrackEffCorr);
  88     hSTrackEffCorr->SetBinEdges(0,zvBinsCorr);
  89     hSTrackEffCorr->SetBinEdges(1,ptBinsEffCorr);
  90     hSTrackEffCorr->SetBinEdges(2,etaBinsCorr);
  91 //     hSTrackEffCorr->SetBinEdges(3,multBinsCorr);
  92     hSTrackEffCorr->GetAxis(0)->SetTitle("Zv (cm)");
  93     hSTrackEffCorr->GetAxis(1)->SetTitle("p_{T} (GeV/c)");
  94     hSTrackEffCorr->GetAxis(2)->SetTitle("#eta");
  95 //     hSTrackEffCorr->GetAxis(3)->SetTitle("mult MB");
  96     hSTrackEffCorr->Sumw2();
  97
  98     THnSparse* hSTrackSecCorr = new THnSparseF("hSTrackSecCorr","Zv:pT:eta",3,binsTrackSecCorr);
  99     hSTrackSecCorr->SetBinEdges(0,zvBinsCorr);
 100     hSTrackSecCorr->SetBinEdges(1,ptBinsSecCorr);
 101     hSTrackSecCorr->SetBinEdges(2,etaBinsCorr);
 102 //     hSTrackSecCorr->SetBinEdges(3,multBinsCorr);
 103     hSTrackSecCorr->GetAxis(0)->SetTitle("Zv (cm)");
 104     hSTrackSecCorr->GetAxis(1)->SetTitle("p_{T} (GeV/c)");
 105     hSTrackSecCorr->GetAxis(2)->SetTitle("#eta");
 106 //     hSTrackSecCorr->GetAxis(3)->SetTitle("mult MB");
 107     hSTrackSecCorr->Sumw2();
 108
 109
 110     // make plots nicer
 111     gROOT->SetStyle("Plain");
 112     gStyle->SetPalette(1);
 113
 114     // array for all correction matrices
 115     TObjArray* CorrMatr = new TObjArray();
 116
 117     // array for all control histograms
 118     TObjArray* ContHist = new TObjArray();
 119
 120
 121     // load mc information
 122     TFile* fmc = TFile::Open(mcfile,"READ");
 123     if (!fmc) return -1;
 124     TList* lmc = dynamic_cast<TList*>(fmc->Get(taskname.Data()));
 125     if (!lmc) return -1;
 126     AlidNdPtAnalysis *obj = dynamic_cast<AlidNdPtAnalysis*>(lmc->FindObject(objname.Data()));
 127     if (!obj) return -1;
 128
 129     //Event statistics
 130     THnSparse *fRecEventMatrix = obj->GetRecEventMatrix(); //all reconstructed events
 131     TH2D* h2RecEventAll = fRecEventMatrix->Projection(0,1)->Clone("h2RecEventAll");
 132     fRecEventMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
 133     TH2D* h2RecEvent = fRecEventMatrix->Projection(0,1)->Clone("h2RecEvent");
 134     Double_t MCReconstructedEvents = h2RecEvent->Integral();
 135     Double_t MCReconstructedEventsAll = h2RecEventAll->Integral();
 136     ContHist->Add(h2RecEvent);
 137     ContHist->Add(h2RecEventAll);
 138
 139     THnSparse *fTriggerEventMatrix = obj->GetTriggerEventMatrix(); //all triggered events
 140     Double_t MCTriggeredEventsAll0mult =  obj->GetTriggerEventMatrix()->Projection(1)->Integral(1,1); // 0 contributers to vertex
 141     TH2D* h2TriggerEventAll = fTriggerEventMatrix->Projection(0,1)->Clone("h2TriggerEventAll");
 142     fTriggerEventMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
 143     Double_t MCTriggeredEvents0mult =  obj->GetTriggerEventMatrix()->Projection(1)->Integral(1,1); // 0 contributers to vertex
 144     TH2D* h2TriggerEvent = fTriggerEventMatrix->Projection(0,1)->Clone("h2TriggerEvent");
 145     Double_t MCTriggeredEvents = h2TriggerEvent->Integral();
 146     Double_t MCTriggeredEventsAll = h2TriggerEventAll->Integral();
 147     ContHist->Add(h2TriggerEvent);
 148     ContHist->Add(h2TriggerEventAll);
 149
 150     THnSparse *fGenEventMatrix = obj->GetGenEventMatrix(); //all generated events
 151     TH2D* h2GenEventAll = fGenEventMatrix->Projection(0,1)->Clone("h2GenEventAll");
 152     fGenEventMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
 153     TH2D* h2GenEvent = fGenEventMatrix->Projection(0,1)->Clone("h2GenEvent");
 154     Double_t MCGeneratedEvents = h2GenEvent->Integral();
 155     Double_t MCGeneratedEventsAll = h2GenEventAll->Integral();
 156     ContHist->Add(h2RecEvent);
 157     ContHist->Add(h2RecEventAll);
 158
 159     THnSparse *fRecEventHist = obj->GetRecEventHist();
 160     THnSparse *fRecEventHist2 = obj->GetRecEventHist2();
 161
 162     printf("=== generated MC events for correction matrices    %lf ===\n",MCGeneratedEvents);
 163     printf("=== triggered MC events for correction matrices    %lf ===\n",MCTriggeredEvents);
 164     printf("=== recontructed MC events for correction matrices %lf ===\n",MCReconstructedEvents);
 165     printf("\n");
 166     printf("=== cut on the zVertex +- %lf ===\n",zVert);
 167     printf("=== generated MC events (before zVertex cut)       %lf ===\n",MCGeneratedEventsAll);
 168     printf("=== triggered MC events (before zVertex cut)       %lf ===\n",MCTriggeredEventsAll);
 169     printf("=== recontructed MC events (before zVertex cut)    %lf ===\n",MCReconstructedEventsAll);
 170
 171     TH1D* h1MCGeneratedEvents = new TH1D("h1MCGeneratedEvents","h1MCGeneratedEvents",1,0,1);
 172     TH1D* h1MCTriggeredEvents = new TH1D("h1MCTriggeredEvents","h1MCTriggeredEvents",1,0,1);
 173     TH1D* h1MCReconstructedEvents = new TH1D("h1MCReconstructedEvents","h1MCReconstructedEvents",1,0,1);
 174     TH1D* h1MCGeneratedEventsAll = new TH1D("h1MCGeneratedEventsAll","h1MCGeneratedEventsAll",1,0,1);
 175     TH1D* h1MCTriggeredEventsAll = new TH1D("h1MCTriggeredEventsAll","h1MCTriggeredEventsAll",1,0,1);
 176     TH1D* h1MCReconstructedEventsAll = new TH1D("h1MCReconstructedEventsAll","h1MCReconstructedEventsAll",1,0,1);
 177     TH1D* h1MCTriggeredEventsAll0mult = new TH1D("h1MCTriggeredEventsAll0mult","h1MCTriggeredEventsAll0mult",1,0,1);
 178     TH1D* h1MCTriggeredEvents0mult = new TH1D("h1MCTriggeredEvents0mult","h1MCTriggeredEvents0mult",1,0,1);
 179
 180     h1MCGeneratedEvents->Fill(0,MCGeneratedEvents);
 181     h1MCGeneratedEvents->SetEntries(MCGeneratedEvents);
 182     h1MCTriggeredEvents->Fill(0,MCTriggeredEvents);
 183     h1MCTriggeredEvents->SetEntries(MCTriggeredEvents);
 184     h1MCReconstructedEvents->Fill(0,MCReconstructedEvents);
 185     h1MCReconstructedEvents->SetEntries(MCReconstructedEvents);
 186     h1MCGeneratedEventsAll->Fill(0,MCGeneratedEventsAll);
 187     h1MCGeneratedEventsAll->SetEntries(MCGeneratedEventsAll);
 188     h1MCTriggeredEventsAll->Fill(0,MCTriggeredEventsAll);
 189     h1MCTriggeredEventsAll->SetEntries(MCTriggeredEventsAll);
 190     h1MCReconstructedEventsAll->Fill(0,MCReconstructedEventsAll);
 191     h1MCReconstructedEventsAll->SetEntries(MCReconstructedEventsAll);
 192     h1MCTriggeredEventsAll0mult->Fill(0,MCTriggeredEventsAll0mult);
 193     h1MCTriggeredEventsAll0mult->SetEntries(MCTriggeredEventsAll0mult);
 194     h1MCTriggeredEvents0mult->Fill(0,MCTriggeredEvents0mult);
 195     h1MCTriggeredEvents0mult->SetEntries(MCTriggeredEvents0mult);
 196
 197
 198     ContHist->Add(h1MCGeneratedEvents);
 199     ContHist->Add(h1MCTriggeredEvents);
 200     ContHist->Add(h1MCReconstructedEvents);
 201     ContHist->Add(h1MCGeneratedEventsAll);
 202     ContHist->Add(h1MCTriggeredEventsAll);
 203     ContHist->Add(h1MCReconstructedEventsAll);
 204     ContHist->Add(h1MCTriggeredEvents0mult);
 205     ContHist->Add(h1MCTriggeredEventsAll0mult);
 206
 207
 208     // efficienfy and correction matrices for tigger and vertex efficiency
 209     TH2D* h2EventTriggerEffAll  = AlidNdPtHelper::GenerateCorrMatrix(h2TriggerEventAll,h2GenEventAll,"h2EventTriggerEffAll");
 210     TH2D* h2EventTriggerCorrAll = AlidNdPtHelper::GenerateCorrMatrix(h2GenEventAll,h2TriggerEventAll,"h2EventTriggerCorrAll");
 211     TH2D* h2EventTriggerEff  = AlidNdPtHelper::GenerateCorrMatrix(h2TriggerEvent,h2GenEvent,"h2EventTriggerEff");
 212     TH2D* h2EventTriggerCorr = AlidNdPtHelper::GenerateCorrMatrix(h2GenEvent,h2TriggerEvent,"h2EventTriggerCorr");
 213
 214     TH2D* h2EventRecEffAll  = AlidNdPtHelper::GenerateCorrMatrix(h2RecEventAll,h2TriggerEventAll,"h2EventRecEffAll");
 215     TH2D* h2EventRecCorrAll = AlidNdPtHelper::GenerateCorrMatrix(h2TriggerEventAll,h2RecEventAll,"h2EventRecCorrAll");
 216     TH2D* h2EventRecEff  = AlidNdPtHelper::GenerateCorrMatrix(h2RecEvent,h2TriggerEvent,"h2EventRecEff");
 217     TH2D* h2EventRecCorr = AlidNdPtHelper::GenerateCorrMatrix(h2TriggerEvent,h2RecEvent,"h2EventRecCorr");
 218
 219     TH2D* h2EventEffAll  = AlidNdPtHelper::GenerateCorrMatrix(h2RecEventAll,h2GenEventAll,"h2EventEffAll");
 220     TH2D* h2EventCorrAll = AlidNdPtHelper::GenerateCorrMatrix(h2GenEventAll,h2RecEventAll,"h2EventCorrAll");
 221     TH2D* h2EventEff  = AlidNdPtHelper::GenerateCorrMatrix(h2RecEvent,h2GenEvent,"h2EventEff");
 222     TH2D* h2EventCorr = AlidNdPtHelper::GenerateCorrMatrix(h2GenEvent,h2RecEvent,"h2EventCorr");
 223
 224     CorrMatr->Add(h2EventTriggerEffAll);
 225     CorrMatr->Add(h2EventTriggerCorrAll);
 226     CorrMatr->Add(h2EventTriggerEff);
 227     CorrMatr->Add(h2EventTriggerCorr);
 228     CorrMatr->Add(h2EventRecEffAll);
 229     CorrMatr->Add(h2EventRecCorrAll);
 230     CorrMatr->Add(h2EventRecEff);
 231     CorrMatr->Add(h2EventRecCorr);
 232     CorrMatr->Add(h2EventEffAll);
 233     CorrMatr->Add(h2EventCorrAll);
 234     CorrMatr->Add(h2EventEff);
 235     CorrMatr->Add(h2EventCorr);
 236
 237     THnSparse* hSMultCorrelation = obj->GetEventMultCorrelationMatrix()->Clone("hSMultCorrelation");
 238     CorrMatr->Add(hSMultCorrelation);
 239     TH2D* h2MultCorrelation = hSMultCorrelation->Projection(0,1)->Clone("h2MultCorrelation");
 240     CorrMatr->Add(h2MultCorrelation);
 241
 242     // all recontructed
 243     THnSparse *fRecTrackMatrix = obj->GetRecTrackMatrix();
 244     fRecTrackMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
 245     fRecTrackMatrix->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
 246     TH3D* h3RecTrack = fRecTrackMatrix->Projection(0,1,2)->Clone("h3RecTrack");
 247     TH2D* h2RecTrack_zv_pt  = fRecTrackMatrix->Projection(0,1)->Clone("h2RecTrack_zv_pt");
 248     TH2D* h2RecTrack_zv_eta = fRecTrackMatrix->Projection(0,2)->Clone("h2RecTrack_zv_eta");
 249     TH2D* h2RecTrack_pt_eta = fRecTrackMatrix->Projection(1,2)->Clone("h2RecTrack_pt_eta");
 250     TH1D* h1RecTrack_zv  = fRecTrackMatrix->Projection(0)->Clone("h1RecTrack_zv");
 251     TH1D* h1RecTrack_pt  = fRecTrackMatrix->Projection(1)->Clone("h1RecTrack_pt");
 252     TH1D* h1RecTrack_eta = fRecTrackMatrix->Projection(2)->Clone("h1RecTrack_eta");
 253     Double_t MCReconstructedTracks = h3RecTrack->Integral();
 254
 255     ContHist->Add(h3RecTrack);
 256     ContHist->Add(h2RecTrack_zv_pt);
 257     ContHist->Add(h2RecTrack_zv_eta);
 258     ContHist->Add(h2RecTrack_pt_eta);
 259     ContHist->Add(h1RecTrack_zv);
 260     ContHist->Add(h1RecTrack_pt);
 261     ContHist->Add(h1RecTrack_eta);
 262
 263      // recontructed primary tracks
 264     THnSparse *fRecPrimTrackMatrix = obj->GetRecPrimTrackMatrix();
 265     THnSparse *fRecTrackMatrixScaled = fRecPrimTrackMatrix->Clone("fRecTrackMatrixScaled"); //used later for secondaries scaling
 266     fRecPrimTrackMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
 267     fRecPrimTrackMatrix->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
 268     TH3D* h3RecPrimTrack = fRecPrimTrackMatrix->Projection(0,1,2)->Clone("h3RecPrimTrack");
 269     TH2D* h2RecPrimTrack_zv_pt  = fRecPrimTrackMatrix->Projection(0,1)->Clone("h2RecPrimTrack_zv_pt");
 270     TH2D* h2RecPrimTrack_zv_eta = fRecPrimTrackMatrix->Projection(0,2)->Clone("h2RecPrimTrack_zv_eta");
 271     TH2D* h2RecPrimTrack_pt_eta = fRecPrimTrackMatrix->Projection(1,2)->Clone("h2RecPrimTrack_pt_eta");
 272     TH1D* h1RecPrimTrack_zv  = fRecPrimTrackMatrix->Projection(0)->Clone("h1RecPrimTrack_zv");
 273     TH1D* h1RecPrimTrack_pt  = fRecPrimTrackMatrix->Projection(1)->Clone("h1RecPrimTrack_pt");
 274     TH1D* h1RecPrimTrack_eta = fRecPrimTrackMatrix->Projection(2)->Clone("h1RecPrimTrack_eta");
 275     Double_t MCReconstructedPrimTracks = h3RecPrimTrack->Integral();
 276
 277     ContHist->Add(h3RecPrimTrack);
 278     ContHist->Add(h2RecPrimTrack_zv_pt);
 279     ContHist->Add(h2RecPrimTrack_zv_eta);
 280     ContHist->Add(h2RecPrimTrack_pt_eta);
 281     ContHist->Add(h1RecPrimTrack_zv);
 282     ContHist->Add(h1RecPrimTrack_pt);
 283     ContHist->Add(h1RecPrimTrack_eta);
 284
 285     // recontructed secondary tracks
 286     THnSparse *fRecSecTrackMatrix = obj->GetRecSecTrackMatrix();
 287     THnSparse *fRecSecTrackMatrixScaled = fRecSecTrackMatrix->Clone("fRecSecTrackMatrixScaled"); //used later for secondaries scaling
 288     fRecSecTrackMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
 289     fRecSecTrackMatrix->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
 290     TH3D* h3RecSecTrack = fRecSecTrackMatrix->Projection(0,1,2)->Clone("h3RecSecTrack");
 291     TH2D* h2RecSecTrack_zv_pt  = fRecSecTrackMatrix->Projection(0,1)->Clone("h2RecSecTrack_zv_pt");
 292     TH2D* h2RecSecTrack_zv_eta = fRecSecTrackMatrix->Projection(0,2)->Clone("h2RecSecTrack_zv_eta");
 293     TH2D* h2RecSecTrack_pt_eta = fRecSecTrackMatrix->Projection(1,2)->Clone("h2RecSecTrack_pt_eta");
 294     TH1D* h1RecSecTrack_zv  = fRecSecTrackMatrix->Projection(0)->Clone("h1RecSecTrack_zv");
 295     TH1D* h1RecSecTrack_pt  = fRecSecTrackMatrix->Projection(1)->Clone("h1RecSecTrack_pt");
 296     TH1D* h1RecSecTrack_eta = fRecSecTrackMatrix->Projection(2)->Clone("h1RecSecTrack_eta");
 297     Double_t MCReconstructedSecTracks = h3RecSecTrack->Integral();
 298
 299     ContHist->Add(h3RecSecTrack);
 300     ContHist->Add(h2RecSecTrack_zv_pt);
 301     ContHist->Add(h2RecSecTrack_zv_eta);
 302     ContHist->Add(h2RecSecTrack_pt_eta);
 303     ContHist->Add(h1RecSecTrack_zv);
 304     ContHist->Add(h1RecSecTrack_pt);
 305     ContHist->Add(h1RecSecTrack_eta);
 306
 307     // generated primary tracks
 308     THnSparse *fGenPrimTrackMatrix = obj->GetGenPrimTrackMatrix();
 309     fGenPrimTrackMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
 310     fGenPrimTrackMatrix->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
 311     TH3D* h3GenPrimTrack = fGenPrimTrackMatrix->Projection(0,1,2)->Clone("h3GenPrimTrack");
 312     TH2D* h2GenPrimTrack_zv_Pt  = fGenPrimTrackMatrix->Projection(0,1)->Clone("h2GenPrimTrack_zv_pt");
 313     TH2D* h2GenPrimTrack_zv_eta = fGenPrimTrackMatrix->Projection(0,2)->Clone("h2GenPrimTrack_zv_eta");
 314     TH2D* h2GenPrimTrack_pt_eta = fGenPrimTrackMatrix->Projection(1,2)->Clone("h2GenPrimTrack_pt_eta");
 315     TH1D* h1GenPrimTrack_zv  = fGenPrimTrackMatrix->Projection(0)->Clone("h1GenPrimTrack_zv");
 316     TH1D* h1GenPrimTrack_pt  = fGenPrimTrackMatrix->Projection(1)->Clone("h1GenPrimTrack_pt");
 317     TH1D* h1GenPrimTrack_eta = fGenPrimTrackMatrix->Projection(2)->Clone("h1GenPrimTrack_eta");
 318     Double_t MCGeneratedPrimTracks = h3GenPrimTrack->Integral();
 319
 320     ContHist->Add(h3GenPrimTrack);
 321     ContHist->Add(h2GenPrimTrack_zv_pt);
 322     ContHist->Add(h2GenPrimTrack_zv_eta);
 323     ContHist->Add(h2GenPrimTrack_pt_eta);
 324     ContHist->Add(h1GenPrimTrack_zv);
 325     ContHist->Add(h1GenPrimTrack_pt);
 326     ContHist->Add(h1GenPrimTrack_eta);
 327     printf("\n");
 328     printf("==============================================================\n");
 329     printf("=== recontructed MC tracks              %lf ===\n",MCReconstructedTracks);
 330     printf("=== recontructed MC secondary tracks    %lf ===\n",MCReconstructedSecTracks);
 331     printf("=== recontructed MC primary tracks      %lf ===\n",MCReconstructedPrimTracks);
 332     printf("=== generated MC primary track          %lf ===\n",MCGeneratedPrimTracks);
 333     printf("==============================================================\n");
 334     printf("\n");
 335
 336     // mc truth histogram (for self-consistency check)
 337     TH1D* dNdPt_MC   = (TH1D*) h1GenPrimTrack_pt->Clone("dNdPt_MC");
 338     TH1D* dNdEta_MC   =(TH1D*) h1GenPrimTrack_eta->Clone("dNdEta_MC");
 339     // normalization and finalization
 340     // 1/N_evt 1/(2 pi pt) 1/width 1/etarange
 341    for (int ii=1; ii <= dNdPt_MC->GetNbinsX() ;ii++) {
 342         Double_t pt = dNdPt_MC->GetBinCenter(ii);
 343         Double_t width = dNdPt_MC->GetBinWidth(ii);
 344         Double_t val = dNdPt_MC->GetBinContent(ii);
 345         Double_t err = dNdPt_MC->GetBinError(ii);
 346         Double_t cval = (val)/(width * 2.0 * TMath::Pi() * 1.6 * MCGeneratedEvents * pt);
 347         Double_t cerr = (err)/(width * 2.0 * TMath::Pi() * 1.6 * MCGeneratedEvents * pt);
 348         dNdPt_MC->SetBinContent(ii,cval);
 349         dNdPt_MC->SetBinError(ii,cerr);
 350     }
 351     dNdPt_MC->SetMarkerStyle(21);
 352     dNdPt_MC->SetTitle("; p_{T} (GeV/c) ; 1/N_{evt} 1/(2#pi p_{T}) (d^{2}N_{ch})/(d#eta dp_{T})^{-2}");
 353     dNdEta_MC->Scale(1,"width");
 354     dNdEta_MC->Scale(1./MCGeneratedEvents);
 355     ContHist->Add(dNdPt_MC);
 356     ContHist->Add(dNdEta_MC);
 357
 358
 359    // Rebin for corrections (pt)
 360    cout << "rebinning...." << endl;
 361
 362    TH1D* h1RecPrimTrack_pt_Rebin = h1RecPrimTrack_pt->Rebin(ptNbins,"h1RecPrimTrack_pt_Rebin",binsPt);
 363    TH1D* h1GenPrimTrack_pt_Rebin = h1GenPrimTrack_pt->Rebin(ptNbins,"h1GenPrimTrack_pt_Rebin",binsPt);
 364    ContHist->Add(h1RecPrimTrack_pt_Rebin);
 365    ContHist->Add(h1GenPrimTrack_pt_Rebin);
 366
 367
 368 //    THnSparse *fSparseTriggerTrackEvent = obj->GetTriggerTrackEventMatrix();//Tracks from triggered events
 369 //    THnSparse *fSparseVtxTrackEvent = obj->GetRecTrackEventMatrix();//Tracks from events with rec. vtx
 370 //    THnSparse *fSparseGenTrackEvent = obj->GetGenTrackEventMatrix();//generated TrackEvent matrix
 371
 372     // tracking efficiencies + corrections
 373    TH2D* h2TrackEff_zv_pt   = AlidNdPtHelper::GenerateCorrMatrix(h2RecPrimTrack_zv_pt,h2GenPrimTrack_zv_pt,"h2TrackEff_zv_pt");
 374    TH2D* h2TrackCorr_zv_pt  = AlidNdPtHelper::GenerateCorrMatrix(h2GenPrimTrack_zv_pt,h2RecPrimTrack_zv_pt,"h2TrackCorr_zv_pt");
 375    TH2D* h2TrackEff_zv_eta  = AlidNdPtHelper::GenerateCorrMatrix(h2RecPrimTrack_zv_eta,h2GenPrimTrack_zv_eta,"h2TrackEff_zv_eta");
 376    TH2D* h2TrackCorr_zv_eta = AlidNdPtHelper::GenerateCorrMatrix(h2GenPrimTrack_zv_eta,h2RecPrimTrack_zv_eta,"h2TrackCorr_zv_eta");
 377    TH2D* h2TrackEff_pt_eta  = AlidNdPtHelper::GenerateCorrMatrix(h2RecPrimTrack_pt_eta,h2GenPrimTrack_pt_eta,"h2TrackEff_pt_eta");
 378    TH2D* h2TrackCorr_pt_eta = AlidNdPtHelper::GenerateCorrMatrix(h2GenPrimTrack_pt_eta,h2RecPrimTrack_pt_eta,"h2TrackCorr_pt_eta");
 379
 380
 381    TH1D* h1TrackEff_zv   = AlidNdPtHelper::GenerateCorrMatrix(h1RecPrimTrack_zv,h1GenPrimTrack_zv,"h1TrackEff_zv");
 382    TH1D* h1TrackCorr_zv  = AlidNdPtHelper::GenerateCorrMatrix(h1GenPrimTrack_zv,h1RecPrimTrack_zv,"h1TrackCorr_zv");
 383    TH1D* h1TrackEff_pt   = AlidNdPtHelper::GenerateCorrMatrix(h1RecPrimTrack_pt_Rebin,h1GenPrimTrack_pt_Rebin,"h1TrackEff_pt");
 384    TH1D* h1TrackCorr_pt  = AlidNdPtHelper::GenerateCorrMatrix(h1GenPrimTrack_pt_Rebin,h1RecPrimTrack_pt_Rebin,"h1TrackCorr_pt");
 385    TH1D* h1TrackEff_eta  = AlidNdPtHelper::GenerateCorrMatrix(h1RecPrimTrack_eta,h1GenPrimTrack_eta,"h1TrackEff_eta");
 386    TH1D* h1TrackCorr_eta = AlidNdPtHelper::GenerateCorrMatrix(h1GenPrimTrack_eta,h1RecPrimTrack_eta,"h1TrackCorr_eta");
 387    CorrMatr->Add(h2TrackEff_zv_pt);
 388    CorrMatr->Add(h2TrackCorr_zv_pt);
 389    CorrMatr->Add(h2TrackEff_zv_eta);
 390    CorrMatr->Add(h2TrackCorr_zv_eta);
 391    CorrMatr->Add(h2TrackEff_pt_eta);
 392    CorrMatr->Add(h2TrackCorr_pt_eta);
 393    CorrMatr->Add(h1TrackEff_zv);
 394    CorrMatr->Add(h1TrackCorr_zv);
 395    CorrMatr->Add(h1TrackEff_pt);
 396    CorrMatr->Add(h1TrackCorr_pt);
 397    CorrMatr->Add(h1TrackEff_eta);
 398    CorrMatr->Add(h1TrackCorr_eta);
 399
 400    // scale the secondaries before calculating correction matrices
 401    for (Long64_t i = 0; i < fRecSecTrackMatrixScaled->GetNbins(); i++) {
 402        Int_t c[3];
 403        Double_t val = fRecSecTrackMatrixScaled->GetBinContent(i,c);
 404        Double_t err = fRecSecTrackMatrixScaled->GetBinError(i);
 405        Double_t pt =  fRecSecTrackMatrixScaled->GetAxis(1)->GetBinCenter(c[1]);
 406        Double_t scale = GetStrangenessCorrFactor(pt,sscale);
 407 //        Double_t scale = AlidNdPtHelper::GetStrangenessCorrFactor(pt);
 408        fRecSecTrackMatrixScaled->SetBinContent(c,val*scale);
 409        fRecSecTrackMatrixScaled->SetBinError(c,err*scale);
 410     }
 411
 412     // for correct determination of secondaries contamination, also the total total tracks have to be scaled
 413     // this is done by taking primaries and adding the scaled secondaries
 414     fRecTrackMatrixScaled->Add(fRecSecTrackMatrixScaled);
 415
 416     fRecSecTrackMatrixScaled->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
 417     fRecSecTrackMatrixScaled->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
 418
 419     TH3D* h3RecSecTrackScaled = fRecSecTrackMatrixScaled->Projection(0,1,2)->Clone("h3RecSecTrackScaled");
 420     TH2D* h2RecSecTrackScaled_zv_pt  = fRecSecTrackMatrixScaled->Projection(0,1)->Clone("h2RecSecTrackScaled_zv_pt");
 421     TH2D* h2RecSecTrackScaled_zv_eta = fRecSecTrackMatrixScaled->Projection(0,2)->Clone("h2RecSecTrackScaled_zv_eta");
 422     TH2D* h2RecSecTrackScaled_pt_eta = fRecSecTrackMatrixScaled->Projection(1,2)->Clone("h2RecSecTrackScaled_pt_eta");
 423     TH1D* h1RecSecTrackScaled_zv  = fRecSecTrackMatrixScaled->Projection(0)->Clone("h1RecSecTrackScaled_zv");
 424     TH1D* h1RecSecTrackScaled_pt  = fRecSecTrackMatrixScaled->Projection(1)->Clone("h1RecSecTrackScaled_pt");
 425     TH1D* h1RecSecTrackScaled_eta = fRecSecTrackMatrixScaled->Projection(2)->Clone("h1RecSecTrackScaled_eta");
 426
 427     ContHist->Add(h3RecSecTrackScaled);
 428     ContHist->Add(h2RecSecTrackScaled_zv_pt);
 429     ContHist->Add(h2RecSecTrackScaled_zv_eta);
 430     ContHist->Add(h2RecSecTrackScaled_pt_eta);
 431     ContHist->Add(h1RecSecTrackScaled_zv);
 432     ContHist->Add(h1RecSecTrackScaled_pt);
 433     ContHist->Add(h1RecSecTrackScaled_eta);
 434
 435     fRecTrackMatrixScaled->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
 436     fRecTrackMatrixScaled->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
 437
 438     TH3D* h3RecTrackScaled = fRecTrackMatrixScaled->Projection(0,1,2)->Clone("h3RecTrackScaled");
 439     TH2D* h2RecTrackScaled_zv_pt  = fRecTrackMatrixScaled->Projection(0,1)->Clone("h2RecTrackScaled_zv_pt");
 440     TH2D* h2RecTrackScaled_zv_eta = fRecTrackMatrixScaled->Projection(0,2)->Clone("h2RecTrackScaled_zv_eta");
 441     TH2D* h2RecTrackScaled_pt_eta = fRecTrackMatrixScaled->Projection(1,2)->Clone("h2RecTrackScaled_pt_eta");
 442     TH1D* h1RecTrackScaled_zv  = fRecTrackMatrixScaled->Projection(0)->Clone("h1RecTrackScaled_zv");
 443     TH1D* h1RecTrackScaled_pt  = fRecTrackMatrixScaled->Projection(1)->Clone("h1RecTrackScaled_pt");
 444     TH1D* h1RecTrackScaled_eta = fRecTrackMatrixScaled->Projection(2)->Clone("h1RecTrackScaled_eta");
 445
 446     ContHist->Add(h3RecTrackScaled);
 447     ContHist->Add(h2RecTrackScaled_zv_pt);
 448     ContHist->Add(h2RecTrackScaled_zv_eta);
 449     ContHist->Add(h2RecTrackScaled_pt_eta);
 450     ContHist->Add(h1RecTrackScaled_zv);
 451     ContHist->Add(h1RecTrackScaled_pt);
 452     ContHist->Add(h1RecTrackScaled_eta);
 453
 454    // Rebin for corrections (pt)
 455    TH1D* h1RecTrackScaled_pt_Rebin = h1RecTrackScaled_pt->Rebin(ptNbins,"h1RecTrackScaled_pt_Rebin",binsPt);
 456    TH1D* h1RecSecTrackScaled_pt_Rebin = h1RecSecTrackScaled_pt->Rebin(ptNbins,"h1RecSecTrackScaled_pt_Rebin",binsPt);
 457    ContHist->Add(h1RecTrackScaled_pt_Rebin);
 458    ContHist->Add(h1RecSecTrackScaled_pt_Rebin);
 459
 460     // create histograms for secondaries contamination and correction
 461
 462     TH2D* h2SecCont_zv_pt  = AlidNdPtHelper::GenerateCorrMatrix(h2RecSecTrackScaled_zv_pt,h2RecTrackScaled_zv_pt,"h2SecCont_zv_pt");
 463     TH2D* h2SecCorr_zv_pt  = AlidNdPtHelper::GenerateContCorrMatrix(h2RecSecTrackScaled_zv_pt,h2RecTrackScaled_zv_pt,"h2SecCorr_zv_pt");
 464     TH2D* h2SecCont_zv_eta = AlidNdPtHelper::GenerateCorrMatrix(h2RecSecTrackScaled_zv_eta,h2RecTrackScaled_zv_eta,"h2SecCont_zv_eta");
 465     TH2D* h2SecCorr_zv_eta = AlidNdPtHelper::GenerateContCorrMatrix(h2RecSecTrackScaled_zv_eta,h2RecTrackScaled_zv_eta,"h2SecCorr_zv_eta");
 466     TH2D* h2SecCont_pt_eta = AlidNdPtHelper::GenerateCorrMatrix(h2RecSecTrackScaled_pt_eta,h2RecTrackScaled_pt_eta,"h2SecCont_pt_eta");
 467     TH2D* h2SecCorr_pt_eta = AlidNdPtHelper::GenerateContCorrMatrix(h2RecSecTrackScaled_pt_eta,h2RecTrackScaled_pt_eta,"h2SecCorr_pt_eta");
 468     TH1D* h1SecCont_zv = AlidNdPtHelper::GenerateCorrMatrix(h1RecSecTrackScaled_zv,h1RecTrackScaled_zv,"h1SecCont_zv");
 469     TH1D* h1SecCorr_zv = AlidNdPtHelper::GenerateContCorrMatrix(h1RecSecTrackScaled_zv,h1RecTrackScaled_zv,"h1SecCorr_zv");
 470     TH1D* h1SecCont_pt = AlidNdPtHelper::GenerateCorrMatrix(h1RecSecTrackScaled_pt_Rebin,h1RecTrackScaled_pt_Rebin,"h1SecCont_pt");
 471     TH1D* h1SecCorr_pt = AlidNdPtHelper::GenerateContCorrMatrix(h1RecSecTrackScaled_pt_Rebin,h1RecTrackScaled_pt_Rebin,"h1SecCorr_pt");
 472     TH1D* h1SecCont_eta = AlidNdPtHelper::GenerateCorrMatrix(h1RecSecTrackScaled_eta,h1RecTrackScaled_eta,"h1SecCont_eta");
 473     TH1D* h1SecCorr_eta = AlidNdPtHelper::GenerateContCorrMatrix(h1RecSecTrackScaled_eta,h1RecTrackScaled_eta,"h1SecCorr_eta");
 474
 475     CorrMatr->Add(h2SecCont_zv_pt);
 476     CorrMatr->Add(h2SecCorr_zv_pt);
 477     CorrMatr->Add(h2SecCont_zv_eta);
 478     CorrMatr->Add(h2SecCorr_zv_eta);
 479     CorrMatr->Add(h2SecCont_pt_eta);
 480     CorrMatr->Add(h2SecCorr_pt_eta);
 481     CorrMatr->Add(h1SecCont_zv);
 482     CorrMatr->Add(h1SecCorr_zv);
 483     CorrMatr->Add(h1SecCont_pt);
 484     CorrMatr->Add(h1SecCorr_pt);
 485     CorrMatr->Add(h1SecCont_eta);
 486     CorrMatr->Add(h1SecCorr_eta);
 487
 488     //vertex distribution for events with no vertex
 489     fGenEventMatrix->GetAxis(0)->SetRange();
 490     fGenEventMatrix->GetAxis(1)->SetRange(1,1);
 491     TH1D* h1GenEventMatrix_bin0_zv = (TH1D*) fGenEventMatrix->Projection(0)->Clone("h1GenEventMatrix_bin0_zv");
 492     ContHist->Add(h1GenEventMatrix_bin0_zv);
 493
 494     fTriggerEventMatrix->GetAxis(0)->SetRange();
 495     fTriggerEventMatrix->GetAxis(1)->SetRange(1,1);
 496     TH1D* h1TriggerEventMatrix_bin0_zv = (TH1D*) fTriggerEventMatrix->Projection(0)->Clone("h1TriggerEventMatrix_bin0_zv");
 497     ContHist->Add(h1TriggerEventMatrix_bin0_zv);
 498
 499
 500     fRecEventMatrix->GetAxis(0)->SetRange();
 501     fRecEventMatrix->GetAxis(1)->SetRange();
 502     TH1D* h1RecEventMatrix_zv = (TH1D*) fRecEventMatrix->Projection(0)->Clone("h1RecEventMatrix_zv");
 503     h1RecEventMatrix_zv->Scale(1./h1RecEventMatrix_zv->Integral());
 504     //ContHist->Add(h1TriggerEventMatrix_zv);
 505
 506     TH1D* h1TriggerEff_bin0_zv   = AlidNdPtHelper::GenerateCorrMatrix(h1TriggerEventMatrix_bin0_zv,h1GenEventMatrix_bin0_zv,"h1TriggerEff_bin0_zv");
 507     TH1D* h1TriggerCorr_bin0_zv  = AlidNdPtHelper::GenerateCorrMatrix(h1GenEventMatrix_bin0_zv,h1TriggerEventMatrix_bin0_zv,"h1TriggerCorr_bin0_zv");
 508     CorrMatr->Add(h1TriggerEff_bin0_zv);
 509     CorrMatr->Add(h1TriggerCorr_bin0_zv);
 510
 511     TH1D* h1Ratio_zv = (TH1D*) h1GenEventMatrix_bin0_zv->Clone("h1Ratio_zv");
 512     h1Ratio_zv->Scale(1./h1Ratio_zv->Integral());
 513     h1Ratio_zv->Divide(h1RecEventMatrix_zv);
 514     CorrMatr->Add(h1Ratio_zv);
 515
 516
 517      // aded to correct sperately for zvert shape of triggered and untriggered bin0 events
 518      fRecEventMatrix->GetAxis(0)->SetRange();
 519      fRecEventMatrix->GetAxis(1)->SetRange();
 520      fTriggerEventMatrix->GetAxis(1)->SetRangeUser(0,0);
 521      fGenEventMatrix->GetAxis(1)->SetRangeUser(0,0);
 522
 523      TH1D* rec_zv = fRecEventMatrix->Projection(0)->Clone("rec_zv");
 524      TH1D* trig0_zv = fTriggerEventMatrix->Projection(0)->Clone("trig0_zv");
 525      TH1D* corr_shape_trig0_zv = AlidNdPtHelper::GenerateCorrMatrix(trig0_zv,rec_zv,"corr_shape_trig0_zv");
 526      CorrMatr->Add(corr_shape_trig0_zv);
 527
 528      TH1D* gen0_zv = fGenEventMatrix->Projection(0)->Clone("gen0_zv");
 529      TH1D* gen0notrig_zv = gen0_zv->Clone("gen0notrig_zv");
 530      gen0notrig_zv->Add(trig0_zv,-1.);
 531      TH1D* corr_shape_notrig0_zv = AlidNdPtHelper::GenerateCorrMatrix(gen0notrig_zv,rec_zv,"corr_shape_notrig0_zv");
 532      CorrMatr->Add(corr_shape_notrig0_zv);
 533
 534
 535    // efficienfy and correction for THnSparse
 536     THnSparse* hSRecSecTrackMatrixScaled = AlidNdPtHelper::RebinTHnSparse(fRecSecTrackMatrixScaled,hSTrackSecCorr,"hSRecSecTrackMatrixScaled");
 537     THnSparse* hSRecTrackMatrixScaled = AlidNdPtHelper::RebinTHnSparse(fRecTrackMatrixScaled,hSTrackSecCorr,"hSRecTrackMatrixScaled");
 538     THnSparse* hSRecPrimTrackMatrix = AlidNdPtHelper::RebinTHnSparse(fRecPrimTrackMatrix,hSTrackEffCorr,"hSRecPrimTrackMatrix");
 539     THnSparse* hSGenPrimTrackMatrix = AlidNdPtHelper::RebinTHnSparse(fGenPrimTrackMatrix,hSTrackEffCorr,"hSGenPrimTrackMatrix");
 540     ContHist->Add(hSRecSecTrackMatrixScaled);
 541     ContHist->Add(hSRecTrackMatrixScaled);
 542     ContHist->Add(hSRecPrimTrackMatrix);
 543     ContHist->Add(hSGenPrimTrackMatrix);
 544
 545     THnSparse* hSTrackEff  = AlidNdPtHelper::GenerateCorrMatrix(hSRecPrimTrackMatrix,hSGenPrimTrackMatrix,"hSTrackEff");
 546     THnSparse* hSTrackCorr = AlidNdPtHelper::GenerateCorrMatrix(hSGenPrimTrackMatrix,hSRecPrimTrackMatrix,"hSTrackCorr");
 547     THnSparse* hSSecCont   = AlidNdPtHelper::GenerateCorrMatrix(hSRecSecTrackMatrixScaled,hSRecTrackMatrixScaled,"hSSecCont");
 548     THnSparse* hSSecCorr   = AlidNdPtHelper::GenerateContCorrMatrix(hSRecSecTrackMatrixScaled,hSRecTrackMatrixScaled,"hSSecCorr");
 549     CorrMatr->Add(hSTrackEff);
 550     CorrMatr->Add(hSTrackCorr);
 551     CorrMatr->Add(hSSecCont);
 552     CorrMatr->Add(hSSecCorr);
 553
 554     // create th3 from thnsparse, used for corrections
 555     TH3D* h3RecPrimTrack_Rebin = hSRecPrimTrackMatrix->Projection(0,1,2)->Clone("h3RecPrimTrack_Rebin");
 556     TH3D* h3GenPrimTrack_Rebin = hSGenPrimTrackMatrix->Projection(0,1,2)->Clone("h3GenPrimTrack_Rebin");
 557     TH3D* h3RecSecTrackScaled_Rebin = hSRecSecTrackMatrixScaled->Projection(0,1,2)->Clone("h3RecSecTrackScaled_Rebin");
 558     TH3D* h3RecTrackScaled_Rebin = hSRecTrackMatrixScaled->Projection(0,1,2)->Clone("h3RecTrackScaled_Rebin");
 559
 560     TH3D* h3TrackEff  = AlidNdPtHelper::GenerateCorrMatrix(h3RecPrimTrack_Rebin,h3GenPrimTrack_Rebin,"h3TrackEff");
 561     TH3D* h3TrackCorr = AlidNdPtHelper::GenerateCorrMatrix(h3GenPrimTrack_Rebin,h3RecPrimTrack_Rebin,"h3TrackCorr");
 562     TH3D* h3SecCont = AlidNdPtHelper::GenerateCorrMatrix(h3RecSecTrackScaled_Rebin,h3RecTrackScaled_Rebin,"h3SecCont");
 563     TH3D* h3SecCorr = AlidNdPtHelper::GenerateContCorrMatrix(h3RecSecTrackScaled_Rebin,h3RecTrackScaled_Rebin,"h3SecCorr");
 564     CorrMatr->Add(h3TrackEff);
 565     CorrMatr->Add(h3TrackCorr);
 566     CorrMatr->Add(h3SecCont);
 567     CorrMatr->Add(h3SecCorr);
 568
 569     // check for trigger/vertex bias on track level
 570     THnSparse* hSRecTrackEventMatrix = (THnSparse*) obj->GetRecTrackEventMatrix()->Clone("hSRecTrackEventMatrix");
 571     THnSparse* hSGenTrackEventMatrix = (THnSparse*) obj->GetGenTrackEventMatrix()->Clone("hSGenTrackEventMatrix");
 572     THnSparse* hSTriggerTrackEventMatrix = (THnSparse*) obj->GetTriggerTrackEventMatrix()->Clone("hSTriggerTrackEventMatrix");
 573
 574     hSRecTrackEventMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
 575     hSRecTrackEventMatrix->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
 576     hSGenTrackEventMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
 577     hSGenTrackEventMatrix->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
 578     hSTriggerTrackEventMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
 579     hSTriggerTrackEventMatrix->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
 580
 581
 582     TH1D* h1TriggerBiasCorr_zv = AlidNdPtHelper::GenerateCorrMatrix(hSTriggerTrackEventMatrix->Projection(0),hSGenTrackEventMatrix->Projection(0),"h1TriggerBiasCorr_zv");
 583     TH1D* h1TriggerBiasCorr_pt = AlidNdPtHelper::GenerateCorrMatrix(hSTriggerTrackEventMatrix->Projection(1),hSGenTrackEventMatrix->Projection(1),"h1TriggerBiasCorr_pt");
 584     TH1D* h1TriggerBiasCorr_eta = AlidNdPtHelper::GenerateCorrMatrix(hSTriggerTrackEventMatrix->Projection(2),hSGenTrackEventMatrix->Projection(2),"h1TriggerBiasCorr_eta");
 585
 586     TH1D* h1VertexBiasCorr_zv  = AlidNdPtHelper::GenerateCorrMatrix(hSRecTrackEventMatrix->Projection(0),hSTriggerTrackEventMatrix->Projection(0),"h1VertexBiasCorr_zv");
 587     TH1D* h1VertexBiasCorr_pt  = AlidNdPtHelper::GenerateCorrMatrix(hSRecTrackEventMatrix->Projection(1),hSTriggerTrackEventMatrix->Projection(1),"h1VertexBiasCorr_pt");
 588     TH1D* h1VertexBiasCorr_eta  = AlidNdPtHelper::GenerateCorrMatrix(hSRecTrackEventMatrix->Projection(2),hSTriggerTrackEventMatrix->Projection(2),"h1VertexBiasCorr_eta");
 589
 590    ContHist->Add(h1TriggerBiasCorr_zv);
 591    ContHist->Add(h1TriggerBiasCorr_pt);
 592    ContHist->Add(h1TriggerBiasCorr_eta);
 593    ContHist->Add(h1VertexBiasCorr_zv);
 594    ContHist->Add(h1VertexBiasCorr_pt);
 595    ContHist->Add(h1VertexBiasCorr_eta);
 596
 597
 598
 599     // plot pictures and save to gifdir
 600     for (i=0; i < CorrMatr->LastIndex(); i++) {
 601         TCanvas* ctmp = PlotHist(CorrMatr->At(i),idstring);
 602         if (gifdir && ctmp) {
 603             TString gif(gifdir);
 604             gif += '/';
 605             gif += ctmp->GetName();
 606             gif += ".gif";
 607             ctmp->SaveAs(gif.Data(),"gif");
 608             delete ctmp;
 609         }
 610     }
 611     for (i=0; i < ContHist->LastIndex(); i++) {
 612         TCanvas* ctmp = PlotHist(ContHist->At(i),idstring);
 613         if (gifdir && ctmp) {
 614             TString gif(gifdir);
 615             gif += '/';
 616             gif += ctmp->GetName();
 617             gif += ".gif";
 618             ctmp->SaveAs(gif.Data(),"gif");
 619             delete ctmp;
 620         }
 621    }
 622
 623     // save all correction matrices and control histograms to file
 624     if (!outfile) { return; }
 625     TFile *out = TFile::Open(outfile,"RECREATE");
 626     CorrMatr->Write();
 627     ContHist->Write();
 628     out->Close();
 629
 630     return MCReconstructedEvents;
 631
 632 }
 633
 634
 635 //_____________________________________________________________________________
 636 Double_t GetStrangenessCorrFactor(Double_t pt, Double_t s)
 637 {
 638     // data driven correction factor for secondaries (PbPb)
 639
 640     if (pt <= 0.17) return 1.0;
 641     if (pt <= 0.4) return GetLinearInterpolationValue(0.17,1.0,0.4,1.07, pt);
 642     if (pt <= 0.6) return GetLinearInterpolationValue(0.4,1.07,0.6,1.25, pt);
 643     if (pt <= 1.2) return GetLinearInterpolationValue(0.6,1.25,1.2,1.5,  pt);
 644     return 1.5;
 645 }
 646
 647
 648 //___________________________________________________________________________
 649 Double_t GetLinearInterpolationValue(const Double_t x1,const  Double_t y1,const  Double_t x2,const  Double_t y2, const Double_t pt)
 650 {
 651     //
 652     // linear interpolation
 653     //
 654     return ((y2-y1)/(x2-x1))*pt+(y2-(((y2-y1)/(x2-x1))*x2));
 655 }
 656
 657 //___________________________________________________________________________
 658 TCanvas* PlotHist(TObject* hobj, const char* label=0)
 659 {
 660     TH1* h = dynamic_cast<TH1*>(hobj);
 661     if (!h) return 0;
 662     if (h->GetDimension() > 2) return 0;
 663     h->SetStats(0);
 664     if ( TString(h->GetName()).Contains("Events")) { h->SetStats(1); }
 665     TString t(label);
 666     if (label) t += "_";
 667     t += h->GetName();
 668     h->SetTitle(t.Data());
 669     TCanvas* c = new TCanvas(t.Data(),t.Data());
 670     if (h->GetDimension() >= 1) {
 671         TString xlabel(h->GetXaxis()->GetTitle());
 672         if (xlabel.Contains("Pt")) { c->SetLogx();  h->GetXaxis()->SetRangeUser(0.1 , 100.); }
 673     }
 674     if (h->GetDimension() == 1) {
 675         if (xlabel.Contains("p_{T}")) { c->SetLogx();  c->SetLogy();  h->GetXaxis()->SetRangeUser(0.1 , 100.); }
 676     }
 677     if (h->GetDimension() == 2) {
 678         TString ylabel(h->GetYaxis()->GetTitle());
 679         if (ylabel.Contains("Pt")) { c->SetLogy(); h->GetYaxis()->SetRangeUser(0.1 , 100.); }
 680         if (ylabel.Contains("p_{T}")) { c->SetLogy(); h->GetYaxis()->SetRangeUser(0.1 , 100.); }
 681         h->Draw("COLZ");
 682     }
 683     if (h->GetDimension() == 1) {
 684         h->Draw();
 685     }
 686     return c;
 687
 688 }
 689
 690 Int_t CheckLoadLibrary(const char* library)
 691 {
 692   // checks if a library is already loaded, if not loads the library
 693
 694   if (strlen(gSystem->GetLibraries(Form("%s.so", library), "", kFALSE)) > 0)
 695     return 1;
 696
 697   return gSystem->Load(library);
 698 }