TENDER becomes Tender
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / ChargedHadrons / 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 }