PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetTriggerQA.cxx

   1 //
   2 // Jet trigger QA analysis task.
   3 //
   4 // Author: M.Verweij
   5
   6 #include <TClonesArray.h>
   7 #include <TH1F.h>
   8 #include <TH2F.h>
   9 #include <TH3F.h>
  10 #include <THnSparse.h>
  11 #include <TList.h>
  12 #include <TLorentzVector.h>
  13
  14 #include "AliVCluster.h"
  15 #include "AliVTrack.h"
  16 #include "AliEmcalJet.h"
  17 #include "AliRhoParameter.h"
  18 #include "AliLog.h"
  19 #include "AliEmcalParticle.h"
  20 #include "AliAODCaloTrigger.h"
  21 #include "AliEMCALGeometry.h"
  22 #include "AliVCaloCells.h"
  23 #include "AliJetContainer.h"
  24 #include "AliClusterContainer.h"
  25 #include "AliParticleContainer.h"
  26 #include "AliEmcalTriggerPatchInfo.h"
  27 #include "AliAODHeader.h"
  28 #include "AliPicoTrack.h"
  29
  30 #include "AliAnalysisTaskEmcalJetTriggerQA.h"
  31
  32 ClassImp(AliAnalysisTaskEmcalJetTriggerQA)
  33
  34 //________________________________________________________________________
  35 AliAnalysisTaskEmcalJetTriggerQA::AliAnalysisTaskEmcalJetTriggerQA() :
  36   AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetTriggerQA", kTRUE),
  37   fDebug(kFALSE),
  38   fTriggerClass(""),
  39   fBitJ1((1<<8)),
  40   fBitJ2((1<<11)),
  41   fContainerFull(0),
  42   fContainerCharged(1),
  43   fMaxPatchEnergy(0),
  44   fTriggerType(-1),
  45   fNFastOR(16),
  46   fhNEvents(0),
  47   fHistRhovsCentFull(0),
  48   fHistRhovsCentCharged(0),
  49   fh3PtEtaPhiTracks(0),
  50   fh3PtEtaPhiTracksOnEmcal(0),
  51   fh3PtEtaPhiTracksProp(0),
  52   fh3PtEtaPhiJetFull(0),
  53   fh3PtEtaPhiJetCharged(0),
  54   fh2NJetsPtFull(0),
  55   fh2NJetsPtCharged(0),
  56   fh3PtEtaAreaJetFull(0),
  57   fh3PtEtaAreaJetCharged(0),
  58   fh2PtNConstituentsCharged(0),
  59   fh2PtNConstituents(0),
  60   fh2PtMeanPtConstituentsCharged(0),
  61   fh2PtMeanPtConstituentsNeutral(0),
  62   fh2PtNEF(0),
  63   fh3NEFEtaPhi(0),
  64   fh2NEFNConstituentsCharged(0),
  65   fh2NEFNConstituentsNeutral(0),
  66   fh2Ptz(0),
  67   fh2PtzCharged(0),
  68   fh2PtLeadJet1VsLeadJet2(0),
  69   fh3EEtaPhiCluster(0),
  70   fh3PtLeadJet1VsPatchEnergy(0),
  71   fh3PtLeadJet2VsPatchEnergy(0),
  72   fh3PatchEnergyEtaPhiCenterJ1(0),
  73   fh3PatchEnergyEtaPhiCenterJ2(0),
  74   fh3PatchEnergyEtaPhiCenterJ1J2(0),
  75   fh3PatchADCEnergyEtaPhiCenterJ1(0),
  76   fh3PatchADCEnergyEtaPhiCenterJ2(0),
  77   fh3PatchADCEnergyEtaPhiCenterJ1J2(0),
  78   fh2CellEnergyVsTime(0),
  79   fh3EClusELeadingCellVsTime(0),
  80   fh3JetReacCent(0),
  81   fh2FullJetCent(0)
  82 {
  83   // Default constructor.
  84
  85   SetMakeGeneralHistograms(kTRUE);
  86 }
  87
  88 //________________________________________________________________________
  89 AliAnalysisTaskEmcalJetTriggerQA::AliAnalysisTaskEmcalJetTriggerQA(const char *name) :
  90   AliAnalysisTaskEmcalJet(name, kTRUE),
  91   fDebug(kFALSE),
  92   fTriggerClass(""),
  93   fBitJ1((1<<8)),
  94   fBitJ2((1<<11)),
  95   fContainerFull(0),
  96   fContainerCharged(1),
  97   fMaxPatchEnergy(0),
  98   fTriggerType(-1),
  99   fNFastOR(16),
 100   fhNEvents(0),
 101   fHistRhovsCentFull(0),
 102   fHistRhovsCentCharged(0),
 103   fh3PtEtaPhiTracks(0),
 104   fh3PtEtaPhiTracksOnEmcal(0),
 105   fh3PtEtaPhiTracksProp(0),
 106   fh3PtEtaPhiJetFull(0),
 107   fh3PtEtaPhiJetCharged(0),
 108   fh2NJetsPtFull(0),
 109   fh2NJetsPtCharged(0),
 110   fh3PtEtaAreaJetFull(0),
 111   fh3PtEtaAreaJetCharged(0),
 112   fh2PtNConstituentsCharged(0),
 113   fh2PtNConstituents(0),
 114   fh2PtMeanPtConstituentsCharged(0),
 115   fh2PtMeanPtConstituentsNeutral(0),
 116   fh2PtNEF(0),
 117   fh3NEFEtaPhi(0),
 118   fh2NEFNConstituentsCharged(0),
 119   fh2NEFNConstituentsNeutral(0),
 120   fh2Ptz(0),
 121   fh2PtzCharged(0),
 122   fh2PtLeadJet1VsLeadJet2(0),
 123   fh3EEtaPhiCluster(0),
 124   fh3PtLeadJet1VsPatchEnergy(0),
 125   fh3PtLeadJet2VsPatchEnergy(0),
 126   fh3PatchEnergyEtaPhiCenterJ1(0),
 127   fh3PatchEnergyEtaPhiCenterJ2(0),
 128   fh3PatchEnergyEtaPhiCenterJ1J2(0),
 129   fh3PatchADCEnergyEtaPhiCenterJ1(0),
 130   fh3PatchADCEnergyEtaPhiCenterJ2(0),
 131   fh3PatchADCEnergyEtaPhiCenterJ1J2(0),
 132   fh2CellEnergyVsTime(0),
 133   fh3EClusELeadingCellVsTime(0),
 134   fh3JetReacCent(0),
 135   fh2FullJetCent(0)
 136 {
 137   // Standard constructor.
 138
 139   SetMakeGeneralHistograms(kTRUE);
 140 }
 141
 142 //________________________________________________________________________
 143 AliAnalysisTaskEmcalJetTriggerQA::~AliAnalysisTaskEmcalJetTriggerQA()
 144 {
 145   // Destructor.
 146   if (fOutput) {
 147     delete fOutput;  // delete output object list
 148     fOutput = 0;
 149   }
 150
 151 }
 152
 153 //________________________________________________________________________
 154 Bool_t AliAnalysisTaskEmcalJetTriggerQA::SelectEvent() {
 155   //
 156   // Decide if event should be selected for analysis
 157   //
 158
 159   fhNEvents->Fill(3.5);
 160
 161   if(!fTriggerClass.IsNull()) {
 162     //Check if requested trigger was fired
 163     TString trigType1 = "J1";
 164     TString trigType2 = "J2";
 165     if(fTriggerClass.Contains("G")) {
 166       trigType1 = "G1";
 167       trigType1 = "G2";
 168     }
 169
 170     TString firedTrigClass = InputEvent()->GetFiredTriggerClasses();
 171
 172     if(fTriggerClass.Contains(trigType1.Data()) && fTriggerClass.Contains(trigType2.Data())) { //if events with J1&&J2 are requested
 173       if(!firedTrigClass.Contains(trigType1.Data()) || !firedTrigClass.Contains(trigType2.Data()) ) //check if both are fired
 174         return kFALSE;
 175     }
 176     else {
 177       if(!firedTrigClass.Contains(fTriggerClass))
 178         return kFALSE;
 179       else if(fTriggerClass.Contains(trigType1.Data()) && firedTrigClass.Contains(trigType2.Data())) //if J2 is requested also add triggers which have J1&&J2. Reject if J1 is requested and J2 is fired
 180         return kFALSE;
 181     }
 182   }
 183
 184   fhNEvents->Fill(1.5);
 185
 186   return kTRUE;
 187
 188 }
 189
 190 //________________________________________________________________________
 191 void AliAnalysisTaskEmcalJetTriggerQA::FindTriggerPatch() {
 192
 193   //Fill trigger patch histos for main trigger
 194
 195   AliEmcalTriggerPatchInfo *patch = GetMainTriggerPatch();
 196   if(patch) {
 197     fMaxPatchEnergy = patch->GetPatchE();
 198     Double_t patchADCGeV = patch->GetADCAmpGeVRough();
 199     if(patch->IsJetLow() && !patch->IsJetHigh()) { //main patch only fired low threshold trigger
 200       fh3PatchEnergyEtaPhiCenterJ2->Fill(patch->GetPatchE(),patch->GetEtaGeo(),patch->GetPhiGeo());
 201       fh3PatchADCEnergyEtaPhiCenterJ2->Fill(patchADCGeV,patch->GetEtaGeo(),patch->GetPhiGeo());
 202     }
 203     else if(patch->IsJetHigh() && !patch->IsJetLow()) { //main patch only fired high threshold trigger - should never happen
 204       fh3PatchEnergyEtaPhiCenterJ1->Fill(patch->GetPatchE(),patch->GetEtaGeo(),patch->GetPhiGeo());
 205       fh3PatchADCEnergyEtaPhiCenterJ1->Fill(patchADCGeV,patch->GetEtaGeo(),patch->GetPhiGeo());
 206     }
 207     else if(patch->IsJetHigh() && patch->IsJetLow()) { //main patch fired both triggers
 208       fh3PatchEnergyEtaPhiCenterJ1J2->Fill(patch->GetPatchE(),patch->GetEtaGeo(),patch->GetPhiGeo());
 209       fh3PatchADCEnergyEtaPhiCenterJ1J2->Fill(patchADCGeV,patch->GetEtaGeo(),patch->GetPhiGeo());
 210     }
 211   }
 212 }
 213
 214 //________________________________________________________________________
 215 void AliAnalysisTaskEmcalJetTriggerQA::UserCreateOutputObjects()
 216 {
 217   // Create user output.
 218
 219   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 220
 221   Bool_t oldStatus = TH1::AddDirectoryStatus();
 222   TH1::AddDirectory(kFALSE);
 223
 224   fhNEvents = new TH1F("fhNEvents","fhNEvents;selection;N_{evt}",5,0,5);
 225   fOutput->Add(fhNEvents);
 226
 227   fHistRhovsCentFull = new TH2F("fHistRhovsCentFull", "fHistRhovsCentFull", 101, -1,  100, 300, 0., 300.);
 228   fHistRhovsCentFull->GetXaxis()->SetTitle("Centrality (%)");
 229   fHistRhovsCentFull->GetYaxis()->SetTitle("s#rho_{ch} (GeV/c * rad^{-1})");
 230   fOutput->Add(fHistRhovsCentFull);
 231
 232   fHistRhovsCentCharged = new TH2F("fHistRhovsCentCharged", "fHistRhovsCentCharged", 101, -1,  100, 300, 0., 300.);
 233   fHistRhovsCentCharged->GetXaxis()->SetTitle("Centrality (%)");
 234   fHistRhovsCentCharged->GetYaxis()->SetTitle("#rho_{ch} (GeV/c * rad^{-1})");
 235   fOutput->Add(fHistRhovsCentCharged);
 236
 237   Int_t fgkNCentBins = 20;
 238   Float_t kMinCent   = 0.;
 239   Float_t kMaxCent   = 100.;
 240   Double_t *binsCent = new Double_t[fgkNCentBins+1];
 241   for(Int_t i=0; i<=fgkNCentBins; i++) binsCent[i]=(Double_t)kMinCent + (kMaxCent-kMinCent)/fgkNCentBins*(Double_t)i ;
 242
 243   Int_t fgkNdEPBins = 18*8;
 244   Float_t kMindEP   = 0.;
 245   Float_t kMaxdEP   = 1.*TMath::Pi()/2.;
 246   Double_t *binsdEP = new Double_t[fgkNdEPBins+1];
 247   for(Int_t i=0; i<=fgkNdEPBins; i++) binsdEP[i]=(Double_t)kMindEP + (kMaxdEP-kMindEP)/fgkNdEPBins*(Double_t)i ;
 248
 249   Int_t fgkNPtBins = 200;
 250   Float_t kMinPt   = -50.;
 251   Float_t kMaxPt   = 150.;
 252   Double_t *binsPt = new Double_t[fgkNPtBins+1];
 253   for(Int_t i=0; i<=fgkNPtBins; i++) binsPt[i]=(Double_t)kMinPt + (kMaxPt-kMinPt)/fgkNPtBins*(Double_t)i ;
 254
 255   Int_t fgkNPhiBins = 18*8;
 256   Float_t kMinPhi   = 0.;
 257   Float_t kMaxPhi   = 2.*TMath::Pi();
 258   Double_t *binsPhi = new Double_t[fgkNPhiBins+1];
 259   for(Int_t i=0; i<=fgkNPhiBins; i++) binsPhi[i]=(Double_t)kMinPhi + (kMaxPhi-kMinPhi)/fgkNPhiBins*(Double_t)i ;
 260
 261   Int_t fgkNEtaBins = 100;
 262   Float_t fgkEtaMin = -1.;
 263   Float_t fgkEtaMax =  1.;
 264   Double_t *binsEta=new Double_t[fgkNEtaBins+1];
 265   for(Int_t i=0; i<=fgkNEtaBins; i++) binsEta[i]=(Double_t)fgkEtaMin + (fgkEtaMax-fgkEtaMin)/fgkNEtaBins*(Double_t)i ;
 266
 267   Int_t fgkNAreaBins = 100;
 268   Float_t kMinArea   = 0.;
 269   Float_t kMaxArea   = 1.;
 270   Double_t *binsArea = new Double_t[fgkNAreaBins+1];
 271   for(Int_t i=0; i<=fgkNAreaBins; i++) binsArea[i]=(Double_t)kMinArea + (kMaxArea-kMinArea)/fgkNAreaBins*(Double_t)i ;
 272
 273   Int_t fgkNConstBins = 100;
 274   Float_t kMinConst   = 0.;
 275   Float_t kMaxConst   = 100.;
 276   Double_t *binsConst = new Double_t[fgkNConstBins+1];
 277   for(Int_t i=0; i<=fgkNConstBins; i++) binsConst[i]=(Double_t)kMinConst + (kMaxConst-kMinConst)/fgkNConstBins*(Double_t)i ;
 278
 279   Int_t fgkNMeanPtBins = 100;
 280   Float_t kMinMeanPt   = 0.;
 281   Float_t kMaxMeanPt   = 20.;
 282   Double_t *binsMeanPt = new Double_t[fgkNMeanPtBins+1];
 283   for(Int_t i=0; i<=fgkNMeanPtBins; i++) binsMeanPt[i]=(Double_t)kMinMeanPt + (kMaxMeanPt-kMinMeanPt)/fgkNMeanPtBins*(Double_t)i ;
 284
 285   Int_t fgkNNEFBins = 101;
 286   Float_t kMinNEF   = 0.;
 287   Float_t kMaxNEF   = 1.01;
 288   Double_t *binsNEF = new Double_t[fgkNNEFBins+1];
 289   for(Int_t i=0; i<=fgkNNEFBins; i++) binsNEF[i]=(Double_t)kMinNEF + (kMaxNEF-kMinNEF)/fgkNNEFBins*(Double_t)i ;
 290
 291   Int_t fgkNzBins = 101;
 292   Float_t kMinz   = 0.;
 293   Float_t kMaxz   = 1.01;
 294   Double_t *binsz = new Double_t[fgkNzBins+1];
 295   for(Int_t i=0; i<=fgkNzBins; i++) binsz[i]=(Double_t)kMinz + (kMaxz-kMinz)/fgkNzBins*(Double_t)i ;
 296
 297   Int_t fgkNJetTypeBins = 2;
 298   Float_t kMinJetType   = -0.5;
 299   Float_t kMaxJetType   = 1.5;
 300   Double_t *binsJetType = new Double_t[fgkNJetTypeBins+1];
 301   for(Int_t i=0; i<=fgkNJetTypeBins; i++) binsJetType[i]=(Double_t)kMinJetType + (kMaxJetType-kMinJetType)/fgkNJetTypeBins*(Double_t)i ;
 302
 303   Int_t fgkNTimeBins = 100;
 304   Float_t kMinTime   = -200.;
 305   Float_t kMaxTime   = 200;
 306   Double_t *binsTime = new Double_t[fgkNTimeBins+1];
 307   for(Int_t i=0; i<=fgkNTimeBins; i++) binsTime[i]=(Double_t)kMinTime + (kMaxTime-kMinTime)/fgkNTimeBins*(Double_t)i ;
 308
 309   Double_t enBinEdges[3][2];
 310   enBinEdges[0][0] = 1.; //10 bins
 311   enBinEdges[0][1] = 0.1;
 312   enBinEdges[1][0] = 5.; //8 bins
 313   enBinEdges[1][1] = 0.5;
 314   enBinEdges[2][0] = 100.;//95 bins
 315   enBinEdges[2][1] = 1.;
 316
 317   const Float_t enmin1 =  0;
 318   const Float_t enmax1 =  enBinEdges[0][0];
 319   const Float_t enmin2 =  enmax1 ;
 320   const Float_t enmax2 =  enBinEdges[1][0];
 321   const Float_t enmin3 =  enmax2 ;
 322   const Float_t enmax3 =  enBinEdges[2][0];//fgkEnMax;
 323   const Int_t nbin11 = (int)((enmax1-enmin1)/enBinEdges[0][1]);
 324   const Int_t nbin12 = (int)((enmax2-enmin2)/enBinEdges[1][1])+nbin11;
 325   const Int_t nbin13 = (int)((enmax3-enmin3)/enBinEdges[2][1])+nbin12;
 326
 327   Int_t fgkNEnBins=nbin13;
 328   Double_t *binsEn=new Double_t[fgkNEnBins+1];
 329   for(Int_t i=0; i<=fgkNEnBins; i++) {
 330     if(i<=nbin11) binsEn[i]=(Double_t)enmin1 + (enmax1-enmin1)/nbin11*(Double_t)i ;
 331     if(i<=nbin12 && i>nbin11) binsEn[i]=(Double_t)enmin2 + (enmax2-enmin2)/(nbin12-nbin11)*((Double_t)i-(Double_t)nbin11) ;
 332     if(i<=nbin13 && i>nbin12) binsEn[i]=(Double_t)enmin3 + (enmax3-enmin3)/(nbin13-nbin12)*((Double_t)i-(Double_t)nbin12) ;
 333   }
 334
 335   fh3PtEtaPhiTracks = new TH3F("fh3PtEtaPhiTracks","fh3PtEtaPhiTracks;#it{p}_{T}^{track};#eta;#varphi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
 336   fOutput->Add(fh3PtEtaPhiTracks);
 337
 338   fh3PtEtaPhiTracksOnEmcal = new TH3F("fh3PtEtaPhiTracksOnEmcal","fh3PtEtaPhiTracksOnEmcal;#it{p}_{T}^{track};#eta;#varphi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
 339   fOutput->Add(fh3PtEtaPhiTracksOnEmcal);
 340
 341   fh3PtEtaPhiTracksProp = new TH3F("fh3PtEtaPhiTracksProp","fh3PtEtaPhiTracksProp;#it{p}_{T}^{track};#eta;#varphi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
 342   fOutput->Add(fh3PtEtaPhiTracksProp);
 343
 344   fh3PtEtaPhiJetFull = new TH3F("fh3PtEtaPhiJetFull","fh3PtEtaPhiJetFull;#it{p}_{T}^{jet};#eta;#varphi",fgkNPtBins,binsPt,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
 345   fOutput->Add(fh3PtEtaPhiJetFull);
 346
 347   fh3PtEtaPhiJetCharged = new TH3F("fh3PtEtaPhiJetCharged","fh3PtEtaPhiJetCharged;#it{p}_{T}^{jet};#eta;#varphi",fgkNPtBins,binsPt,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
 348   fOutput->Add(fh3PtEtaPhiJetCharged);
 349
 350   fh2NJetsPtFull = new TH2F("fh2NJetsPtFull","fh2NJetsPtFull;N_{jets};#it{p}_{T}^{jet}",20,-0.5,19.5,fgkNPtBins,binsPt);
 351   fOutput->Add(fh2NJetsPtFull);
 352
 353   fh2NJetsPtCharged = new TH2F("fh2NJetsPtCharged","fh2NJetsPtCharged;N_{jets};#it{p}_{T}^{jet}",20,-0.5,19.5,fgkNPtBins,binsPt);
 354   fOutput->Add(fh2NJetsPtCharged);
 355
 356   fh3PtEtaAreaJetFull = new TH3F("fh3PtEtaAreaJetFull","fh3PtEtaAreaJetFull;#it{p}_{T}^{jet};#eta;A",fgkNPtBins,binsPt,fgkNEtaBins,binsEta,fgkNAreaBins,binsArea);
 357   fOutput->Add(fh3PtEtaAreaJetFull);
 358
 359   fh3PtEtaAreaJetCharged = new TH3F("fh3PtEtaAreaJetCharged","fh3PtEtaAreaJetCharged;#it{p}_{T}^{jet};#eta;A",fgkNPtBins,binsPt,fgkNEtaBins,binsEta,fgkNAreaBins,binsArea);
 360   fOutput->Add(fh3PtEtaAreaJetCharged);
 361
 362   fh2PtNConstituentsCharged = new TH2F("fh2PtNConstituentsCharged","fh2PtNConstituentsCharged;#it{p}_{T}^{jet};N_{charged constituents}",fgkNPtBins,binsPt,fgkNConstBins,binsConst);
 363   fOutput->Add(fh2PtNConstituentsCharged);
 364
 365   fh2PtNConstituents = new TH2F("fh2PtNConstituents","fh2PtNConstituents;#it{p}_{T}^{jet};N_{constituents}",fgkNPtBins,binsPt,fgkNConstBins,binsConst);
 366   fOutput->Add(fh2PtNConstituents);
 367
 368   fh2PtMeanPtConstituentsCharged = new TH2F("fh2PtMeanPtConstituentsCharged","fh2PtMeanPtConstituentsCharged;#it{p}_{T}^{jet};charged #langle #it{p}_{T} #rangle",fgkNPtBins,binsPt,fgkNMeanPtBins,binsMeanPt);
 369   fOutput->Add(fh2PtMeanPtConstituentsCharged);
 370
 371   fh2PtMeanPtConstituentsNeutral = new TH2F("fh2PtMeanPtConstituentsNeutral","fh2PtMeanPtConstituentsNeutral;#it{p}_{T}^{jet};neutral langle #it{p}_{T} #rangle",fgkNPtBins,binsPt,fgkNMeanPtBins,binsMeanPt);
 372   fOutput->Add(fh2PtMeanPtConstituentsNeutral);
 373
 374   fh2PtNEF = new TH2F("fh2PtNEF","fh2PtNEF;#it{p}_{T}^{jet};NEF",fgkNPtBins,binsPt,fgkNNEFBins,binsNEF);
 375   fOutput->Add(fh2PtNEF);
 376
 377   fh3NEFEtaPhi = new TH3F("fh3NEFEtaPhi","fh3NEFEtaPhi;NEF;#eta;#varphi",fgkNNEFBins,binsNEF,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
 378   fOutput->Add(fh3NEFEtaPhi);
 379
 380   fh2NEFNConstituentsCharged = new TH2F("fh2NEFNConstituentsCharged","fh2NEFNConstituentsCharged;NEF;N_{charged constituents}",fgkNNEFBins,binsNEF,fgkNConstBins,binsConst);
 381   fOutput->Add(fh2NEFNConstituentsCharged);
 382
 383   fh2NEFNConstituentsNeutral = new TH2F("fh2NEFNConstituentsNeutral","fh2NEFNConstituentsNeutral;NEF;N_{clusters}",fgkNNEFBins,binsNEF,fgkNConstBins,binsConst);
 384   fOutput->Add(fh2NEFNConstituentsNeutral);
 385
 386   fh2Ptz = new TH2F("fh2Ptz","fh2Ptz;#it{p}_{T}^{jet};z=p_{t,trk}^{proj}/p_{jet}",fgkNPtBins,binsPt,fgkNzBins,binsz);
 387   fOutput->Add(fh2Ptz);
 388
 389   fh2PtzCharged = new TH2F("fh2PtzCharged","fh2Ptz;#it{p}_{T}^{ch jet};z=p_{t,trk}^{proj}/p_{ch jet}",fgkNPtBins,binsPt,fgkNzBins,binsz);
 390   fOutput->Add(fh2PtzCharged);
 391
 392   fh2PtLeadJet1VsLeadJet2 = new TH2F("fh2PtLeadJet1VsLeadJet2","fh2PtLeadJet1VsLeadJet2;#it{p}_{T}^{jet 1};#it{p}_{T}^{jet 2}",fgkNPtBins,binsPt,fgkNPtBins,binsPt);
 393   fOutput->Add(fh2PtLeadJet1VsLeadJet2);
 394
 395   fh3EEtaPhiCluster = new TH3F("fh3EEtaPhiCluster","fh3EEtaPhiCluster;E_{clus};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
 396   fOutput->Add(fh3EEtaPhiCluster);
 397
 398   fh3PtLeadJet1VsPatchEnergy = new TH3F("fh3PtLeadJet1VsPatchEnergy","fh3PtLeadJet1VsPatchEnergy;#it{p}_{T}^{jet 1};Amplitude_{patch};trig type",fgkNPtBins,binsPt,fgkNPtBins,binsPt,fgkNJetTypeBins,binsJetType);
 399   fOutput->Add(fh3PtLeadJet1VsPatchEnergy);
 400   fh3PtLeadJet2VsPatchEnergy = new TH3F("fh3PtLeadJet2VsPatchEnergy","fh3PtLeadJet2VsPatchEnergy;#it{p}_{T}^{jet 1};Amplitude_{patch};trig type",fgkNPtBins,binsPt,fgkNPtBins,binsPt,fgkNJetTypeBins,binsJetType);
 401   fOutput->Add(fh3PtLeadJet2VsPatchEnergy);
 402
 403   fh3PatchEnergyEtaPhiCenterJ1 = new TH3F("fh3PatchEnergyEtaPhiCenterJ1","fh3PatchEnergyEtaPhiCenterJ1;E_{patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
 404   fOutput->Add(fh3PatchEnergyEtaPhiCenterJ1);
 405
 406   fh3PatchEnergyEtaPhiCenterJ2 = new TH3F("fh3PatchEnergyEtaPhiCenterJ2","fh3PatchEnergyEtaPhiCenterJ2;E_{patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
 407   fOutput->Add(fh3PatchEnergyEtaPhiCenterJ2);
 408
 409   fh3PatchEnergyEtaPhiCenterJ1J2 = new TH3F("fh3PatchEnergyEtaPhiCenterJ1J2","fh3PatchEnergyEtaPhiCenterJ1J2;E_{patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
 410   fOutput->Add(fh3PatchEnergyEtaPhiCenterJ1J2);
 411
 412   fh3PatchADCEnergyEtaPhiCenterJ1 = new TH3F("fh3PatchADCEnergyEtaPhiCenterJ1","fh3PatchADCEnergyEtaPhiCenterJ1;E_{ADC,patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
 413   fOutput->Add(fh3PatchADCEnergyEtaPhiCenterJ1);
 414
 415   fh3PatchADCEnergyEtaPhiCenterJ2 = new TH3F("fh3PatchADCEnergyEtaPhiCenterJ2","fh3PatchADCEnergyEtaPhiCenterJ2;E_{ADC,patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
 416   fOutput->Add(fh3PatchADCEnergyEtaPhiCenterJ2);
 417
 418   fh3PatchADCEnergyEtaPhiCenterJ1J2 = new TH3F("fh3PatchADCEnergyEtaPhiCenterJ1J2","fh3PatchADCEnergyEtaPhiCenterJ1J2;E_{ADC,patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
 419   fOutput->Add(fh3PatchADCEnergyEtaPhiCenterJ1J2);
 420
 421   fh2CellEnergyVsTime = new TH2F("fh2CellEnergyVsTime","fh2CellEnergyVsTime;E_{cell};time",fgkNEnBins,binsEn,fgkNTimeBins,binsTime);
 422   fOutput->Add(fh2CellEnergyVsTime);
 423
 424   fh3EClusELeadingCellVsTime = new TH3F("fh3EClusELeadingCellVsTime","fh3EClusELeadingCellVsTime;E_{cluster};E_{leading cell};time_{leading cell}",fgkNEnBins,binsEn,fgkNEnBins,binsEn,fgkNTimeBins,binsTime);
 425   fOutput->Add(fh3EClusELeadingCellVsTime);
 426
 427   fh3JetReacCent = new TH3F("fh3JetReacCent","fh3JetReacCent;E_{Jet};Centrality;dEP",fgkNEnBins,binsEn,fgkNCentBins,binsCent,fgkNdEPBins,binsdEP);
 428   fOutput->Add(fh3JetReacCent);
 429
 430   fh2FullJetCent = new TH2F("fh2FullJetCent","fh2FullJetCent;Centrality;dEP",fgkNCentBins,binsCent,fgkNdEPBins,binsdEP);
 431   fOutput->Add(fh2FullJetCent);
 432
 433
 434   // =========== Switch on Sumw2 for all histos ===========
 435   for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
 436     TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
 437     if (h1){
 438       h1->Sumw2();
 439       continue;
 440     }
 441     TH2 *h2 = dynamic_cast<TH2*>(fOutput->At(i));
 442     if (h2){
 443       h2->Sumw2();
 444       continue;
 445     }
 446     TH3 *h3 = dynamic_cast<TH3*>(fOutput->At(i));
 447     if (h3){
 448       h3->Sumw2();
 449       continue;
 450     }
 451     THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
 452     if(hn)hn->Sumw2();
 453   }
 454
 455   TH1::AddDirectory(oldStatus);
 456
 457   PostData(1, fOutput); // Post data for ALL output slots > 0 here.
 458
 459   if(binsCent)              delete [] binsCent;
 460   if(binsdEP)               delete [] binsdEP;
 461   if(binsEn)                delete [] binsEn;
 462   if(binsPt)                delete [] binsPt;
 463   if(binsPhi)               delete [] binsPhi;
 464   if(binsEta)               delete [] binsEta;
 465   if(binsArea)              delete [] binsArea;
 466   if(binsConst)             delete [] binsConst;
 467   if(binsMeanPt)            delete [] binsMeanPt;
 468   if(binsNEF)               delete [] binsNEF;
 469   if(binsz)                 delete [] binsz;
 470   if(binsJetType)           delete [] binsJetType;
 471   if(binsTime)              delete [] binsTime;
 472
 473 }
 474
 475 //________________________________________________________________________
 476 Bool_t AliAnalysisTaskEmcalJetTriggerQA::FillHistograms()
 477 {
 478   // Fill histograms.
 479
 480   //Tracks
 481   AliParticleContainer *partCont = GetParticleContainer(0);
 482   if (partCont) {
 483     Int_t i = 0;
 484     AliPicoTrack *track = dynamic_cast<AliPicoTrack*>(partCont->GetNextAcceptParticle(0));
 485     while(track) {
 486       Double_t trkphi = track->Phi()*TMath::RadToDeg();
 487       fh3PtEtaPhiTracks->Fill(track->Pt(),track->Eta(),track->Phi());
 488       fh3PtEtaPhiTracksOnEmcal->Fill(track->GetTrackPtOnEMCal(),track->GetTrackEtaOnEMCal(),track->GetTrackPhiOnEMCal());
 489       if(track->IsEMCAL()) {
 490         i++;
 491         if(TMath::Abs(track->Eta())<0.9 && trkphi > 10 && trkphi < 250 )
 492           fh3PtEtaPhiTracksProp->Fill(track->Pt(),track->Eta(),track->Phi());
 493       }
 494       track = dynamic_cast<AliPicoTrack*>(partCont->GetNextAcceptParticle());
 495     }
 496   }
 497
 498   //Clusters
 499   AliClusterContainer  *clusCont = GetClusterContainer(0);
 500   if (clusCont) {
 501     Int_t nclusters = clusCont->GetNClusters();
 502     TString arrName = clusCont->GetArrayName();
 503     for (Int_t ic = 0; ic < nclusters; ic++) {
 504       AliVCluster *cluster = static_cast<AliVCluster*>(clusCont->GetCluster(ic));
 505       if (!cluster) {
 506         AliDebug(2,Form("Could not receive cluster %d", ic));
 507         continue;
 508       }
 509       if (!cluster->IsEMCAL()) {
 510         AliDebug(2,Form("%s: Cluster is not emcal",GetName()));
 511         continue;
 512       }
 513
 514       TLorentzVector lp;
 515       cluster->GetMomentum(lp, const_cast<Double_t*>(fVertex));
 516       fh3EEtaPhiCluster->Fill(lp.E(),lp.Eta(),lp.Phi());
 517       if(fCaloCells) {
 518         Double_t leadCellE = GetEnergyLeadingCell(cluster);
 519         Double_t leadCellT = cluster->GetTOF();
 520         fh3EClusELeadingCellVsTime->Fill(lp.E(),leadCellE,leadCellT*1e9);
 521       }
 522     }
 523   }
 524
 525   //cells
 526   if(fCaloCells) {
 527     const Short_t nCells   = fCaloCells->GetNumberOfCells();
 528
 529     for(Int_t iCell=0; iCell<nCells; ++iCell) {
 530       Short_t cellId = fCaloCells->GetCellNumber(iCell);
 531       Double_t cellE = fCaloCells->GetCellAmplitude(cellId);
 532       Double_t cellT = fCaloCells->GetCellTime(cellId);
 533
 534       AliDebug(2,Form("cell energy = %f  time = %f",cellE,cellT*1e9));
 535       fh2CellEnergyVsTime->Fill(cellE,cellT*1e9);
 536     }
 537   }
 538
 539   //Jets
 540   Double_t ptLeadJet1 = 0.;
 541   Double_t ptLeadJet2 = 0.;
 542
 543   fHistRhovsCentFull->Fill(fCent,GetRhoVal(fContainerFull));
 544   fHistRhovsCentCharged->Fill(fCent,GetRhoVal(fContainerCharged));
 545
 546   TArrayI *nJetsArr = new TArrayI(fh2NJetsPtFull->GetNbinsY()+1);
 547   nJetsArr->Reset(0);
 548   nJetsArr->Set(fh2NJetsPtFull->GetNbinsY()+1);
 549
 550   if (GetJetContainer(fContainerFull)) {
 551     const Int_t njets = GetNJets(fContainerFull);
 552     for (Int_t ij = 0; ij < njets; ij++) {
 553
 554       AliEmcalJet* jet = GetAcceptJetFromArray(ij,fContainerFull);
 555       if (!jet)
 556         continue; //jet not selected
 557
 558       Double_t jetPt = jet->Pt() - GetRhoVal(fContainerFull)*jet->Area();
 559       Double_t jetPhi = jet->Phi();
 560       Double_t dEPJetFull = -500.0;
 561       dEPJetFull = RelativeEP(jetPhi , fEPV0);
 562
 563       fh3JetReacCent->Fill(jet->E(),fCent,dEPJetFull);
 564       fh2FullJetCent->Fill(fCent,dEPJetFull);
 565
 566       if(jetPt>ptLeadJet1) ptLeadJet1=jetPt;
 567       fh3PtEtaPhiJetFull->Fill(jetPt,jet->Eta(),jet->Phi());
 568       fh3PtEtaAreaJetFull->Fill(jetPt,jet->Eta(),jet->Area());
 569
 570       //count jets above certain pT threshold
 571       Int_t ptbin = fh2NJetsPtFull->GetYaxis()->FindBin(jetPt);
 572       for(Int_t iptbin = ptbin; iptbin<=fh2NJetsPtFull->GetNbinsY(); iptbin++)
 573         nJetsArr->AddAt(nJetsArr->At(iptbin)+1,iptbin);
 574
 575       fh2PtNConstituentsCharged->Fill(jetPt,jet->GetNumberOfTracks());
 576       fh2PtNConstituents->Fill(jetPt,jet->GetNumberOfConstituents());
 577
 578       fh2PtNEF->Fill(jetPt,jet->NEF());
 579       fh3NEFEtaPhi->Fill(jet->NEF(),jet->Eta(),jet->Phi());
 580       fh2NEFNConstituentsCharged->Fill(jet->NEF(),jet->GetNumberOfTracks());
 581       fh2NEFNConstituentsNeutral->Fill(jet->NEF(),jet->GetNumberOfClusters());
 582
 583       AliVParticle *vp;
 584       Double_t sumPtCh = 0.;
 585       for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
 586         vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
 587         if(!vp) continue;
 588         fh2Ptz->Fill(jetPt,GetZ(vp,jet));
 589         sumPtCh+=vp->Pt();
 590       }
 591
 592       if(jet->GetNumberOfTracks()>0)
 593         fh2PtMeanPtConstituentsCharged->Fill(jetPt,sumPtCh/(double)(jet->GetNumberOfTracks()) );
 594
 595       AliVCluster *vc = 0x0;
 596       Double_t sumPtNe = 0.;
 597       if (clusCont) {
 598         for(Int_t icc=0; icc<jet->GetNumberOfClusters(); icc++) {
 599           vc = static_cast<AliVCluster*>(clusCont->GetCluster(icc));
 600           if(!vc) continue;
 601           TLorentzVector lp;
 602           vc->GetMomentum(lp, const_cast<Double_t*>(fVertex));
 603           sumPtNe+=lp.Pt();
 604         }
 605
 606         if(jet->GetNumberOfClusters()>0)
 607           fh2PtMeanPtConstituentsNeutral->Fill(jetPt,sumPtNe/(double)(jet->GetNumberOfClusters()) );
 608       }
 609     } //full jet loop
 610
 611     for(Int_t i=1; i<=fh2NJetsPtFull->GetNbinsY(); i++) {
 612       Int_t nJetsInEvent = nJetsArr->At(i);
 613       fh2NJetsPtFull->Fill(nJetsInEvent,fh2NJetsPtFull->GetYaxis()->GetBinCenter(i));
 614     }
 615   }
 616
 617   //Reset array to zero to also count charged jets
 618   nJetsArr->Reset(0);
 619
 620   //Loop over charged jets
 621   if (GetJetContainer(fContainerCharged)) {
 622     const Int_t njets = GetNJets(fContainerCharged);
 623     for (Int_t ij = 0; ij < njets; ij++) {
 624
 625       AliEmcalJet* jet = GetAcceptJetFromArray(ij,fContainerCharged);
 626       if (!jet)
 627         continue; //jet not selected
 628
 629       Double_t jetPt = jet->Pt() - GetRhoVal(fContainerCharged)*jet->Area();
 630       if(jetPt>ptLeadJet2) ptLeadJet2=jetPt;
 631       fh3PtEtaPhiJetCharged->Fill(jetPt,jet->Eta(),jet->Phi());
 632       fh3PtEtaAreaJetCharged->Fill(jetPt,jet->Eta(),jet->Area());
 633
 634       AliVParticle *vp;
 635       for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
 636         vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
 637         if(!vp) continue;
 638         fh2PtzCharged->Fill(jetPt,GetZ(vp,jet));
 639       }
 640
 641       //count jets above certain pT threshold
 642       Int_t ptbin = fh2NJetsPtCharged->GetYaxis()->FindBin(jetPt);
 643       for(Int_t iptbin = ptbin; iptbin<=fh2NJetsPtCharged->GetNbinsY(); iptbin++)
 644         nJetsArr->AddAt(nJetsArr->At(iptbin)+1,iptbin);
 645
 646     }//ch jet loop
 647     for(Int_t i=1; i<=fh2NJetsPtCharged->GetNbinsY(); i++) {
 648       Int_t nJetsInEvent = nJetsArr->At(i);
 649       fh2NJetsPtCharged->Fill(nJetsInEvent,fh2NJetsPtCharged->GetYaxis()->GetBinCenter(i));
 650     }
 651   }
 652
 653   if(GetJetContainer(fContainerFull) && GetJetContainer(fContainerCharged)) {
 654     fh2PtLeadJet1VsLeadJet2->Fill(ptLeadJet1,ptLeadJet2);
 655   }
 656
 657   fh3PtLeadJet1VsPatchEnergy->Fill(ptLeadJet1,fMaxPatchEnergy,fTriggerType);
 658   fh3PtLeadJet2VsPatchEnergy->Fill(ptLeadJet2,fMaxPatchEnergy,fTriggerType);
 659
 660   if(nJetsArr) delete nJetsArr;
 661
 662   return kTRUE;
 663 }
 664
 665 //________________________________________________________________________
 666 Bool_t AliAnalysisTaskEmcalJetTriggerQA::Run()
 667 {
 668   // Run analysis code here, if needed. It will be executed before FillHistograms().
 669
 670   //Check if event is selected (vertex & pile-up)
 671   if(!SelectEvent())
 672     return kFALSE;
 673
 674   if(fTriggerPatchInfo)
 675     FindTriggerPatch();
 676
 677   return kTRUE;  // If return kFALSE FillHistogram() will NOT be executed.
 678 }
 679
 680 //_______________________________________________________________________
 681 void AliAnalysisTaskEmcalJetTriggerQA::Terminate(Option_t *)
 682 {
 683   // Called once at the end of the analysis.
 684 }
 685 //________________________________________________________________________
 686 Double_t AliAnalysisTaskEmcalJetTriggerQA::GetZ(const AliVParticle *trk, const AliEmcalJet *jet)          const
 687 {
 688   // Get Z of constituent trk
 689   return GetZ(trk->Px(),trk->Py(),trk->Pz(),jet->Px(),jet->Py(),jet->Pz());
 690 }
 691
 692 //________________________________________________________________________
 693 Double_t AliAnalysisTaskEmcalJetTriggerQA::GetZ(const Double_t trkPx, const Double_t trkPy, const Double_t trkPz, const Double_t jetPx, const Double_t jetPy, const Double_t jetPz) const
 694 {
 695   //
 696   // Get the z of a constituent inside of a jet
 697   //
 698   Double_t pJetSq = jetPx*jetPx+jetPy*jetPy+jetPz*jetPz;
 699   if(pJetSq>0.)
 700     return (trkPx*jetPx+trkPy*jetPy+trkPz*jetPz)/pJetSq;
 701   else {
 702     AliWarning(Form("%s: strange, pjet*pjet seems to be zero pJetSq: %f",GetName(), pJetSq));
 703     return 0;
 704   }
 705 }
 706
 707 //________________________________________________________________________
 708 Int_t AliAnalysisTaskEmcalJetTriggerQA::GetLeadingCellId(const AliVCluster *clus) const
 709 {
 710   //Get energy of leading cell in cluster
 711
 712   if(!fCaloCells)
 713     return -1;
 714
 715   Double_t emax = -1.;
 716   Int_t iCellAbsIdMax = -1;
 717   Int_t nCells = clus->GetNCells();
 718   for(Int_t i = 0; i<nCells; i++) {
 719     Int_t absId = clus->GetCellAbsId(i);
 720     Double_t cellE = fCaloCells->GetCellAmplitude(absId);
 721     if(cellE>emax) {
 722       emax = cellE;
 723       iCellAbsIdMax = absId;
 724     }
 725   }
 726
 727   return iCellAbsIdMax;
 728 }
 729
 730 //________________________________________________________________________
 731 Double_t AliAnalysisTaskEmcalJetTriggerQA::GetEnergyLeadingCell(const AliVCluster *clus) const
 732 {
 733   //Get energy of leading cell in cluster
 734   if(!fCaloCells)
 735     return -1.;
 736
 737   Int_t absID = GetLeadingCellId(clus);
 738   if(absID>-1)
 739     return fCaloCells->GetCellAmplitude(absID);
 740   else
 741     return -1.;
 742 }
 743
 744 //________________________________________________________________________
 745 Double_t AliAnalysisTaskEmcalJetTriggerQA::GetECross(Int_t absID) const {
 746
 747   //Get Ecross = sum of energy of neighbouring cells (using uncalibrated energy)
 748
 749   if(!fCaloCells)
 750     return -1.;
 751
 752   Double_t ecross = -1.;
 753
 754   Int_t absID1 = -1;
 755   Int_t absID2 = -1;
 756   Int_t absID3 = -1;
 757   Int_t absID4 = -1;
 758
 759   Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
 760   fGeom->GetCellIndex(absID,imod,iTower,iIphi,iIeta);
 761   fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
 762
 763   if( iphi < AliEMCALGeoParams::fgkEMCALRows-1)
 764     absID1 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta);
 765   if( iphi > 0 )
 766     absID2 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta);
 767
 768   if( ieta == AliEMCALGeoParams::fgkEMCALCols-1 && !(imod%2) ) {
 769     absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod+1, iphi, 0);
 770     absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod,   iphi, ieta-1);
 771   }
 772   else if( ieta == 0 && imod%2 ) {
 773     absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod,   iphi, ieta+1);
 774     absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod-1, iphi, AliEMCALGeoParams::fgkEMCALCols-1);
 775   }
 776   else  {
 777     if( ieta < AliEMCALGeoParams::fgkEMCALCols-1 )
 778       absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
 779     if( ieta > 0 )
 780       absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
 781   }
 782
 783   Double_t ecell1 = fCaloCells->GetCellAmplitude(absID1);
 784   Double_t ecell2 = fCaloCells->GetCellAmplitude(absID2);
 785   Double_t ecell3 = fCaloCells->GetCellAmplitude(absID3);
 786   Double_t ecell4 = fCaloCells->GetCellAmplitude(absID4);
 787
 788   ecross = ecell1+ecell2+ecell3+ecell4;
 789
 790   return ecross;
 791 }
 792
 793 //_________________________________________________________________________
 794 Float_t AliAnalysisTaskEmcalJetTriggerQA:: RelativeEP(Double_t objAng, Double_t EPAng) const
 795 { // function to calculate angle between object and EP in the 1st quadrant (0,Pi/2)
 796     Double_t dphi = (EPAng - objAng);
 797
 798     // ran into trouble with a few dEP<-Pi so trying this...
 799     if( dphi<-1*TMath::Pi() ){
 800         dphi = dphi + 1*TMath::Pi();
 801     }
 802
 803     if( dphi>1*TMath::Pi()){
 804         dphi = dphi - 1*TMath::Pi();
 805     }
 806
 807
 808
 809     if( (dphi>0) && (dphi<1*TMath::Pi()/2) ){
 810         // Do nothing! we are in quadrant 1
 811     }else if( (dphi>1*TMath::Pi()/2) && (dphi<1*TMath::Pi()) ){
 812         dphi = 1*TMath::Pi() - dphi;
 813     }else if( (dphi<0) && (dphi>-1*TMath::Pi()/2) ){
 814         dphi = fabs(dphi);
 815     }else if( (dphi<-1*TMath::Pi()/2) && (dphi>-1*TMath::Pi()) ){
 816         dphi = dphi + 1*TMath::Pi();
 817     }
 818
 819     // test
 820     if( dphi < 0 || dphi > TMath::Pi()/2 )
 821         cout<<"dphi:  "<<dphi<<endl;
 822
 823     return dphi;   // dphi in [0, Pi/2]
 824 }
 825