PWGUD/dNdPt/macros/GenerateCorrMatr_PbPb.C

   1 &
   2 Double_t GenerateCorrMatr_PbPb(const char* mcfile, const char* mctask, const char* idstring ,const char* outfile, const char* gifdir = 0)
   3 //Double_t GenerateCorrMatr_PbPb_TPCIT()
   4 {
   5
   6 Int_t c_first = 1;
   7 Int_t c_last = 11;
   8
   9 TString id = TString(idstring);
  10 if ( id.Contains("c0-5") ) { c_first = c_last = 1; }
  11 if ( id.Contains("c5-10") ) { c_first = c_last = 2; }
  12 if ( id.Contains("c10-20") ) { c_first = c_last = 3; }
  13 if ( id.Contains("c20-30") ) { c_first = c_last = 4; }
  14 if ( id.Contains("c30-40") ) { c_first = c_last = 5; }
  15 if ( id.Contains("c40-50") ) { c_first = c_last = 6; }
  16 if ( id.Contains("c50-60") ) { c_first = c_last = 7; }
  17 if ( id.Contains("c60-70") ) { c_first = c_last = 8; }
  18 if ( id.Contains("c70-80") ) { c_first = c_last = 9; }
  19 if ( id.Contains("c80-90") ) { c_first = c_last = 10; }
  20 if ( id.Contains("c90-100") ) { c_first = c_last = 11; }
  21
  22 //tmp setting
  23    //const char* mcfile = "/lustre/alice/train/V006.MC_PbPb/2011-03-15_0037.5917/mergedPeriods/MC_PbPb/2.76ATeV/LHC10h8/mknichel_dNdPtPbPb_TPCITS_VZERO1.root";
  24
  25    //const char* mcfile = "/lustre/alice/train/V006.MC_PbPb/2011-03-13_0236.5847/mergedRuns/MC_PbPb/2.76ATeV/LHC11a7/137161.ana/mknichel_dNdPtPbPb_TPCITS_VZERO1.root";
  26      //const char* mctask =  "mknichel_dNdPtPbPb_TPCITS_VZERO";
  27      //const char* idstring = "c70";
  28      //const char* outfile = "/u/mknichel/alice/dNdPt_PbPb/2011-03-15/corrMatr_LHC10h8_TPCITS_c70.root";
  29      //const char* gifdir ="/u/mknichel/alice/dNdPt_PbPb/2011-03-15/LHC10h8";
  30
  31     // settings vor zVertex cut (event and track level)
  32     Double_t zVert = 10.0;
  33
  34     // setting on eta cut (track level)
  35     Double_t eta = 0.8;
  36
  37     // strangeness scaling factor (for secondaries from strange decays)
  38     Double_t sscale = 2.0;
  39
  40     //load required libraries
  41     //load required libraries
  42     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");
  43
  44   gSystem->Load("libCore");
  45   gSystem->Load("libPhysics");
  46   gSystem->Load("libMinuit");
  47   gSystem->Load("libGui");
  48   gSystem->Load("libXMLParser");
  49
  50   gSystem->Load("libGeom");
  51   gSystem->Load("libVMC");
  52
  53   gSystem->Load("libNet");
  54
  55
  56   gSystem->Load("libSTEERBase");
  57   gSystem->Load("libESD");
  58   gSystem->Load("libCDB");
  59   gSystem->Load("libRAWDatabase");
  60   gSystem->Load("libRAWDatarec");
  61   gSystem->Load("libANALYSIS");
  62
  63
  64
  65     gSystem->Load("libANALYSIS.so");
  66     gSystem->Load("libANALYSISalice.so");
  67     gSystem->Load("libTENDER.so");
  68     gSystem->Load("libCORRFW.so");
  69     gSystem->Load("libPWG0base.so");
  70     gSystem->Load("libPWG0dep");
  71     gSystem->Load("libPWG0selectors.so");
  72
  73
  74
  75
  76
  77     // make plots nicer
  78     gROOT->SetStyle("Plain");
  79     gStyle->SetPalette(1);
  80
  81     // array for all correction matrices
  82     TObjArray* CorrMatr = new TObjArray();
  83
  84     // array for all control histograms
  85     TObjArray* ContHist = new TObjArray();
  86
  87
  88     // load mc information
  89     TFile* fmc = TFile::Open(mcfile,"READ");
  90     if (!fmc) return -1;
  91     TList* lmc = dynamic_cast<TList*>(fmc->Get(mctask));
  92     if (!lmc) return -1;
  93     AlidNdPtAnalysisPbPb *obj = dynamic_cast<AlidNdPtAnalysisPbPb*>(lmc->FindObject("dNdPtAnalysisPbPb"));
  94     if (!obj) return -1;
  95
  96 //
  97 // create rebinned thnsparse
  98 //
  99   //ptbins before rebinning: 73;
 100   const Int_t etaNbins = 30;
 101   const Int_t zvNbins = 12;
 102   const Int_t centralityBins = 11;
 103
 104   //Double_t binsMult[multNbins+1] = {-0.5,  10000.5 }; // forPbPb --reduced!
 105
 106    // Set pt binning for corrections
 107     const Int_t ptNbins = 30;
 108     Double_t binsPt[31] = {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.4, 1.6, 1.8, 2.0, 3.0, 4.0, 50.0, 100.0 };
 109
 110   Double_t binsEta[etaNbins+1] = {
 111         -1.5, -1.4, -1.3, -1.2, -1.1, -1.0, -0.9, -0.8, -0.7, -0.6,
 112         -0.5, -0.4, -0.3, -0.2, -0.1,  0.0,  0.1,  0.2,  0.3,  0.4,
 113          0.5,  0.6,  0.7,  0.8,  0.9,  1.0,  1.1,  1.2,  1.3,  1.4,
 114          1.5};
 115
 116   Double_t binsZv[zvNbins+1] = {-30.,-25.,-20.,-15.,-10.,-5.,0.,5.,10.,15.,20.,25.,30.};
 117
 118   Double_t binsCentrality[centralityBins+1] = {0., 5., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100.};
 119
 120   Int_t binsTrackEventCorrMatrix[4]={zvNbins,ptNbins,etaNbins,centralityBins};
 121   //Int_t binsTrackEvent[4]={zvNbins,ptNbins,etaNbins,centralityBins};
 122   //Int_t binsTrackPtCorrelationMatrix[4]={ptNbins,ptNbins,etaNbins,centralityBins};
 123 /*
 124   THnSparseF* fGenPrimTrackMatrix = new THnSparseF("fGenPrimTrackMatrix","mcZv:mcPt:mcEta:Centrality",4,binsTrackEventCorrMatrix);
 125   fGenPrimTrackMatrix->SetBinEdges(0,binsZv);
 126   fGenPrimTrackMatrix->SetBinEdges(1,binsPt);
 127   fGenPrimTrackMatrix->SetBinEdges(2,binsEta);
 128   fGenPrimTrackMatrix->SetBinEdges(3,binsCentrality);
 129   fGenPrimTrackMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 130   fGenPrimTrackMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
 131   fGenPrimTrackMatrix->GetAxis(2)->SetTitle("mcEta");
 132   fGenPrimTrackMatrix->GetAxis(3)->SetTitle("Centrality");
 133   fGenPrimTrackMatrix->Sumw2();
 134
 135   THnSparseF* fRecPrimTrackMatrix = new THnSparseF("fRecPrimTrackMatrix","mcZv:mcPt:mcEta:Centrality",4,binsTrackEventCorrMatrix);
 136   fRecPrimTrackMatrix->SetBinEdges(0,binsZv);
 137   fRecPrimTrackMatrix->SetBinEdges(1,binsPt);
 138   fRecPrimTrackMatrix->SetBinEdges(2,binsEta);
 139   fRecPrimTrackMatrix->SetBinEdges(3,binsCentrality);
 140   fRecPrimTrackMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 141   fRecPrimTrackMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
 142   fRecPrimTrackMatrix->GetAxis(2)->SetTitle("mcEta");
 143   fRecPrimTrackMatrix->GetAxis(3)->SetTitle("Centrality");
 144   fRecPrimTrackMatrix->Sumw2();
 145
 146   THnSparseF* fRecTrackMatrix = new THnSparseF("fRecTrackMatrix","mcZv:mcPt:mcEta:Centrality",4,binsTrackEventCorrMatrix);
 147   fRecTrackMatrix->SetBinEdges(0,binsZv);
 148   fRecTrackMatrix->SetBinEdges(1,binsPt);
 149   fRecTrackMatrix->SetBinEdges(2,binsEta);
 150   fRecTrackMatrix->SetBinEdges(3,binsCentrality);
 151   fRecTrackMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 152   fRecTrackMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
 153   fRecTrackMatrix->GetAxis(2)->SetTitle("mcEta");
 154   fRecTrackMatrix->GetAxis(3)->SetTitle("Centrality");
 155   fRecTrackMatrix->Sumw2();
 156
 157   THnSparseF* fRecSecTrackMatrix = new THnSparseF("fRecSecTrackMatrix","mcZv:mcPt:mcEta:Centrality",4,binsTrackEventCorrMatrix);
 158   fRecSecTrackMatrix->SetBinEdges(0,binsZv);
 159   fRecSecTrackMatrix->SetBinEdges(1,binsPt);
 160   fRecSecTrackMatrix->SetBinEdges(2,binsEta);
 161   fRecSecTrackMatrix->SetBinEdges(3,binsCentrality);
 162   fRecSecTrackMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 163   fRecSecTrackMatrix->GetAxis(1)->SetTitle("Pt (GeV/c)");
 164   fRecSecTrackMatrix->GetAxis(2)->SetTitle("Eta");
 165   fRecSecTrackMatrix->GetAxis(3)->SetTitle("Centrality");
 166   fRecSecTrackMatrix->Sumw2();
 167 */
 168
 169     //Event statistics
 170     THnSparse *fRecEventMatrix = obj->GetRecEventMatrix(); //all reconstructed events
 171     fRecEventMatrix->GetAxis(2)->SetRange(c_first,c_last); // select centrality
 172     TH2D* h2RecEventAll = fRecEventMatrix->Projection(0,1)->Clone("h2RecEventAll");
 173     fRecEventMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
 174     TH2D* h2RecEvent = fRecEventMatrix->Projection(0,1)->Clone("h2RecEvent");
 175     Double_t MCReconstructedEvents = h2RecEvent->Integral();
 176     Double_t MCReconstructedEventsAll = h2RecEventAll->Integral();
 177     ContHist->Add(h2RecEvent);
 178     ContHist->Add(h2RecEventAll);
 179
 180     THnSparse *fTriggerEventMatrix = obj->GetTriggerEventMatrix(); //all triggered events
 181     fTriggerEventMatrix->GetAxis(2)->SetRange(c_first,c_last); // select centrality
 182     TH2D* h2TriggerEventAll = fTriggerEventMatrix->Projection(0,1)->Clone("h2TriggerEventAll");
 183     fTriggerEventMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
 184     TH2D* h2TriggerEvent = fTriggerEventMatrix->Projection(0,1)->Clone("h2TriggerEvent");
 185     Double_t MCTriggeredEvents = h2TriggerEvent->Integral();
 186     Double_t MCTriggeredEventsAll = h2TriggerEventAll->Integral();
 187     ContHist->Add(h2TriggerEvent);
 188     ContHist->Add(h2TriggerEventAll);
 189
 190     THnSparse *fGenEventMatrix = obj->GetGenEventMatrix(); //all generated events
 191     fGenEventMatrix->GetAxis(2)->SetRange(c_first,c_last); // select centrality
 192     TH2D* h2GenEventAll = fGenEventMatrix->Projection(0,1)->Clone("h2GenEventAll");
 193     fGenEventMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
 194     TH2D* h2GenEvent = fGenEventMatrix->Projection(0,1)->Clone("h2GenEvent");
 195     Double_t MCGeneratedEvents = h2GenEvent->Integral();
 196     Double_t MCGeneratedEventsAll = h2GenEventAll->Integral();
 197     ContHist->Add(h2RecEvent);
 198     ContHist->Add(h2RecEventAll);
 199
 200     printf("=== generated MC events for correction matrices    %lf ===\n",MCGeneratedEvents);
 201     printf("=== triggered MC events for correction matrices    %lf ===\n",MCTriggeredEvents);
 202     printf("=== recontructed MC events for correction matrices %lf ===\n",MCReconstructedEvents);
 203     printf("\n");
 204     printf("=== cut on the zVertex +- %lf ===\n",zVert);
 205     printf("=== generated MC events (before zVertex cut)       %lf ===\n",MCGeneratedEventsAll);
 206     printf("=== triggered MC events (before zVertex cut)       %lf ===\n",MCTriggeredEventsAll);
 207     printf("=== recontructed MC events (before zVertex cut)    %lf ===\n",MCReconstructedEventsAll);
 208
 209     TH1D* h1MCGeneratedEvents = new TH1D("h1MCGeneratedEvents","h1MCGeneratedEvents",1,0,1);
 210     TH1D* h1MCTriggeredEvents = new TH1D("h1MCTriggeredEvents","h1MCTriggeredEvents",1,0,1);
 211     TH1D* h1MCReconstructedEvents = new TH1D("h1MCReconstructedEvents","h1MCReconstructedEvents",1,0,1);
 212     TH1D* h1MCGeneratedEventsAll = new TH1D("h1MCGeneratedEventsAll","h1MCGeneratedEventsAll",1,0,1);
 213     TH1D* h1MCTriggeredEventsAll = new TH1D("h1MCTriggeredEventsAll","h1MCTriggeredEventsAll",1,0,1);
 214     TH1D* h1MCReconstructedEventsAll = new TH1D("h1MCReconstructedEventsAll","h1MCReconstructedEventsAll",1,0,1);
 215
 216     h1MCGeneratedEvents->Fill(0,MCGeneratedEvents);
 217     h1MCGeneratedEvents->SetEntries(MCGeneratedEvents);
 218     h1MCTriggeredEvents->Fill(0,MCTriggeredEvents);
 219     h1MCTriggeredEvents->SetEntries(MCTriggeredEvents);
 220     h1MCReconstructedEvents->Fill(0,MCReconstructedEvents);
 221     h1MCReconstructedEvents->SetEntries(MCReconstructedEvents);
 222     h1MCGeneratedEventsAll->Fill(0,MCGeneratedEventsAll);
 223     h1MCGeneratedEventsAll->SetEntries(MCGeneratedEventsAll);
 224     h1MCTriggeredEventsAll->Fill(0,MCTriggeredEventsAll);
 225     h1MCTriggeredEventsAll->SetEntries(MCTriggeredEventsAll);
 226     h1MCReconstructedEventsAll->Fill(0,MCReconstructedEventsAll);
 227     h1MCReconstructedEventsAll->SetEntries(MCReconstructedEventsAll);
 228
 229     ContHist->Add(h1MCGeneratedEvents);
 230     ContHist->Add(h1MCTriggeredEvents);
 231     ContHist->Add(h1MCReconstructedEvents);
 232     ContHist->Add(h1MCGeneratedEventsAll);
 233     ContHist->Add(h1MCTriggeredEventsAll);
 234     ContHist->Add(h1MCReconstructedEventsAll);
 235
 236
 237     // efficienfy and correction matrices for tigger and vertex efficiency
 238     TH2D* h2EventTriggerEffAll  = AlidNdPtHelper::GenerateCorrMatrix(h2TriggerEventAll,h2GenEventAll,"h2EventTriggerEffAll");
 239     TH2D* h2EventTriggerCorrAll = AlidNdPtHelper::GenerateCorrMatrix(h2GenEventAll,h2TriggerEventAll,"h2EventTriggerCorrAll");
 240     TH2D* h2EventTriggerEff  = AlidNdPtHelper::GenerateCorrMatrix(h2TriggerEvent,h2GenEvent,"h2EventTriggerEff");
 241     TH2D* h2EventTriggerCorr = AlidNdPtHelper::GenerateCorrMatrix(h2GenEvent,h2TriggerEvent,"h2EventTriggerCorr");
 242
 243     TH2D* h2EventRecEffAll  = AlidNdPtHelper::GenerateCorrMatrix(h2RecEventAll,h2TriggerEventAll,"h2EventRecEffAll");
 244     TH2D* h2EventRecCorrAll = AlidNdPtHelper::GenerateCorrMatrix(h2TriggerEventAll,h2RecEventAll,"h2EventRecCorrAll");
 245     TH2D* h2EventRecEff  = AlidNdPtHelper::GenerateCorrMatrix(h2RecEvent,h2TriggerEvent,"h2EventRecEff");
 246     TH2D* h2EventRecCorr = AlidNdPtHelper::GenerateCorrMatrix(h2TriggerEvent,h2RecEvent,"h2EventRecCorr");
 247
 248     TH2D* h2EventEffAll  = AlidNdPtHelper::GenerateCorrMatrix(h2RecEventAll,h2GenEventAll,"h2EventEffAll");
 249     TH2D* h2EventCorrAll = AlidNdPtHelper::GenerateCorrMatrix(h2GenEventAll,h2RecEventAll,"h2EventCorrAll");
 250     TH2D* h2EventEff  = AlidNdPtHelper::GenerateCorrMatrix(h2RecEvent,h2GenEvent,"h2EventEff");
 251     TH2D* h2EventCorr = AlidNdPtHelper::GenerateCorrMatrix(h2GenEvent,h2RecEvent,"h2EventCorr");
 252
 253     CorrMatr->Add(h2EventTriggerEffAll);
 254     CorrMatr->Add(h2EventTriggerCorrAll);
 255     CorrMatr->Add(h2EventTriggerEff);
 256     CorrMatr->Add(h2EventTriggerCorr);
 257     CorrMatr->Add(h2EventRecEffAll);
 258     CorrMatr->Add(h2EventRecCorrAll);
 259     CorrMatr->Add(h2EventRecEff);
 260     CorrMatr->Add(h2EventRecCorr);
 261     CorrMatr->Add(h2EventEffAll);
 262     CorrMatr->Add(h2EventCorrAll);
 263     CorrMatr->Add(h2EventEff);
 264     CorrMatr->Add(h2EventCorr);
 265
 266
 267
 268     // all recontructed
 269     THnSparse *fRecTrackMatrix = obj->GetRecTrackMatrix(); //all reconstructed tracks
 270     fRecTrackMatrix->GetAxis(3)->SetRange(c_first,c_last); // select centrality
 271     fRecTrackMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
 272     fRecTrackMatrix->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
 273
 274     TH3D* h3RecTrack = fRecTrackMatrix->Projection(0,1,2)->Clone("h3RecTrack");
 275     TH2D* h2RecTrack_zv_pt  = fRecTrackMatrix->Projection(0,1)->Clone("h2RecTrack_zv_pt");
 276     TH2D* h2RecTrack_zv_eta = fRecTrackMatrix->Projection(0,2)->Clone("h2RecTrack_zv_eta");
 277     TH2D* h2RecTrack_pt_eta = fRecTrackMatrix->Projection(1,2)->Clone("h2RecTrack_pt_eta");
 278     TH1D* h1RecTrack_zv  = fRecTrackMatrix->Projection(0)->Clone("h1RecTrack_zv");
 279     TH1D* h1RecTrack_pt  = fRecTrackMatrix->Projection(1)->Clone("h1RecTrack_pt");
 280     TH1D* h1RecTrack_eta = fRecTrackMatrix->Projection(2)->Clone("h1RecTrack_eta");
 281     Double_t MCReconstructedTracks = h3RecTrack->Integral();
 282
 283     h1RecTrack_pt = (TH1D*) h1RecTrack_pt->Rebin(ptNbins,"",binsPt); // rebin hist
 284
 285     ContHist->Add(h3RecTrack);
 286     ContHist->Add(h2RecTrack_zv_pt);
 287     ContHist->Add(h2RecTrack_zv_eta);
 288     ContHist->Add(h2RecTrack_pt_eta);
 289     ContHist->Add(h1RecTrack_zv);
 290     ContHist->Add(h1RecTrack_pt);
 291     ContHist->Add(h1RecTrack_eta);
 292
 293      // recontructed primary tracks
 294     THnSparse *fRecPrimTrackMatrix = obj->GetRecPrimTrackMatrix();
 295     fRecPrimTrackMatrix->GetAxis(3)->SetRange(c_first,c_last); // select centrality
 296     THnSparse *fRecTrackMatrixScaled = fRecPrimTrackMatrix->Clone("fRecTrackMatrixScaled"); //used later for secondaries scaling
 297     fRecPrimTrackMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
 298     fRecPrimTrackMatrix->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
 299     TH3D* h3RecPrimTrack = fRecPrimTrackMatrix->Projection(0,1,2)->Clone("h3RecPrimTrack");
 300     TH2D* h2RecPrimTrack_zv_pt  = fRecPrimTrackMatrix->Projection(0,1)->Clone("h2RecPrimTrack_zv_pt");
 301     TH2D* h2RecPrimTrack_zv_eta = fRecPrimTrackMatrix->Projection(0,2)->Clone("h2RecPrimTrack_zv_eta");
 302     TH2D* h2RecPrimTrack_pt_eta = fRecPrimTrackMatrix->Projection(1,2)->Clone("h2RecPrimTrack_pt_eta");
 303     TH1D* h1RecPrimTrack_zv  = fRecPrimTrackMatrix->Projection(0)->Clone("h1RecPrimTrack_zv");
 304     TH1D* h1RecPrimTrack_pt  = fRecPrimTrackMatrix->Projection(1)->Clone("h1RecPrimTrack_pt");
 305     TH1D* h1RecPrimTrack_eta = fRecPrimTrackMatrix->Projection(2)->Clone("h1RecPrimTrack_eta");
 306     Double_t MCReconstructedPrimTracks = h3RecPrimTrack->Integral();
 307
 308     h1RecPrimTrack_pt = (TH1D*) h1RecPrimTrack_pt->Rebin(ptNbins,"",binsPt); // rebin hist
 309
 310     ContHist->Add(h3RecPrimTrack);
 311     ContHist->Add(h2RecPrimTrack_zv_pt);
 312     ContHist->Add(h2RecPrimTrack_zv_eta);
 313     ContHist->Add(h2RecPrimTrack_pt_eta);
 314     ContHist->Add(h1RecPrimTrack_zv);
 315     ContHist->Add(h1RecPrimTrack_pt);
 316     ContHist->Add(h1RecPrimTrack_eta);
 317
 318      // recontructed secondary tracks
 319     THnSparse *fRecSecTrackMatrix = obj->GetRecSecTrackMatrix();
 320     fRecSecTrackMatrix->GetAxis(3)->SetRange(c_first,c_last); // select centrality
 321     THnSparse *fRecSecTrackMatrixScaled = fRecSecTrackMatrix->Clone("fRecSecTrackMatrixScaled"); //used later for secondaries scaling
 322     fRecSecTrackMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
 323     fRecSecTrackMatrix->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
 324     TH3D* h3RecSecTrack = fRecSecTrackMatrix->Projection(0,1,2)->Clone("h3RecSecTrack");
 325     TH2D* h2RecSecTrack_zv_pt  = fRecSecTrackMatrix->Projection(0,1)->Clone("h2RecSecTrack_zv_pt");
 326     TH2D* h2RecSecTrack_zv_eta = fRecSecTrackMatrix->Projection(0,2)->Clone("h2RecSecTrack_zv_eta");
 327     TH2D* h2RecSecTrack_pt_eta = fRecSecTrackMatrix->Projection(1,2)->Clone("h2RecSecTrack_pt_eta");
 328     TH1D* h1RecSecTrack_zv  = fRecSecTrackMatrix->Projection(0)->Clone("h1RecSecTrack_zv");
 329     TH1D* h1RecSecTrack_pt  = fRecSecTrackMatrix->Projection(1)->Clone("h1RecSecTrack_pt");
 330     TH1D* h1RecSecTrack_eta = fRecSecTrackMatrix->Projection(2)->Clone("h1RecSecTrack_eta");
 331     Double_t MCReconstructedSecTracks = h3RecSecTrack->Integral();
 332
 333     h1RecSecTrack_pt = (TH1D*) h1RecSecTrack_pt->Rebin(ptNbins,"",binsPt); // rebin hist
 334
 335     ContHist->Add(h3RecSecTrack);
 336     ContHist->Add(h2RecSecTrack_zv_pt);
 337     ContHist->Add(h2RecSecTrack_zv_eta);
 338     ContHist->Add(h2RecSecTrack_pt_eta);
 339     ContHist->Add(h1RecSecTrack_zv);
 340     ContHist->Add(h1RecSecTrack_pt);
 341     ContHist->Add(h1RecSecTrack_eta);
 342
 343     // generated primary tracks
 344     THnSparse *fGenPrimTrackMatrix =obj->GetGenPrimTrackMatrix();
 345     fGenPrimTrackMatrix->GetAxis(3)->SetRange(c_first,c_last); // select centrality
 346     fGenPrimTrackMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
 347     fGenPrimTrackMatrix->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
 348     TH3D* h3GenPrimTrack = fGenPrimTrackMatrix->Projection(0,1,2)->Clone("h3GenPrimTrack");
 349     TH2D* h2GenPrimTrack_zv_Pt  = fGenPrimTrackMatrix->Projection(0,1)->Clone("h2GenPrimTrack_zv_pt");
 350     TH2D* h2GenPrimTrack_zv_eta = fGenPrimTrackMatrix->Projection(0,2)->Clone("h2GenPrimTrack_zv_eta");
 351     TH2D* h2GenPrimTrack_pt_eta = fGenPrimTrackMatrix->Projection(1,2)->Clone("h2GenPrimTrack_pt_eta");
 352     TH1D* h1GenPrimTrack_zv  = fGenPrimTrackMatrix->Projection(0)->Clone("h1GenPrimTrack_zv");
 353     TH1D* h1GenPrimTrack_pt  = fGenPrimTrackMatrix->Projection(1)->Clone("h1GenPrimTrack_pt");
 354     TH1D* h1GenPrimTrack_eta = fGenPrimTrackMatrix->Projection(2)->Clone("h1GenPrimTrack_eta");
 355     Double_t MCGeneratedPrimTracks = h3GenPrimTrack->Integral();
 356
 357
 358     // mc truth histogram (for self-consistency check)
 359     TH1D* dNdPt_MC   = h1GenPrimTrack_pt->Clone("dNdPt_MC");
 360     // normalization and finalization
 361     // 1/N_evt 1/(2 pi pt) 1/width 1/etarange
 362     for (int ii=1; ii <= dNdPt_MC->GetNbinsX() ;ii++) {
 363         Double_t pt = dNdPt_MC->GetBinCenter(ii);
 364         Double_t width = dNdPt_MC->GetBinWidth(ii);
 365         Double_t val = dNdPt_MC->GetBinContent(ii);
 366         Double_t err = dNdPt_MC->GetBinError(ii);
 367         Double_t cval = (val)/(width * 2.0 * TMath::Pi() * 1.6 * MCGeneratedEvents * pt);
 368         Double_t cerr = (err)/(width * 2.0 * TMath::Pi() * 1.6 * MCGeneratedEvents * pt);
 369         dNdPt_MC->SetBinContent(ii,cval);
 370         dNdPt_MC->SetBinError(ii,cerr);
 371     }
 372     dNdPt_MC->SetMarkerStyle(21);
 373     dNdPt_MC->SetTitle("; p_{T} (GeV/c) ; 1/N_{evt} 1/(2#pi p_{T}) (d^{2}N_{ch})/(d#eta dp_{T})^{-2}");
 374     ContHist->Add(dNdPt_MC);
 375
 376
 377
 378     h1GenPrimTrack_pt = (TH1D*) h1GenPrimTrack_pt->Rebin(ptNbins,"",binsPt); // rebin hist
 379
 380     ContHist->Add(h3GenPrimTrack);
 381     ContHist->Add(h2GenPrimTrack_zv_pt);
 382     ContHist->Add(h2GenPrimTrack_zv_eta);
 383     ContHist->Add(h2GenPrimTrack_pt_eta);
 384     ContHist->Add(h1GenPrimTrack_zv);
 385     ContHist->Add(h1GenPrimTrack_pt);
 386     ContHist->Add(h1GenPrimTrack_eta);
 387     printf("\n");
 388     printf("==============================================================\n");
 389     printf("=== recontructed MC tracks              %lf ===\n",MCReconstructedTracks);
 390     printf("=== recontructed MC secondary tracks    %lf ===\n",MCReconstructedSecTracks);
 391     printf("=== recontructed MC primary tracks      %lf ===\n",MCReconstructedPrimTracks);
 392     printf("=== generated MC primary track          %lf ===\n",MCGeneratedPrimTracks);
 393     printf("==============================================================\n");
 394     printf("\n");
 395
 396 //    THnSparse *fSparseTriggerTrackEvent = obj->GetTriggerTrackEventMatrix();//Tracks from triggered events
 397 //    THnSparse *fSparseVtxTrackEvent = obj->GetRecTrackEventMatrix();//Tracks from events with rec. vtx
 398 //    THnSparse *fSparseGenTrackEvent = obj->GetGenTrackEventMatrix();//generated TrackEvent matrix
 399
 400     // tracking efficiencies + corrections
 401    TH3D* h3TrackEff  = AlidNdPtHelper::GenerateCorrMatrix(h3RecPrimTrack,h3GenPrimTrack,"h3TrackEff");
 402    TH3D* h3TrackCorr = AlidNdPtHelper::GenerateCorrMatrix(h3GenPrimTrack,h3RecPrimTrack,"h3TrackCorr");
 403
 404    TH2D* h2TrackEff_zv_pt   = AlidNdPtHelper::GenerateCorrMatrix(h2RecPrimTrack_zv_pt,h2GenPrimTrack_zv_pt,"h2TrackEff_zv_pt");
 405    TH2D* h2TrackCorr_zv_pt  = AlidNdPtHelper::GenerateCorrMatrix(h2GenPrimTrack_zv_pt,h2RecPrimTrack_zv_pt,"h2TrackCorr_zv_pt");
 406    TH2D* h2TrackEff_zv_eta  = AlidNdPtHelper::GenerateCorrMatrix(h2RecPrimTrack_zv_eta,h2GenPrimTrack_zv_eta,"h2TrackEff_zv_eta");
 407    TH2D* h2TrackCorr_zv_eta = AlidNdPtHelper::GenerateCorrMatrix(h2GenPrimTrack_zv_eta,h2RecPrimTrack_zv_eta,"h2TrackCorr_zv_eta");
 408    TH2D* h2TrackEff_pt_eta  = AlidNdPtHelper::GenerateCorrMatrix(h2RecPrimTrack_pt_eta,h2GenPrimTrack_pt_eta,"h2TrackEff_pt_eta");
 409    TH2D* h2TrackCorr_pt_eta = AlidNdPtHelper::GenerateCorrMatrix(h2GenPrimTrack_pt_eta,h2RecPrimTrack_pt_eta,"h2TrackCorr_pt_eta");
 410
 411    TH1D* h1TrackEff_zv   = AlidNdPtHelper::GenerateCorrMatrix(h1RecPrimTrack_zv,h1GenPrimTrack_zv,"h1TrackEff_zv");
 412    TH1D* h1TrackCorr_zv  = AlidNdPtHelper::GenerateCorrMatrix(h1GenPrimTrack_zv,h1RecPrimTrack_zv,"h1TrackCorr_zv");
 413    TH1D* h1TrackEff_pt   = AlidNdPtHelper::GenerateCorrMatrix(h1RecPrimTrack_pt,h1GenPrimTrack_pt,"h1TrackEff_pt");
 414    TH1D* h1TrackCorr_pt  = AlidNdPtHelper::GenerateCorrMatrix(h1GenPrimTrack_pt,h1RecPrimTrack_pt,"h1TrackCorr_pt");
 415    TH1D* h1TrackEff_eta  = AlidNdPtHelper::GenerateCorrMatrix(h1RecPrimTrack_eta,h1GenPrimTrack_eta,"h1TrackEff_eta");
 416    TH1D* h1TrackCorr_eta = AlidNdPtHelper::GenerateCorrMatrix(h1GenPrimTrack_eta,h1RecPrimTrack_eta,"h1TrackCorr_eta");
 417
 418    CorrMatr->Add(h3TrackEff);
 419    CorrMatr->Add(h3TrackCorr);
 420    CorrMatr->Add(h2TrackEff_zv_pt);
 421    CorrMatr->Add(h2TrackCorr_zv_pt);
 422    CorrMatr->Add(h2TrackEff_zv_eta);
 423    CorrMatr->Add(h2TrackCorr_zv_eta);
 424    CorrMatr->Add(h2TrackEff_pt_eta);
 425    CorrMatr->Add(h2TrackCorr_pt_eta);
 426    CorrMatr->Add(h1TrackEff_zv);
 427    CorrMatr->Add(h1TrackCorr_zv);
 428    CorrMatr->Add(h1TrackEff_pt);
 429    CorrMatr->Add(h1TrackCorr_pt);
 430    CorrMatr->Add(h1TrackEff_eta);
 431    CorrMatr->Add(h1TrackCorr_eta);
 432
 433    // scale the secondaries before calculating correction matrices
 434    for (Long64_t i = 0; i < fRecSecTrackMatrixScaled->GetNbins(); i++) {
 435        Int_t c[3];
 436        Double_t val = fRecSecTrackMatrixScaled->GetBinContent(i,c);
 437        Double_t err = fRecSecTrackMatrixScaled->GetBinError(i);
 438        Double_t pt =  fRecSecTrackMatrixScaled->GetAxis(1)->GetBinCenter(c[1]);
 439        Double_t scale = GetStrangenessCorrFactorPbPb(pt, sscale);
 440        fRecSecTrackMatrixScaled->SetBinContent(c,val*scale);
 441        fRecSecTrackMatrixScaled->SetBinError(c,err*scale);
 442     }
 443
 444     // for correct determination of secondaries contamination, also the total total tracks have to be scaled
 445     // this is done by taking primaries and adding the scaled secondaries
 446     fRecTrackMatrixScaled->Add(fRecSecTrackMatrixScaled);
 447
 448     fRecSecTrackMatrixScaled->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
 449     fRecSecTrackMatrixScaled->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
 450     fRecSecTrackMatrixScaled->GetAxis(3)->SetRange(c_first,c_last); // select centrality
 451
 452     TH3D* h3RecSecTrackScaled = fRecSecTrackMatrixScaled->Projection(0,1,2)->Clone("h3RecSecTrackScaled");
 453     TH2D* h2RecSecTrackScaled_zv_pt  = fRecSecTrackMatrixScaled->Projection(0,1)->Clone("h2RecSecTrackScaled_zv_pt");
 454     TH2D* h2RecSecTrackScaled_zv_eta = fRecSecTrackMatrixScaled->Projection(0,2)->Clone("h2RecSecTrackScaled_zv_eta");
 455     TH2D* h2RecSecTrackScaled_pt_eta = fRecSecTrackMatrixScaled->Projection(1,2)->Clone("h2RecSecTrackScaled_pt_eta");
 456     TH1D* h1RecSecTrackScaled_zv  = fRecSecTrackMatrixScaled->Projection(0)->Clone("h1RecSecTrackScaled_zv");
 457     TH1D* h1RecSecTrackScaled_pt  = fRecSecTrackMatrixScaled->Projection(1)->Clone("h1RecSecTrackScaled_pt");
 458     TH1D* h1RecSecTrackScaled_eta = fRecSecTrackMatrixScaled->Projection(2)->Clone("h1RecSecTrackScaled_eta");
 459
 460     h1RecSecTrackScaled_pt = (TH1D*) h1RecSecTrackScaled_pt->Rebin(ptNbins,"",binsPt); // rebin hist
 461
 462     ContHist->Add(h3RecSecTrackScaled);
 463     ContHist->Add(h2RecSecTrackScaled_zv_pt);
 464     ContHist->Add(h2RecSecTrackScaled_zv_eta);
 465     ContHist->Add(h2RecSecTrackScaled_pt_eta);
 466     ContHist->Add(h1RecSecTrackScaled_zv);
 467     ContHist->Add(h1RecSecTrackScaled_pt);
 468     ContHist->Add(h1RecSecTrackScaled_eta);
 469
 470     fRecTrackMatrixScaled->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
 471     fRecTrackMatrixScaled->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
 472
 473     TH3D* h3RecTrackScaled = fRecTrackMatrixScaled->Projection(0,1,2)->Clone("h3RecTrackScaled");
 474     TH2D* h2RecTrackScaled_zv_pt  = fRecTrackMatrixScaled->Projection(0,1)->Clone("h2RecTrackScaled_zv_pt");
 475     TH2D* h2RecTrackScaled_zv_eta = fRecTrackMatrixScaled->Projection(0,2)->Clone("h2RecTrackScaled_zv_eta");
 476     TH2D* h2RecTrackScaled_pt_eta = fRecTrackMatrixScaled->Projection(1,2)->Clone("h2RecTrackScaled_pt_eta");
 477     TH1D* h1RecTrackScaled_zv  = fRecTrackMatrixScaled->Projection(0)->Clone("h1RecTrackScaled_zv");
 478     TH1D* h1RecTrackScaled_pt  = fRecTrackMatrixScaled->Projection(1)->Clone("h1RecTrackScaled_pt");
 479     TH1D* h1RecTrackScaled_eta = fRecTrackMatrixScaled->Projection(2)->Clone("h1RecTrackScaled_eta");
 480
 481
 482     h1RecTrackScaled_pt = (TH1D*) h1RecTrackScaled_pt->Rebin(ptNbins,"",binsPt); // rebin hist
 483
 484
 485     ContHist->Add(h3RecTrackScaled);
 486     ContHist->Add(h2RecTrackScaled_zv_pt);
 487     ContHist->Add(h2RecTrackScaled_zv_eta);
 488     ContHist->Add(h2RecTrackScaled_pt_eta);
 489     ContHist->Add(h1RecTrackScaled_zv);
 490     ContHist->Add(h1RecTrackScaled_pt);
 491     ContHist->Add(h1RecTrackScaled_eta);
 492
 493     // create histograms for secondaries contamination and correction
 494
 495     TH3D* h3SecCont = AlidNdPtHelper::GenerateCorrMatrix(h3RecSecTrackScaled,h3RecTrackScaled,"h3SecCont");
 496     TH3D* h3SecCorr = AlidNdPtHelper::GenerateContCorrMatrix(h3RecSecTrackScaled,h3RecTrackScaled,"h3SecCorr");
 497     TH2D* h2SecCont_zv_pt  = AlidNdPtHelper::GenerateCorrMatrix(h2RecSecTrackScaled_zv_pt,h2RecTrackScaled_zv_pt,"h2SecCont_zv_pt");
 498     TH2D* h2SecCorr_zv_pt  = AlidNdPtHelper::GenerateContCorrMatrix(h2RecSecTrackScaled_zv_pt,h2RecTrackScaled_zv_pt,"h2SecCorr_zv_pt");
 499     TH2D* h2SecCont_zv_eta = AlidNdPtHelper::GenerateCorrMatrix(h2RecSecTrackScaled_zv_eta,h2RecTrackScaled_zv_eta,"h2SecCont_zv_eta");
 500     TH2D* h2SecCorr_zv_eta = AlidNdPtHelper::GenerateContCorrMatrix(h2RecSecTrackScaled_zv_eta,h2RecTrackScaled_zv_eta,"h2SecCorr_zv_eta");
 501     TH2D* h2SecCont_pt_eta = AlidNdPtHelper::GenerateCorrMatrix(h2RecSecTrackScaled_pt_eta,h2RecTrackScaled_pt_eta,"h2SecCont_pt_eta");
 502     TH2D* h2SecCorr_pt_eta = AlidNdPtHelper::GenerateContCorrMatrix(h2RecSecTrackScaled_pt_eta,h2RecTrackScaled_pt_eta,"h2SecCorr_pt_eta");
 503     TH1D* h1SecCont_zv = AlidNdPtHelper::GenerateCorrMatrix(h1RecSecTrackScaled_zv,h1RecTrackScaled_zv,"h1SecCont_zv");
 504     TH1D* h1SecCorr_zv = AlidNdPtHelper::GenerateContCorrMatrix(h1RecSecTrackScaled_zv,h1RecTrackScaled_zv,"h1SecCorr_zv");
 505     TH1D* h1SecCont_pt = AlidNdPtHelper::GenerateCorrMatrix(h1RecSecTrackScaled_pt,h1RecTrackScaled_pt,"h1SecCont_pt");
 506     TH1D* h1SecCorr_pt = AlidNdPtHelper::GenerateContCorrMatrix(h1RecSecTrackScaled_pt,h1RecTrackScaled_pt,"h1SecCorr_pt");
 507     TH1D* h1SecCont_eta = AlidNdPtHelper::GenerateCorrMatrix(h1RecSecTrackScaled_eta,h1RecTrackScaled_eta,"h1SecCont_eta");
 508     TH1D* h1SecCorr_eta = AlidNdPtHelper::GenerateContCorrMatrix(h1RecSecTrackScaled_eta,h1RecTrackScaled_eta,"h1SecCorr_eta");
 509
 510     CorrMatr->Add(h3SecCont);
 511     CorrMatr->Add(h3SecCorr);
 512     CorrMatr->Add(h2SecCont_zv_pt);
 513     CorrMatr->Add(h2SecCorr_zv_pt);
 514     CorrMatr->Add(h2SecCont_zv_eta);
 515     CorrMatr->Add(h2SecCorr_zv_eta);
 516     CorrMatr->Add(h2SecCont_pt_eta);
 517     CorrMatr->Add(h2SecCorr_pt_eta);
 518     CorrMatr->Add(h1SecCont_zv);
 519     CorrMatr->Add(h1SecCorr_zv);
 520     CorrMatr->Add(h1SecCont_pt);
 521     CorrMatr->Add(h1SecCorr_pt);
 522     CorrMatr->Add(h1SecCont_eta);
 523     CorrMatr->Add(h1SecCorr_eta);
 524
 525     // plot pictures and save to gifdir
 526     for (i=0; i < CorrMatr->LastIndex(); i++) {
 527         TCanvas* ctmp = PlotHist(CorrMatr->At(i),idstring);
 528         if (gifdir && ctmp) {
 529             TString gif(gifdir);
 530             gif += '/';
 531             gif += ctmp->GetName();
 532             gif += ".gif";
 533             ctmp->SaveAs(gif.Data(),"gif");
 534             delete ctmp;
 535         }
 536     }
 537     for (i=0; i < ContHist->LastIndex(); i++) {
 538         TCanvas* ctmp = PlotHist(ContHist->At(i),idstring);
 539         if (gifdir && ctmp) {
 540             TString gif(gifdir);
 541             gif += '/';
 542             gif += ctmp->GetName();
 543             gif += ".gif";
 544             ctmp->SaveAs(gif.Data(),"gif");
 545             delete ctmp;
 546         }
 547    }
 548
 549     // save all correction matrices and control histograms to file
 550     if (!outfile) { return; }
 551     TFile *out = TFile::Open(outfile,"RECREATE");
 552     CorrMatr->Write();
 553     ContHist->Write();
 554     out->Close();
 555
 556     return MCReconstructedEvents;
 557
 558 }
 559
 560
 561 //_____________________________________________________________________________
 562 Double_t GetStrangenessCorrFactorPbPb(Double_t pt, Double_t s)
 563 {
 564     // data driven correction factor for secondaries (PbPb)
 565
 566     if (pt <= 0.25) return 1.0;
 567     if (pt <= 0.5) return GetLinearInterpolationValue(0.25,1.0,0.5,0.60+0.40*s, pt);
 568     if (pt <= 1.0) return GetLinearInterpolationValue(0.5,0.60+0.40*s,1.0,0.53+0.47*s, pt);
 569     if (pt <= 2.0) return GetLinearInterpolationValue(1.0,0.53+0.47*s,2.0,0.44+0.56*s, pt);
 570     if (pt <= 5.0) return GetLinearInterpolationValue(2.0,0.44+0.56*s,5.0,0.33+0.67*s, pt);
 571     return 0.33+0.67*s;
 572 }
 573
 574
 575 //___________________________________________________________________________
 576 Double_t GetLinearInterpolationValue(const Double_t x1,const  Double_t y1,const  Double_t x2,const  Double_t y2, const Double_t pt)
 577 {
 578     //
 579     // linear interpolation
 580     //
 581     return ((y2-y1)/(x2-x1))*pt+(y2-(((y2-y1)/(x2-x1))*x2));
 582 }
 583
 584 //___________________________________________________________________________
 585 TCanvas* PlotHist(TObject* hobj, const char* label=0)
 586 {
 587     TH1* h = dynamic_cast<TH1*>(hobj);
 588     if (!h) return 0;
 589     if (h->GetDimension() > 2) return 0;
 590     h->SetStats(0);
 591     if ( TString(h->GetName()).Contains("Events")) { h->SetStats(1); }
 592     TString t(label);
 593     if (label) t += "_";
 594     t += h->GetName();
 595     h->SetTitle(t.Data());
 596     TCanvas* c = new TCanvas(t.Data(),t.Data());
 597     if (h->GetDimension() >= 1) {
 598         TString xlabel(h->GetXaxis()->GetTitle());
 599         if (xlabel.Contains("Pt")) { c->SetLogx();  h->GetXaxis()->SetRangeUser(0.1 , 100.); }
 600     }
 601     if (h->GetDimension() == 1) {
 602         if (xlabel.Contains("p_{T}")) { c->SetLogx();  c->SetLogy();  h->GetXaxis()->SetRangeUser(0.1 , 100.); }
 603     }
 604     if (h->GetDimension() == 2) {
 605         TString ylabel(h->GetYaxis()->GetTitle());
 606         if (ylabel.Contains("Pt")) { c->SetLogy(); h->GetYaxis()->SetRangeUser(0.1 , 100.); }
 607         if (ylabel.Contains("p_{T}")) { c->SetLogy(); h->GetYaxis()->SetRangeUser(0.1 , 100.); }
 608         h->Draw("COLZ");
 609     }
 610     if (h->GetDimension() == 1) {
 611         h->Draw();
 612     }
 613     return c;
 614
 615 }
 616
 617 Int_t CheckLoadLibrary(const char* library)
 618 {
 619   // checks if a library is already loaded, if not loads the library
 620
 621   if (strlen(gSystem->GetLibraries(Form("%s.so", library), "", kFALSE)) > 0)
 622     return 1;
 623
 624   return gSystem->Load(library);
 625 }