Scripts to make FMD analysis correction objects
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis / AliFMDAnalysisTaskBFCorrelation.cxx
1 #include <TROOT.h>
2 #include <TSystem.h>
3 #include <TInterpreter.h>
4 #include <TChain.h>
5 #include <TFile.h>
6 #include <TList.h>
7 #include <iostream>
8 #include <string>
9 #include <TCanvas.h>
10 #include "TH1F.h"
11 #include "TH2F.h"
12 #include "TProfile.h"
13 #include "AliFMDAnalysisTaskBFCorrelation.h"
14 #include "AliAnalysisManager.h"
15 #include "AliESDFMD.h"
16 #include "AliESDEvent.h"
17 #include "AliAODEvent.h"
18 #include "AliAODHandler.h"
19 #include "AliMCEventHandler.h"
20 #include "AliStack.h"
21 #include "AliLog.h"
22 #include "AliESDVertex.h"
23 #include "TMath.h"
24 #include "AliFMDAnaParameters.h"
25 //#include "AliFMDGeometry.h"
26 #include "AliGenEventHeader.h"
27 #include "AliGenPythiaEventHeader.h"
28 #include "AliHeader.h"
29 //#include "TDatabasePDG.h"
30 //#include "TParticlePDG.h"
31 #include "AliESDInputHandler.h"
32
33 ClassImp(AliFMDAnalysisTaskBFCorrelation)
34
35 AliFMDAnalysisTaskBFCorrelation::AliFMDAnalysisTaskBFCorrelation()
36 : fDebug(0),
37   fOutputList(0),
38   fInputList(0),
39   fInternalList(0),
40   fVertexString(0x0),
41   fStandalone(kTRUE),
42   fEvent(0),
43   fnBinsX(200),
44   fXmin(-6),
45   fXmax(6),
46   fnBinsY(20),
47   fYmin(0),
48   fYmax(2 * TMath::Pi()),
49   c(0),
50   debug0(0),
51   debug1(0)
52 {
53   // Default constructor
54   DefineInput (0, TList::Class());
55   DefineOutput(0, TList::Class());
56 }
57
58 //_____________________________________________________________________
59 AliFMDAnalysisTaskBFCorrelation::AliFMDAnalysisTaskBFCorrelation(const char* name, Bool_t SE):
60   AliAnalysisTask(name,name),
61   fDebug(0),
62   fOutputList(0),
63   fInputList(0),
64   fInternalList(0),
65   fVertexString(0x0),
66   fStandalone(kTRUE),
67   fEvent(0),
68   fnBinsX(200),
69   fXmin(-6),
70   fXmax(6),
71   fnBinsY(20),
72   fYmin(0),
73   fYmax(2 * TMath::Pi()),
74   c(0),
75   debug0(0),
76   debug1(0)
77 {
78   fStandalone = SE;
79   if(fStandalone) {
80     DefineInput (0, TList::Class());
81     DefineInput(1, TObjString::Class());
82     DefineOutput(0, TList::Class());
83   }
84 }
85
86 //_____________________________________________________________________
87 void AliFMDAnalysisTaskBFCorrelation::CreateOutputObjects()
88 {
89   // Setup the list for storing results, if it does not exist
90
91   if(!fOutputList) {
92     fOutputList = new TList();
93     fOutputList->SetName("BackgroundCorrected");
94   }
95
96   // Setup the list for temporary storage during the analysis
97   
98   if(!fInternalList) {
99     fInternalList = new TList();
100     fInternalList->SetName("InternalBFList");
101   }
102   
103   // Set up histograms for analysis and storage. 4 different binnings
104   
105   for (Int_t i = 1; i <= 4; i++) { 
106     
107     // Temporary histograms for storing hist per eta summed over phi
108     
109     TH1D *hESDMult = new TH1D(Form("hESDMult_binning%d", i), 
110                               Form("Multiplicity per event vs eta ESD (%d bins)", 30*i), 
111                               18*i, -9, 9);
112     TH1D *hMCMult  = new TH1D(Form("hMCMult_binning%d", i), 
113                               Form("Multiplicity per event vs eta MC-truth (%d bins)", 30*i), 
114                               18*i, -9, 9);
115     fInternalList->Add(hESDMult);
116     fInternalList->Add(hMCMult);
117
118     
119     // Histograms for storing the 5 values to calculate b (ESD & MC)
120     
121     TH1D *pnFnB = new TH1D(Form("pnFnB_binning%d", i),
122                            Form("Correlation ESD (%d bins)", i*9), 
123                            9*i, 0, 9);
124     TH1D *pnF2  = new TH1D(Form("pnF2_binning%d", i),
125                            Form("Fwd^2 ESD (%d bins)", i*9),
126                            9*i, 0, 9);
127     TH1D *pnB2  = new TH1D(Form("pnB2_binning%d", i), 
128                            Form("Bwd^2 ESD (%d bins)", i*9), 
129                            9*i, 0, 9);
130     TH1D *pnF   = new TH1D(Form("pnF_binning%d", i), 
131                            Form("Fwd ESD (%d bins)", i*9), 
132                            9*i, 0, 9);
133     TH1D *pnB   = new TH1D(Form("pnB_binning%d", i),
134                            Form("Bwd ESD (%d bins)", i*9),
135                            9*i, 0, 9);    
136     
137     TH1D *pnFnB_MC = new TH1D(Form("pnFnB_MC_binning%d", i),
138                               Form("Correlation MC (%d bins)", i*9), 
139                               9*i, 0, 9);
140     TH1D *pnF2_MC  = new TH1D(Form("pnF2_MC_binning%d", i),
141                               Form("Fwd^2 MC (%d bins)", i*9), 
142                               9*i, 0, 9);
143     TH1D *pnB2_MC  = new TH1D(Form("pnB2_MC_binning%d", i),
144                               Form("Bwd^2 MC (%d bins)", i*9),
145                               9*i, 0, 9);
146     TH1D *pnF_MC   = new TH1D(Form("pnF_MC_binning%d", i),
147                               Form("Fwd MC (%d bins)", i*9), 
148                               9*i, 0, 9);
149     TH1D *pnB_MC   = new TH1D(Form("pnB_MC_binning%d", i),
150                               Form("Backward MC (%d bins)" ,i*9), 
151                               9*i, 0, 9);
152     fOutputList->Add(pnFnB);
153     fOutputList->Add(pnF2);
154     fOutputList->Add(pnB2);
155     fOutputList->Add(pnF);
156     fOutputList->Add(pnB);
157     
158     fOutputList->Add(pnFnB_MC);
159     fOutputList->Add(pnF2_MC);
160     fOutputList->Add(pnB2_MC);
161     fOutputList->Add(pnF_MC);
162     fOutputList->Add(pnB_MC);
163     
164     // Debugging histograms : "dNdEta"
165     
166     //    TH1D *hESDMultvsEta = new TH1D(Form("hESDMultvsEta_binning%d" ,i), 
167     //                     Form("Mult vs Eta (%d bins, ESD)", 18*i),
168     //                     18*i, -9, 9);
169     TH1D *hESDMultvsEta = new TH1D(Form("hESDMultvsEta_binning%d" ,i), 
170                                    Form("Mult vs Eta (%d bins, ESD)", i*18),
171                                    200, -6, 6);
172     TH1D *hMCMultvsEta = new TH1D(Form("hMCMultvsEta_binning%d", i),
173                                   Form("Mult vs Eta (%d bins, MC)", 18*i), 
174                                   200, -6, 6);
175     //TH1D *hMCMultvsEta = new TH1D(Form("hMCMultvsEta_binning%d", i),
176     //                    Form("Mult vs Eta (%d bins, MC)", 18*i), 
177     //                    18*i, -9, 9);
178     fOutputList->Add(hESDMultvsEta);
179     fOutputList->Add(hMCMultvsEta);
180
181     TH1D *hESDMultvsEtaC = new TH1D(Form("hESDMultvsEtaC_binning%d" ,i), 
182                                    Form("Mult vs Eta C (%d bins, ESD)", 18*i),
183                                    18*i, -9, 9);
184     TH1D *hMCMultvsEtaC  = new TH1D(Form("hMCMultvsEtaC_binning%d", i),
185                                   Form("Mult vs Eta C (%d bins, MC)", 18*i), 
186                                   18*i, -9, 9);
187     fOutputList->Add(hESDMultvsEtaC);
188     fOutputList->Add(hMCMultvsEtaC);
189
190     // Debugging histograms : Different distributions
191
192     TH2D *hESDBwdDist = new TH2D(Form("hESDBwdDist_binning%d", i),
193                                  Form("Distribution of Bwd values (%d bins)", i*9),
194                                  9*i, 0, 9, 20, 0, 20);
195     
196     TH2D *hESDBwd2Dist = new TH2D(Form("hESDBwd2Dist_binning%d", i),
197                                   Form("Distribution of Bwd^2 values (%d bins)", i*9),
198                                   9*i, 0, 9, 400, 0, 400);
199     
200     TH2D *hMCBwdDist = new TH2D(Form("hMCBwdDist_binning%d", i),
201                                 Form("Distribution of Bwd values (%d bins)", i*9),
202                                 9*i, 0, 9, 20, 0, 20);
203     
204     TH2D *hMCBwd2Dist = new TH2D(Form("hMCBwd2Dist_binning%d", i),
205                                  Form("Distribution of Bwd^2 values (%d bins)", i*9),
206                                  9*i, 0, 9, 400, 0, 400);
207     fOutputList->Add(hESDBwdDist);
208     fOutputList->Add(hESDBwd2Dist);
209     fOutputList->Add(hMCBwdDist);
210     fOutputList->Add(hMCBwd2Dist);
211   }
212 }
213
214 //_____________________________________________________________________
215 void AliFMDAnalysisTaskBFCorrelation::ConnectInputData(Option_t */*option*/)
216 {
217   if(fStandalone) {
218     fInputList   = (TList*)GetInputData(0);
219     fVertexString = (TObjString*)GetInputData(1);
220   }
221 }
222
223 //_____________________________________________________________________
224 void AliFMDAnalysisTaskBFCorrelation::Exec(Option_t */*option*/)
225 {
226   fEvent++;
227   if (fEvent % 1000 == 0) 
228     std::cout << "Event # " << fEvent << std::endl;
229   
230   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
231   
232   fVertexString = (TObjString*)fInputList->At(0);
233   Int_t vtxbin   = fVertexString->GetString().Atoi();
234   if(vtxbin !=4) return;
235   SetBounds();
236   CountESDHits();
237   CalculateValues("ESD");
238   
239   if(pars->GetProcessPrimary())
240     ProcessPrimary();
241   
242   if(fStandalone) {
243     PostData(0, fOutputList); 
244   }
245 }
246 //_____________________________________________________________________
247 void AliFMDAnalysisTaskBFCorrelation::CountESDHits() {
248
249   TH2F *dNdEtadPhiHist = (TH2F*)fInputList->FindObject("dNdetadphiHistogramSPDTrVtx");
250
251   for (Int_t i = 1; i<=4; i++) {
252     
253     TH1D *hESDMult = (TH1D*)fInternalList->FindObject(Form("hESDMult_binning%d", i));
254     hESDMult->Reset();
255   }
256   
257   Int_t etamax = dNdEtadPhiHist->GetNbinsX();
258   Int_t phimax = dNdEtadPhiHist->GetNbinsY();
259   
260   TH1D *hprod = dNdEtadPhiHist->ProjectionX();
261   
262   TH1D *hESDMultvsEta = (TH1D*)fOutputList->FindObject(Form("hESDMultvsEta_binning%d" ,4));
263   hESDMultvsEta->Add(hprod);
264   
265   for (Int_t etabin = 1; etabin <= etamax; etabin++) {
266     Float_t val = 0;
267     Float_t eta = dNdEtadPhiHist->GetXaxis()->GetBinCenter(etabin);
268     
269     
270     for (Int_t i = 1; i <= 4; i++) {
271        TH1D *hESDMult = (TH1D*)fInternalList->FindObject(Form("hESDMult_binning%d", i));
272        hESDMult->Fill(eta, hprod->GetBinContent(etabin));
273     }
274   }
275   
276   //hESDMult->Add(hprod);
277   
278 }
279 //_____________________________________________________________________
280 void AliFMDAnalysisTaskBFCorrelation::CalculateValues(TString type) {
281
282   type.ToUpper();
283
284   const char *stype = type.Data(); 
285
286   TString sinput;
287   TString soutput;
288   TString soutputC;
289   TString snFnB;
290   TString snF2;
291   TString snB2;
292   TString snF;
293   TString snB;
294  
295   if (type == "ESD") {
296     snFnB  = "pnFnB";
297     snF2   = "pnF2";
298     snB2   = "pnB2";
299     snF    = "pnF";
300     snB    = "pnB";
301   } else {
302     snFnB  = "pnFnB_MC";
303     snF2   = "pnF2_MC";
304     snB2   = "pnB2_MC";
305     snF    = "pnF_MC";
306     snB    = "pnB_MC";
307   }
308
309   for (Int_t i = 1; i <= 4; i++) {
310     
311     sinput   = Form("h%sMult_binning%d", stype, i);
312     soutput  = Form("h%sMultvsEta_binning%d", stype, i);
313     soutputC = Form("h%sMultvsEtaC_binning%d", stype, i);
314     
315     TH1D *input  = (TH1D*)fInternalList->FindObject(sinput); 
316     TH1D *output = (TH1D*)fOutputList->FindObject(soutput);
317     TH1D *outputC = (TH1D*)fOutputList->FindObject(soutputC);
318
319     //output->Add(input);
320     
321     for (Int_t bin = 1; bin <= input->GetNbinsX()/2; bin++) {
322       
323       Double_t eta    = input->GetBinCenter(bin);
324       Double_t bwd    = input->GetBinContent(bin);
325       Double_t bwd2   = TMath::Power(bwd, 2);
326       Double_t fwd    = input->GetBinContent(input->FindBin(TMath::Abs(eta)));
327       Double_t fwd2   = TMath::Power(fwd, 2);
328       Double_t fwdbwd = fwd*bwd; 
329
330       snFnB += Form("_binning%d", i);
331       snF2  += Form("_binning%d", i);
332       snB2  += Form("_binning%d", i);
333       snF   += Form("_binning%d", i);
334       snB   += Form("_binning%d", i);
335       
336       TH1D *nFnB  = static_cast<TH1D*>(fOutputList->FindObject(snFnB)); 
337       TH1D *nF2   = static_cast<TH1D*>(fOutputList->FindObject(snF2)); 
338       TH1D *nB2   = static_cast<TH1D*>(fOutputList->FindObject(snB2)); 
339       TH1D *nF    = static_cast<TH1D*>(fOutputList->FindObject(snF)); 
340       TH1D *nB    = static_cast<TH1D*>(fOutputList->FindObject(snB)); 
341       
342       nFnB->Fill(TMath::Abs(eta), fwdbwd);
343       nF2->Fill(TMath::Abs(eta),  fwd2);
344       nB2->Fill(TMath::Abs(eta),  bwd2);
345       nF->Fill(TMath::Abs(eta),   fwd);
346       nB->Fill(TMath::Abs(eta),   bwd);
347       
348       outputC->Fill(TMath::Abs(eta), fwd);
349       outputC->Fill(eta, bwd);
350
351       if (type == "ESD") {
352         
353         snFnB.Remove(5,9);
354         snF2.Remove(4,9);
355         snB2.Remove(4,9);
356         snF.Remove(3,9);
357         snB.Remove(3,9);
358         TH2D *hESDBwdDist = (TH2D*)fOutputList->FindObject(Form("hESDBwdDist_binning%d", i));
359         TH2D *hESDBwd2Dist = (TH2D*)fOutputList->FindObject(Form("hESDBwd2Dist_binning%d", i));
360         hESDBwdDist->Fill(TMath::Abs(eta), bwd);
361         hESDBwd2Dist->Fill(TMath::Abs(eta), bwd2);
362
363       } else {
364         
365         snFnB.Remove(8,9);
366         snF2.Remove(7,9);
367         snB2.Remove(7,9);
368         snF.Remove(6,9);
369         snB.Remove(6,9);
370         TH2D *hMCBwdDist = (TH2D*)fOutputList->FindObject(Form("hMCBwdDist_binning%d", i));
371         TH2D *hMCBwd2Dist = (TH2D*)fOutputList->FindObject(Form("hMCBwd2Dist_binning%d", i));
372         hMCBwdDist->Fill(TMath::Abs(eta), bwd);
373         hMCBwd2Dist->Fill(TMath::Abs(eta), bwd2);
374       }
375     }
376   }
377 }
378
379 //_____________________________________________________________________
380 void AliFMDAnalysisTaskBFCorrelation::SetBounds() {
381   
382   TH2D *hTemp = (TH2D*)fInputList->FindObject("multTrVtx_FMD1I_vtxbin0");
383   
384   fnBinsX = hTemp->GetNbinsX();
385   fnBinsY = hTemp->GetNbinsY();
386   fXmin   = hTemp->GetXaxis()->GetXmin();
387   fXmax   = hTemp->GetXaxis()->GetXmax();
388   fYmin   = hTemp->GetYaxis()->GetXmin();
389   fYmax   = hTemp->GetYaxis()->GetXmax();
390 }  
391
392 //_____________________________________________________________________
393 void AliFMDAnalysisTaskBFCorrelation::Terminate(Option_t */*option*/) {
394   
395 }
396 //_____________________________________________________________________
397 void AliFMDAnalysisTaskBFCorrelation::ProcessPrimary() {
398   
399   AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
400   AliMCEvent* mcEvent = eventHandler->MCEvent();
401   if(!mcEvent)
402     return;
403     
404   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
405   
406   AliMCParticle* particle = 0;
407   AliStack* stack = mcEvent->Stack();
408   
409   TH1F* hPrimary = (TH1F*)fOutputList->FindObject("hMultvsEta");
410   AliHeader* header            = mcEvent->Header();
411   AliGenEventHeader* genHeader = header->GenEventHeader();
412   
413   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
414   
415   if (!pythiaGenHeader) {
416     std::cout<<" no pythia header!"<<std::endl;
417     return; 
418   }
419
420   Int_t pythiaType = pythiaGenHeader->ProcessType();
421   
422   if(pythiaType==92||pythiaType==93){
423     std::cout<<"single diffractive"<<std::endl;
424     return;
425   }
426   if(pythiaType==94){
427     std::cout<<"double diffractive"<<std::endl;
428     return;
429   }
430
431
432   TArrayF vertex;
433   genHeader->PrimaryVertex(vertex);
434   if(TMath::Abs(vertex.At(2)) > pars->GetVtxCutZ())
435     return;
436   //Double_t delta           = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
437   //Double_t vertexBinDouble = (vertex.At(2) + pars->GetVtxCutZ()) / delta;
438   //Int_t    vertexBin       = (Int_t)vertexBinDouble;
439     
440   //Bool_t firstTrack = kTRUE;
441   
442   // we loop over the primaries only unless we need the hits (diagnostics running slowly)
443   Int_t nTracks = stack->GetNprimary();
444   //  if(pars->GetProcessHits())
445   //  nTracks = stack->GetNtrack();
446   TH1D *hMCMultvsEta = (TH1D*)fOutputList->FindObject(Form("hMCMultvsEta_binning%d" ,4));  
447   
448   for (Int_t i = 1; i <= 4; i++) {
449     
450     TH1D *hMCMult = (TH1D*)fInternalList->FindObject(Form("hMCMult_binning%d", i));
451     hMCMult->Reset();
452
453     for(Int_t ii = 0 ;ii<nTracks;ii++) {
454       particle = (AliMCParticle*) mcEvent->GetTrack(ii);
455       if(!particle)
456       continue;
457       
458       if(stack->IsPhysicalPrimary(ii) && particle->Charge() != 0) {
459         if (i == 1) 
460           hPrimary->Fill(particle->Eta());
461         if(i==4)
462           hMCMultvsEta->Fill(particle->Eta());
463         hMCMult->Fill(particle->Eta());
464       }
465     }
466   }
467   CalculateValues("MC");
468 }
469 //_____________________________________________________________________
470 //
471 //
472 // EOF