don't lie in the log!
[u/mrichter/AliRoot.git] / PWGPP / VZERO / AliAnaVZEROEPFlatenning.cxx
CommitLineData
345e2385 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
6b4456fd 13#include "AliVEvent.h"
345e2385 14#include "AliESDEvent.h"
6b4456fd 15#include "AliAODEvent.h"
16#include "AliInputEventHandler.h"
345e2385 17#include "AliAnaVZEROEPFlatenning.h"
18#include "AliCentrality.h"
19#include "AliEventplane.h"
20
21// VZERO includes
6b4456fd 22#include "AliVVZERO.h"
345e2385 23
24ClassImp(AliAnaVZEROEPFlatenning)
25
26AliAnaVZEROEPFlatenning::AliAnaVZEROEPFlatenning()
6b4456fd 27 : AliAnalysisTaskSE("AliAnaVZEROEPFlatenning"), fEvent(0), fOutputList(0),
345e2385 28 fMBTrigName("CPBI"),
fe509094 29 fUsePhysSel(kFALSE),
ed806549 30 fPsiARawCentr(NULL),
31 fPsiAFlatCentr(NULL),
32 fPsiCRawCentr(NULL),
33 fPsiCFlatCentr(NULL),
34 fPsiACRawCentr(NULL),
35 fPsiACFlatCentr(NULL)
345e2385 36{
37 // Default constructor
38 // Init pointers
ed806549 39 for(Int_t i = 0; i < 11; ++i) fX2[i] = fY2[i] = fX2Y2[i] = fCos8Psi[i] = NULL;
345e2385 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;
fe509094 42 for(Int_t i = 0; i < 8; ++i) fC2[i] = fS2[i] = fC4[i] = fS4[i] = NULL;
ed806549 43 for(Int_t i = 0; i < 11; ++i) fX2Corr[i] = fY2Corr[i] = fX2Y2Corr[i] = NULL;
345e2385 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//________________________________________________________________________
53AliAnaVZEROEPFlatenning::AliAnaVZEROEPFlatenning(const char *name)
6b4456fd 54 : AliAnalysisTaskSE(name), fEvent(0), fOutputList(0),
345e2385 55 fMBTrigName("CPBI"),
fe509094 56 fUsePhysSel(kFALSE),
ed806549 57 fPsiARawCentr(NULL),
58 fPsiAFlatCentr(NULL),
59 fPsiCRawCentr(NULL),
60 fPsiCFlatCentr(NULL),
61 fPsiACRawCentr(NULL),
62 fPsiACFlatCentr(NULL)
345e2385 63{
64 // Constructor
65 // Init pointers
ed806549 66 for(Int_t i = 0; i < 11; ++i) fX2[i] = fY2[i] = fX2Y2[i] = fCos8Psi[i] = NULL;
345e2385 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;
fe509094 69 for(Int_t i = 0; i < 8; ++i) fC2[i] = fS2[i] = fC4[i] = fS4[i] = NULL;
ed806549 70 for(Int_t i = 0; i < 11; ++i) fX2Corr[i] = fY2Corr[i] = fX2Y2Corr[i] = NULL;
345e2385 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//________________________________________________________________________
80void 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//________________________________________________________________________
103void AliAnaVZEROEPFlatenning::UserCreateOutputObjects()
104{
105 // Create histograms
106 // Called once
107
108 fOutputList = new TList();
109 fOutputList->SetOwner(kTRUE);
110
ed806549 111 for(Int_t i = 0; i < 11; ++i) {
345e2385 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 }
ed806549 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 }
345e2385 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 }
ed806549 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);
345e2385 150
fe509094 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
345e2385 162 PostData(1, fOutputList);
163}
164
165//________________________________________________________________________
166void AliAnaVZEROEPFlatenning::Init()
167{
168 // Nothing here
169 // ....
170}
171
172//________________________________________________________________________
173void AliAnaVZEROEPFlatenning::UserExec(Option_t *)
174{
175 // Main loop
176 // Called for each event
177
178
6b4456fd 179 fEvent = InputEvent();
180 if (!fEvent) {
181 printf("ERROR: fEvent not available\n");
345e2385 182 return;
183 }
184
6b4456fd 185 AliVVZERO* esdV0 = fEvent->GetVZEROData();
345e2385 186 if (!esdV0) {
6b4456fd 187 Printf("ERROR: esd/aod V0 not available");
345e2385 188 return;
189 }
190
345e2385 191 // Phys sel
192 Bool_t goodEvent = kTRUE;
193 Bool_t isSelected;
194 if (fUsePhysSel)
0bd9e4f8 195 isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kMB | AliVEvent::kSemiCentral | AliVEvent::kCentral));
345e2385 196 else
197 isSelected = ((esdV0->GetV0ADecision()==1) && (esdV0->GetV0CDecision()==1));
198
199 if (!isSelected) goodEvent = kFALSE;
200
6b4456fd 201 AliCentrality *centrality = fEvent->GetCentrality();
0bd9e4f8 202 // Float_t spdPercentile = centrality->GetCentralityPercentile("V0M");
345e2385 203 Float_t spdPercentile = centrality->GetCentralityPercentile("CL1");
204
6b4456fd 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 }
345e2385 238
239 if (goodEvent) {
ed806549 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;
345e2385 244 for(Int_t iring = 0; iring < 8; ++iring) {
ed806549 245 Double_t ringMult = 0;
246 Double_t c2 = 0, s2 = 0;
345e2385 247 for(Int_t iCh = iring*8; iCh < (iring+1)*8; ++iCh) {
248 Double_t phi = TMath::Pi()/8. + TMath::Pi()/4.*(iCh%8);
6b4456fd 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));
fe509094 256 }
ed806549 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;
345e2385 270 }
345e2385 271 fX2[iring]->Fill(spdPercentile,c2);
272 fY2[iring]->Fill(spdPercentile,s2);
273 fX2Y2[iring]->Fill(spdPercentile,c2*s2);
274
fe509094 275 Double_t qxOut = 0, qyOut = 0;
6b4456fd 276 Double_t psiRingRaw = fEvent->GetEventplane()->CalculateVZEROEventPlane(fEvent,iring,iring,2,qxOut,qyOut);
345e2385 277 fPsiRingRawCentr[iring]->Fill(spdPercentile,psiRingRaw);
ed806549 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);
345e2385 282 Int_t ibin = fCos8PsiIn[iring]->FindBin(spdPercentile);
ed806549 283 Double_t psiRingFlatFinal = psiRingFlat + (fCos8PsiIn[iring]->GetBinContent(ibin)*TMath::Sin(2.*4.*psiRingFlat))/2.;
345e2385 284 if (psiRingFlatFinal > TMath::Pi()/2) psiRingFlatFinal -= TMath::Pi();
285 if (psiRingFlatFinal <-TMath::Pi()/2) psiRingFlatFinal += TMath::Pi();
286 fPsiRingFlatFinalCentr[iring]->Fill(spdPercentile,psiRingFlatFinal);
ed806549 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);
ed806549 304 Double_t psiA = TMath::ATan2(qyA,qxA)/2.;
305 fPsiAFlatCentr->Fill(spdPercentile,psiA);
0bd9e4f8 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);
ed806549 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);
ed806549 319 Double_t psiC = TMath::ATan2(qyC,qxC)/2.;
320 fPsiCFlatCentr->Fill(spdPercentile,psiC);
0bd9e4f8 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);
ed806549 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));
ed806549 334 Double_t psiAC = TMath::ATan2(qyA+qyC,qxA+qxC)/2.;
335 fPsiACFlatCentr->Fill(spdPercentile,psiAC);
0bd9e4f8 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);
ed806549 342 fPsiACRawCentr->Fill(spdPercentile,psiACOrg);
343 fCos8Psi[10]->Fill(spdPercentile, 2./4.*TMath::Cos(2.*4.*psiACOrg));
345e2385 344 }
345 }
346
347 PostData(1, fOutputList);
348}
349
350//________________________________________________________________________
351void 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
fe509094 363Double_t AliAnaVZEROEPFlatenning::CalculateVZEROEventPlane(const AliVEvent * event, Int_t ring, Float_t centrality, Double_t &qxTierce, Double_t &qyTierce) const
345e2385 364{
365 // Calculate the VZERO event plane
366 // taking into account the recentering/twist and rescaling
367 // of the cumulants
fe509094 368 qxTierce = qyTierce = 0.;
345e2385 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
345e2385 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
ed806549 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);
345e2385 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
fe509094 415 qxTierce = qxSeconde/aPlus;
416 qyTierce = qySeconde/aMinus;
345e2385 417
418 return (TMath::ATan2(qyTierce,qxTierce)/2.);
419}