PWGUD/dNdPt/macros/GenerateCorrMatr.C

   1 Double_t GenerateCorrMatr(const char* mcfile, const char* mctask, const char* idstring ,const char* outfile, const char* gifdir = 0)
   2 //Double_t GenerateCorrMatr()
   3 {
   4
   5 //tmp setting
   6 //   const char* mcfile = "/lustre/alice/train/V005.PbPb_MC/2010-11-29_0108.4251/mergedPeriods/MC/PbPb/LHC10h8/mknichel_dNdPtPbPb_gt_v0_c0_syst4.root";
   7 //     const char* mctask = "mknichel_dNdPtPbPb_gt_v0_c0_syst4";
   8 //     const char* idstring = "gt_v0_c0_syst4";
   9 //     const char* outfile = "/u/mknichel/alice/dNdPt/2010-11-29_0318/corrMatr_gt_v0_c0_syst4.root";
  10 //     const char* gifdir = "/u/mknichel/alice/dNdPt/2010-11-29";
  11
  12     // settings vor zVertex cut (event and track level)
  13     Double_t zVert = 10.0;
  14
  15     // setting on eta cut (track level)
  16     Double_t eta = 0.8;
  17
  18     // strangeness scaling factor (for secondaries from strange decays)
  19     Double_t sscale = 1.4;
  20
  21     //load required libraries
  22     //load required libraries
  23     gSystem->AddIncludePath("-I$ROOTSYS/include -I$ALICE_ROOT/ -I$ALICE_ROOT/include -I$ALICE_ROOT/STEER  -I$ALICE_ROOT/ANALYSIS -I$ALICE_ROOT/PWG0 -I$ALICE_ROOT/PWGPP -I$ALICE_ROOT/PWG2 -I$ALICE_ROOT/PWG3 -I$ALICE_ROOT/PWG3/vertexingHF -I$ALICE_ROOT/PWG4 -I$ALICE_ROOT/CORRFW -I$ALICE_ROOT/TPC -I$ALICE_ROOT/TRD -I$ALICE_ROOT/PWG3/muon -I$ALICE_ROOT/JETAN -I$ALICE_ROOT/ANALYSIS/Tender");
  24
  25   gSystem->Load("libCore");
  26   gSystem->Load("libPhysics");
  27   gSystem->Load("libMinuit");
  28   gSystem->Load("libGui");
  29   gSystem->Load("libXMLParser");
  30
  31   gSystem->Load("libGeom");
  32   gSystem->Load("libVMC");
  33
  34   gSystem->Load("libNet");
  35
  36
  37   gSystem->Load("libSTEERBase");
  38   gSystem->Load("libESD");
  39   gSystem->Load("libCDB");
  40   gSystem->Load("libRAWDatabase");
  41   gSystem->Load("libRAWDatarec");
  42   gSystem->Load("libANALYSIS");
  43
  44
  45
  46     gSystem->Load("libANALYSIS.so");
  47     gSystem->Load("libANALYSISalice.so");
  48     gSystem->Load("libTENDER.so");
  49     gSystem->Load("libCORRFW.so");
  50     gSystem->Load("libPWG0base.so");
  51     gSystem->Load("libPWG0dep");
  52     gSystem->Load("libPWG0selectors.so");
  53
  54
  55
  56
  57
  58     // make plots nicer
  59     gROOT->SetStyle("Plain");
  60     gStyle->SetPalette(1);
  61
  62     // array for all correction matrices
  63     TObjArray* CorrMatr = new TObjArray();
  64
  65     // array for all control histograms
  66     TObjArray* ContHist = new TObjArray();
  67
  68
  69     // load mc information
  70     TFile* fmc = TFile::Open(mcfile,"READ");
  71     if (!fmc) return -1;
  72     TList* lmc = dynamic_cast<TList*>(fmc->Get(mctask));
  73     if (!lmc) return -1;
  74     AlidNdPtAnalysisPbPb *obj = dynamic_cast<AlidNdPtAnalysisPbPb*>(lmc->FindObject("dNdPtAnalysisPbPb"));
  75     if (!obj) return -1;
  76
  77     //Event statistics
  78     THnSparse *fRecEventMatrix = obj->GetRecEventMatrix(); //all reconstructed events
  79     TH2D* h2RecEventAll = fRecEventMatrix->Projection(0,1)->Clone("h2RecEventAll");
  80     fRecEventMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
  81     TH2D* h2RecEvent = fRecEventMatrix->Projection(0,1)->Clone("h2RecEvent");
  82     Double_t MCReconstructedEvents = h2RecEvent->Integral();
  83     Double_t MCReconstructedEventsAll = h2RecEventAll->Integral();
  84     ContHist->Add(h2RecEvent);
  85     ContHist->Add(h2RecEventAll);
  86
  87     THnSparse *fTriggerEventMatrix = obj->GetTriggerEventMatrix(); //all triggered events
  88     TH2D* h2TriggerEventAll = fTriggerEventMatrix->Projection(0,1)->Clone("h2TriggerEventAll");
  89     fTriggerEventMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
  90     TH2D* h2TriggerEvent = fTriggerEventMatrix->Projection(0,1)->Clone("h2TriggerEvent");
  91     Double_t MCTriggeredEvents = h2TriggerEvent->Integral();
  92     Double_t MCTriggeredEventsAll = h2TriggerEventAll->Integral();
  93     ContHist->Add(h2TriggerEvent);
  94     ContHist->Add(h2TriggerEventAll);
  95
  96     THnSparse *fGenEventMatrix = obj->GetGenEventMatrix(); //all generated events
  97     TH2D* h2GenEventAll = fGenEventMatrix->Projection(0,1)->Clone("h2GenEventAll");
  98     fGenEventMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
  99     TH2D* h2GenEvent = fGenEventMatrix->Projection(0,1)->Clone("h2GenEvent");
 100     Double_t MCGeneratedEvents = h2GenEvent->Integral();
 101     Double_t MCGeneratedEventsAll = h2GenEventAll->Integral();
 102     ContHist->Add(h2RecEvent);
 103     ContHist->Add(h2RecEventAll);
 104
 105     printf("=== generated MC events for correction matrices    %lf ===\n",MCGeneratedEvents);
 106     printf("=== triggered MC events for correction matrices    %lf ===\n",MCTriggeredEvents);
 107     printf("=== recontructed MC events for correction matrices %lf ===\n",MCReconstructedEvents);
 108     printf("\n");
 109     printf("=== cut on the zVertex +- %lf ===\n",zVert);
 110     printf("=== generated MC events (before zVertex cut)       %lf ===\n",MCGeneratedEventsAll);
 111     printf("=== triggered MC events (before zVertex cut)       %lf ===\n",MCTriggeredEventsAll);
 112     printf("=== recontructed MC events (before zVertex cut)    %lf ===\n",MCReconstructedEventsAll);
 113
 114     TH1D* h1MCGeneratedEvents = new TH1D("h1MCGeneratedEvents","h1MCGeneratedEvents",1,0,1);
 115     TH1D* h1MCTriggeredEvents = new TH1D("h1MCTriggeredEvents","h1MCTriggeredEvents",1,0,1);
 116     TH1D* h1MCReconstructedEvents = new TH1D("h1MCReconstructedEvents","h1MCReconstructedEvents",1,0,1);
 117     TH1D* h1MCGeneratedEventsAll = new TH1D("h1MCGeneratedEventsAll","h1MCGeneratedEventsAll",1,0,1);
 118     TH1D* h1MCTriggeredEventsAll = new TH1D("h1MCTriggeredEventsAll","h1MCTriggeredEventsAll",1,0,1);
 119     TH1D* h1MCReconstructedEventsAll = new TH1D("h1MCReconstructedEventsAll","h1MCReconstructedEventsAll",1,0,1);
 120
 121     h1MCGeneratedEvents->Fill(0,MCGeneratedEvents);
 122     h1MCGeneratedEvents->SetEntries(MCGeneratedEvents);
 123     h1MCTriggeredEvents->Fill(0,MCTriggeredEvents);
 124     h1MCTriggeredEvents->SetEntries(MCTriggeredEvents);
 125     h1MCReconstructedEvents->Fill(0,MCReconstructedEvents);
 126     h1MCReconstructedEvents->SetEntries(MCReconstructedEvents);
 127     h1MCGeneratedEventsAll->Fill(0,MCGeneratedEventsAll);
 128     h1MCGeneratedEventsAll->SetEntries(MCGeneratedEventsAll);
 129     h1MCTriggeredEventsAll->Fill(0,MCTriggeredEventsAll);
 130     h1MCTriggeredEventsAll->SetEntries(MCTriggeredEventsAll);
 131     h1MCReconstructedEventsAll->Fill(0,MCReconstructedEventsAll);
 132     h1MCReconstructedEventsAll->SetEntries(MCReconstructedEventsAll);
 133
 134     ContHist->Add(h1MCGeneratedEvents);
 135     ContHist->Add(h1MCTriggeredEvents);
 136     ContHist->Add(h1MCReconstructedEvents);
 137     ContHist->Add(h1MCGeneratedEventsAll);
 138     ContHist->Add(h1MCTriggeredEventsAll);
 139     ContHist->Add(h1MCReconstructedEventsAll);
 140
 141     // efficienfy and correction matrices for tigger and vertex efficiency
 142     TH2D* h2EventTriggerEffAll  = AlidNdPtHelper::GenerateCorrMatrix(h2TriggerEventAll,h2GenEventAll,"h2EventTriggerEffAll");
 143     TH2D* h2EventTriggerCorrAll = AlidNdPtHelper::GenerateCorrMatrix(h2GenEventAll,h2TriggerEventAll,"h2EventTriggerCorrAll");
 144     TH2D* h2EventTriggerEff  = AlidNdPtHelper::GenerateCorrMatrix(h2TriggerEvent,h2GenEvent,"h2EventTriggerEff");
 145     TH2D* h2EventTriggerCorr = AlidNdPtHelper::GenerateCorrMatrix(h2GenEvent,h2TriggerEvent,"h2EventTriggerCorr");
 146
 147     TH2D* h2EventRecEffAll  = AlidNdPtHelper::GenerateCorrMatrix(h2RecEventAll,h2TriggerEventAll,"h2EventRecEffAll");
 148     TH2D* h2EventRecCorrAll = AlidNdPtHelper::GenerateCorrMatrix(h2TriggerEventAll,h2RecEventAll,"h2EventRecCorrAll");
 149     TH2D* h2EventRecEff  = AlidNdPtHelper::GenerateCorrMatrix(h2RecEvent,h2TriggerEvent,"h2EventRecEff");
 150     TH2D* h2EventRecCorr = AlidNdPtHelper::GenerateCorrMatrix(h2TriggerEvent,h2RecEvent,"h2EventRecCorr");
 151
 152     TH2D* h2EventEffAll  = AlidNdPtHelper::GenerateCorrMatrix(h2RecEventAll,h2GenEventAll,"h2EventEffAll");
 153     TH2D* h2EventCorrAll = AlidNdPtHelper::GenerateCorrMatrix(h2GenEventAll,h2RecEventAll,"h2EventCorrAll");
 154     TH2D* h2EventEff  = AlidNdPtHelper::GenerateCorrMatrix(h2RecEvent,h2GenEvent,"h2EventEff");
 155     TH2D* h2EventCorr = AlidNdPtHelper::GenerateCorrMatrix(h2GenEvent,h2RecEvent,"h2EventCorr");
 156
 157     CorrMatr->Add(h2EventTriggerEffAll);
 158     CorrMatr->Add(h2EventTriggerCorrAll);
 159     CorrMatr->Add(h2EventTriggerEff);
 160     CorrMatr->Add(h2EventTriggerCorr);
 161     CorrMatr->Add(h2EventRecEffAll);
 162     CorrMatr->Add(h2EventRecCorrAll);
 163     CorrMatr->Add(h2EventRecEff);
 164     CorrMatr->Add(h2EventRecCorr);
 165     CorrMatr->Add(h2EventEffAll);
 166     CorrMatr->Add(h2EventCorrAll);
 167     CorrMatr->Add(h2EventEff);
 168     CorrMatr->Add(h2EventCorr);
 169
 170
 171
 172     // all recontructed
 173     THnSparse *fRecTrackMatrix = obj->GetRecTrackMatrix();
 174     fRecTrackMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
 175     fRecTrackMatrix->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
 176     TH3D* h3RecTrack = fRecTrackMatrix->Projection(0,1,2)->Clone("h3RecTrack");
 177     TH2D* h2RecTrack_zv_pt  = fRecTrackMatrix->Projection(0,1)->Clone("h2RecTrack_zv_pt");
 178     TH2D* h2RecTrack_zv_eta = fRecTrackMatrix->Projection(0,2)->Clone("h2RecTrack_zv_eta");
 179     TH2D* h2RecTrack_pt_eta = fRecTrackMatrix->Projection(1,2)->Clone("h2RecTrack_pt_eta");
 180     TH1D* h1RecTrack_zv  = fRecTrackMatrix->Projection(0)->Clone("h1RecTrack_zv");
 181     TH1D* h1RecTrack_pt  = fRecTrackMatrix->Projection(1)->Clone("h1RecTrack_pt");
 182     TH1D* h1RecTrack_eta = fRecTrackMatrix->Projection(2)->Clone("h1RecTrack_eta");
 183     Double_t MCReconstructedTracks = h3RecTrack->Integral();
 184
 185     ContHist->Add(h3RecTrack);
 186     ContHist->Add(h2RecTrack_zv_pt);
 187     ContHist->Add(h2RecTrack_zv_eta);
 188     ContHist->Add(h2RecTrack_pt_eta);
 189     ContHist->Add(h1RecTrack_zv);
 190     ContHist->Add(h1RecTrack_pt);
 191     ContHist->Add(h1RecTrack_eta);
 192
 193      // recontructed primary tracks
 194     THnSparse *fRecPrimTrackMatrix = obj->GetRecPrimTrackMatrix();
 195     THnSparse *fRecTrackMatrixScaled = fRecPrimTrackMatrix->Clone("fRecTrackMatrixScaled"); //used later for secondaries scaling
 196     fRecPrimTrackMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
 197     fRecPrimTrackMatrix->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
 198     TH3D* h3RecPrimTrack = fRecPrimTrackMatrix->Projection(0,1,2)->Clone("h3RecPrimTrack");
 199     TH2D* h2RecPrimTrack_zv_pt  = fRecPrimTrackMatrix->Projection(0,1)->Clone("h2RecPrimTrack_zv_pt");
 200     TH2D* h2RecPrimTrack_zv_eta = fRecPrimTrackMatrix->Projection(0,2)->Clone("h2RecPrimTrack_zv_eta");
 201     TH2D* h2RecPrimTrack_pt_eta = fRecPrimTrackMatrix->Projection(1,2)->Clone("h2RecPrimTrack_pt_eta");
 202     TH1D* h1RecPrimTrack_zv  = fRecPrimTrackMatrix->Projection(0)->Clone("h1RecPrimTrack_zv");
 203     TH1D* h1RecPrimTrack_pt  = fRecPrimTrackMatrix->Projection(1)->Clone("h1RecPrimTrack_pt");
 204     TH1D* h1RecPrimTrack_eta = fRecPrimTrackMatrix->Projection(2)->Clone("h1RecPrimTrack_eta");
 205     Double_t MCReconstructedPrimTracks = h3RecPrimTrack->Integral();
 206
 207     ContHist->Add(h3RecPrimTrack);
 208     ContHist->Add(h2RecPrimTrack_zv_pt);
 209     ContHist->Add(h2RecPrimTrack_zv_eta);
 210     ContHist->Add(h2RecPrimTrack_pt_eta);
 211     ContHist->Add(h1RecPrimTrack_zv);
 212     ContHist->Add(h1RecPrimTrack_pt);
 213     ContHist->Add(h1RecPrimTrack_eta);
 214
 215      // recontructed secondary tracks
 216     THnSparse *fRecSecTrackMatrix = obj->GetRecSecTrackMatrix();
 217     THnSparse *fRecSecTrackMatrixScaled = fRecSecTrackMatrix->Clone("fRecSecTrackMatrixScaled"); //used later for secondaries scaling
 218     fRecSecTrackMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
 219     fRecSecTrackMatrix->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
 220     TH3D* h3RecSecTrack = fRecSecTrackMatrix->Projection(0,1,2)->Clone("h3RecSecTrack");
 221     TH2D* h2RecSecTrack_zv_pt  = fRecSecTrackMatrix->Projection(0,1)->Clone("h2RecSecTrack_zv_pt");
 222     TH2D* h2RecSecTrack_zv_eta = fRecSecTrackMatrix->Projection(0,2)->Clone("h2RecSecTrack_zv_eta");
 223     TH2D* h2RecSecTrack_pt_eta = fRecSecTrackMatrix->Projection(1,2)->Clone("h2RecSecTrack_pt_eta");
 224     TH1D* h1RecSecTrack_zv  = fRecSecTrackMatrix->Projection(0)->Clone("h1RecSecTrack_zv");
 225     TH1D* h1RecSecTrack_pt  = fRecSecTrackMatrix->Projection(1)->Clone("h1RecSecTrack_pt");
 226     TH1D* h1RecSecTrack_eta = fRecSecTrackMatrix->Projection(2)->Clone("h1RecSecTrack_eta");
 227     Double_t MCReconstructedSecTracks = h3RecSecTrack->Integral();
 228
 229     ContHist->Add(h3RecSecTrack);
 230     ContHist->Add(h2RecSecTrack_zv_pt);
 231     ContHist->Add(h2RecSecTrack_zv_eta);
 232     ContHist->Add(h2RecSecTrack_pt_eta);
 233     ContHist->Add(h1RecSecTrack_zv);
 234     ContHist->Add(h1RecSecTrack_pt);
 235     ContHist->Add(h1RecSecTrack_eta);
 236
 237     // generated primary tracks
 238     THnSparse *fGenPrimTrackMatrix = obj->GetGenPrimTrackMatrix();
 239     fGenPrimTrackMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
 240     fGenPrimTrackMatrix->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
 241     TH3D* h3GenPrimTrack = fGenPrimTrackMatrix->Projection(0,1,2)->Clone("h3GenPrimTrack");
 242     TH2D* h2GenPrimTrack_zv_Pt  = fGenPrimTrackMatrix->Projection(0,1)->Clone("h2GenPrimTrack_zv_pt");
 243     TH2D* h2GenPrimTrack_zv_eta = fGenPrimTrackMatrix->Projection(0,2)->Clone("h2GenPrimTrack_zv_eta");
 244     TH2D* h2GenPrimTrack_pt_eta = fGenPrimTrackMatrix->Projection(1,2)->Clone("h2GenPrimTrack_pt_eta");
 245     TH1D* h1GenPrimTrack_zv  = fGenPrimTrackMatrix->Projection(0)->Clone("h1GenPrimTrack_zv");
 246     TH1D* h1GenPrimTrack_pt  = fGenPrimTrackMatrix->Projection(1)->Clone("h1GenPrimTrack_pt");
 247     TH1D* h1GenPrimTrack_eta = fGenPrimTrackMatrix->Projection(2)->Clone("h1GenPrimTrack_eta");
 248     Double_t MCGeneratedPrimTracks = h3GenPrimTrack->Integral();
 249
 250     ContHist->Add(h3GenPrimTrack);
 251     ContHist->Add(h2GenPrimTrack_zv_pt);
 252     ContHist->Add(h2GenPrimTrack_zv_eta);
 253     ContHist->Add(h2GenPrimTrack_pt_eta);
 254     ContHist->Add(h1GenPrimTrack_zv);
 255     ContHist->Add(h1GenPrimTrack_pt);
 256     ContHist->Add(h1GenPrimTrack_eta);
 257     printf("\n");
 258     printf("==============================================================\n");
 259     printf("=== recontructed MC tracks              %lf ===\n",MCReconstructedTracks);
 260     printf("=== recontructed MC secondary tracks    %lf ===\n",MCReconstructedSecTracks);
 261     printf("=== recontructed MC primary tracks      %lf ===\n",MCReconstructedPrimTracks);
 262     printf("=== generated MC primary track          %lf ===\n",MCGeneratedPrimTracks);
 263     printf("==============================================================\n");
 264     printf("\n");
 265
 266 //    THnSparse *fSparseTriggerTrackEvent = obj->GetTriggerTrackEventMatrix();//Tracks from triggered events
 267 //    THnSparse *fSparseVtxTrackEvent = obj->GetRecTrackEventMatrix();//Tracks from events with rec. vtx
 268 //    THnSparse *fSparseGenTrackEvent = obj->GetGenTrackEventMatrix();//generated TrackEvent matrix
 269
 270     // tracking efficiencies + corrections
 271    TH3D* h3TrackEff  = AlidNdPtHelper::GenerateCorrMatrix(h3RecPrimTrack,h3GenPrimTrack,"h3TrackEff");
 272    TH3D* h3TrackCorr = AlidNdPtHelper::GenerateCorrMatrix(h3GenPrimTrack,h3RecPrimTrack,"h3TrackCorr");
 273
 274    TH2D* h2TrackEff_zv_pt   = AlidNdPtHelper::GenerateCorrMatrix(h2RecPrimTrack_zv_pt,h2GenPrimTrack_zv_pt,"h2TrackEff_zv_pt");
 275    TH2D* h2TrackCorr_zv_pt  = AlidNdPtHelper::GenerateCorrMatrix(h2GenPrimTrack_zv_pt,h2RecPrimTrack_zv_pt,"h2TrackCorr_zv_pt");
 276    TH2D* h2TrackEff_zv_eta  = AlidNdPtHelper::GenerateCorrMatrix(h2RecPrimTrack_zv_eta,h2GenPrimTrack_zv_eta,"h2TrackEff_zv_eta");
 277    TH2D* h2TrackCorr_zv_eta = AlidNdPtHelper::GenerateCorrMatrix(h2GenPrimTrack_zv_eta,h2RecPrimTrack_zv_eta,"h2TrackCorr_zv_eta");
 278    TH2D* h2TrackEff_pt_eta  = AlidNdPtHelper::GenerateCorrMatrix(h2RecPrimTrack_pt_eta,h2GenPrimTrack_pt_eta,"h2TrackEff_pt_eta");
 279    TH2D* h2TrackCorr_pt_eta = AlidNdPtHelper::GenerateCorrMatrix(h2GenPrimTrack_pt_eta,h2RecPrimTrack_pt_eta,"h2TrackCorr_pt_eta");
 280
 281    TH1D* h1TrackEff_zv   = AlidNdPtHelper::GenerateCorrMatrix(h1RecPrimTrack_zv,h1GenPrimTrack_zv,"h1TrackEff_zv");
 282    TH1D* h1TrackCorr_zv  = AlidNdPtHelper::GenerateCorrMatrix(h1GenPrimTrack_zv,h1RecPrimTrack_zv,"h1TrackCorr_zv");
 283    TH1D* h1TrackEff_pt   = AlidNdPtHelper::GenerateCorrMatrix(h1RecPrimTrack_pt,h1GenPrimTrack_pt,"h1TrackEff_pt");
 284    TH1D* h1TrackCorr_pt  = AlidNdPtHelper::GenerateCorrMatrix(h1GenPrimTrack_pt,h1RecPrimTrack_pt,"h1TrackCorr_pt");
 285    TH1D* h1TrackEff_eta  = AlidNdPtHelper::GenerateCorrMatrix(h1RecPrimTrack_eta,h1GenPrimTrack_eta,"h1TrackEff_eta");
 286    TH1D* h1TrackCorr_eta = AlidNdPtHelper::GenerateCorrMatrix(h1GenPrimTrack_eta,h1RecPrimTrack_eta,"h1TrackCorr_eta");
 287
 288    CorrMatr->Add(h3TrackEff);
 289    CorrMatr->Add(h3TrackCorr);
 290    CorrMatr->Add(h2TrackEff_zv_pt);
 291    CorrMatr->Add(h2TrackCorr_zv_pt);
 292    CorrMatr->Add(h2TrackEff_zv_eta);
 293    CorrMatr->Add(h2TrackCorr_zv_eta);
 294    CorrMatr->Add(h2TrackEff_pt_eta);
 295    CorrMatr->Add(h2TrackCorr_pt_eta);
 296    CorrMatr->Add(h1TrackEff_zv);
 297    CorrMatr->Add(h1TrackCorr_zv);
 298    CorrMatr->Add(h1TrackEff_pt);
 299    CorrMatr->Add(h1TrackCorr_pt);
 300    CorrMatr->Add(h1TrackEff_eta);
 301    CorrMatr->Add(h1TrackCorr_eta);
 302
 303    // scale the secondaries before calculating correction matrices
 304    for (Long64_t i = 0; i < fRecSecTrackMatrixScaled->GetNbins(); i++) {
 305        Int_t c[3];
 306        Double_t val = fRecSecTrackMatrixScaled->GetBinContent(i,c);
 307        Double_t err = fRecSecTrackMatrixScaled->GetBinError(i);
 308        Double_t pt =  fRecSecTrackMatrixScaled->GetAxis(1)->GetBinCenter(c[1]);
 309        Double_t scale = GetStrangenessCorrFactorPbPb(pt, sscale);
 310        fRecSecTrackMatrixScaled->SetBinContent(c,val*scale);
 311        fRecSecTrackMatrixScaled->SetBinError(c,err*scale);
 312     }
 313
 314     // for correct determination of secondaries contamination, also the total total tracks have to be scaled
 315     // this is done by taking primaries and adding the scaled secondaries
 316     fRecTrackMatrixScaled->Add(fRecSecTrackMatrixScaled);
 317
 318     fRecSecTrackMatrixScaled->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
 319     fRecSecTrackMatrixScaled->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
 320
 321     TH3D* h3RecSecTrackScaled = fRecSecTrackMatrixScaled->Projection(0,1,2)->Clone("h3RecSecTrackScaled");
 322     TH2D* h2RecSecTrackScaled_zv_pt  = fRecSecTrackMatrixScaled->Projection(0,1)->Clone("h2RecSecTrackScaled_zv_pt");
 323     TH2D* h2RecSecTrackScaled_zv_eta = fRecSecTrackMatrixScaled->Projection(0,2)->Clone("h2RecSecTrackScaled_zv_eta");
 324     TH2D* h2RecSecTrackScaled_pt_eta = fRecSecTrackMatrixScaled->Projection(1,2)->Clone("h2RecSecTrackScaled_pt_eta");
 325     TH1D* h1RecSecTrackScaled_zv  = fRecSecTrackMatrixScaled->Projection(0)->Clone("h1RecSecTrackScaled_zv");
 326     TH1D* h1RecSecTrackScaled_pt  = fRecSecTrackMatrixScaled->Projection(1)->Clone("h1RecSecTrackScaled_pt");
 327     TH1D* h1RecSecTrackScaled_eta = fRecSecTrackMatrixScaled->Projection(2)->Clone("h1RecSecTrackScaled_eta");
 328
 329     ContHist->Add(h3RecSecTrackScaled);
 330     ContHist->Add(h2RecSecTrackScaled_zv_pt);
 331     ContHist->Add(h2RecSecTrackScaled_zv_eta);
 332     ContHist->Add(h2RecSecTrackScaled_pt_eta);
 333     ContHist->Add(h1RecSecTrackScaled_zv);
 334     ContHist->Add(h1RecSecTrackScaled_pt);
 335     ContHist->Add(h1RecSecTrackScaled_eta);
 336
 337     fRecTrackMatrixScaled->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
 338     fRecTrackMatrixScaled->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
 339
 340     TH3D* h3RecTrackScaled = fRecTrackMatrixScaled->Projection(0,1,2)->Clone("h3RecTrackScaled");
 341     TH2D* h2RecTrackScaled_zv_pt  = fRecTrackMatrixScaled->Projection(0,1)->Clone("h2RecTrackScaled_zv_pt");
 342     TH2D* h2RecTrackScaled_zv_eta = fRecTrackMatrixScaled->Projection(0,2)->Clone("h2RecTrackScaled_zv_eta");
 343     TH2D* h2RecTrackScaled_pt_eta = fRecTrackMatrixScaled->Projection(1,2)->Clone("h2RecTrackScaled_pt_eta");
 344     TH1D* h1RecTrackScaled_zv  = fRecTrackMatrixScaled->Projection(0)->Clone("h1RecTrackScaled_zv");
 345     TH1D* h1RecTrackScaled_pt  = fRecTrackMatrixScaled->Projection(1)->Clone("h1RecTrackScaled_pt");
 346     TH1D* h1RecTrackScaled_eta = fRecTrackMatrixScaled->Projection(2)->Clone("h1RecTrackScaled_eta");
 347
 348     ContHist->Add(h3RecTrackScaled);
 349     ContHist->Add(h2RecTrackScaled_zv_pt);
 350     ContHist->Add(h2RecTrackScaled_zv_eta);
 351     ContHist->Add(h2RecTrackScaled_pt_eta);
 352     ContHist->Add(h1RecTrackScaled_zv);
 353     ContHist->Add(h1RecTrackScaled_pt);
 354     ContHist->Add(h1RecTrackScaled_eta);
 355
 356     // create histograms for secondaries contamination and correction
 357
 358     TH3D* h3SecCont = AlidNdPtHelper::GenerateCorrMatrix(h3RecSecTrackScaled,h3RecTrackScaled,"h3SecCont");
 359     TH3D* h3SecCorr = AlidNdPtHelper::GenerateContCorrMatrix(h3RecSecTrackScaled,h3RecTrackScaled,"h3SecCorr");
 360     TH2D* h2SecCont_zv_pt  = AlidNdPtHelper::GenerateCorrMatrix(h2RecSecTrackScaled_zv_pt,h2RecTrackScaled_zv_pt,"h2SecCont_zv_pt");
 361     TH2D* h2SecCorr_zv_pt  = AlidNdPtHelper::GenerateContCorrMatrix(h2RecSecTrackScaled_zv_pt,h2RecTrackScaled_zv_pt,"h2SecCorr_zv_pt");
 362     TH2D* h2SecCont_zv_eta = AlidNdPtHelper::GenerateCorrMatrix(h2RecSecTrackScaled_zv_eta,h2RecTrackScaled_zv_eta,"h2SecCont_zv_eta");
 363     TH2D* h2SecCorr_zv_eta = AlidNdPtHelper::GenerateContCorrMatrix(h2RecSecTrackScaled_zv_eta,h2RecTrackScaled_zv_eta,"h2SecCorr_zv_eta");
 364     TH2D* h2SecCont_pt_eta = AlidNdPtHelper::GenerateCorrMatrix(h2RecSecTrackScaled_pt_eta,h2RecTrackScaled_pt_eta,"h2SecCont_pt_eta");
 365     TH2D* h2SecCorr_pt_eta = AlidNdPtHelper::GenerateContCorrMatrix(h2RecSecTrackScaled_pt_eta,h2RecTrackScaled_pt_eta,"h2SecCorr_pt_eta");
 366     TH1D* h1SecCont_zv = AlidNdPtHelper::GenerateCorrMatrix(h1RecSecTrackScaled_zv,h1RecTrackScaled_zv,"h1SecCont_zv");
 367     TH1D* h1SecCorr_zv = AlidNdPtHelper::GenerateContCorrMatrix(h1RecSecTrackScaled_zv,h1RecTrackScaled_zv,"h1SecCorr_zv");
 368     TH1D* h1SecCont_pt = AlidNdPtHelper::GenerateCorrMatrix(h1RecSecTrackScaled_pt,h1RecTrackScaled_pt,"h1SecCont_pt");
 369     TH1D* h1SecCorr_pt = AlidNdPtHelper::GenerateContCorrMatrix(h1RecSecTrackScaled_pt,h1RecTrackScaled_pt,"h1SecCorr_pt");
 370     TH1D* h1SecCont_eta = AlidNdPtHelper::GenerateCorrMatrix(h1RecSecTrackScaled_eta,h1RecTrackScaled_eta,"h1SecCont_eta");
 371     TH1D* h1SecCorr_eta = AlidNdPtHelper::GenerateContCorrMatrix(h1RecSecTrackScaled_eta,h1RecTrackScaled_eta,"h1SecCorr_eta");
 372
 373     CorrMatr->Add(h3SecCont);
 374     CorrMatr->Add(h3SecCorr);
 375     CorrMatr->Add(h2SecCont_zv_pt);
 376     CorrMatr->Add(h2SecCorr_zv_pt);
 377     CorrMatr->Add(h2SecCont_zv_eta);
 378     CorrMatr->Add(h2SecCorr_zv_eta);
 379     CorrMatr->Add(h2SecCont_pt_eta);
 380     CorrMatr->Add(h2SecCorr_pt_eta);
 381     CorrMatr->Add(h1SecCont_zv);
 382     CorrMatr->Add(h1SecCorr_zv);
 383     CorrMatr->Add(h1SecCont_pt);
 384     CorrMatr->Add(h1SecCorr_pt);
 385     CorrMatr->Add(h1SecCont_eta);
 386     CorrMatr->Add(h1SecCorr_eta);
 387
 388     // plot pictures and save to gifdir
 389     for (i=0; i < CorrMatr->LastIndex(); i++) {
 390         TCanvas* ctmp = PlotHist(CorrMatr->At(i),idstring);
 391         if (gifdir && ctmp) {
 392             TString gif(gifdir);
 393             gif += '/';
 394             gif += ctmp->GetName();
 395             gif += ".gif";
 396             ctmp->SaveAs(gif.Data(),"gif");
 397             delete ctmp;
 398         }
 399     }
 400     for (i=0; i < ContHist->LastIndex(); i++) {
 401         TCanvas* ctmp = PlotHist(ContHist->At(i),idstring);
 402         if (gifdir && ctmp) {
 403             TString gif(gifdir);
 404             gif += '/';
 405             gif += ctmp->GetName();
 406             gif += ".gif";
 407             ctmp->SaveAs(gif.Data(),"gif");
 408             delete ctmp;
 409         }
 410    }
 411
 412     // save all correction matrices and control histograms to file
 413     if (!outfile) { return; }
 414     TFile *out = TFile::Open(outfile,"RECREATE");
 415     CorrMatr->Write();
 416     ContHist->Write();
 417     out->Close();
 418
 419     return MCReconstructedEvents;
 420
 421 }
 422
 423
 424 //_____________________________________________________________________________
 425 Double_t GetStrangenessCorrFactorPbPb(Double_t pt, Double_t s)
 426 {
 427     // data driven correction factor for secondaries (PbPb)
 428
 429     if (pt <= 0.25) return 1.0;
 430     if (pt <= 0.5) return GetLinearInterpolationValue(0.25,1.0,0.5,0.60+0.40*s, pt);
 431     if (pt <= 1.0) return GetLinearInterpolationValue(0.5,0.60+0.40*s,1.0,0.53+0.47*s, pt);
 432     if (pt <= 2.0) return GetLinearInterpolationValue(1.0,0.53+0.47*s,2.0,0.44+0.56*s, pt);
 433     if (pt <= 5.0) return GetLinearInterpolationValue(2.0,0.44+0.56*s,5.0,0.33+0.67*s, pt);
 434     return 0.33+0.67*s;
 435 }
 436
 437
 438 //___________________________________________________________________________
 439 Double_t GetLinearInterpolationValue(const Double_t x1,const  Double_t y1,const  Double_t x2,const  Double_t y2, const Double_t pt)
 440 {
 441     //
 442     // linear interpolation
 443     //
 444     return ((y2-y1)/(x2-x1))*pt+(y2-(((y2-y1)/(x2-x1))*x2));
 445 }
 446
 447 //___________________________________________________________________________
 448 TCanvas* PlotHist(TObject* hobj, const char* label=0)
 449 {
 450     TH1* h = dynamic_cast<TH1*>(hobj);
 451     if (!h) return 0;
 452     if (h->GetDimension() > 2) return 0;
 453     h->SetStats(0);
 454     if ( TString(h->GetName()).Contains("Events")) { h->SetStats(1); }
 455     TString t(label);
 456     if (label) t += "_";
 457     t += h->GetName();
 458     h->SetTitle(t.Data());
 459     TCanvas* c = new TCanvas(t.Data(),t.Data());
 460     if (h->GetDimension() >= 1) {
 461         TString xlabel(h->GetXaxis()->GetTitle());
 462         if (xlabel.Contains("Pt")) { c->SetLogx();  h->GetXaxis()->SetRangeUser(0.1 , 50.); }
 463     }
 464     if (h->GetDimension() == 1) {
 465         if (xlabel.Contains("p_{T}")) { c->SetLogx();  c->SetLogy();  h->GetXaxis()->SetRangeUser(0.1 , 50.); }
 466     }
 467     if (h->GetDimension() == 2) {
 468         TString ylabel(h->GetYaxis()->GetTitle());
 469         if (ylabel.Contains("Pt")) { c->SetLogy(); h->GetYaxis()->SetRangeUser(0.1 , 50.); }
 470         if (ylabel.Contains("p_{T}")) { c->SetLogy(); h->GetYaxis()->SetRangeUser(0.1 , 50.); }
 471         h->Draw("COLZ");
 472     }
 473     if (h->GetDimension() == 1) {
 474         h->Draw();
 475     }
 476     return c;
 477
 478 }
 479
 480 Int_t CheckLoadLibrary(const char* library)
 481 {
 482   // checks if a library is already loaded, if not loads the library
 483
 484   if (strlen(gSystem->GetLibraries(Form("%s.so", library), "", kFALSE)) > 0)
 485     return 1;
 486
 487   return gSystem->Load(library);
 488 }