don't lie in the log!
[u/mrichter/AliRoot.git] / PWGPP / VZERO / AliAnaVZEROEPFlatenning.cxx
1 #include <stdio.h>
2 #include <stdlib.h>
3
4 #include "TH2F.h"
5 #include "TFile.h"
6 #include "TProfile.h"
7 #include "TList.h"
8 #include "TChain.h"
9
10 #include "AliAnalysisTask.h"
11 #include "AliAnalysisManager.h"
12
13 #include "AliVEvent.h"
14 #include "AliESDEvent.h"
15 #include "AliAODEvent.h"
16 #include "AliInputEventHandler.h"
17 #include "AliAnaVZEROEPFlatenning.h"
18 #include "AliCentrality.h"
19 #include "AliEventplane.h"
20
21 // VZERO includes
22 #include "AliVVZERO.h"
23
24 ClassImp(AliAnaVZEROEPFlatenning)
25
26 AliAnaVZEROEPFlatenning::AliAnaVZEROEPFlatenning() 
27   : AliAnalysisTaskSE("AliAnaVZEROEPFlatenning"), fEvent(0), fOutputList(0),
28   fMBTrigName("CPBI"),
29   fUsePhysSel(kFALSE),
30   fPsiARawCentr(NULL),
31   fPsiAFlatCentr(NULL),
32   fPsiCRawCentr(NULL),
33   fPsiCFlatCentr(NULL),
34   fPsiACRawCentr(NULL),
35   fPsiACFlatCentr(NULL)
36 {
37   // Default constructor
38   // Init pointers
39   for(Int_t i = 0; i < 11; ++i) fX2[i] = fY2[i] = fX2Y2[i] = fCos8Psi[i] = NULL;
40   for(Int_t i = 0; i < 8; ++i) fX2In[i] = fY2In[i] = fX2Y2In[i] = fCos8PsiIn[i] = NULL;
41   for(Int_t i = 0; i < 8; ++i) fPsiRingRawCentr[i] = fPsiRingFlatCentr[i] = fPsiRingFlatFinalCentr[i] = NULL;
42   for(Int_t i = 0; i < 8; ++i) fC2[i] = fS2[i] = fC4[i] = fS4[i] = NULL;
43   for(Int_t i = 0; i < 11; ++i) fX2Corr[i] = fY2Corr[i] = fX2Y2Corr[i] = NULL;
44   // Define input and output slots here
45   // Input slot #0 works with a TChain
46   DefineInput(0, TChain::Class());
47   // Output slot #0 id reserved by the base class for AOD
48   // Output slot #1 writes into a TH1 container
49   DefineOutput(1, TList::Class());
50 }
51
52 //________________________________________________________________________
53 AliAnaVZEROEPFlatenning::AliAnaVZEROEPFlatenning(const char *name) 
54   : AliAnalysisTaskSE(name), fEvent(0), fOutputList(0),
55   fMBTrigName("CPBI"),
56   fUsePhysSel(kFALSE),
57   fPsiARawCentr(NULL),
58   fPsiAFlatCentr(NULL),
59   fPsiCRawCentr(NULL),
60   fPsiCFlatCentr(NULL),
61   fPsiACRawCentr(NULL),
62   fPsiACFlatCentr(NULL)
63 {
64   // Constructor
65   // Init pointers
66   for(Int_t i = 0; i < 11; ++i) fX2[i] = fY2[i] = fX2Y2[i] = fCos8Psi[i] = NULL;
67   for(Int_t i = 0; i < 8; ++i) fX2In[i] = fY2In[i] = fX2Y2In[i] = fCos8PsiIn[i] = NULL;
68   for(Int_t i = 0; i < 8; ++i) fPsiRingRawCentr[i] = fPsiRingFlatCentr[i] = fPsiRingFlatFinalCentr[i] = NULL;
69   for(Int_t i = 0; i < 8; ++i) fC2[i] = fS2[i] = fC4[i] = fS4[i] = NULL;
70   for(Int_t i = 0; i < 11; ++i) fX2Corr[i] = fY2Corr[i] = fX2Y2Corr[i] = NULL;
71   // Define input and output slots here
72   // Input slot #0 works with a TChain
73   DefineInput(0, TChain::Class());
74   // Output slot #0 id reserved by the base class for AOD
75   // Output slot #1 writes into a TH1 container
76   DefineOutput(1, TList::Class());
77 }
78
79 //________________________________________________________________________
80 void AliAnaVZEROEPFlatenning::SetInput(const char *filename)
81 {
82   // Read recentering
83   // parameters from an input file
84
85   AliInfo(Form("Reading input histograms from %s",filename));
86   TFile *f = TFile::Open(filename);
87   TList *list = (TList*)f->Get("coutput");
88
89   for(Int_t i = 0; i < 8; ++i) {
90     fX2In[i] = (TProfile*)list->FindObject(Form("fX2_%d",i))->Clone(Form("fX2In_%d",i));
91     fX2In[i]->SetDirectory(0);
92     fY2In[i] = (TProfile*)list->FindObject(Form("fY2_%d",i))->Clone(Form("fY2In_%d",i));
93     fY2In[i]->SetDirectory(0);
94     fX2Y2In[i] = (TProfile*)list->FindObject(Form("fX2Y2_%d",i))->Clone(Form("fX2Y2In_%d",i));
95     fX2Y2In[i]->SetDirectory(0);
96     fCos8PsiIn[i] = (TProfile*)list->FindObject(Form("fCos8Psi_%d",i))->Clone(Form("fCos8PsiIn_%d",i));
97     fCos8PsiIn[i]->SetDirectory(0);
98   }
99   f->Close();
100 }
101
102 //________________________________________________________________________
103 void AliAnaVZEROEPFlatenning::UserCreateOutputObjects()
104 {
105   // Create histograms
106   // Called once
107
108   fOutputList = new TList();
109   fOutputList->SetOwner(kTRUE);
110
111   for(Int_t i = 0; i < 11; ++i) {
112     fX2[i] = new TProfile(Form("fX2_%d",i),"",21,0,105,"s");
113     fOutputList->Add(fX2[i]);
114     fY2[i] = new TProfile(Form("fY2_%d",i),"",21,0,105,"s");
115     fOutputList->Add(fY2[i]);
116     fX2Y2[i] = new TProfile(Form("fX2Y2_%d",i),"",21,0,105,"s");
117     fOutputList->Add(fX2Y2[i]);
118     fCos8Psi[i] = new TProfile(Form("fCos8Psi_%d",i),"",21,0,105,"s");
119     fOutputList->Add(fCos8Psi[i]);
120   }
121   for(Int_t i = 0; i < 11; ++i) {
122     fX2Corr[i] = new TProfile(Form("fX2Corr_%d",i),"",21,0,105,"s");
123     fOutputList->Add(fX2Corr[i]);
124     fY2Corr[i] = new TProfile(Form("fY2Corr_%d",i),"",21,0,105,"s");
125     fOutputList->Add(fY2Corr[i]);
126     fX2Y2Corr[i] = new TProfile(Form("fX2Y2Corr_%d",i),"",21,0,105,"s");
127     fOutputList->Add(fX2Y2Corr[i]);
128   }
129
130   for(Int_t i = 0; i < 8; ++i) {
131     fPsiRingRawCentr[i] = new TH2F(Form("fPsiRingRawCentr_%d",i),"",105,0,105,100,-TMath::Pi()/2,TMath::Pi()/2);
132     fOutputList->Add(fPsiRingRawCentr[i]);
133     fPsiRingFlatCentr[i] = new TH2F(Form("fPsiRingFlatCentr_%d",i),"",105,0,105,100,-TMath::Pi()/2,TMath::Pi()/2);
134     fOutputList->Add(fPsiRingFlatCentr[i]);
135     fPsiRingFlatFinalCentr[i] = new TH2F(Form("fPsiRingFlatFinalCentr_%d",i),"",105,0,105,100,-TMath::Pi()/2,TMath::Pi()/2);
136     fOutputList->Add(fPsiRingFlatFinalCentr[i]);
137   }
138   fPsiARawCentr = new TH2F("fPsiARawCentr","",105,0,105,100,-TMath::Pi()/2,TMath::Pi()/2);
139   fOutputList->Add(fPsiARawCentr);
140   fPsiAFlatCentr = new TH2F("fPsiAFlatCentr","",105,0,105,100,-TMath::Pi()/2,TMath::Pi()/2);
141   fOutputList->Add(fPsiAFlatCentr);
142   fPsiCRawCentr = new TH2F("fPsiCRawCentr","",105,0,105,100,-TMath::Pi()/2,TMath::Pi()/2);
143   fOutputList->Add(fPsiCRawCentr);
144   fPsiCFlatCentr = new TH2F("fPsiCFlatCentr","",105,0,105,100,-TMath::Pi()/2,TMath::Pi()/2);
145   fOutputList->Add(fPsiCFlatCentr);
146   fPsiACRawCentr = new TH2F("fPsiACRawCentr","",105,0,105,100,-TMath::Pi()/2,TMath::Pi()/2);
147   fOutputList->Add(fPsiACRawCentr);
148   fPsiACFlatCentr = new TH2F("fPsiACFlatCentr","",105,0,105,100,-TMath::Pi()/2,TMath::Pi()/2);
149   fOutputList->Add(fPsiACFlatCentr);
150
151   for(Int_t i = 0; i < 8; ++i) {
152     fC2[i] = new TProfile(Form("fC2_%d",i),"",21,0,105,"s");
153     fOutputList->Add(fC2[i]);
154     fS2[i] = new TProfile(Form("fS2_%d",i),"",21,0,105,"s");
155     fOutputList->Add(fS2[i]);
156     fC4[i] = new TProfile(Form("fC4_%d",i),"",21,0,105,"s");
157     fOutputList->Add(fC4[i]);
158     fS4[i] = new TProfile(Form("fS4_%d",i),"",21,0,105,"s");
159     fOutputList->Add(fS4[i]);
160   }
161
162   PostData(1, fOutputList);
163 }
164
165 //________________________________________________________________________
166 void AliAnaVZEROEPFlatenning::Init()
167 {
168   // Nothing here
169   // ....
170 }
171
172 //________________________________________________________________________
173 void AliAnaVZEROEPFlatenning::UserExec(Option_t *) 
174 {
175   // Main loop
176   // Called for each event
177
178
179   fEvent = InputEvent();
180   if (!fEvent) {
181     printf("ERROR: fEvent not available\n");
182     return;
183   }
184
185   AliVVZERO* esdV0 = fEvent->GetVZEROData();
186   if (!esdV0) {
187     Printf("ERROR: esd/aod V0  not available");
188     return;
189   }
190
191   // Phys sel
192   Bool_t goodEvent = kTRUE;
193   Bool_t isSelected;
194   if (fUsePhysSel)
195     isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kMB | AliVEvent::kSemiCentral | AliVEvent::kCentral));
196   else
197     isSelected = ((esdV0->GetV0ADecision()==1) && (esdV0->GetV0CDecision()==1));
198
199   if (!isSelected) goodEvent = kFALSE;
200
201   AliCentrality *centrality = fEvent->GetCentrality();
202   //  Float_t spdPercentile = centrality->GetCentralityPercentile("V0M");
203   Float_t spdPercentile = centrality->GetCentralityPercentile("CL1");
204
205   TString inputDataType = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->GetDataType(); // can be "ESD" or "AOD"
206   if (inputDataType == "ESD") {
207     AliESDEvent* esdEvent = static_cast<AliESDEvent*>(fEvent);
208     // Trigger
209     TString trigStr(esdEvent->GetFiredTriggerClasses());
210     if (!trigStr.Contains("-B-")) return;
211     if (!trigStr.Contains(fMBTrigName.Data())) return;
212
213     const AliESDVertex *primaryVtx = esdEvent->GetPrimaryVertexSPD();
214     if (!primaryVtx) goodEvent = kFALSE;
215     if (!primaryVtx->GetStatus()) goodEvent = kFALSE;
216     Double_t tPrimaryVtxPosition[3];
217     primaryVtx->GetXYZ(tPrimaryVtxPosition);
218     if (TMath::Abs(tPrimaryVtxPosition[2]) > 10.0) goodEvent = kFALSE;
219   }
220   else if (inputDataType == "AOD") {
221     AliAODEvent* aodEvent = static_cast<AliAODEvent*>(fEvent);
222     // Trigger
223     TString trigStr(aodEvent->GetHeader()->GetFiredTriggerClasses());
224     if (!trigStr.Contains("-B-")) return;
225     if (!trigStr.Contains(fMBTrigName.Data())) return;
226
227     const AliAODVertex *primaryVtx = aodEvent->GetPrimaryVertexSPD();
228     if (!primaryVtx) goodEvent = kFALSE;
229     if (primaryVtx->GetNContributors()==0) goodEvent = kFALSE;
230     Double_t tPrimaryVtxPosition[3];
231     primaryVtx->GetXYZ(tPrimaryVtxPosition);
232     if (TMath::Abs(tPrimaryVtxPosition[2]) > 10.0) goodEvent = kFALSE;
233   }
234   else {
235     AliFatal("Trying to get the vertex from neither ESD nor AOD event!");
236     return;
237   }
238
239   if (goodEvent) {
240     Double_t totMult = 0, totMultA = 0, totMultC = 0;
241     Double_t c2A = 0, c2C = 0, s2A = 0, s2C = 0;
242     Double_t qxA = 0, qyA = 0;
243     Double_t qxC = 0, qyC = 0;
244     for(Int_t iring = 0; iring < 8; ++iring) {
245       Double_t ringMult = 0;
246       Double_t c2 = 0, s2 = 0;
247       for(Int_t iCh = iring*8; iCh < (iring+1)*8; ++iCh) {
248         Double_t phi = TMath::Pi()/8. + TMath::Pi()/4.*(iCh%8);
249         c2 += fEvent->GetVZEROEqMultiplicity(iCh)*TMath::Cos(2.*phi);
250         s2 += fEvent->GetVZEROEqMultiplicity(iCh)*TMath::Sin(2.*phi);
251         if (fEvent->GetVZEROEqMultiplicity(iCh) > 1e-6) {
252           fC2[iring]->Fill(spdPercentile,TMath::Cos(2.*phi),fEvent->GetVZEROEqMultiplicity(iCh));
253           fS2[iring]->Fill(spdPercentile,TMath::Sin(2.*phi),fEvent->GetVZEROEqMultiplicity(iCh));
254           fC4[iring]->Fill(spdPercentile,TMath::Cos(4.*phi),fEvent->GetVZEROEqMultiplicity(iCh));
255           fS4[iring]->Fill(spdPercentile,TMath::Sin(4.*phi),fEvent->GetVZEROEqMultiplicity(iCh));
256         }
257         ringMult += fEvent->GetVZEROEqMultiplicity(iCh);
258       }
259       if (ringMult < 1e-6) continue;
260       totMult += ringMult;
261       if (iring >= 4) {
262         totMultA += ringMult;
263         c2A += c2;
264         s2A += s2;
265       }
266       else {
267         totMultC += ringMult;
268         c2C += c2;
269         s2C += s2;
270       }
271       fX2[iring]->Fill(spdPercentile,c2);
272       fY2[iring]->Fill(spdPercentile,s2);
273       fX2Y2[iring]->Fill(spdPercentile,c2*s2);
274
275       Double_t qxOut = 0, qyOut = 0;
276       Double_t psiRingRaw = fEvent->GetEventplane()->CalculateVZEROEventPlane(fEvent,iring,iring,2,qxOut,qyOut);
277       fPsiRingRawCentr[iring]->Fill(spdPercentile,psiRingRaw);
278       fCos8Psi[iring]->Fill(spdPercentile, 2./4.*TMath::Cos(2.*4.*psiRingRaw));
279       Double_t qxTierce = 0, qyTierce = 0;
280       Double_t psiRingFlat = CalculateVZEROEventPlane(fEvent,iring,spdPercentile,qxTierce,qyTierce);
281       fPsiRingFlatCentr[iring]->Fill(spdPercentile,psiRingFlat);
282       Int_t ibin = fCos8PsiIn[iring]->FindBin(spdPercentile);
283       Double_t psiRingFlatFinal = psiRingFlat + (fCos8PsiIn[iring]->GetBinContent(ibin)*TMath::Sin(2.*4.*psiRingFlat))/2.;
284       if (psiRingFlatFinal > TMath::Pi()/2) psiRingFlatFinal -= TMath::Pi();
285       if (psiRingFlatFinal <-TMath::Pi()/2) psiRingFlatFinal += TMath::Pi();
286       fPsiRingFlatFinalCentr[iring]->Fill(spdPercentile,psiRingFlatFinal);
287       fX2Corr[iring]->Fill(spdPercentile,qxTierce);
288       fY2Corr[iring]->Fill(spdPercentile,qyTierce);
289       fX2Y2Corr[iring]->Fill(spdPercentile,qxTierce*qyTierce);
290       if (iring >= 4) {
291         qxA += qxTierce;
292         qyA += qyTierce;
293       }
294       else {
295         qxC += qxTierce;
296         qyC += qyTierce;
297      }
298     }
299
300     if (totMultA > 1e-6) {
301       fX2[8]->Fill(spdPercentile,c2A);
302       fY2[8]->Fill(spdPercentile,s2A);
303       fX2Y2[8]->Fill(spdPercentile,c2A*s2A);
304       Double_t psiA = TMath::ATan2(qyA,qxA)/2.;
305       fPsiAFlatCentr->Fill(spdPercentile,psiA);
306       //      Double_t psiAOrg = fEvent->GetEventplane()->GetEventplane("V0A",fEvent,2);
307       Double_t qxEP = 0, qyEP = 0;
308       Double_t psiAOrg = fEvent->GetEventplane()->CalculateVZEROEventPlane(fEvent,8,2,qxEP,qyEP);
309       fX2Corr[8]->Fill(spdPercentile,qxEP);
310       fY2Corr[8]->Fill(spdPercentile,qyEP);
311       fX2Y2Corr[8]->Fill(spdPercentile,qxEP*qyEP);
312       fPsiARawCentr->Fill(spdPercentile,psiAOrg);
313       fCos8Psi[8]->Fill(spdPercentile, 2./4.*TMath::Cos(2.*4.*psiAOrg));
314     }
315     if (totMultC > 1e-6) {
316       fX2[9]->Fill(spdPercentile,c2C);
317       fY2[9]->Fill(spdPercentile,s2C);
318       fX2Y2[9]->Fill(spdPercentile,c2C*s2C);
319       Double_t psiC = TMath::ATan2(qyC,qxC)/2.;
320       fPsiCFlatCentr->Fill(spdPercentile,psiC);
321       //      Double_t psiCOrg = fEvent->GetEventplane()->GetEventplane("V0C",fEvent,2);
322       Double_t qxEP = 0, qyEP = 0;
323       Double_t psiCOrg = fEvent->GetEventplane()->CalculateVZEROEventPlane(fEvent,9,2,qxEP,qyEP);
324       fX2Corr[9]->Fill(spdPercentile,qxEP);
325       fY2Corr[9]->Fill(spdPercentile,qyEP);
326       fX2Y2Corr[9]->Fill(spdPercentile,qxEP*qyEP);
327       fPsiCRawCentr->Fill(spdPercentile,psiCOrg);
328       fCos8Psi[9]->Fill(spdPercentile, 2./4.*TMath::Cos(2.*4.*psiCOrg));
329    }
330     if (totMult > 1e-6) {
331       fX2[10]->Fill(spdPercentile,c2A+c2C);
332       fY2[10]->Fill(spdPercentile,s2A+s2C);
333       fX2Y2[10]->Fill(spdPercentile,(c2A+c2C)*(s2A+s2C));
334       Double_t psiAC = TMath::ATan2(qyA+qyC,qxA+qxC)/2.;
335       fPsiACFlatCentr->Fill(spdPercentile,psiAC);
336       //      Double_t psiACOrg = fEvent->GetEventplane()->GetEventplane("V0",fEvent,2);
337       Double_t qxEP = 0, qyEP = 0;
338       Double_t psiACOrg = fEvent->GetEventplane()->CalculateVZEROEventPlane(fEvent,10,2,qxEP,qyEP);
339       fX2Corr[10]->Fill(spdPercentile,qxEP);
340       fY2Corr[10]->Fill(spdPercentile,qyEP);
341       fX2Y2Corr[10]->Fill(spdPercentile,qxEP*qyEP);
342       fPsiACRawCentr->Fill(spdPercentile,psiACOrg);
343       fCos8Psi[10]->Fill(spdPercentile, 2./4.*TMath::Cos(2.*4.*psiACOrg));
344     }
345   }
346
347   PostData(1, fOutputList);
348 }      
349
350 //________________________________________________________________________
351 void AliAnaVZEROEPFlatenning::Terminate(Option_t *) 
352 {
353   // Check the output list and store output config file
354   // Called once at the end of the query
355
356   fOutputList = dynamic_cast<TList*> (GetOutputData(1));
357   if (!fOutputList) {
358     printf("ERROR: Output list not available\n");
359     return;
360   }
361 }
362
363 Double_t AliAnaVZEROEPFlatenning::CalculateVZEROEventPlane(const AliVEvent *  event, Int_t ring, Float_t centrality, Double_t &qxTierce, Double_t &qyTierce) const
364 {
365   // Calculate the VZERO event plane
366   // taking into account the recentering/twist and rescaling
367   // of the cumulants
368   qxTierce = qyTierce = 0.;
369   if(!event) {
370     AliError("No Event received");
371     return -1000.;
372   }
373   AliVVZERO *vzeroData = event->GetVZEROData();
374   if(!vzeroData) {
375     AliError("Enable to get VZERO Data");
376     return -1000.;
377   }
378
379   Int_t ibin = fX2In[ring]->FindBin(centrality);
380   Double_t meanX2 = fX2In[ring]->GetBinContent(ibin);
381   Double_t meanY2 = fY2In[ring]->GetBinContent(ibin);
382   Double_t sigmaX2 = fX2In[ring]->GetBinError(ibin);
383   Double_t sigmaY2 = fY2In[ring]->GetBinError(ibin);
384   Double_t rho = (fX2Y2In[ring]->GetBinContent(ibin)-meanX2*meanY2)/sigmaX2/sigmaY2;
385
386   Double_t b = rho*sigmaX2*sigmaY2*
387     TMath::Sqrt(2.*(sigmaX2*sigmaX2+sigmaY2*sigmaY2-2.*sigmaX2*sigmaY2*TMath::Sqrt(1.-rho*rho))/
388                 ((sigmaX2*sigmaX2-sigmaY2*sigmaY2)*(sigmaX2*sigmaX2-sigmaY2*sigmaY2)+
389                  4.*sigmaX2*sigmaX2*sigmaY2*sigmaY2*rho*rho));
390   Double_t aPlus = TMath::Sqrt(2.*sigmaX2*sigmaX2-b*b);
391   Double_t aMinus= TMath::Sqrt(2.*sigmaY2*sigmaY2-b*b);
392
393   Double_t lambdaPlus = b/aPlus;
394   Double_t lambdaMinus = b/aMinus;
395
396   Double_t qx=0., qy=0.;
397   for(Int_t iCh = ring*8; iCh < (ring+1)*8; ++iCh) {
398     Double_t phi = TMath::Pi()/8. + TMath::Pi()/4.*(iCh%8);
399     Double_t mult = event->GetVZEROEqMultiplicity(iCh);
400     qx += mult*TMath::Cos(2.*phi);
401     qy += mult*TMath::Sin(2.*phi);
402   }
403   // Recentering
404   Double_t qxPrime = qx - meanX2;
405   Double_t qyPrime = qy - meanY2;
406   // Twist
407   Double_t trans[2][2];
408   trans[0][0] = 1./(1.-lambdaPlus*lambdaMinus);
409   trans[0][1] = -lambdaMinus/(1.-lambdaPlus*lambdaMinus);
410   trans[1][0] = -lambdaPlus/(1.-lambdaPlus*lambdaMinus);
411   trans[1][1] = 1./(1.-lambdaPlus*lambdaMinus);
412   Double_t qxSeconde = trans[0][0]*qxPrime + trans[0][1]*qyPrime;
413   Double_t qySeconde = trans[1][0]*qxPrime + trans[1][1]*qyPrime;
414   // Rescaling
415   qxTierce = qxSeconde/aPlus;
416   qyTierce = qySeconde/aMinus;
417
418   return (TMath::ATan2(qyTierce,qxTierce)/2.);
419 }