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 }