10 #include "AliAnalysisTask.h"
11 #include "AliAnalysisManager.h"
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"
22 #include "AliVVZERO.h"
24 ClassImp(AliAnaVZEROEPFlatenning)
26 AliAnaVZEROEPFlatenning::AliAnaVZEROEPFlatenning()
27 : AliAnalysisTaskSE("AliAnaVZEROEPFlatenning"), fEvent(0), fOutputList(0),
37 // Default constructor
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());
52 //________________________________________________________________________
53 AliAnaVZEROEPFlatenning::AliAnaVZEROEPFlatenning(const char *name)
54 : AliAnalysisTaskSE(name), fEvent(0), fOutputList(0),
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());
79 //________________________________________________________________________
80 void AliAnaVZEROEPFlatenning::SetInput(const char *filename)
83 // parameters from an input file
85 AliInfo(Form("Reading input histograms from %s",filename));
86 TFile *f = TFile::Open(filename);
87 TList *list = (TList*)f->Get("coutput");
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);
102 //________________________________________________________________________
103 void AliAnaVZEROEPFlatenning::UserCreateOutputObjects()
108 fOutputList = new TList();
109 fOutputList->SetOwner(kTRUE);
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]);
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]);
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]);
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);
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]);
162 PostData(1, fOutputList);
165 //________________________________________________________________________
166 void AliAnaVZEROEPFlatenning::Init()
172 //________________________________________________________________________
173 void AliAnaVZEROEPFlatenning::UserExec(Option_t *)
176 // Called for each event
179 fEvent = InputEvent();
181 printf("ERROR: fEvent not available\n");
185 AliVVZERO* esdV0 = fEvent->GetVZEROData();
187 Printf("ERROR: esd/aod V0 not available");
192 Bool_t goodEvent = kTRUE;
195 isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kMB | AliVEvent::kSemiCentral | AliVEvent::kCentral));
197 isSelected = ((esdV0->GetV0ADecision()==1) && (esdV0->GetV0CDecision()==1));
199 if (!isSelected) goodEvent = kFALSE;
201 AliCentrality *centrality = fEvent->GetCentrality();
202 // Float_t spdPercentile = centrality->GetCentralityPercentile("V0M");
203 Float_t spdPercentile = centrality->GetCentralityPercentile("CL1");
205 TString inputDataType = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->GetDataType(); // can be "ESD" or "AOD"
206 if (inputDataType == "ESD") {
207 AliESDEvent* esdEvent = static_cast<AliESDEvent*>(fEvent);
209 TString trigStr(esdEvent->GetFiredTriggerClasses());
210 if (!trigStr.Contains("-B-")) return;
211 if (!trigStr.Contains(fMBTrigName.Data())) return;
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;
220 else if (inputDataType == "AOD") {
221 AliAODEvent* aodEvent = static_cast<AliAODEvent*>(fEvent);
223 TString trigStr(aodEvent->GetHeader()->GetFiredTriggerClasses());
224 if (!trigStr.Contains("-B-")) return;
225 if (!trigStr.Contains(fMBTrigName.Data())) return;
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;
235 AliFatal("Trying to get the vertex from neither ESD nor AOD event!");
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));
257 ringMult += fEvent->GetVZEROEqMultiplicity(iCh);
259 if (ringMult < 1e-6) continue;
262 totMultA += ringMult;
267 totMultC += ringMult;
271 fX2[iring]->Fill(spdPercentile,c2);
272 fY2[iring]->Fill(spdPercentile,s2);
273 fX2Y2[iring]->Fill(spdPercentile,c2*s2);
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);
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));
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));
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));
347 PostData(1, fOutputList);
350 //________________________________________________________________________
351 void AliAnaVZEROEPFlatenning::Terminate(Option_t *)
353 // Check the output list and store output config file
354 // Called once at the end of the query
356 fOutputList = dynamic_cast<TList*> (GetOutputData(1));
358 printf("ERROR: Output list not available\n");
363 Double_t AliAnaVZEROEPFlatenning::CalculateVZEROEventPlane(const AliVEvent * event, Int_t ring, Float_t centrality, Double_t &qxTierce, Double_t &qyTierce) const
365 // Calculate the VZERO event plane
366 // taking into account the recentering/twist and rescaling
368 qxTierce = qyTierce = 0.;
370 AliError("No Event received");
373 AliVVZERO *vzeroData = event->GetVZEROData();
375 AliError("Enable to get VZERO Data");
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;
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);
393 Double_t lambdaPlus = b/aPlus;
394 Double_t lambdaMinus = b/aMinus;
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);
404 Double_t qxPrime = qx - meanX2;
405 Double_t qyPrime = qy - meanY2;
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;
415 qxTierce = qxSeconde/aPlus;
416 qyTierce = qySeconde/aMinus;
418 return (TMath::ATan2(qyTierce,qxTierce)/2.);