Place the config and root file at the right place
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFEVZEROEventPlane.cxx
1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 *                                                                        *
4 * Author: The ALICE Off-line Project.                                    *
5 * Contributors are mentioned in the code where appropriate.              *
6 *                                                                        *
7 * Permission to use, copy, modify and distribute this software and its   *
8 * documentation strictly for non-commercial purposes is hereby granted   *
9 * without fee, provided that the above copyright notice appears in all   *
10 * copies and that both the copyright notice and this permission notice   *
11 * appear in the supporting documentation. The authors make no claims     *
12 * about the suitability of this software for any purpose. It is          *
13 * provided "as is" without express or implied warranty.                  *
14 **************************************************************************/
15 //
16 //
17 // VZERO event plane task for 2010
18 // Gain equalization + Recentering 
19 // Need a root file
20 // 
21 // Author:
22 //   Raphaelle Bailhache <R.Bailhache@gsi.de>
23 //
24 //
25
26 #include "AliHFEVZEROEventPlane.h"
27
28 // AliRoot includes
29 #include "AliAnalysisManager.h"
30 #include "AliInputEventHandler.h"
31 #include "AliVEvent.h"
32 //#include "AliESDEvent.h"
33 #include "AliCentrality.h"
34 //#include "AliESDVZERO.h"
35 #include "AliVVZERO.h"
36 #include "TFile.h"
37 #include "AliOADBContainer.h"
38 #include "TH2F.h"
39 #include "TF1.h"
40 #include "TList.h"
41 #include "TChain.h"
42 #include "THnSparse.h"
43 #include "TString.h"
44 #include "TVector2.h"
45 #include "AliEventplane.h"
46
47
48 // STL includes
49 #include <iostream>
50 using namespace std;
51
52
53 //_____________________________________________________________________________
54 AliHFEVZEROEventPlane::AliHFEVZEROEventPlane():
55   TNamed(),
56   fEventPlaneV0A(-100.0),
57   fEventPlaneV0C(-100.0),
58   fEventPlaneV0(-100.0),
59   fESD(0),
60   fRun(-1),
61   fMultV0(0x0),
62   fV0Cpol(100),
63   fV0Apol(100),
64   fnamefile(TString("")),
65   fOutputList(0x0),
66   fMultV0Before(0x0),
67   fMultV0After(0x0)
68 {
69   //
70   // Default constructor (should not be used)
71   //
72
73   for(Int_t k = 0; k < fgknCentrBin; k++) {
74     for(Int_t iside = 0; iside < 2; iside++) {
75       for(Int_t icoord = 0; icoord < 2; icoord++) {
76         fQBefore[k][iside][icoord] = 0x0;
77         fQAfter[k][iside][icoord] = 0x0;
78       }
79     }
80   }
81   
82 }
83
84 //______________________________________________________________________________
85 AliHFEVZEROEventPlane::AliHFEVZEROEventPlane(const char *name, const Char_t *title):
86   TNamed(name,title),
87   fEventPlaneV0A(-100.0),
88   fEventPlaneV0C(-100.0),
89   fEventPlaneV0(-100.0),
90   fESD(0),
91   fRun(-1),
92   fMultV0(0x0),
93   fV0Cpol(100),
94   fV0Apol(100),
95   fnamefile(TString("")),
96   fOutputList(0x0),
97   fMultV0Before(0x0),
98   fMultV0After(0x0)
99 {
100   //
101   // Constructor
102   //
103
104   TString namelist("QA_hist_");
105   namelist += name;
106   fOutputList = new TList();
107   fOutputList->SetName((const char*)namelist);
108   fOutputList->SetOwner(kTRUE);
109
110   // Multiplicity
111   fMultV0Before = new TProfile("MultV0Before","",64,0,64);
112   fMultV0Before->Sumw2();
113   fMultV0After = new TProfile("MultV0After","",64,0,64);
114   fMultV0After->Sumw2();
115   fOutputList->Add(fMultV0Before);
116   fOutputList->Add(fMultV0After);
117
118   // Recentering
119   for(Int_t k = 0; k < fgknCentrBin; k++) {
120     for(Int_t iside = 0; iside < 2; iside++) {
121       for(Int_t icoord = 0; icoord < 2; icoord++) {
122         
123         TString namecontbefore("");
124         TString namecontafter("");
125         
126         if(iside==0 && icoord==0) {
127           namecontbefore += "hQxc2_";
128           namecontbefore += k;
129         }
130         else if(iside==1 && icoord==0){
131           namecontbefore += "hQxa2_";
132           namecontbefore += k;
133         }
134         else if(iside==0 && icoord==1){
135           namecontbefore += "hQyc2_";
136           namecontbefore += k;    
137         }
138         else if(iside==1 && icoord==1) {
139           namecontbefore += "hQya2_";
140           namecontbefore += k;
141         }
142         //
143         if(iside==0 && icoord==0){
144           namecontafter += "hQxc2_";
145           namecontafter += k;
146           namecontafter += "_after";
147         }
148         else if(iside==1 && icoord==0){
149           namecontafter += "hQxa2_";
150           namecontafter += k;
151           namecontafter += "_after";
152         }
153         else if(iside==0 && icoord==1) {
154           namecontafter += "hQyc2_";
155           namecontafter += k;
156           namecontafter += "_after";
157         }
158         else if(iside==1 && icoord==1) {
159           namecontafter += "hQya2_";
160           namecontafter += k;
161           namecontafter += "_after";    
162         }
163         //
164         fQBefore[k][iside][icoord] = new TH1F(((const char*)namecontbefore),"",1600,-800.0,800.0);
165         fQBefore[k][iside][icoord]->Sumw2();
166         fOutputList->Add(fQBefore[k][iside][icoord]);
167         fQAfter[k][iside][icoord] = new TH1F(((const char*)namecontafter),"",1600,-800.0,800.0);
168         fQAfter[k][iside][icoord]->Sumw2();
169         fOutputList->Add(fQAfter[k][iside][icoord]);
170         //
171       }
172     }
173   }
174   
175 }
176 //____________________________________________________________
177 AliHFEVZEROEventPlane::AliHFEVZEROEventPlane(const AliHFEVZEROEventPlane &ref):
178   TNamed(ref),
179   fEventPlaneV0A(ref.fEventPlaneV0A),
180   fEventPlaneV0C(ref.fEventPlaneV0C),
181   fEventPlaneV0(ref.fEventPlaneV0),
182   fESD(0),
183   fRun(-1),
184   fMultV0(0x0),
185   fV0Cpol(100),
186   fV0Apol(100),
187   fnamefile(ref.fnamefile),
188   fOutputList(0x0),
189   fMultV0Before(0x0),
190   fMultV0After(0x0)
191 {
192   //
193   // Copy Constructor
194   //
195   ref.Copy(*this);
196 }
197
198 //____________________________________________________________
199 AliHFEVZEROEventPlane &AliHFEVZEROEventPlane::operator=(const AliHFEVZEROEventPlane &ref){
200   //
201   // Assignment operator
202   //
203   if(this == &ref) 
204     ref.Copy(*this);
205   return *this;
206 }
207
208 //____________________________________________________________
209 void AliHFEVZEROEventPlane::Copy(TObject &o) const {
210   // 
211   // Copy into object o
212   //
213   AliHFEVZEROEventPlane &target = dynamic_cast<AliHFEVZEROEventPlane &>(o);
214
215   target.fEventPlaneV0A = fEventPlaneV0A;
216   target.fEventPlaneV0C = fEventPlaneV0C;
217   target.fEventPlaneV0 = fEventPlaneV0;
218   target.fnamefile = fnamefile;
219   
220 }
221
222 //__________________________________________________________________
223 AliHFEVZEROEventPlane::~AliHFEVZEROEventPlane(){
224   //
225   // Destruktor
226   //
227   if(fOutputList){
228     fOutputList->Delete();
229     delete fOutputList;
230   }
231   fOutputList = 0x0;
232   
233 }
234 //______________________________________________________________________________
235 void AliHFEVZEROEventPlane::ProcessEvent(AliVEvent *event) 
236 {
237   //
238   // Process the event
239   // 
240   
241   // Reset
242   fEventPlaneV0A = -100.0;
243   fEventPlaneV0C = -100.0;
244   fEventPlaneV0  = -100.0;
245   
246   //
247   Int_t run = event->GetRunNumber();
248   Bool_t ok = kTRUE;
249   if(run != fRun){
250     if(fMultV0) delete fMultV0;
251     fMultV0 = 0x0;
252     // Load the calibrations run dependent
253     if(!OpenInfoCalbration(run)) ok = kFALSE;
254     fRun=run;
255   }
256
257   if(ok) Analyze(event);
258
259 }
260
261 //________________________________________________________________________
262 void AliHFEVZEROEventPlane::Analyze(AliVEvent* esdEvent)
263 {  
264   //
265   // Do VZERO calibration + centering
266   //
267
268   if(!fMultV0) {
269     //printf("Did not find calibration VZERO\n");
270     return;
271   }
272   //printf("HERE!!!\n");
273
274   //Centrality
275   Float_t v0Centr  = -10.;
276   AliCentrality* centrality = esdEvent->GetCentrality();
277   if (centrality){
278     v0Centr  = centrality->GetCentralityPercentile("V0M");
279   }
280   
281   // Analyse only for 0-80% PbPb collisions
282   Int_t iC = -1;    
283   if (v0Centr >0 && v0Centr < 80){
284     if(v0Centr < 5) iC = 0;
285     else if(v0Centr < 10) iC = 1;
286     else if(v0Centr < 20) iC = 2;
287     else if(v0Centr < 30) iC = 3;
288     else if(v0Centr < 40) iC = 4;
289     else if(v0Centr < 50) iC = 5;
290     else if(v0Centr < 60) iC = 6;
291     else if(v0Centr < 70) iC = 7;
292     else iC = 8;
293     
294     //general   
295     Double_t qxa2 = 0, qya2 = 0;
296     Double_t qxc2 = 0, qyc2 = 0;
297     
298     //V0 info    
299     AliVVZERO* esdV0 = esdEvent->GetVZEROData();
300     
301     for (Int_t iv0 = 0; iv0 < 64; iv0++) {
302       Double_t phiV0 = TMath::PiOver4()*(0.5 + iv0 % 8);
303       Float_t multv0 = esdV0->GetMultiplicity(iv0);
304       if (iv0 < 32){
305         //printf("Value %f\n",fMultV0->GetBinContent(iv0+1));
306         if(fMultV0->GetBinContent(iv0+1)>0.0) {
307           qxc2 += TMath::Cos(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
308           qyc2 += TMath::Sin(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
309           fMultV0After->Fill(iv0,multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1));
310         }
311       } else {
312         if(fMultV0->GetBinContent(iv0+1)>0.0) {
313           qxa2 += TMath::Cos(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
314           qya2 += TMath::Sin(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
315           fMultV0After->Fill(iv0,multv0*fV0Apol/fMultV0->GetBinContent(iv0+1));
316         }
317       }
318       fMultV0Before->Fill(iv0,multv0);
319     }
320
321     // Fill histos
322     fQBefore[iC][1][0] -> Fill(qxa2);
323     fQBefore[iC][1][1] -> Fill(qya2);
324     fQBefore[iC][0][0] -> Fill(qxc2);
325     fQBefore[iC][0][1] -> Fill(qyc2);
326     
327     //grab for each centrality the proper histo with the Qx and Qy to do the recentering
328     Double_t qxamean2 = fMeanQ[iC][1][0];
329     Double_t qxarms2  = fWidthQ[iC][1][0];
330     Double_t qyamean2 = fMeanQ[iC][1][1];
331     Double_t qyarms2  = fWidthQ[iC][1][1];
332     
333     Double_t qxcmean2 = fMeanQ[iC][0][0];
334     Double_t qxcrms2  = fWidthQ[iC][0][0];
335     Double_t qycmean2 = fMeanQ[iC][0][1];
336     Double_t qycrms2  = fWidthQ[iC][0][1];      
337
338     if((TMath::Abs(qxarms2) < 0.00000001) || (TMath::Abs(qyarms2) < 0.00000001) || (TMath::Abs(qxcrms2) < 0.00000001) || (TMath::Abs(qycrms2) < 0.00000001)) return;    
339
340     Double_t qxaCor2 = (qxa2 - qxamean2)/qxarms2;
341     Double_t qyaCor2 = (qya2 - qyamean2)/qyarms2;
342     Double_t qxcCor2 = (qxc2 - qxcmean2)/qxcrms2;
343     Double_t qycCor2 = (qyc2 - qycmean2)/qycrms2;
344
345     Double_t qxCor2 = qxaCor2 + qxcCor2;
346     Double_t qyCor2 = qyaCor2 + qycCor2;
347
348     // Fill histos
349     fQAfter[iC][1][0] -> Fill(qxaCor2);
350     fQAfter[iC][1][1] -> Fill(qyaCor2);
351     fQAfter[iC][0][0] -> Fill(qxcCor2);
352     fQAfter[iC][0][1] -> Fill(qycCor2);
353     
354     Double_t evPlAngV0ACor2 = TVector2::Phi_0_2pi(TMath::ATan2(qyaCor2, qxaCor2))/2.;
355     Double_t evPlAngV0CCor2 = TVector2::Phi_0_2pi(TMath::ATan2(qycCor2, qxcCor2))/2.;
356     Double_t evPlAngV0Cor2  = TVector2::Phi_0_2pi(TMath::ATan2(qyCor2, qxCor2))/2.;
357
358     fEventPlaneV0A = evPlAngV0ACor2;
359     fEventPlaneV0C = evPlAngV0CCor2;
360     fEventPlaneV0  = evPlAngV0Cor2;
361     
362     //printf("Eventplane V0A %f, V0C %f\n",fEventPlaneV0A,fEventPlaneV0C);
363
364   }
365 }
366 //_____________________________________________________________________________
367 Bool_t AliHFEVZEROEventPlane::OpenInfoCalbration(Int_t run)
368 {
369   //
370   // Take the calibration coefficients
371   //  
372
373   //printf("Name of the file %s\n",(const char*)fnamefile);
374  
375   TFile *foadb = TFile::Open(fnamefile.Data());
376   if(!foadb){
377     printf("OADB file %s cannot be opened\n",fnamefile.Data());
378     return kFALSE;
379   }
380
381   //printf("test\n");
382
383   AliOADBContainer *cont = (AliOADBContainer*) foadb->Get("hMultV0BefCorr");
384   if(!cont){
385     printf("OADB object hMultV0BefCorr is not available in the file\n");
386     return kFALSE;      
387   }
388
389   if(!(cont->GetObject(run))){
390     printf("OADB object hMultV0BefCorr is not available for run %i\n",run);
391     return kFALSE;      
392   }
393   //TProfile *multV0 = ((TH2F *) cont->GetObject(run))->ProfileX();
394   TProfile *multV0 = ((TProfile *) cont->GetObject(run));
395   if(!multV0) return kFALSE;
396   fMultV0 = (TProfile *) multV0->Clone();
397   fMultV0->SetDirectory(0);
398   
399   TF1 *fpol0 = new TF1("fpol0","pol0"); 
400   fMultV0->Fit(fpol0,"Q0","",0,31);
401   fV0Cpol = fpol0->GetParameter(0);
402   fMultV0->Fit(fpol0,"Q0","",32,64);
403   fV0Apol = fpol0->GetParameter(0);
404   
405   for(Int_t iside=0;iside<2;iside++){
406     for(Int_t icoord=0;icoord<2;icoord++){
407       for(Int_t i=0;i  < fgknCentrBin;i++){
408
409         TString namecont("");
410         if(iside==0 && icoord==0) {
411           namecont += "hQxc2_";
412           namecont += i;
413         }
414         else if(iside==1 && icoord==0){
415           namecont += "hQxa2_";
416           namecont += i;
417         }
418         else if(iside==0 && icoord==1){
419           namecont += "hQyc2_";
420           namecont += i;
421         }
422         else if(iside==1 && icoord==1) {
423           namecont += "hQya2_";
424           namecont += i;
425         }
426         
427         cont = (AliOADBContainer*) foadb->Get(namecont);
428         if(!cont){
429           printf("OADB object %s is not available in the file\n",(const char*)namecont);
430           delete fpol0;
431           return kFALSE;        
432         }
433         
434         if(!(cont->GetObject(run))){
435           printf("OADB object %s is not available for run %i\n",(const char*)namecont,run);
436           delete fpol0;
437           return kFALSE;        
438         }
439         fMeanQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
440         fWidthQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
441       }
442     }
443   }
444   
445   delete fpol0;
446   foadb->Close();
447   
448   return kTRUE;
449
450 }