PWGCF/FLOW/CME/AliAnalysisTaskCMEv2A.cxx

   1 #include <iostream>
   2 #include <cstdlib>
   3 #include <sys/time.h>
   4
   5 // ROOT classes
   6 #include "TChain.h"
   7 #include "TTree.h"
   8 #include "TList.h"
   9 #include "TH1F.h"
  10 #include "TH2F.h"
  11 #include "TProfile.h"
  12 #include "TExMap.h"
  13 #include "TRandom3.h"
  14 #include "TMath.h"
  15
  16 // Alice analysis base class
  17 #include "AliAnalysisTaskSE.h"
  18
  19 // Alice analysis additional classes
  20 #include "AliAnalysisManager.h"
  21 //#include "AliInputEventHandler.h"
  22
  23 // Alice AOD classes
  24 #include "AliAODInputHandler.h"
  25 //#include "AliAODHandler.h"
  26 #include "AliAODEvent.h"
  27 #include "AliAODVertex.h"
  28 #include "AliAODVZERO.h"
  29
  30 // Alice classes
  31 #include "AliCentrality.h"
  32 #include "AliEventplane.h"
  33 #include "AliAnalysisUtils.h"
  34
  35 // Alice MC classes
  36 #include "AliMCEvent.h"
  37 //#include "AliMCEventHandler.h"
  38 #include "AliAODMCParticle.h"
  39
  40 // Alice "V" classes
  41 #include "AliVParticle.h"
  42 //#include "AliVEvent.h"
  43 //#include "AliVVertex.h"
  44 //#include "AliVVZERO.h"
  45
  46 // Alice PID classes
  47 //#include "AliAODPid.h"
  48 //#include "AliAODpidUtil.h"
  49 #include "AliPID.h"
  50 //#include "AliPIDCombined.h"
  51 #include "AliPIDResponse.h"
  52
  53
  54 // This class
  55 #include "AliAnalysisTaskCMEv2A.h"
  56
  57 ClassImp(AliAnalysisTaskCMEv2A); // import class inheriting from TObject
  58
  59 using std::cout;
  60 using std::endl;
  61
  62 const float pi = 3.141592653589793;
  63
  64 // function prototype...
  65 float GetDPhiStar(float phi1, float pt1, float charge1, float phi2, float pt2, float charge2, float radius, float bSign);
  66
  67
  68 //------------------------------------------------------------------
  69 AliAnalysisTaskCMEv2A::AliAnalysisTaskCMEv2A() : AliAnalysisTaskSE()
  70 {
  71   // Default class constructor
  72
  73   cout<<"Default class constructor called.  Prepare for fun analysis time!"<<endl;
  74   this->SetParameters();
  75
  76 }
  77
  78
  79
  80 //--------------------------------------------------------------------------------------
  81 AliAnalysisTaskCMEv2A::AliAnalysisTaskCMEv2A(const char *name) : AliAnalysisTaskSE(name)
  82 {
  83   // Class Constructor with name
  84
  85   cout<<"User defined class constructor called.  Prepare for fun analysis time!"<<endl;
  86   this->SetParameters();
  87
  88   // Define input and output slots here
  89   // Input slot #0 works with a TChain
  90   DefineInput(0,TChain::Class());
  91
  92   // Output slot #0 is reserved by the base class for AOD
  93   // Output slot #1 writes into a TH1 container
  94   //DefineOutput(0,TTree::Class());
  95   DefineOutput(1,TList::Class());
  96
  97   // user can change as needed
  98   debug = 0;
  99   doMC = false;
 100   trigger = AliVEvent::kMB;
 101   dopupcut = true;
 102   doeffcorr = true;
 103   centhandle = 1;
 104   fbit = 128;
 105   zvtxcut = 10.0;
 106   centcut = 5.0;
 107   nclscut = 60;
 108   dcacutz = 5.0;
 109   dcacutxy = 5.0;
 110   dodcacuts = false;
 111   outeta = 0.8;
 112   ineta = 0.5;
 113   excleta = 0.0;
 114   ptmin = 0.2;
 115   ptmax = 5.0;
 116   doacuts = false;
 117   nspid = 2.0;
 118   cbinlo = 3;
 119   cbinhi = 4;
 120   donested = true;
 121   dopaircut = true;
 122   centlo = 20.0;
 123   centhi = 60.0;
 124
 125 }
 126
 127
 128
 129 //---------------------------------------------------
 130 AliAnalysisTaskCMEv2A::~AliAnalysisTaskCMEv2A()
 131 {
 132   // Default class destructor
 133
 134   cout<<"Default class destructor called.  Analysis fun time has ended"<<endl;
 135
 136 }
 137
 138
 139
 140 //---------------------------------------------------
 141 void AliAnalysisTaskCMEv2A::UserCreateOutputObjects()
 142 {
 143
 144   // Create output objects, like histograms - called once
 145
 146   cout<<"UserCreateOutputObjects called, now making histrograms and things"<<endl;
 147
 148
 149
 150   // -------------------------- //
 151   // --- create output list --- //
 152   // -------------------------- //
 153
 154   fOutputList = new TList();
 155   fOutputList->SetName(GetName());
 156   fOutputList->SetOwner(kTRUE);
 157
 158
 159
 160   // ------------------------- //
 161   // --- create histograms --- //
 162   // ------------------------- //
 163
 164   fHistPt = new TH1F("fHistPt","",50,0.0,5.0);
 165   fHistPhi = new TH1F("fHistPhi","",63,0.0,6.3);
 166   fHistEta = new TH1F("fHistEta","",300,-1.5,1.5);
 167   fHistCharge = new TH1F("fHistCharge","",3,-1.5,1.5);
 168   fHistTPCncls = new TH1F("fHistTPCncls","",161,-0.5,160.5);
 169   fHistDedx = new TH1F("fHistDedx","",1000,-20,980);
 170   fHistDCAxy = new TH1F("fHistDCAxy","",350,-3.5,3.5);
 171   fHistDCAz = new TH1F("fHistDCAz","",350,-3.5,3.5);
 172   fHistDCAxyAfter = new TH1F("fHistDCAxyAfter","",350,-3.5,3.5);
 173   fHistDCAzAfter = new TH1F("fHistDCAzAfter","",350,-3.5,3.5);
 174   fOutputList->Add(fHistPt);
 175   fOutputList->Add(fHistPhi);
 176   fOutputList->Add(fHistEta);
 177   fOutputList->Add(fHistCharge);
 178   fOutputList->Add(fHistTPCncls);
 179   fOutputList->Add(fHistDedx);
 180   fOutputList->Add(fHistDCAxy);
 181   fOutputList->Add(fHistDCAz);
 182   fOutputList->Add(fHistDCAxyAfter);
 183   fOutputList->Add(fHistDCAzAfter);
 184
 185   fHistPosPt = new TH1F("fHistPosPt","",50,0.0,5.0);
 186   fHistPosPhi = new TH1F("fHistPosPhi","",63,0.0,6.3);
 187   fHistPosEta = new TH1F("fHistPosEta","",300,-1.5,1.5);
 188   fOutputList->Add(fHistPosPt);
 189   fOutputList->Add(fHistPosPhi);
 190   fOutputList->Add(fHistPosEta);
 191
 192   fHistNegPt = new TH1F("fHistNegPt","",50,0.0,5.0);
 193   fHistNegPhi = new TH1F("fHistNegPhi","",63,0.0,6.3);
 194   fHistNegEta = new TH1F("fHistNegEta","",300,-1.5,1.5);
 195   fOutputList->Add(fHistNegPt);
 196   fOutputList->Add(fHistNegPhi);
 197   fOutputList->Add(fHistNegEta);
 198
 199
 200
 201   fHistPtPion = new TH1F("fHistPtPion","",50,0.0,5.0);
 202   fHistPtPionH = new TH1F("fHistPtPionH","",50,0.0,5.0);
 203   fHistPtPionL = new TH1F("fHistPtPionL","",50,0.0,5.0);
 204   fHistNsigmaPion = new TH1F("fHistNsigmaPion","",100,-5.0,5.0);
 205   fHistPtProt = new TH1F("fHistPtProt","",50,0.0,5.0);
 206   fHistPtProtH = new TH1F("fHistPtProtH","",50,0.0,5.0);
 207   fHistPtProtL = new TH1F("fHistPtProtL","",50,0.0,5.0);
 208   fHistNsigmaProt = new TH1F("fHistNsigmaProt","",100,-5.0,5.0);
 209   fOutputList->Add(fHistPtPion);
 210   fOutputList->Add(fHistPtPionH);
 211   fOutputList->Add(fHistPtPionL);
 212   fOutputList->Add(fHistNsigmaPion);
 213   fOutputList->Add(fHistPtProt);
 214   fOutputList->Add(fHistPtProtH);
 215   fOutputList->Add(fHistPtProtL);
 216   fOutputList->Add(fHistNsigmaProt);
 217
 218
 219   fHistPtMC = new TH1F("fHistPtMC","",50,0.0,5.0);
 220   fHistPhiMC = new TH1F("fHistPhiMC","",63,0.0,6.3);
 221   fHistEtaMC = new TH1F("fHistEtaMC","",110,-1.1,1.1);
 222   fHistChargeMC = new TH1F("fHistChargeMC","",3,-1.5,1.5);
 223   fHistTPCnclsMC = new TH1F("fHistTPCnclsMC","",161,-0.5,160.5);
 224   fHistDedxMC = new TH1F("fHistDedxMC","",1000,-20,980);
 225   fHistDCAxyMC = new TH1F("fHistDCAxyMC","",350,-3.5,3.5);
 226   fHistDCAzMC = new TH1F("fHistDCAzMC","",350,-3.5,3.5);
 227   fOutputList->Add(fHistPtMC);
 228   fOutputList->Add(fHistPhiMC);
 229   fOutputList->Add(fHistEtaMC);
 230   fOutputList->Add(fHistChargeMC);
 231   fOutputList->Add(fHistTPCnclsMC);
 232   fOutputList->Add(fHistDedxMC);
 233   fOutputList->Add(fHistDCAxyMC);
 234   fOutputList->Add(fHistDCAzMC);
 235
 236
 237
 238
 239   fHistPlaneV0h2 = new TH1F("fHistPlaneV0h2","",640,-3.2,3.2);
 240   fHistPlaneV0Ah2 = new TH1F("fHistPlaneV0Ah2","",640,-3.2,3.2);
 241   fHistPlaneV0Ch2 = new TH1F("fHistPlaneV0Ch2","",640,-3.2,3.2);
 242   fHistPlaneV0ACDCh2 = new TH1F("fHistPlaneV0ACDCh2","",640,-3.2,3.2);
 243   fOutputList->Add(fHistPlaneV0h2);
 244   fOutputList->Add(fHistPlaneV0Ah2);
 245   fOutputList->Add(fHistPlaneV0Ch2);
 246   fOutputList->Add(fHistPlaneV0ACDCh2);
 247
 248
 249   fHistZVtx = new TH1F("fHistZVtx","",80,-20,20);
 250   fHistZVtxD = new TH1F("fHistZVtxD","",100,-5,5);
 251   fHistCentTRK = new TH1F("fHistCentTRK","",100,0,100);
 252   fHistCentV0M = new TH1F("fHistCentV0M","",100,0,100);
 253   fHistCentDIFF = new TH1F("fHistCentDIFF","",162,-40.5,40.5);
 254   fHistCentDIAG = new TH1F("fHistCentDIAG","",6,-3.5,2.5);
 255   fOutputList->Add(fHistZVtx);
 256   fOutputList->Add(fHistZVtxD);
 257   fOutputList->Add(fHistCentTRK);
 258   fOutputList->Add(fHistCentV0M);
 259   fOutputList->Add(fHistCentDIFF);
 260   fOutputList->Add(fHistCentDIAG);
 261
 262   fHistZVtxMC = new TH1F("fHistZVtxMC","",80,-20,20);
 263   fHistZVtxDiffMC = new TH1F("fHistZVtxDiffMC","",100,-5,5);
 264   fOutputList->Add(fHistZVtxMC);
 265   fOutputList->Add(fHistZVtxDiffMC);
 266
 267   fHistCtrkDIAG = new TH1F("fHistCtrkDIAG","",1001,-0.5,1000.5);
 268   fHistVtxRDIAG = new TH1F("fHistVtxRDIAG","",100,-5,5);
 269   fHistVtxSDIAG = new TH1F("fHistVtxSDIAG","",100,-5,5);
 270   fHistVtxRDIBG = new TH1F("fHistVtxRDIBG","",100,-5,5);
 271   fHistVtxSDIBG = new TH1F("fHistVtxSDIBG","",100,-5,5);
 272   fOutputList->Add(fHistCtrkDIAG);
 273   fOutputList->Add(fHistVtxRDIAG);
 274   fOutputList->Add(fHistVtxSDIAG);
 275   fOutputList->Add(fHistVtxRDIBG);
 276   fOutputList->Add(fHistVtxSDIBG);
 277
 278
 279   fHistCentTRKAVEkMB = new TH1F("fHistCentTRKAVEkMB","",100,0,100);
 280   fHistCentV0MAVEkMB = new TH1F("fHistCentV0MAVEkMB","",100,0,100);
 281   fHistCentTRKAVEkCentral = new TH1F("fHistCentTRKAVEkCentral","",100,0,100);
 282   fHistCentV0MAVEkCentral = new TH1F("fHistCentV0MAVEkCentral","",100,0,100);
 283   fHistCentTRKAVEkSemiCentral = new TH1F("fHistCentTRKAVEkSemiCentral","",100,0,100);
 284   fHistCentV0MAVEkSemiCentral = new TH1F("fHistCentV0MAVEkSemiCentral","",100,0,100);
 285   fHistCentTRKAVEkA3 = new TH1F("fHistCentTRKAVEkA3","",100,0,100);
 286   fHistCentV0MAVEkA3 = new TH1F("fHistCentV0MAVEkA3","",100,0,100);
 287   fHistCentTRKAVEkSel = new TH1F("fHistCentTRKAVEkSel","",100,0,100);
 288   fHistCentV0MAVEkSel = new TH1F("fHistCentV0MAVEkSel","",100,0,100);
 289   fOutputList->Add(fHistCentTRKAVEkMB);
 290   fOutputList->Add(fHistCentV0MAVEkMB);
 291   fOutputList->Add(fHistCentTRKAVEkCentral);
 292   fOutputList->Add(fHistCentV0MAVEkCentral);
 293   fOutputList->Add(fHistCentTRKAVEkSemiCentral);
 294   fOutputList->Add(fHistCentV0MAVEkSemiCentral);
 295   fOutputList->Add(fHistCentTRKAVEkA3);
 296   fOutputList->Add(fHistCentV0MAVEkA3);
 297   fOutputList->Add(fHistCentTRKAVEkSel);
 298   fOutputList->Add(fHistCentV0MAVEkSel);
 299
 300
 301   // -------------------------- //
 302   // --- flow and asymmetry --- //
 303   // -------------------------- //
 304
 305   float tpmax = 1e10;
 306   float tpmin = -1e10;
 307
 308   fProfMeanChargePt = new TProfile("fProfMeanChargePt","",50,0,5.0,tpmin,tpmax,"");
 309   fProfMeanChargeEta = new TProfile("fProfMeanChargeEta","",300,-1.5,1.5,tpmin,tpmax,"");
 310
 311   h_eta_pos_F1 = new TH1F("h_eta_pos_F1","",300,-1.5,1.5);
 312   h_eta_neg_F1 = new TH1F("h_eta_neg_F1","",300,-1.5,1.5);
 313   h_eta_pos_F3 = new TH1F("h_eta_pos_F3","",300,-1.5,1.5);
 314   h_eta_neg_F3 = new TH1F("h_eta_neg_F3","",300,-1.5,1.5);
 315   fOutputList->Add(h_eta_pos_F1);
 316   fOutputList->Add(h_eta_neg_F1);
 317   fOutputList->Add(h_eta_pos_F3);
 318   fOutputList->Add(h_eta_neg_F3);
 319   h_cut_eta_pos_F1 = new TH1F("h_cut_eta_pos_F1","",300,-1.5,1.5);
 320   h_cut_eta_neg_F1 = new TH1F("h_cut_eta_neg_F1","",300,-1.5,1.5);
 321   h_cut_eta_pos_F3 = new TH1F("h_cut_eta_pos_F3","",300,-1.5,1.5);
 322   h_cut_eta_neg_F3 = new TH1F("h_cut_eta_neg_F3","",300,-1.5,1.5);
 323   fOutputList->Add(h_cut_eta_pos_F1);
 324   fOutputList->Add(h_cut_eta_neg_F1);
 325   fOutputList->Add(h_cut_eta_pos_F3);
 326   fOutputList->Add(h_cut_eta_neg_F3);
 327
 328   // --- charge asymmetry (average charge in event)
 329   h_AT_eta = new TProfile("h_AT_eta","",150,-1.5,1.5,tpmin,tpmax,"");
 330   h_AT_etaF1 = new TProfile("h_AT_etaF1","",150,-1.5,1.5,tpmin,tpmax,"");
 331   h_AT_etaF3 = new TProfile("h_AT_etaF3","",150,-1.5,1.5,tpmin,tpmax,"");
 332   h_AT_etaMC = new TProfile("h_AT_etaMC","",150,-1.5,1.5,tpmin,tpmax,"");
 333   fOutputList->Add(h_AT_eta);
 334   fOutputList->Add(h_AT_etaF1);
 335   fOutputList->Add(h_AT_etaF3);
 336   fOutputList->Add(h_AT_etaMC);
 337
 338   // --- charge asymmetry (average charge in event)
 339   h2_AT_eta = new TH2F("h2_AT_eta","",150,-1.5,1.5,100,-1.0,1.0);
 340   h2_AT_etaF1 = new TH2F("h2_AT_etaF1","",150,-1.5,1.5,100,-1.0,1.0);
 341   h2_AT_etaF3 = new TH2F("h2_AT_etaF3","",150,-1.5,1.5,100,-1.0,1.0);
 342   h2_AT_etaMC = new TH2F("h2_AT_etaMC","",150,-1.5,1.5,100,-1.0,1.0);
 343   fOutputList->Add(h2_AT_eta);
 344   fOutputList->Add(h2_AT_etaF1);
 345   fOutputList->Add(h2_AT_etaF3);
 346   fOutputList->Add(h2_AT_etaMC);
 347
 348   // --- charge asymmetry (average charge in event)
 349   h_AT_cut_eta = new TProfile("h_AT_cut_eta","",150,-1.5,1.5,tpmin,tpmax,"");
 350   h_AT_cut_etaF1 = new TProfile("h_AT_cut_etaF1","",150,-1.5,1.5,tpmin,tpmax,"");
 351   h_AT_cut_etaF3 = new TProfile("h_AT_cut_etaF3","",150,-1.5,1.5,tpmin,tpmax,"");
 352   h_AT_cut_etaMC = new TProfile("h_AT_cut_etaMC","",150,-1.5,1.5,tpmin,tpmax,"");
 353   fOutputList->Add(h_AT_cut_eta);
 354   fOutputList->Add(h_AT_cut_etaF1);
 355   fOutputList->Add(h_AT_cut_etaF3);
 356   fOutputList->Add(h_AT_cut_etaMC);
 357
 358   // --- charge asymmetry (average charge in event)
 359   h2_AT_cut_eta = new TH2F("h2_AT_cut_eta","",150,-1.5,1.5,100,-1.0,1.0);
 360   h2_AT_cut_etaF1 = new TH2F("h2_AT_cut_etaF1","",150,-1.5,1.5,100,-1.0,1.0);
 361   h2_AT_cut_etaF3 = new TH2F("h2_AT_cut_etaF3","",150,-1.5,1.5,100,-1.0,1.0);
 362   h2_AT_cut_etaMC = new TH2F("h2_AT_cut_etaMC","",150,-1.5,1.5,100,-1.0,1.0);
 363   fOutputList->Add(h2_AT_cut_eta);
 364   fOutputList->Add(h2_AT_cut_etaF1);
 365   fOutputList->Add(h2_AT_cut_etaF3);
 366   fOutputList->Add(h2_AT_cut_etaMC);
 367
 368   // --- charge asymmetry (average charge in event)
 369   h_A_cent = new TProfile("h_A_cent","",100,0,100,tpmin,tpmax,"");
 370   h_rA_cent = new TProfile("h_rA_cent","",100,0,100,tpmin,tpmax,"");
 371   h_r1A_cent = new TProfile("h_r1A_cent","",100,0,100,tpmin,tpmax,"");
 372   h_r2A_cent = new TProfile("h_r2A_cent","",100,0,100,tpmin,tpmax,"");
 373   h_A_centF1 = new TProfile("h_A_centF1","",100,0,100,tpmin,tpmax,"");
 374   h_A_centF3 = new TProfile("h_A_centF3","",100,0,100,tpmin,tpmax,"");
 375   h_A_centMC = new TProfile("h_A_centMC","",100,0,100,tpmin,tpmax,"");
 376   fOutputList->Add(h_A_cent);
 377   fOutputList->Add(h_rA_cent);
 378   fOutputList->Add(h_r1A_cent);
 379   fOutputList->Add(h_r2A_cent);
 380   fOutputList->Add(h_A_centF1);
 381   fOutputList->Add(h_A_centF3);
 382   fOutputList->Add(h_A_centMC);
 383
 384   // --- charge asymmetry (average charge in event)
 385   h2_A_cent = new TH2F("h2_A_cent","",100,0,100,100,-1.0,1.0);
 386   h2_rA_cent = new TH2F("h2_rA_cent","",100,0,100,100,-1.0,1.0);
 387   h2_r1A_cent = new TH2F("h2_r1A_cent","",100,0,100,100,-1.0,1.0);
 388   h2_r2A_cent = new TH2F("h2_r2A_cent","",100,0,100,100,-1.0,1.0);
 389   h2_A_centF1 = new TH2F("h2_A_centF1","",100,0,100,100,-1.0,1.0);
 390   h2_A_centF3 = new TH2F("h2_A_centF3","",100,0,100,100,-1.0,1.0);
 391   h2_A_centMC = new TH2F("h2_A_centMC","",100,0,100,100,-1.0,1.0);
 392   fOutputList->Add(h2_A_cent);
 393   fOutputList->Add(h2_rA_cent);
 394   fOutputList->Add(h2_r1A_cent);
 395   fOutputList->Add(h2_r2A_cent);
 396   fOutputList->Add(h2_A_centF1);
 397   fOutputList->Add(h2_A_centF3);
 398   fOutputList->Add(h2_A_centMC);
 399
 400   // --- MONTE CARLO second harmonic vs centrality
 401   h_MCq22_cent = new TProfile("h_MCq22_cent","",100,0,100,tpmin,tpmax,"");
 402   h_MCq22_centP = new TProfile("h_MCq22_centP","",100,0,100,tpmin,tpmax,"");
 403   h_MCq22_centN = new TProfile("h_MCq22_centN","",100,0,100,tpmin,tpmax,"");
 404   fOutputList->Add(h_MCq22_cent);
 405   fOutputList->Add(h_MCq22_centP);
 406   fOutputList->Add(h_MCq22_centN);
 407   h_MCAq22_cent = new TProfile("h_MCAq22_cent","",100,0,100,tpmin,tpmax,"");
 408   h_MCAq22_centP = new TProfile("h_MCAq22_centP","",100,0,100,tpmin,tpmax,"");
 409   h_MCAq22_centN = new TProfile("h_MCAq22_centN","",100,0,100,tpmin,tpmax,"");
 410   fOutputList->Add(h_MCAq22_cent);
 411   fOutputList->Add(h_MCAq22_centP);
 412   fOutputList->Add(h_MCAq22_centN);
 413   h_MCq22gap0_cent = new TProfile("h_MCq22gap0_cent","",100,0,100,tpmin,tpmax,"");
 414   h_MCq22gap0_centP = new TProfile("h_MCq22gap0_centP","",100,0,100,tpmin,tpmax,"");
 415   h_MCq22gap0_centN = new TProfile("h_MCq22gap0_centN","",100,0,100,tpmin,tpmax,"");
 416   fOutputList->Add(h_MCq22gap0_cent);
 417   fOutputList->Add(h_MCq22gap0_centP);
 418   fOutputList->Add(h_MCq22gap0_centN);
 419   h_MCq22gap1_cent = new TProfile("h_MCq22gap1_cent","",100,0,100,tpmin,tpmax,"");
 420   h_MCq22gap1_centP = new TProfile("h_MCq22gap1_centP","",100,0,100,tpmin,tpmax,"");
 421   h_MCq22gap1_centN = new TProfile("h_MCq22gap1_centN","",100,0,100,tpmin,tpmax,"");
 422   fOutputList->Add(h_MCq22gap1_cent);
 423   fOutputList->Add(h_MCq22gap1_centP);
 424   fOutputList->Add(h_MCq22gap1_centN);
 425
 426   hMC_AT_X_deta = new TProfile("hMC_AT_X_deta","",32,-1.6,1.6,tpmin,tpmax,"");
 427   hMC_AT_X_detaP = new TProfile("hMC_AT_X_detaP","",32,-1.6,1.6,tpmin,tpmax,"");
 428   hMC_AT_X_detaN = new TProfile("hMC_AT_X_detaN","",32,-1.6,1.6,tpmin,tpmax,"");
 429   fOutputList->Add(hMC_AT_X_deta);
 430   fOutputList->Add(hMC_AT_X_detaP);
 431   fOutputList->Add(hMC_AT_X_detaN);
 432   //hMC_diffq22_X_deta = new TProfile("hMC_diffq22_X_deta","",32,-1.6,1.6,tpmin,tpmax,"");
 433   hMC_diffq22_X_detaP = new TProfile("hMC_diffq22_X_detaP","",32,-1.6,1.6,tpmin,tpmax,"");
 434   hMC_diffq22_X_detaN = new TProfile("hMC_diffq22_X_detaN","",32,-1.6,1.6,tpmin,tpmax,"");
 435   //hMC_ATdiffq22_X_deta = new TProfile("hMC_ATdiffq22_X_deta","",32,-1.6,1.6,tpmin,tpmax,"");
 436   hMC_ATdiffq22_X_detaP = new TProfile("hMC_ATdiffq22_X_detaP","",32,-1.6,1.6,tpmin,tpmax,"");
 437   hMC_ATdiffq22_X_detaN = new TProfile("hMC_ATdiffq22_X_detaN","",32,-1.6,1.6,tpmin,tpmax,"");
 438   fOutputList->Add(hMC_diffq22_X_detaP);
 439   fOutputList->Add(hMC_diffq22_X_detaN);
 440   fOutputList->Add(hMC_ATdiffq22_X_detaP);
 441   fOutputList->Add(hMC_ATdiffq22_X_detaN);
 442   //hMC_diffq32_X_deta = new TProfile("hMC_diffq32_X_deta","",32,-1.6,1.6,tpmin,tpmax,"");
 443   hMC_diffq32_X_detaP = new TProfile("hMC_diffq32_X_detaP","",32,-1.6,1.6,tpmin,tpmax,"");
 444   hMC_diffq32_X_detaN = new TProfile("hMC_diffq32_X_detaN","",32,-1.6,1.6,tpmin,tpmax,"");
 445   //hMC_ATdiffq32_X_deta = new TProfile("hMC_ATdiffq32_X_deta","",32,-1.6,1.6,tpmin,tpmax,"");
 446   hMC_ATdiffq32_X_detaP = new TProfile("hMC_ATdiffq32_X_detaP","",32,-1.6,1.6,tpmin,tpmax,"");
 447   hMC_ATdiffq32_X_detaN = new TProfile("hMC_ATdiffq32_X_detaN","",32,-1.6,1.6,tpmin,tpmax,"");
 448   fOutputList->Add(hMC_diffq32_X_detaP);
 449   fOutputList->Add(hMC_diffq32_X_detaN);
 450   fOutputList->Add(hMC_ATdiffq32_X_detaP);
 451   fOutputList->Add(hMC_ATdiffq32_X_detaN);
 452   //hMC_diffq42_X_deta = new TProfile("hMC_diffq42_X_deta","",32,-1.6,1.6,tpmin,tpmax,"");
 453   hMC_diffq42_X_detaP = new TProfile("hMC_diffq42_X_detaP","",32,-1.6,1.6,tpmin,tpmax,"");
 454   hMC_diffq42_X_detaN = new TProfile("hMC_diffq42_X_detaN","",32,-1.6,1.6,tpmin,tpmax,"");
 455   //hMC_ATdiffq42_X_deta = new TProfile("hMC_ATdiffq42_X_deta","",32,-1.6,1.6,tpmin,tpmax,"");
 456   hMC_ATdiffq42_X_detaP = new TProfile("hMC_ATdiffq42_X_detaP","",32,-1.6,1.6,tpmin,tpmax,"");
 457   hMC_ATdiffq42_X_detaN = new TProfile("hMC_ATdiffq42_X_detaN","",32,-1.6,1.6,tpmin,tpmax,"");
 458   fOutputList->Add(hMC_diffq42_X_detaP);
 459   fOutputList->Add(hMC_diffq42_X_detaN);
 460   fOutputList->Add(hMC_ATdiffq42_X_detaP);
 461   fOutputList->Add(hMC_ATdiffq42_X_detaN);
 462
 463
 464   // --- q-vector components to check recentering
 465   h_X22_cent = new TProfile("h_X22_cent","",100,0,100,tpmin,tpmax,"");
 466   h_X22_centP = new TProfile("h_X22_centP","",100,0,100,tpmin,tpmax,"");
 467   h_X22_centN = new TProfile("h_X22_centN","",100,0,100,tpmin,tpmax,"");
 468   h_Y22_cent = new TProfile("h_Y22_cent","",100,0,100,tpmin,tpmax,"");
 469   h_Y22_centP = new TProfile("h_Y22_centP","",100,0,100,tpmin,tpmax,"");
 470   h_Y22_centN = new TProfile("h_Y22_centN","",100,0,100,tpmin,tpmax,"");
 471   fOutputList->Add(h_X22_cent);
 472   fOutputList->Add(h_X22_centP);
 473   fOutputList->Add(h_X22_centN);
 474   fOutputList->Add(h_Y22_cent);
 475   fOutputList->Add(h_Y22_centP);
 476   fOutputList->Add(h_Y22_centN);
 477   // --- outer left
 478   h_X22_lo_cent = new TProfile("h_X22_lo_cent","",100,0,100,tpmin,tpmax,"");
 479   h_X22_lo_centP = new TProfile("h_X22_lo_centP","",100,0,100,tpmin,tpmax,"");
 480   h_X22_lo_centN = new TProfile("h_X22_lo_centN","",100,0,100,tpmin,tpmax,"");
 481   h_Y22_lo_cent = new TProfile("h_Y22_lo_cent","",100,0,100,tpmin,tpmax,"");
 482   h_Y22_lo_centP = new TProfile("h_Y22_lo_centP","",100,0,100,tpmin,tpmax,"");
 483   h_Y22_lo_centN = new TProfile("h_Y22_lo_centN","",100,0,100,tpmin,tpmax,"");
 484   fOutputList->Add(h_X22_lo_cent);
 485   fOutputList->Add(h_X22_lo_centP);
 486   fOutputList->Add(h_X22_lo_centN);
 487   fOutputList->Add(h_Y22_lo_cent);
 488   fOutputList->Add(h_Y22_lo_centP);
 489   fOutputList->Add(h_Y22_lo_centN);
 490   // --- inner left
 491   h_X22_li_cent = new TProfile("h_X22_li_cent","",100,0,100,tpmin,tpmax,"");
 492   h_X22_li_centP = new TProfile("h_X22_li_centP","",100,0,100,tpmin,tpmax,"");
 493   h_X22_li_centN = new TProfile("h_X22_li_centN","",100,0,100,tpmin,tpmax,"");
 494   h_Y22_li_cent = new TProfile("h_Y22_li_cent","",100,0,100,tpmin,tpmax,"");
 495   h_Y22_li_centP = new TProfile("h_Y22_li_centP","",100,0,100,tpmin,tpmax,"");
 496   h_Y22_li_centN = new TProfile("h_Y22_li_centN","",100,0,100,tpmin,tpmax,"");
 497   fOutputList->Add(h_X22_li_cent);
 498   fOutputList->Add(h_X22_li_centP);
 499   fOutputList->Add(h_X22_li_centN);
 500   fOutputList->Add(h_Y22_li_cent);
 501   fOutputList->Add(h_Y22_li_centP);
 502   fOutputList->Add(h_Y22_li_centN);
 503   // --- inner right
 504   h_X22_ri_cent = new TProfile("h_X22_ri_cent","",100,0,100,tpmin,tpmax,"");
 505   h_X22_ri_centP = new TProfile("h_X22_ri_centP","",100,0,100,tpmin,tpmax,"");
 506   h_X22_ri_centN = new TProfile("h_X22_ri_centN","",100,0,100,tpmin,tpmax,"");
 507   h_Y22_ri_cent = new TProfile("h_Y22_ri_cent","",100,0,100,tpmin,tpmax,"");
 508   h_Y22_ri_centP = new TProfile("h_Y22_ri_centP","",100,0,100,tpmin,tpmax,"");
 509   h_Y22_ri_centN = new TProfile("h_Y22_ri_centN","",100,0,100,tpmin,tpmax,"");
 510   fOutputList->Add(h_X22_ri_cent);
 511   fOutputList->Add(h_X22_ri_centP);
 512   fOutputList->Add(h_X22_ri_centN);
 513   fOutputList->Add(h_Y22_ri_cent);
 514   fOutputList->Add(h_Y22_ri_centP);
 515   fOutputList->Add(h_Y22_ri_centN);
 516   // --- outer right
 517   h_X22_ro_cent = new TProfile("h_X22_ro_cent","",100,0,100,tpmin,tpmax,"");
 518   h_X22_ro_centP = new TProfile("h_X22_ro_centP","",100,0,100,tpmin,tpmax,"");
 519   h_X22_ro_centN = new TProfile("h_X22_ro_centN","",100,0,100,tpmin,tpmax,"");
 520   h_Y22_ro_cent = new TProfile("h_Y22_ro_cent","",100,0,100,tpmin,tpmax,"");
 521   h_Y22_ro_centP = new TProfile("h_Y22_ro_centP","",100,0,100,tpmin,tpmax,"");
 522   h_Y22_ro_centN = new TProfile("h_Y22_ro_centN","",100,0,100,tpmin,tpmax,"");
 523   fOutputList->Add(h_X22_ro_cent);
 524   fOutputList->Add(h_X22_ro_centP);
 525   fOutputList->Add(h_X22_ro_centN);
 526   fOutputList->Add(h_Y22_ro_cent);
 527   fOutputList->Add(h_Y22_ro_centP);
 528   fOutputList->Add(h_Y22_ro_centN);
 529
 530   // --- second harmonic vs centrality
 531   h_q22_cent = new TProfile("h_q22_cent","",100,0,100,tpmin,tpmax,"");
 532   h_q22_centP = new TProfile("h_q22_centP","",100,0,100,tpmin,tpmax,"");
 533   h_q22_centN = new TProfile("h_q22_centN","",100,0,100,tpmin,tpmax,"");
 534   fOutputList->Add(h_q22_cent);
 535   fOutputList->Add(h_q22_centP);
 536   fOutputList->Add(h_q22_centN);
 537   h_q23_cent = new TProfile("h_q23_cent","",100,0,100,tpmin,tpmax,"");
 538   h_q23_centP = new TProfile("h_q23_centP","",100,0,100,tpmin,tpmax,"");
 539   h_q23_centN = new TProfile("h_q23_centN","",100,0,100,tpmin,tpmax,"");
 540   fOutputList->Add(h_q23_cent);
 541   fOutputList->Add(h_q23_centP);
 542   fOutputList->Add(h_q23_centN);
 543   h_q24_cent = new TProfile("h_q24_cent","",100,0,100,tpmin,tpmax,"");
 544   h_q24_centP = new TProfile("h_q24_centP","",100,0,100,tpmin,tpmax,"");
 545   h_q24_centN = new TProfile("h_q24_centN","",100,0,100,tpmin,tpmax,"");
 546   fOutputList->Add(h_q24_cent);
 547   fOutputList->Add(h_q24_centP);
 548   fOutputList->Add(h_q24_centN);
 549   h_q22gap0_cent = new TProfile("h_q22gap0_cent","",100,0,100,tpmin,tpmax,"");
 550   h_q22gap0_centP = new TProfile("h_q22gap0_centP","",100,0,100,tpmin,tpmax,"");
 551   h_q22gap0_centN = new TProfile("h_q22gap0_centN","",100,0,100,tpmin,tpmax,"");
 552   fOutputList->Add(h_q22gap0_cent);
 553   fOutputList->Add(h_q22gap0_centP);
 554   fOutputList->Add(h_q22gap0_centN);
 555   h_q22gap1_cent = new TProfile("h_q22gap1_cent","",100,0,100,tpmin,tpmax,"");
 556   h_q22gap1_centP = new TProfile("h_q22gap1_centP","",100,0,100,tpmin,tpmax,"");
 557   h_q22gap1_centN = new TProfile("h_q22gap1_centN","",100,0,100,tpmin,tpmax,"");
 558   fOutputList->Add(h_q22gap1_cent);
 559   fOutputList->Add(h_q22gap1_centP);
 560   fOutputList->Add(h_q22gap1_centN);
 561
 562   // --- third harmonic vs centrality
 563   h_q32_cent = new TProfile("h_q32_cent","",100,0,100,tpmin,tpmax,"");
 564   h_q32_centP = new TProfile("h_q32_centP","",100,0,100,tpmin,tpmax,"");
 565   h_q32_centN = new TProfile("h_q32_centN","",100,0,100,tpmin,tpmax,"");
 566   fOutputList->Add(h_q32_cent);
 567   fOutputList->Add(h_q32_centP);
 568   fOutputList->Add(h_q32_centN);
 569   h_q32gap0_cent = new TProfile("h_q32gap0_cent","",100,0,100,tpmin,tpmax,"");
 570   h_q32gap0_centP = new TProfile("h_q32gap0_centP","",100,0,100,tpmin,tpmax,"");
 571   h_q32gap0_centN = new TProfile("h_q32gap0_centN","",100,0,100,tpmin,tpmax,"");
 572   fOutputList->Add(h_q32gap0_cent);
 573   fOutputList->Add(h_q32gap0_centP);
 574   fOutputList->Add(h_q32gap0_centN);
 575   h_q32gap1_cent = new TProfile("h_q32gap1_cent","",100,0,100,tpmin,tpmax,"");
 576   h_q32gap1_centP = new TProfile("h_q32gap1_centP","",100,0,100,tpmin,tpmax,"");
 577   h_q32gap1_centN = new TProfile("h_q32gap1_centN","",100,0,100,tpmin,tpmax,"");
 578   fOutputList->Add(h_q32gap1_cent);
 579   fOutputList->Add(h_q32gap1_centP);
 580   fOutputList->Add(h_q32gap1_centN);
 581
 582   // --- fourth harmonic vs centrality
 583   h_q42_cent = new TProfile("h_q42_cent","",100,0,100,tpmin,tpmax,"");
 584   h_q42_centP = new TProfile("h_q42_centP","",100,0,100,tpmin,tpmax,"");
 585   h_q42_centN = new TProfile("h_q42_centN","",100,0,100,tpmin,tpmax,"");
 586   fOutputList->Add(h_q42_cent);
 587   fOutputList->Add(h_q42_centP);
 588   fOutputList->Add(h_q42_centN);
 589   h_q42gap0_cent = new TProfile("h_q42gap0_cent","",100,0,100,tpmin,tpmax,"");
 590   h_q42gap0_centP = new TProfile("h_q42gap0_centP","",100,0,100,tpmin,tpmax,"");
 591   h_q42gap0_centN = new TProfile("h_q42gap0_centN","",100,0,100,tpmin,tpmax,"");
 592   fOutputList->Add(h_q42gap0_cent);
 593   fOutputList->Add(h_q42gap0_centP);
 594   fOutputList->Add(h_q42gap0_centN);
 595   h_q42gap1_cent = new TProfile("h_q42gap1_cent","",100,0,100,tpmin,tpmax,"");
 596   h_q42gap1_centP = new TProfile("h_q42gap1_centP","",100,0,100,tpmin,tpmax,"");
 597   h_q42gap1_centN = new TProfile("h_q42gap1_centN","",100,0,100,tpmin,tpmax,"");
 598   fOutputList->Add(h_q42gap1_cent);
 599   fOutputList->Add(h_q42gap1_centP);
 600   fOutputList->Add(h_q42gap1_centN);
 601
 602   // harmonics vs charge asymmetry in centrality bins
 603   for(int i=0; i<10; i++)
 604     {
 605       // --- second harmonic
 606       h_q22qasym_cent[i] = new TProfile(Form("h_q22qasym_cent_%d",i),Form("h_q22qasym_cent_%d",i),220,-1.1,1.1,tpmin,tpmax,"");
 607       h_q22qasymP_cent[i] = new TProfile(Form("h_q22qasymP_cent_%d",i),Form("h_q22qasymP_cent_%d",i),220,-1.1,1.1,tpmin,tpmax,"");
 608       h_q22qasymN_cent[i] = new TProfile(Form("h_q22qasymN_cent_%d",i),Form("h_q22qasymN_cent_%d",i),220,-1.1,1.1,tpmin,tpmax,"");
 609       h_q22qasym_gap0_cent[i] = new TProfile(Form("h_q22qasym_gap0_cent_%d",i),Form("h_q22qasym_gap0_cent_%d",i),220,-1.1,1.1,tpmin,tpmax,"");
 610       h_q22qasymP_gap0_cent[i] = new TProfile(Form("h_q22qasymP_gap0_cent_%d",i),Form("h_q22qasymP_gap0_cent_%d",i),220,-1.1,1.1,tpmin,tpmax,"");
 611       h_q22qasymN_gap0_cent[i] = new TProfile(Form("h_q22qasymN_gap0_cent_%d",i),Form("h_q22qasymN_gap0_cent_%d",i),220,-1.1,1.1,tpmin,tpmax,"");
 612       h_q22qasym_gap1_cent[i] = new TProfile(Form("h_q22qasym_gap1_cent_%d",i),Form("h_q22qasym_gap1_cent_%d",i),220,-1.1,1.1,tpmin,tpmax,"");
 613       h_q22qasymP_gap1_cent[i] = new TProfile(Form("h_q22qasymP_gap1_cent_%d",i),Form("h_q22qasymP_gap1_cent_%d",i),220,-1.1,1.1,tpmin,tpmax,"");
 614       h_q22qasymN_gap1_cent[i] = new TProfile(Form("h_q22qasymN_gap1_cent_%d",i),Form("h_q22qasymN_gap1_cent_%d",i),220,-1.1,1.1,tpmin,tpmax,"");
 615       h_q24qasym_cent[i] = new TProfile(Form("h_q24qasym_cent_%d",i),Form("h_q24qasym_cent_%d",i),220,-1.1,1.1,tpmin,tpmax,"");
 616       h_q24qasymP_cent[i] = new TProfile(Form("h_q24qasymP_cent_%d",i),Form("h_q24qasymP_cent_%d",i),220,-1.1,1.1,tpmin,tpmax,"");
 617       h_q24qasymN_cent[i] = new TProfile(Form("h_q24qasymN_cent_%d",i),Form("h_q24qasymN_cent_%d",i),220,-1.1,1.1,tpmin,tpmax,"");
 618
 619       // --- third harmonic
 620       h_q32qasym_cent[i] = new TProfile(Form("h_q32qasym_cent_%d",i),Form("h_q32qasym_cent_%d",i),220,-1.1,1.1,tpmin,tpmax,"");
 621       h_q32qasymP_cent[i] = new TProfile(Form("h_q32qasymP_cent_%d",i),Form("h_q32qasymP_cent_%d",i),220,-1.1,1.1,tpmin,tpmax,"");
 622       h_q32qasymN_cent[i] = new TProfile(Form("h_q32qasymN_cent_%d",i),Form("h_q32qasymN_cent_%d",i),220,-1.1,1.1,tpmin,tpmax,"");
 623       h_q32qasym_gap0_cent[i] = new TProfile(Form("h_q32qasym_gap0_cent_%d",i),Form("h_q32qasym_gap0_cent_%d",i),220,-1.1,1.1,tpmin,tpmax,"");
 624       h_q32qasymP_gap0_cent[i] = new TProfile(Form("h_q32qasymP_gap0_cent_%d",i),Form("h_q32qasymP_gap0_cent_%d",i),220,-1.1,1.1,tpmin,tpmax,"");
 625       h_q32qasymN_gap0_cent[i] = new TProfile(Form("h_q32qasymN_gap0_cent_%d",i),Form("h_q32qasymN_gap0_cent_%d",i),220,-1.1,1.1,tpmin,tpmax,"");
 626       h_q32qasym_gap1_cent[i] = new TProfile(Form("h_q32qasym_gap1_cent_%d",i),Form("h_q32qasym_gap1_cent_%d",i),220,-1.1,1.1,tpmin,tpmax,"");
 627       h_q32qasymP_gap1_cent[i] = new TProfile(Form("h_q32qasymP_gap1_cent_%d",i),Form("h_q32qasymP_gap1_cent_%d",i),220,-1.1,1.1,tpmin,tpmax,"");
 628       h_q32qasymN_gap1_cent[i] = new TProfile(Form("h_q32qasymN_gap1_cent_%d",i),Form("h_q32qasymN_gap1_cent_%d",i),220,-1.1,1.1,tpmin,tpmax,"");
 629       h_q34qasym_cent[i] = new TProfile(Form("h_q34qasym_cent_%d",i),Form("h_q34qasym_cent_%d",i),220,-1.1,1.1,tpmin,tpmax,"");
 630       h_q34qasymP_cent[i] = new TProfile(Form("h_q34qasymP_cent_%d",i),Form("h_q34qasymP_cent_%d",i),220,-1.1,1.1,tpmin,tpmax,"");
 631       h_q34qasymN_cent[i] = new TProfile(Form("h_q34qasymN_cent_%d",i),Form("h_q34qasymN_cent_%d",i),220,-1.1,1.1,tpmin,tpmax,"");
 632
 633       // --- fourth harmonic
 634       h_q42qasym_cent[i] = new TProfile(Form("h_q42qasym_cent_%d",i),Form("h_q42qasym_cent_%d",i),220,-1.1,1.1,tpmin,tpmax,"");
 635       h_q42qasymP_cent[i] = new TProfile(Form("h_q42qasymP_cent_%d",i),Form("h_q42qasymP_cent_%d",i),220,-1.1,1.1,tpmin,tpmax,"");
 636       h_q42qasymN_cent[i] = new TProfile(Form("h_q42qasymN_cent_%d",i),Form("h_q42qasymN_cent_%d",i),220,-1.1,1.1,tpmin,tpmax,"");
 637       h_q42qasym_gap0_cent[i] = new TProfile(Form("h_q42qasym_gap0_cent_%d",i),Form("h_q42qasym_gap0_cent_%d",i),220,-1.1,1.1,tpmin,tpmax,"");
 638       h_q42qasymP_gap0_cent[i] = new TProfile(Form("h_q42qasymP_gap0_cent_%d",i),Form("h_q42qasymP_gap0_cent_%d",i),220,-1.1,1.1,tpmin,tpmax,"");
 639       h_q42qasymN_gap0_cent[i] = new TProfile(Form("h_q42qasymN_gap0_cent_%d",i),Form("h_q42qasymN_gap0_cent_%d",i),220,-1.1,1.1,tpmin,tpmax,"");
 640       h_q42qasym_gap1_cent[i] = new TProfile(Form("h_q42qasym_gap1_cent_%d",i),Form("h_q42qasym_gap1_cent_%d",i),220,-1.1,1.1,tpmin,tpmax,"");
 641       h_q42qasymP_gap1_cent[i] = new TProfile(Form("h_q42qasymP_gap1_cent_%d",i),Form("h_q42qasymP_gap1_cent_%d",i),220,-1.1,1.1,tpmin,tpmax,"");
 642       h_q42qasymN_gap1_cent[i] = new TProfile(Form("h_q42qasymN_gap1_cent_%d",i),Form("h_q42qasymN_gap1_cent_%d",i),220,-1.1,1.1,tpmin,tpmax,"");
 643       h_q44qasym_cent[i] = new TProfile(Form("h_q44qasym_cent_%d",i),Form("h_q44qasym_cent_%d",i),220,-1.1,1.1,tpmin,tpmax,"");
 644       h_q44qasymP_cent[i] = new TProfile(Form("h_q44qasymP_cent_%d",i),Form("h_q44qasymP_cent_%d",i),220,-1.1,1.1,tpmin,tpmax,"");
 645       h_q44qasymN_cent[i] = new TProfile(Form("h_q44qasymN_cent_%d",i),Form("h_q44qasymN_cent_%d",i),220,-1.1,1.1,tpmin,tpmax,"");
 646     }
 647
 648
 649   // harmonics vs charge asymmetry in centrality bins
 650   for(int i=cbinlo; i<cbinhi; i++)
 651     {
 652       // --- second harmonic
 653       fOutputList->Add(h_q22qasym_cent[i]);
 654       fOutputList->Add(h_q22qasymP_cent[i]);
 655       fOutputList->Add(h_q22qasymN_cent[i]);
 656       fOutputList->Add(h_q22qasym_gap0_cent[i]);
 657       fOutputList->Add(h_q22qasymP_gap0_cent[i]);
 658       fOutputList->Add(h_q22qasymN_gap0_cent[i]);
 659       fOutputList->Add(h_q22qasym_gap1_cent[i]);
 660       fOutputList->Add(h_q22qasymP_gap1_cent[i]);
 661       fOutputList->Add(h_q22qasymN_gap1_cent[i]);
 662       fOutputList->Add(h_q24qasym_cent[i]);
 663       fOutputList->Add(h_q24qasymP_cent[i]);
 664       fOutputList->Add(h_q24qasymN_cent[i]);
 665
 666       // --- third harmonic
 667       fOutputList->Add(h_q32qasym_cent[i]);
 668       fOutputList->Add(h_q32qasymP_cent[i]);
 669       fOutputList->Add(h_q32qasymN_cent[i]);
 670       fOutputList->Add(h_q32qasym_gap0_cent[i]);
 671       fOutputList->Add(h_q32qasymP_gap0_cent[i]);
 672       fOutputList->Add(h_q32qasymN_gap0_cent[i]);
 673       fOutputList->Add(h_q32qasym_gap1_cent[i]);
 674       fOutputList->Add(h_q32qasymP_gap1_cent[i]);
 675       fOutputList->Add(h_q32qasymN_gap1_cent[i]);
 676       fOutputList->Add(h_q34qasym_cent[i]);
 677       fOutputList->Add(h_q34qasymP_cent[i]);
 678       fOutputList->Add(h_q34qasymN_cent[i]);
 679
 680       // --- fourth harmonic
 681       fOutputList->Add(h_q42qasym_cent[i]);
 682       fOutputList->Add(h_q42qasymP_cent[i]);
 683       fOutputList->Add(h_q42qasymN_cent[i]);
 684       fOutputList->Add(h_q42qasym_gap0_cent[i]);
 685       fOutputList->Add(h_q42qasymP_gap0_cent[i]);
 686       fOutputList->Add(h_q42qasymN_gap0_cent[i]);
 687       fOutputList->Add(h_q42qasym_gap1_cent[i]);
 688       fOutputList->Add(h_q42qasymP_gap1_cent[i]);
 689       fOutputList->Add(h_q42qasymN_gap1_cent[i]);
 690       fOutputList->Add(h_q44qasym_cent[i]);
 691       fOutputList->Add(h_q44qasymP_cent[i]);
 692       fOutputList->Add(h_q44qasymN_cent[i]);
 693     }
 694
 695
 696
 697   // --------------------------------------- //
 698   // --- three-particle-like correlators --- //
 699   // --------------------------------------- //
 700
 701
 702   // --- second harmonic
 703   h_Aq22_cent = new TProfile("h_Aq22_cent","",100,0,100,tpmin,tpmax,"");
 704   h_Aq22_centP = new TProfile("h_Aq22_centP","",100,0,100,tpmin,tpmax,"");
 705   h_Aq22_centN = new TProfile("h_Aq22_centN","",100,0,100,tpmin,tpmax,"");
 706   fOutputList->Add(h_Aq22_cent);
 707   fOutputList->Add(h_Aq22_centP);
 708   fOutputList->Add(h_Aq22_centN);
 709   h_Aq22_SL_cent = new TProfile("h_Aq22_SL_cent","",100,0,100,tpmin,tpmax,"");
 710   h_Aq22_SL_centP = new TProfile("h_Aq22_SL_centP","",100,0,100,tpmin,tpmax,"");
 711   h_Aq22_SL_centN = new TProfile("h_Aq22_SL_centN","",100,0,100,tpmin,tpmax,"");
 712   fOutputList->Add(h_Aq22_SL_cent);
 713   fOutputList->Add(h_Aq22_SL_centP);
 714   fOutputList->Add(h_Aq22_SL_centN);
 715   h_Aq22_SLL_cent = new TProfile("h_Aq22_SLL_cent","",100,0,100,tpmin,tpmax,"");
 716   h_Aq22_SLL_centP = new TProfile("h_Aq22_SLL_centP","",100,0,100,tpmin,tpmax,"");
 717   h_Aq22_SLL_centN = new TProfile("h_Aq22_SLL_centN","",100,0,100,tpmin,tpmax,"");
 718   fOutputList->Add(h_Aq22_SLL_cent);
 719   fOutputList->Add(h_Aq22_SLL_centP);
 720   fOutputList->Add(h_Aq22_SLL_centN);
 721   h_Aq22_SLR_cent = new TProfile("h_Aq22_SLR_cent","",100,0,100,tpmin,tpmax,"");
 722   h_Aq22_SLR_centP = new TProfile("h_Aq22_SLR_centP","",100,0,100,tpmin,tpmax,"");
 723   h_Aq22_SLR_centN = new TProfile("h_Aq22_SLR_centN","",100,0,100,tpmin,tpmax,"");
 724   fOutputList->Add(h_Aq22_SLR_cent);
 725   fOutputList->Add(h_Aq22_SLR_centP);
 726   fOutputList->Add(h_Aq22_SLR_centN);
 727
 728   // --- third harmonic
 729   h_Aq32_cent = new TProfile("h_Aq32_cent","",100,0,100,tpmin,tpmax,"");
 730   h_Aq32_centP = new TProfile("h_Aq32_centP","",100,0,100,tpmin,tpmax,"");
 731   h_Aq32_centN = new TProfile("h_Aq32_centN","",100,0,100,tpmin,tpmax,"");
 732   fOutputList->Add(h_Aq32_cent);
 733   fOutputList->Add(h_Aq32_centP);
 734   fOutputList->Add(h_Aq32_centN);
 735   h_Aq32_SL_cent = new TProfile("h_Aq32_SL_cent","",100,0,100,tpmin,tpmax,"");
 736   h_Aq32_SL_centP = new TProfile("h_Aq32_SL_centP","",100,0,100,tpmin,tpmax,"");
 737   h_Aq32_SL_centN = new TProfile("h_Aq32_SL_centN","",100,0,100,tpmin,tpmax,"");
 738   fOutputList->Add(h_Aq32_SL_cent);
 739   fOutputList->Add(h_Aq32_SL_centP);
 740   fOutputList->Add(h_Aq32_SL_centN);
 741   h_Aq32_SLL_cent = new TProfile("h_Aq32_SLL_cent","",100,0,100,tpmin,tpmax,"");
 742   h_Aq32_SLL_centP = new TProfile("h_Aq32_SLL_centP","",100,0,100,tpmin,tpmax,"");
 743   h_Aq32_SLL_centN = new TProfile("h_Aq32_SLL_centN","",100,0,100,tpmin,tpmax,"");
 744   fOutputList->Add(h_Aq32_SLL_cent);
 745   fOutputList->Add(h_Aq32_SLL_centP);
 746   fOutputList->Add(h_Aq32_SLL_centN);
 747   h_Aq32_SLR_cent = new TProfile("h_Aq32_SLR_cent","",100,0,100,tpmin,tpmax,"");
 748   h_Aq32_SLR_centP = new TProfile("h_Aq32_SLR_centP","",100,0,100,tpmin,tpmax,"");
 749   h_Aq32_SLR_centN = new TProfile("h_Aq32_SLR_centN","",100,0,100,tpmin,tpmax,"");
 750   fOutputList->Add(h_Aq32_SLR_cent);
 751   fOutputList->Add(h_Aq32_SLR_centP);
 752   fOutputList->Add(h_Aq32_SLR_centN);
 753
 754   // --- fourth harmonic
 755   h_Aq42_cent = new TProfile("h_Aq42_cent","",100,0,100,tpmin,tpmax,"");
 756   h_Aq42_centP = new TProfile("h_Aq42_centP","",100,0,100,tpmin,tpmax,"");
 757   h_Aq42_centN = new TProfile("h_Aq42_centN","",100,0,100,tpmin,tpmax,"");
 758   fOutputList->Add(h_Aq42_cent);
 759   fOutputList->Add(h_Aq42_centP);
 760   fOutputList->Add(h_Aq42_centN);
 761   h_Aq42_SL_cent = new TProfile("h_Aq42_SL_cent","",100,0,100,tpmin,tpmax,"");
 762   h_Aq42_SL_centP = new TProfile("h_Aq42_SL_centP","",100,0,100,tpmin,tpmax,"");
 763   h_Aq42_SL_centN = new TProfile("h_Aq42_SL_centN","",100,0,100,tpmin,tpmax,"");
 764   fOutputList->Add(h_Aq42_SL_cent);
 765   fOutputList->Add(h_Aq42_SL_centP);
 766   fOutputList->Add(h_Aq42_SL_centN);
 767   h_Aq42_SLL_cent = new TProfile("h_Aq42_SLL_cent","",100,0,100,tpmin,tpmax,"");
 768   h_Aq42_SLL_centP = new TProfile("h_Aq42_SLL_centP","",100,0,100,tpmin,tpmax,"");
 769   h_Aq42_SLL_centN = new TProfile("h_Aq42_SLL_centN","",100,0,100,tpmin,tpmax,"");
 770   fOutputList->Add(h_Aq42_SLL_cent);
 771   fOutputList->Add(h_Aq42_SLL_centP);
 772   fOutputList->Add(h_Aq42_SLL_centN);
 773   h_Aq42_SLR_cent = new TProfile("h_Aq42_SLR_cent","",100,0,100,tpmin,tpmax,"");
 774   h_Aq42_SLR_centP = new TProfile("h_Aq42_SLR_centP","",100,0,100,tpmin,tpmax,"");
 775   h_Aq42_SLR_centN = new TProfile("h_Aq42_SLR_centN","",100,0,100,tpmin,tpmax,"");
 776   fOutputList->Add(h_Aq42_SLR_cent);
 777   fOutputList->Add(h_Aq42_SLR_centP);
 778   fOutputList->Add(h_Aq42_SLR_centN);
 779
 780
 781
 782
 783   // ------------------------------------
 784   // --- differential cumulant histograms
 785   // ------------------------------------
 786
 787   h_diffq22_pt = new TProfile("h_diffq22_pt","",50,0,5,tpmin,tpmax,"");
 788   h_diffq22_ptP = new TProfile("h_diffq22_ptP","",50,0,5,tpmin,tpmax,"");
 789   h_diffq22_ptN = new TProfile("h_diffq22_ptN","",50,0,5,tpmin,tpmax,"");
 790   h_diffAq22_pt = new TProfile("h_diffAq22_pt","",50,0,5,tpmin,tpmax,"");
 791   h_diffAq22_ptP = new TProfile("h_diffAq22_ptP","",50,0,5,tpmin,tpmax,"");
 792   h_diffAq22_ptN = new TProfile("h_diffAq22_ptN","",50,0,5,tpmin,tpmax,"");
 793   fOutputList->Add(h_diffq22_pt);
 794   fOutputList->Add(h_diffq22_ptP);
 795   fOutputList->Add(h_diffq22_ptN);
 796   fOutputList->Add(h_diffAq22_pt);
 797   fOutputList->Add(h_diffAq22_ptP);
 798   fOutputList->Add(h_diffAq22_ptN);
 799   h_diffq22_eta = new TProfile("h_diffq22_eta","",32,-0.8,0.8,tpmin,tpmax,"");
 800   h_diffq22_etaP = new TProfile("h_diffq22_etaP","",32,-0.8,0.8,tpmin,tpmax,"");
 801   h_diffq22_etaN = new TProfile("h_diffq22_etaN","",32,-0.8,0.8,tpmin,tpmax,"");
 802   h_diffAq22_eta = new TProfile("h_diffAq22_eta","",32,-0.8,0.8,tpmin,tpmax,"");
 803   h_diffAq22_etaP = new TProfile("h_diffAq22_etaP","",32,-0.8,0.8,tpmin,tpmax,"");
 804   h_diffAq22_etaN = new TProfile("h_diffAq22_etaN","",32,-0.8,0.8,tpmin,tpmax,"");
 805   fOutputList->Add(h_diffq22_eta);
 806   fOutputList->Add(h_diffq22_etaP);
 807   fOutputList->Add(h_diffq22_etaN);
 808   fOutputList->Add(h_diffAq22_eta);
 809   fOutputList->Add(h_diffAq22_etaP);
 810   fOutputList->Add(h_diffAq22_etaN);
 811
 812   h_diffq32_pt = new TProfile("h_diffq32_pt","",50,0,5,tpmin,tpmax,"");
 813   h_diffq32_ptP = new TProfile("h_diffq32_ptP","",50,0,5,tpmin,tpmax,"");
 814   h_diffq32_ptN = new TProfile("h_diffq32_ptN","",50,0,5,tpmin,tpmax,"");
 815   h_diffAq32_pt = new TProfile("h_diffAq32_pt","",50,0,5,tpmin,tpmax,"");
 816   h_diffAq32_ptP = new TProfile("h_diffAq32_ptP","",50,0,5,tpmin,tpmax,"");
 817   h_diffAq32_ptN = new TProfile("h_diffAq32_ptN","",50,0,5,tpmin,tpmax,"");
 818   fOutputList->Add(h_diffq32_pt);
 819   fOutputList->Add(h_diffq32_ptP);
 820   fOutputList->Add(h_diffq32_ptN);
 821   fOutputList->Add(h_diffAq32_pt);
 822   fOutputList->Add(h_diffAq32_ptP);
 823   fOutputList->Add(h_diffAq32_ptN);
 824   h_diffq32_eta = new TProfile("h_diffq32_eta","",32,-0.8,0.8,tpmin,tpmax,"");
 825   h_diffq32_etaP = new TProfile("h_diffq32_etaP","",32,-0.8,0.8,tpmin,tpmax,"");
 826   h_diffq32_etaN = new TProfile("h_diffq32_etaN","",32,-0.8,0.8,tpmin,tpmax,"");
 827   h_diffAq32_eta = new TProfile("h_diffAq32_eta","",32,-0.8,0.8,tpmin,tpmax,"");
 828   h_diffAq32_etaP = new TProfile("h_diffAq32_etaP","",32,-0.8,0.8,tpmin,tpmax,"");
 829   h_diffAq32_etaN = new TProfile("h_diffAq32_etaN","",32,-0.8,0.8,tpmin,tpmax,"");
 830   fOutputList->Add(h_diffq32_eta);
 831   fOutputList->Add(h_diffq32_etaP);
 832   fOutputList->Add(h_diffq32_etaN);
 833   fOutputList->Add(h_diffAq32_eta);
 834   fOutputList->Add(h_diffAq32_etaP);
 835   fOutputList->Add(h_diffAq32_etaN);
 836
 837   h_diffq42_pt = new TProfile("h_diffq42_pt","",50,0,5,tpmin,tpmax,"");
 838   h_diffq42_ptP = new TProfile("h_diffq42_ptP","",50,0,5,tpmin,tpmax,"");
 839   h_diffq42_ptN = new TProfile("h_diffq42_ptN","",50,0,5,tpmin,tpmax,"");
 840   h_diffAq42_pt = new TProfile("h_diffAq42_pt","",50,0,5,tpmin,tpmax,"");
 841   h_diffAq42_ptP = new TProfile("h_diffAq42_ptP","",50,0,5,tpmin,tpmax,"");
 842   h_diffAq42_ptN = new TProfile("h_diffAq42_ptN","",50,0,5,tpmin,tpmax,"");
 843   fOutputList->Add(h_diffq42_pt);
 844   fOutputList->Add(h_diffq42_ptP);
 845   fOutputList->Add(h_diffq42_ptN);
 846   fOutputList->Add(h_diffAq42_pt);
 847   fOutputList->Add(h_diffAq42_ptP);
 848   fOutputList->Add(h_diffAq42_ptN);
 849   h_diffq42_eta = new TProfile("h_diffq42_eta","",32,-0.8,0.8,tpmin,tpmax,"");
 850   h_diffq42_etaP = new TProfile("h_diffq42_etaP","",32,-0.8,0.8,tpmin,tpmax,"");
 851   h_diffq42_etaN = new TProfile("h_diffq42_etaN","",32,-0.8,0.8,tpmin,tpmax,"");
 852   h_diffAq42_eta = new TProfile("h_diffAq42_eta","",32,-0.8,0.8,tpmin,tpmax,"");
 853   h_diffAq42_etaP = new TProfile("h_diffAq42_etaP","",32,-0.8,0.8,tpmin,tpmax,"");
 854   h_diffAq42_etaN = new TProfile("h_diffAq42_etaN","",32,-0.8,0.8,tpmin,tpmax,"");
 855   fOutputList->Add(h_diffq42_eta);
 856   fOutputList->Add(h_diffq42_etaP);
 857   fOutputList->Add(h_diffq42_etaN);
 858   fOutputList->Add(h_diffAq42_eta);
 859   fOutputList->Add(h_diffAq42_etaP);
 860   fOutputList->Add(h_diffAq42_etaN);
 861
 862
 863
 864   h_AT_X_deta = new TProfile("h_AT_X_deta","",32,-1.6,1.6,tpmin,tpmax,"");
 865   h_AT_X_detaP = new TProfile("h_AT_X_detaP","",32,-1.6,1.6,tpmin,tpmax,"");
 866   h_AT_X_detaN = new TProfile("h_AT_X_detaN","",32,-1.6,1.6,tpmin,tpmax,"");
 867   fOutputList->Add(h_AT_X_deta);
 868   fOutputList->Add(h_AT_X_detaP);
 869   fOutputList->Add(h_AT_X_detaN);
 870   //h_diffq22_X_deta = new TProfile("h_diffq22_X_deta","",32,-1.6,1.6,tpmin,tpmax,"");
 871   h_diffq22_X_detaP = new TProfile("h_diffq22_X_detaP","",32,-1.6,1.6,tpmin,tpmax,"");
 872   h_diffq22_X_detaN = new TProfile("h_diffq22_X_detaN","",32,-1.6,1.6,tpmin,tpmax,"");
 873   //h_ATdiffq22_X_deta = new TProfile("h_ATdiffq22_X_deta","",32,-1.6,1.6,tpmin,tpmax,"");
 874   h_ATdiffq22_X_detaP = new TProfile("h_ATdiffq22_X_detaP","",32,-1.6,1.6,tpmin,tpmax,"");
 875   h_ATdiffq22_X_detaN = new TProfile("h_ATdiffq22_X_detaN","",32,-1.6,1.6,tpmin,tpmax,"");
 876   fOutputList->Add(h_diffq22_X_detaP);
 877   fOutputList->Add(h_diffq22_X_detaN);
 878   fOutputList->Add(h_ATdiffq22_X_detaP);
 879   fOutputList->Add(h_ATdiffq22_X_detaN);
 880   //h_diffq32_X_deta = new TProfile("h_diffq32_X_deta","",32,-1.6,1.6,tpmin,tpmax,"");
 881   h_diffq32_X_detaP = new TProfile("h_diffq32_X_detaP","",32,-1.6,1.6,tpmin,tpmax,"");
 882   h_diffq32_X_detaN = new TProfile("h_diffq32_X_detaN","",32,-1.6,1.6,tpmin,tpmax,"");
 883   //h_ATdiffq32_X_deta = new TProfile("h_ATdiffq32_X_deta","",32,-1.6,1.6,tpmin,tpmax,"");
 884   h_ATdiffq32_X_detaP = new TProfile("h_ATdiffq32_X_detaP","",32,-1.6,1.6,tpmin,tpmax,"");
 885   h_ATdiffq32_X_detaN = new TProfile("h_ATdiffq32_X_detaN","",32,-1.6,1.6,tpmin,tpmax,"");
 886   fOutputList->Add(h_diffq32_X_detaP);
 887   fOutputList->Add(h_diffq32_X_detaN);
 888   fOutputList->Add(h_ATdiffq32_X_detaP);
 889   fOutputList->Add(h_ATdiffq32_X_detaN);
 890   //h_diffq42_X_deta = new TProfile("h_diffq42_X_deta","",32,-1.6,1.6,tpmin,tpmax,"");
 891   h_diffq42_X_detaP = new TProfile("h_diffq42_X_detaP","",32,-1.6,1.6,tpmin,tpmax,"");
 892   h_diffq42_X_detaN = new TProfile("h_diffq42_X_detaN","",32,-1.6,1.6,tpmin,tpmax,"");
 893   //h_ATdiffq42_X_deta = new TProfile("h_ATdiffq42_X_deta","",32,-1.6,1.6,tpmin,tpmax,"");
 894   h_ATdiffq42_X_detaP = new TProfile("h_ATdiffq42_X_detaP","",32,-1.6,1.6,tpmin,tpmax,"");
 895   h_ATdiffq42_X_detaN = new TProfile("h_ATdiffq42_X_detaN","",32,-1.6,1.6,tpmin,tpmax,"");
 896   fOutputList->Add(h_diffq42_X_detaP);
 897   fOutputList->Add(h_diffq42_X_detaN);
 898   fOutputList->Add(h_ATdiffq42_X_detaP);
 899   fOutputList->Add(h_ATdiffq42_X_detaN);
 900
 901
 902   // --- now field selection
 903   // --- field 1
 904   h_AT_X_deta_F1 = new TProfile("h_AT_X_deta_F1","",32,-1.6,1.6,tpmin,tpmax,"");
 905   h_AT_X_detaP_F1 = new TProfile("h_AT_X_detaP_F1","",32,-1.6,1.6,tpmin,tpmax,"");
 906   h_AT_X_detaN_F1 = new TProfile("h_AT_X_detaN_F1","",32,-1.6,1.6,tpmin,tpmax,"");
 907   fOutputList->Add(h_AT_X_deta_F1);
 908   fOutputList->Add(h_AT_X_detaP_F1);
 909   fOutputList->Add(h_AT_X_detaN_F1);
 910   //h_diffq22_X_deta_F1 = new TProfile("h_diffq22_X_deta_F1","",32,-1.6,1.6,tpmin,tpmax,"");
 911   h_diffq22_X_detaP_F1 = new TProfile("h_diffq22_X_detaP_F1","",32,-1.6,1.6,tpmin,tpmax,"");
 912   h_diffq22_X_detaN_F1 = new TProfile("h_diffq22_X_detaN_F1","",32,-1.6,1.6,tpmin,tpmax,"");
 913   //h_ATdiffq22_X_deta_F1 = new TProfile("h_ATdiffq22_X_deta_F1","",32,-1.6,1.6,tpmin,tpmax,"");
 914   h_ATdiffq22_X_detaP_F1 = new TProfile("h_ATdiffq22_X_detaP_F1","",32,-1.6,1.6,tpmin,tpmax,"");
 915   h_ATdiffq22_X_detaN_F1 = new TProfile("h_ATdiffq22_X_detaN_F1","",32,-1.6,1.6,tpmin,tpmax,"");
 916   fOutputList->Add(h_diffq22_X_detaP_F1);
 917   fOutputList->Add(h_diffq22_X_detaN_F1);
 918   fOutputList->Add(h_ATdiffq22_X_detaP_F1);
 919   fOutputList->Add(h_ATdiffq22_X_detaN_F1);
 920   //h_diffq32_X_deta_F1 = new TProfile("h_diffq32_X_deta_F1","",32,-1.6,1.6,tpmin,tpmax,"");
 921   h_diffq32_X_detaP_F1 = new TProfile("h_diffq32_X_detaP_F1","",32,-1.6,1.6,tpmin,tpmax,"");
 922   h_diffq32_X_detaN_F1 = new TProfile("h_diffq32_X_detaN_F1","",32,-1.6,1.6,tpmin,tpmax,"");
 923   //h_ATdiffq32_X_deta_F1 = new TProfile("h_ATdiffq32_X_deta_F1","",32,-1.6,1.6,tpmin,tpmax,"");
 924   h_ATdiffq32_X_detaP_F1 = new TProfile("h_ATdiffq32_X_detaP_F1","",32,-1.6,1.6,tpmin,tpmax,"");
 925   h_ATdiffq32_X_detaN_F1 = new TProfile("h_ATdiffq32_X_detaN_F1","",32,-1.6,1.6,tpmin,tpmax,"");
 926   fOutputList->Add(h_diffq32_X_detaP_F1);
 927   fOutputList->Add(h_diffq32_X_detaN_F1);
 928   fOutputList->Add(h_ATdiffq32_X_detaP_F1);
 929   fOutputList->Add(h_ATdiffq32_X_detaN_F1);
 930   //h_diffq42_X_deta_F1 = new TProfile("h_diffq42_X_deta_F1","",32,-1.6,1.6,tpmin,tpmax,"");
 931   h_diffq42_X_detaP_F1 = new TProfile("h_diffq42_X_detaP_F1","",32,-1.6,1.6,tpmin,tpmax,"");
 932   h_diffq42_X_detaN_F1 = new TProfile("h_diffq42_X_detaN_F1","",32,-1.6,1.6,tpmin,tpmax,"");
 933   //h_ATdiffq42_X_deta_F1 = new TProfile("h_ATdiffq42_X_deta_F1","",32,-1.6,1.6,tpmin,tpmax,"");
 934   h_ATdiffq42_X_detaP_F1 = new TProfile("h_ATdiffq42_X_detaP_F1","",32,-1.6,1.6,tpmin,tpmax,"");
 935   h_ATdiffq42_X_detaN_F1 = new TProfile("h_ATdiffq42_X_detaN_F1","",32,-1.6,1.6,tpmin,tpmax,"");
 936   fOutputList->Add(h_diffq42_X_detaP_F1);
 937   fOutputList->Add(h_diffq42_X_detaN_F1);
 938   fOutputList->Add(h_ATdiffq42_X_detaP_F1);
 939   fOutputList->Add(h_ATdiffq42_X_detaN_F1);
 940   // --- field 3
 941   h_AT_X_deta_F3 = new TProfile("h_AT_X_deta_F3","",32,-1.6,1.6,tpmin,tpmax,"");
 942   h_AT_X_detaP_F3 = new TProfile("h_AT_X_detaP_F3","",32,-1.6,1.6,tpmin,tpmax,"");
 943   h_AT_X_detaN_F3 = new TProfile("h_AT_X_detaN_F3","",32,-1.6,1.6,tpmin,tpmax,"");
 944   fOutputList->Add(h_AT_X_deta_F3);
 945   fOutputList->Add(h_AT_X_detaP_F3);
 946   fOutputList->Add(h_AT_X_detaN_F3);
 947   //h_diffq22_X_deta_F3 = new TProfile("h_diffq22_X_deta_F3","",32,-1.6,1.6,tpmin,tpmax,"");
 948   h_diffq22_X_detaP_F3 = new TProfile("h_diffq22_X_detaP_F3","",32,-1.6,1.6,tpmin,tpmax,"");
 949   h_diffq22_X_detaN_F3 = new TProfile("h_diffq22_X_detaN_F3","",32,-1.6,1.6,tpmin,tpmax,"");
 950   //h_ATdiffq22_X_deta_F3 = new TProfile("h_ATdiffq22_X_deta_F3","",32,-1.6,1.6,tpmin,tpmax,"");
 951   h_ATdiffq22_X_detaP_F3 = new TProfile("h_ATdiffq22_X_detaP_F3","",32,-1.6,1.6,tpmin,tpmax,"");
 952   h_ATdiffq22_X_detaN_F3 = new TProfile("h_ATdiffq22_X_detaN_F3","",32,-1.6,1.6,tpmin,tpmax,"");
 953   fOutputList->Add(h_diffq22_X_detaP_F3);
 954   fOutputList->Add(h_diffq22_X_detaN_F3);
 955   fOutputList->Add(h_ATdiffq22_X_detaP_F3);
 956   fOutputList->Add(h_ATdiffq22_X_detaN_F3);
 957   //h_diffq32_X_deta_F3 = new TProfile("h_diffq32_X_deta_F3","",32,-1.6,1.6,tpmin,tpmax,"");
 958   h_diffq32_X_detaP_F3 = new TProfile("h_diffq32_X_detaP_F3","",32,-1.6,1.6,tpmin,tpmax,"");
 959   h_diffq32_X_detaN_F3 = new TProfile("h_diffq32_X_detaN_F3","",32,-1.6,1.6,tpmin,tpmax,"");
 960   //h_ATdiffq32_X_deta_F3 = new TProfile("h_ATdiffq32_X_deta_F3","",32,-1.6,1.6,tpmin,tpmax,"");
 961   h_ATdiffq32_X_detaP_F3 = new TProfile("h_ATdiffq32_X_detaP_F3","",32,-1.6,1.6,tpmin,tpmax,"");
 962   h_ATdiffq32_X_detaN_F3 = new TProfile("h_ATdiffq32_X_detaN_F3","",32,-1.6,1.6,tpmin,tpmax,"");
 963   fOutputList->Add(h_diffq32_X_detaP_F3);
 964   fOutputList->Add(h_diffq32_X_detaN_F3);
 965   fOutputList->Add(h_ATdiffq32_X_detaP_F3);
 966   fOutputList->Add(h_ATdiffq32_X_detaN_F3);
 967   //h_diffq42_X_deta_F3 = new TProfile("h_diffq42_X_deta_F3","",32,-1.6,1.6,tpmin,tpmax,"");
 968   h_diffq42_X_detaP_F3 = new TProfile("h_diffq42_X_detaP_F3","",32,-1.6,1.6,tpmin,tpmax,"");
 969   h_diffq42_X_detaN_F3 = new TProfile("h_diffq42_X_detaN_F3","",32,-1.6,1.6,tpmin,tpmax,"");
 970   //h_ATdiffq42_X_deta_F3 = new TProfile("h_ATdiffq42_X_deta_F3","",32,-1.6,1.6,tpmin,tpmax,"");
 971   h_ATdiffq42_X_detaP_F3 = new TProfile("h_ATdiffq42_X_detaP_F3","",32,-1.6,1.6,tpmin,tpmax,"");
 972   h_ATdiffq42_X_detaN_F3 = new TProfile("h_ATdiffq42_X_detaN_F3","",32,-1.6,1.6,tpmin,tpmax,"");
 973   fOutputList->Add(h_diffq42_X_detaP_F3);
 974   fOutputList->Add(h_diffq42_X_detaN_F3);
 975   fOutputList->Add(h_ATdiffq42_X_detaP_F3);
 976   fOutputList->Add(h_ATdiffq42_X_detaN_F3);
 977
 978
 979
 980   h_AT_X_dpt = new TProfile("h_AT_X_dpt","",32,-1.6,1.6,tpmin,tpmax,"");
 981   h_AT_X_dptP = new TProfile("h_AT_X_dptP","",32,-1.6,1.6,tpmin,tpmax,"");
 982   h_AT_X_dptN = new TProfile("h_AT_X_dptN","",32,-1.6,1.6,tpmin,tpmax,"");
 983   fOutputList->Add(h_AT_X_dpt);
 984   fOutputList->Add(h_AT_X_dptP);
 985   fOutputList->Add(h_AT_X_dptN);
 986   //h_diffq22_X_dpt = new TProfile("h_diffq22_X_dpt","",32,-1.6,1.6,tpmin,tpmax,"");
 987   h_diffq22_X_dptP = new TProfile("h_diffq22_X_dptP","",32,-1.6,1.6,tpmin,tpmax,"");
 988   h_diffq22_X_dptN = new TProfile("h_diffq22_X_dptN","",32,-1.6,1.6,tpmin,tpmax,"");
 989   //h_ATdiffq22_X_dpt = new TProfile("h_ATdiffq22_X_dpt","",32,-1.6,1.6,tpmin,tpmax,"");
 990   h_ATdiffq22_X_dptP = new TProfile("h_ATdiffq22_X_dptP","",32,-1.6,1.6,tpmin,tpmax,"");
 991   h_ATdiffq22_X_dptN = new TProfile("h_ATdiffq22_X_dptN","",32,-1.6,1.6,tpmin,tpmax,"");
 992   fOutputList->Add(h_diffq22_X_dptP);
 993   fOutputList->Add(h_diffq22_X_dptN);
 994   fOutputList->Add(h_ATdiffq22_X_dptP);
 995   fOutputList->Add(h_ATdiffq22_X_dptN);
 996   //h_diffq32_X_dpt = new TProfile("h_diffq32_X_dpt","",32,-1.6,1.6,tpmin,tpmax,"");
 997   h_diffq32_X_dptP = new TProfile("h_diffq32_X_dptP","",32,-1.6,1.6,tpmin,tpmax,"");
 998   h_diffq32_X_dptN = new TProfile("h_diffq32_X_dptN","",32,-1.6,1.6,tpmin,tpmax,"");
 999   //h_ATdiffq32_X_dpt = new TProfile("h_ATdiffq32_X_dpt","",32,-1.6,1.6,tpmin,tpmax,"");
1000   h_ATdiffq32_X_dptP = new TProfile("h_ATdiffq32_X_dptP","",32,-1.6,1.6,tpmin,tpmax,"");
1001   h_ATdiffq32_X_dptN = new TProfile("h_ATdiffq32_X_dptN","",32,-1.6,1.6,tpmin,tpmax,"");
1002   fOutputList->Add(h_diffq32_X_dptP);
1003   fOutputList->Add(h_diffq32_X_dptN);
1004   fOutputList->Add(h_ATdiffq32_X_dptP);
1005   fOutputList->Add(h_ATdiffq32_X_dptN);
1006   //h_diffq42_X_dpt = new TProfile("h_diffq42_X_dpt","",32,-1.6,1.6,tpmin,tpmax,"");
1007   h_diffq42_X_dptP = new TProfile("h_diffq42_X_dptP","",32,-1.6,1.6,tpmin,tpmax,"");
1008   h_diffq42_X_dptN = new TProfile("h_diffq42_X_dptN","",32,-1.6,1.6,tpmin,tpmax,"");
1009   //h_ATdiffq42_X_dpt = new TProfile("h_ATdiffq42_X_dpt","",32,-1.6,1.6,tpmin,tpmax,"");
1010   h_ATdiffq42_X_dptP = new TProfile("h_ATdiffq42_X_dptP","",32,-1.6,1.6,tpmin,tpmax,"");
1011   h_ATdiffq42_X_dptN = new TProfile("h_ATdiffq42_X_dptN","",32,-1.6,1.6,tpmin,tpmax,"");
1012   fOutputList->Add(h_diffq42_X_dptP);
1013   fOutputList->Add(h_diffq42_X_dptN);
1014   fOutputList->Add(h_ATdiffq42_X_dptP);
1015   fOutputList->Add(h_ATdiffq42_X_dptN);
1016
1017
1018
1019   // --- now Xsub...
1020
1021   h_AT_S_deta = new TProfile("h_AT_S_deta","",32,-1.6,1.6,tpmin,tpmax,"");
1022   h_AT_S_detaP = new TProfile("h_AT_S_detaP","",32,-1.6,1.6,tpmin,tpmax,"");
1023   h_AT_S_detaN = new TProfile("h_AT_S_detaN","",32,-1.6,1.6,tpmin,tpmax,"");
1024   fOutputList->Add(h_AT_S_deta);
1025   fOutputList->Add(h_AT_S_detaP);
1026   fOutputList->Add(h_AT_S_detaN);
1027   //h_diffq22_S_deta = new TProfile("h_diffq22_S_deta","",32,-1.6,1.6,tpmin,tpmax,"");
1028   h_diffq22_S_detaP = new TProfile("h_diffq22_S_detaP","",32,-1.6,1.6,tpmin,tpmax,"");
1029   h_diffq22_S_detaN = new TProfile("h_diffq22_S_detaN","",32,-1.6,1.6,tpmin,tpmax,"");
1030   //h_ATdiffq22_S_deta = new TProfile("h_ATdiffq22_S_deta","",32,-1.6,1.6,tpmin,tpmax,"");
1031   h_ATdiffq22_S_detaP = new TProfile("h_ATdiffq22_S_detaP","",32,-1.6,1.6,tpmin,tpmax,"");
1032   h_ATdiffq22_S_detaN = new TProfile("h_ATdiffq22_S_detaN","",32,-1.6,1.6,tpmin,tpmax,"");
1033   fOutputList->Add(h_diffq22_S_detaP);
1034   fOutputList->Add(h_diffq22_S_detaN);
1035   fOutputList->Add(h_ATdiffq22_S_detaP);
1036   fOutputList->Add(h_ATdiffq22_S_detaN);
1037   //h_diffq32_S_deta = new TProfile("h_diffq32_S_deta","",32,-1.6,1.6,tpmin,tpmax,"");
1038   h_diffq32_S_detaP = new TProfile("h_diffq32_S_detaP","",32,-1.6,1.6,tpmin,tpmax,"");
1039   h_diffq32_S_detaN = new TProfile("h_diffq32_S_detaN","",32,-1.6,1.6,tpmin,tpmax,"");
1040   //h_ATdiffq32_S_deta = new TProfile("h_ATdiffq32_S_deta","",32,-1.6,1.6,tpmin,tpmax,"");
1041   h_ATdiffq32_S_detaP = new TProfile("h_ATdiffq32_S_detaP","",32,-1.6,1.6,tpmin,tpmax,"");
1042   h_ATdiffq32_S_detaN = new TProfile("h_ATdiffq32_S_detaN","",32,-1.6,1.6,tpmin,tpmax,"");
1043   fOutputList->Add(h_diffq32_S_detaP);
1044   fOutputList->Add(h_diffq32_S_detaN);
1045   fOutputList->Add(h_ATdiffq32_S_detaP);
1046   fOutputList->Add(h_ATdiffq32_S_detaN);
1047   //h_diffq42_S_deta = new TProfile("h_diffq42_S_deta","",32,-1.6,1.6,tpmin,tpmax,"");
1048   h_diffq42_S_detaP = new TProfile("h_diffq42_S_detaP","",32,-1.6,1.6,tpmin,tpmax,"");
1049   h_diffq42_S_detaN = new TProfile("h_diffq42_S_detaN","",32,-1.6,1.6,tpmin,tpmax,"");
1050   //h_ATdiffq42_S_deta = new TProfile("h_ATdiffq42_S_deta","",32,-1.6,1.6,tpmin,tpmax,"");
1051   h_ATdiffq42_S_detaP = new TProfile("h_ATdiffq42_S_detaP","",32,-1.6,1.6,tpmin,tpmax,"");
1052   h_ATdiffq42_S_detaN = new TProfile("h_ATdiffq42_S_detaN","",32,-1.6,1.6,tpmin,tpmax,"");
1053   fOutputList->Add(h_diffq42_S_detaP);
1054   fOutputList->Add(h_diffq42_S_detaN);
1055   fOutputList->Add(h_ATdiffq42_S_detaP);
1056   fOutputList->Add(h_ATdiffq42_S_detaN);
1057
1058   // ---
1059
1060   // --- random fun
1061
1062   // ---
1063
1064   h_a1q22_cent = new TProfile("h_a1q22_cent","q22 vs centrality",100,0,100,tpmin,tpmax,"");
1065   h_a1q22_centP = new TProfile("h_a1q22_centP","q22(pos) vs centrality",100,0,100,tpmin,tpmax,"");
1066   h_a1q22_centN = new TProfile("h_a1q22_centN","q22(neg) vs centrality",100,0,100,tpmin,tpmax,"");
1067   fOutputList->Add(h_a1q22_cent);
1068   fOutputList->Add(h_a1q22_centP);
1069   fOutputList->Add(h_a1q22_centN);
1070
1071   h_a2q22_cent = new TProfile("h_a2q22_cent","q22 vs centrality",100,0,100,tpmin,tpmax,"");
1072   h_a2q22_centP = new TProfile("h_a2q22_centP","q22(pos) vs centrality",100,0,100,tpmin,tpmax,"");
1073   h_a2q22_centN = new TProfile("h_a2q22_centN","q22(neg) vs centrality",100,0,100,tpmin,tpmax,"");
1074   fOutputList->Add(h_a2q22_cent);
1075   fOutputList->Add(h_a2q22_centP);
1076   fOutputList->Add(h_a2q22_centN);
1077
1078   h_rAq22_X1_cent = new TProfile("h_rAq22_X1_cent","q22 vs centrality",100,0,100,tpmin,tpmax,"");
1079   h_rAq22_X1_centP = new TProfile("h_rAq22_X1_centP","q22(pos) vs centrality",100,0,100,tpmin,tpmax,"");
1080   h_rAq22_X1_centN = new TProfile("h_rAq22_X1_centN","q22(neg) vs centrality",100,0,100,tpmin,tpmax,"");
1081   h_rAq22_X2_cent = new TProfile("h_rAq22_X2_cent","q22 vs centrality",100,0,100,tpmin,tpmax,"");
1082   h_rAq22_X2_centP = new TProfile("h_rAq22_X2_centP","q22(pos) vs centrality",100,0,100,tpmin,tpmax,"");
1083   h_rAq22_X2_centN = new TProfile("h_rAq22_X2_centN","q22(neg) vs centrality",100,0,100,tpmin,tpmax,"");
1084   h_rAq22_X3_cent = new TProfile("h_rAq22_X3_cent","q22 vs centrality",100,0,100,tpmin,tpmax,"");
1085   h_rAq22_X3_centP = new TProfile("h_rAq22_X3_centP","q22(pos) vs centrality",100,0,100,tpmin,tpmax,"");
1086   h_rAq22_X3_centN = new TProfile("h_rAq22_X3_centN","q22(neg) vs centrality",100,0,100,tpmin,tpmax,"");
1087   h_rAq22_X4_cent = new TProfile("h_rAq22_X4_cent","q22 vs centrality",100,0,100,tpmin,tpmax,"");
1088   h_rAq22_X4_centP = new TProfile("h_rAq22_X4_centP","q22(pos) vs centrality",100,0,100,tpmin,tpmax,"");
1089   h_rAq22_X4_centN = new TProfile("h_rAq22_X4_centN","q22(neg) vs centrality",100,0,100,tpmin,tpmax,"");
1090   fOutputList->Add(h_rAq22_X1_cent);
1091   fOutputList->Add(h_rAq22_X1_centP);
1092   fOutputList->Add(h_rAq22_X1_centN);
1093   fOutputList->Add(h_rAq22_X2_cent);
1094   fOutputList->Add(h_rAq22_X2_centP);
1095   fOutputList->Add(h_rAq22_X2_centN);
1096   fOutputList->Add(h_rAq22_X3_cent);
1097   fOutputList->Add(h_rAq22_X3_centP);
1098   fOutputList->Add(h_rAq22_X3_centN);
1099   fOutputList->Add(h_rAq22_X4_cent);
1100   fOutputList->Add(h_rAq22_X4_centP);
1101   fOutputList->Add(h_rAq22_X4_centN);
1102
1103
1104
1105   // ---------------------------- //
1106   // --- done with histograms --- //
1107   // ---------------------------- //
1108
1109
1110   // ------------------------------------ //
1111   // --- send data to the output list --- //
1112   // ------------------------------------ //
1113   PostData(1,fOutputList);
1114 }
1115
1116
1117
1118 //----------------------------------------------
1119 void AliAnalysisTaskCMEv2A::UserExec(Option_t *)
1120 {
1121
1122   // Main analyis loop, called for each event
1123
1124   if(debug>0) cout<<"Processing event with debug = "<<debug<<endl;
1125
1126   // Analysis manager
1127   AliAnalysisManager *manager = AliAnalysisManager::GetAnalysisManager();
1128   if(!manager)
1129     {
1130       cout<<"FATAL: could not get Analysis Manager."<<endl;
1131       return;
1132     }
1133
1134   // Input handler
1135   AliAODInputHandler *handler = (AliAODInputHandler *)manager->GetInputEventHandler();
1136   if(!handler)
1137     {
1138       cout<<"FATAL: could not get Input Handler."<<endl;
1139       return;
1140     }
1141
1142
1143   // AOD Event object from tree
1144   AliAODEvent *fAOD = (AliAODEvent *)InputEvent();
1145   if(!fAOD)
1146     {
1147       if(debug>-1) cout<<"ERROR: AOD event object not available. Discarding event..."<<endl;
1148       return;
1149     }
1150
1151
1152   int runnumber = fAOD->GetRunNumber();
1153   float mag = fAOD->GetMagneticField();
1154
1155   if(debug>0)
1156     {
1157       cout<<"runnumber is "<<runnumber<<endl;
1158       cout<<"magnetic field is "<<mag<<endl;
1159     }
1160
1161   // MC Event object from tree
1162   AliMCEvent *fMC = MCEvent();
1163   if(!fMC&&doMC)
1164     {
1165       if(debug>-1) cout<<"ERROR: MC event object not available. Discarding event..."<<endl;
1166       return;
1167     }
1168
1169   // get pid
1170   AliPIDResponse *fPID = handler->GetPIDResponse();
1171   if(!fPID && !doMC)    // use PID only if no MC
1172     {
1173       if(debug>-1) cout<<"ERROR: PIDResponse object not available. Discarding event..."<<endl;
1174       return;
1175     }
1176
1177   // // Eventplane object (from AOD)
1178   // AliEventplane *fEventplane = fAOD->GetEventplane();
1179   // if(!fEventplane)
1180   //   {
1181   //     if(debug>-1) cout<<"ERROR: Eventplane object not available. Discarding event..."<<endl;
1182   //     return;
1183   //   }
1184
1185   // float psi_V0_h2 = fEventplane->GetEventplane("V0",fAOD,2);
1186   // float psi_V0A_h2 = fEventplane->GetEventplane("V0A",fAOD,2);
1187   // float psi_V0C_h2 = fEventplane->GetEventplane("V0C",fAOD,2);
1188
1189   // fHistPlaneV0h2->Fill(psi_V0_h2);
1190   // fHistPlaneV0Ah2->Fill(psi_V0A_h2);
1191   // fHistPlaneV0Ch2->Fill(psi_V0C_h2);
1192   // fHistPlaneV0ACDCh2->Fill(psi_V0A_h2-psi_V0C_h2);
1193
1194
1195
1196   // Centrality object (from AOD)
1197   AliCentrality *fCentrality = fAOD->GetCentrality();
1198   if(!fCentrality)
1199     {
1200       if(debug>-1) cout<<"ERROR: Centrality object not available. Discarding event..."<<endl;
1201       return;
1202     }
1203
1204   float centTRK = fCentrality->GetCentralityPercentile("TRK");
1205   float centV0M = fCentrality->GetCentralityPercentile("V0M");
1206   float centSPD = fCentrality->GetCentralityPercentile("CL1");//outer SPD?
1207   float cent = centTRK;
1208   if(centhandle==2) cent = centV0M;
1209   if(centhandle==3) cent = centSPD;
1210
1211   int icent = int(cent)/10;
1212
1213   int centstatus = 0;
1214   if(centTRK<0.0||centV0M<0.0) centstatus = -1;
1215   if(centTRK<0.0&&centV0M<0.0) centstatus = -2;
1216   if(centTRK>0.0||centV0M>0.0) centstatus = 1;
1217   if(centTRK>0.0&&centV0M>0.0) centstatus = 2;
1218   if(centTRK==0.0&&centV0M==0.0) centstatus = 0;
1219
1220   if(debug>0)
1221     {
1222       cout<<"centTRK "<<centTRK<<endl;
1223       cout<<"centV0M "<<centV0M<<endl;
1224       cout<<"centSPD "<<centSPD<<endl;
1225       cout<<"centrality selection is "<<centhandle<<endl;
1226       cout<<"cent is "<<cent<<endl;
1227     }
1228   fHistCentTRK->Fill(centTRK);
1229   fHistCentV0M->Fill(centV0M);
1230   fHistCentDIFF->Fill(centTRK-centV0M);
1231
1232   //ULong64_t mask = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
1233   ULong64_t mask = handler->IsEventSelected();
1234   ULong64_t amb = AliVEvent::kMB;
1235   ULong64_t acn = AliVEvent::kCentral;
1236   ULong64_t asc = AliVEvent::kSemiCentral;
1237
1238   if(debug>0) cout<<"trigger selection is "<<trigger<<endl;
1239   if(debug>0) cout<<"trigger mask is "<<mask<<endl;
1240
1241   if(mask&amb)
1242     {
1243       fHistCentTRKAVEkMB->Fill(centTRK);
1244       fHistCentV0MAVEkMB->Fill(centV0M);
1245     }
1246   if(mask&acn)
1247     {
1248       fHistCentTRKAVEkCentral->Fill(centTRK);
1249       fHistCentV0MAVEkCentral->Fill(centV0M);
1250     }
1251   if(mask&asc)
1252     {
1253       fHistCentTRKAVEkSemiCentral->Fill(centTRK);
1254       fHistCentV0MAVEkSemiCentral->Fill(centV0M);
1255     }
1256   if(mask&(amb|acn|asc))
1257     {
1258       fHistCentTRKAVEkA3->Fill(centTRK);
1259       fHistCentV0MAVEkA3->Fill(centV0M);
1260     }
1261
1262   if(mask&trigger)
1263     {
1264       fHistCentTRKAVEkSel->Fill(centTRK);
1265       fHistCentV0MAVEkSel->Fill(centV0M);
1266     }
1267   else
1268     {
1269       if(debug>0) cout<<"wrong trigger, rejecting event"<<endl;
1270       return;
1271     }
1272
1273   if(fabs(centTRK-centV0M)>centcut)
1274     {
1275       if(debug>0) cout<<"centrality difference outside cut, rejecting event"<<endl;
1276       return;
1277     }
1278
1279
1280
1281   // AOD vertex objects
1282   AliAODVertex *fVtx = fAOD->GetPrimaryVertex();
1283   AliAODVertex *fVtxSPD = fAOD->GetPrimaryVertexSPD();
1284   if(!fVtx)
1285     {
1286       if(debug>-1) cout<<"ERROR: Vertex object not available. Discarding event..."<<endl;
1287       return;
1288     }
1289   float zvtxV0 = fVtx->GetZ();
1290   float zvtxSPD = fVtxSPD->GetZ();
1291   if(debug>0) cout<<"zvtxV0 is "<<zvtxV0<<endl;
1292   if(debug>0) cout<<"zvtxSPD is "<<zvtxSPD<<endl;
1293   if(centstatus==2)
1294     {
1295       fHistZVtx->Fill(zvtxV0);
1296       fHistZVtxD->Fill(zvtxV0-zvtxSPD);
1297     }
1298   if(fabs(zvtxV0)>zvtxcut)
1299     {
1300       if(debug>0) cout<<"vertex outside cut, rejecting event"<<endl;
1301       return;
1302     }
1303   float eventX = fVtx->GetX();
1304   float eventY = fVtx->GetY();
1305   float eventZ = fVtx->GetZ();
1306
1307
1308
1309   // // get V0 object
1310   // AliAODVZERO *fV0 = fAOD->GetVZEROData();
1311   // for(int i=0; i<64; i++)
1312   //   {
1313   //     float phiV0 = pi/4.0*(0.5+i%8);
1314   //     float multV0 = fV0->GetMultiplicity(i);
1315   //   }
1316
1317
1318   int d_ntrk = fAOD->GetNumberOfTracks();
1319   int d_ntrkMC = 0;
1320   if(fMC) d_ntrkMC = fMC->GetNumberOfTracks();
1321
1322   if(centstatus==-2&&d_ntrk>0) centstatus = -3;
1323   if(debug>0) cout<<"there are "<<d_ntrk<<" tracks in this event"<<endl;
1324   if(debug>0&&fMC) cout<<"there are "<<d_ntrkMC<<" Monte Carlo tracks in this event"<<endl;
1325   if(debug>0) cout<<"centrality diagnostic is "<<centstatus<<endl;
1326
1327   fHistCentDIAG->Fill(centstatus);
1328   if(centstatus==-3)
1329     {
1330       fHistCtrkDIAG->Fill(d_ntrk);
1331       fHistVtxRDIAG->Fill(zvtxV0);
1332       fHistVtxSDIAG->Fill(zvtxSPD);
1333     }
1334   if(centstatus==-2)
1335     {
1336       fHistVtxRDIBG->Fill(zvtxV0);
1337       fHistVtxSDIBG->Fill(zvtxSPD);
1338     }
1339
1340   // --- AliAnalysisUtils for pileup cuts
1341   if(dopupcut)
1342     {
1343       AliAnalysisUtils *fUtils = new AliAnalysisUtils();
1344       if(!fUtils)
1345         {
1346           if(debug>-1) cout<<"ERROR: cannot find AliAnalysisUtils..."<<endl;
1347           return;
1348         }
1349       fUtils->SetUseOutOfBunchPileUp(true);
1350       bool pileup = fUtils->IsPileUpEvent(fAOD);
1351       //bool pileup = fAOD->IsPileUpFromSPD(3,0.8,3.0,2.0,5.0); // default parameters in AliAODEvent.h
1352       if(pileup)
1353         {
1354           if(debug>0) cout<<"Rejecting event for pileup (AliAnalysisUtils)"<<endl;
1355           return;
1356         }
1357       if(fUtils) delete fUtils;
1358       //
1359       // ---
1360       //
1361       /*
1362       const AliAODVertex* vtxSPD = fAOD->GetPrimaryVertexSPD();
1363       if(!vtxSPD||vtxSPD->GetNContributors()<=0)
1364         {
1365           if(debug>0) cout<<"rejecting pileup event (no vtxSPD or zero contributors) "<<vtxSPD<<endl;
1366           return;
1367         }
1368
1369       const AliAODVertex* vtxTPC = 0;
1370       int nVertices = fAOD->GetNumberOfVertices();
1371       if(debug>0) cout<<"number of vertices is "<<nVertices<<endl;
1372       for(int iVertices = 0; iVertices < nVertices; iVertices++)
1373         {
1374           const AliAODVertex* vertex = fAOD->GetVertex(iVertices);
1375           if(debug>0) cout<<"vertex type is "<<vertex->GetType()<<endl;
1376           if(vertex->GetType()!=AliAODVertex::kMainTPC) continue;
1377           vtxTPC = vertex;
1378         }
1379       if(!vtxTPC||vtxTPC->GetNContributors()<=0)
1380         {
1381           if(debug>0) cout<<"rejecting pileup event (no vtxTPC or zero contributors)"<<vtxTPC<<endl;
1382           return;
1383         }
1384
1385       float diffZ = vtxSPD->GetZ() - vtxTPC->GetZ();
1386       if(fabs(diffZ)>2.0)
1387         {
1388           if(debug>0) cout<<"rejecting pileup event with vtxTPC "<<vtxTPC->GetZ()<<" vtxSPD "<<vtxSPD->GetZ()<<endl;
1389           return;
1390         }
1391       */
1392       //
1393       // ---
1394       //
1395     }
1396
1397
1398   int ntrk = 0;
1399   int ntrkpos = 0;
1400   int ntrkneg = 0;
1401   int ntrkL = 0;
1402   int ntrkposL = 0;
1403   int ntrknegL = 0;
1404   int ntrkR = 0;
1405   int ntrkposR = 0;
1406   int ntrknegR = 0;
1407   int ntrkA1 = 0;
1408   int ntrkposA1 = 0;
1409   int ntrknegA1 = 0;
1410   int ntrkA2 = 0;
1411   int ntrkposA2 = 0;
1412   int ntrknegA2 = 0;
1413   int ntrkMC = 0;
1414   int ntrkposMC = 0;
1415   int ntrknegMC = 0;
1416
1417   int cutntrk = 0;
1418   int cutntrkpos = 0;
1419   int cutntrkneg = 0;
1420   int cutntrkL = 0;
1421   int cutntrkposL = 0;
1422   int cutntrknegL = 0;
1423   int cutntrkR = 0;
1424   int cutntrkposR = 0;
1425   int cutntrknegR = 0;
1426   int cutntrkMC = 0;
1427   int cutntrkposMC = 0;
1428   int cutntrknegMC = 0;
1429
1430   const int hmax = 9; // number of harmonics including 0 for multiplicity counting
1431   float tpcXplo[hmax], tpcYplo[hmax], tpcXpro[hmax], tpcYpro[hmax];
1432   float tpcXnlo[hmax], tpcYnlo[hmax], tpcXnro[hmax], tpcYnro[hmax];
1433   float tpcXpli[hmax], tpcYpli[hmax], tpcXpri[hmax], tpcYpri[hmax];
1434   float tpcXnli[hmax], tpcYnli[hmax], tpcXnri[hmax], tpcYnri[hmax];
1435   float tpcXpl[hmax], tpcYpl[hmax], tpcXpr[hmax], tpcYpr[hmax];
1436   float tpcXnl[hmax], tpcYnl[hmax], tpcXnr[hmax], tpcYnr[hmax];
1437   for(int i=0; i<hmax;i++)
1438     {
1439       tpcXplo[i] = 0.0;
1440       tpcYplo[i] = 0.0;
1441       tpcXpro[i] = 0.0;
1442       tpcYpro[i] = 0.0;
1443       tpcXnlo[i] = 0.0;
1444       tpcYnlo[i] = 0.0;
1445       tpcXnro[i] = 0.0;
1446       tpcYnro[i] = 0.0;
1447       //
1448       tpcXpli[i] = 0.0;
1449       tpcYpli[i] = 0.0;
1450       tpcXpri[i] = 0.0;
1451       tpcYpri[i] = 0.0;
1452       tpcXnli[i] = 0.0;
1453       tpcYnli[i] = 0.0;
1454       tpcXnri[i] = 0.0;
1455       tpcYnri[i] = 0.0;
1456       //
1457       tpcXpl[i] = 0.0;
1458       tpcYpl[i] = 0.0;
1459       tpcXpr[i] = 0.0;
1460       tpcYpr[i] = 0.0;
1461       tpcXnl[i] = 0.0;
1462       tpcYnl[i] = 0.0;
1463       tpcXnr[i] = 0.0;
1464       tpcYnr[i] = 0.0;
1465     }
1466
1467   float MCtpcXplo[hmax], MCtpcYplo[hmax], MCtpcXpro[hmax], MCtpcYpro[hmax];
1468   float MCtpcXnlo[hmax], MCtpcYnlo[hmax], MCtpcXnro[hmax], MCtpcYnro[hmax];
1469   float MCtpcXpli[hmax], MCtpcYpli[hmax], MCtpcXpri[hmax], MCtpcYpri[hmax];
1470   float MCtpcXnli[hmax], MCtpcYnli[hmax], MCtpcXnri[hmax], MCtpcYnri[hmax];
1471   float MCtpcXpl[hmax], MCtpcYpl[hmax], MCtpcXpr[hmax], MCtpcYpr[hmax];
1472   float MCtpcXnl[hmax], MCtpcYnl[hmax], MCtpcXnr[hmax], MCtpcYnr[hmax];
1473   for(int i=0; i<hmax;i++)
1474     {
1475       MCtpcXplo[i] = 0.0;
1476       MCtpcYplo[i] = 0.0;
1477       MCtpcXpro[i] = 0.0;
1478       MCtpcYpro[i] = 0.0;
1479       MCtpcXnlo[i] = 0.0;
1480       MCtpcYnlo[i] = 0.0;
1481       MCtpcXnro[i] = 0.0;
1482       MCtpcYnro[i] = 0.0;
1483       //
1484       MCtpcXpli[i] = 0.0;
1485       MCtpcYpli[i] = 0.0;
1486       MCtpcXpri[i] = 0.0;
1487       MCtpcYpri[i] = 0.0;
1488       MCtpcXnli[i] = 0.0;
1489       MCtpcYnli[i] = 0.0;
1490       MCtpcXnri[i] = 0.0;
1491       MCtpcYnri[i] = 0.0;
1492       //
1493       MCtpcXpl[i] = 0.0;
1494       MCtpcYpl[i] = 0.0;
1495       MCtpcXpr[i] = 0.0;
1496       MCtpcYpr[i] = 0.0;
1497       MCtpcXnl[i] = 0.0;
1498       MCtpcYnl[i] = 0.0;
1499       MCtpcXnr[i] = 0.0;
1500       MCtpcYnr[i] = 0.0;
1501     }
1502
1503   // ---------------------------------- //
1504   // --- Now looping over MC tracks --- //
1505   // ---------------------------------- //
1506   if(fMC)
1507     {
1508       // variables for 3rd particle
1509       float MCpt3[d_ntrkMC];
1510       float MCeta3[d_ntrkMC];
1511       //float MCphi3[d_ntrkMC];//new...
1512       int MCcharge3[d_ntrkMC];
1513       // initial variables to out of range values
1514       // for cases when the track loop skips certain tracks
1515       for(int i=0; i<d_ntrk; i++)
1516         {
1517           MCpt3[i] = -99;
1518           MCeta3[i] = -99;
1519           //MCphi3[i] = -99;
1520           MCcharge3[i] = 0;
1521         }
1522       // --- first MC track loop
1523       for(int itrkMC = 0; itrkMC<d_ntrkMC; itrkMC++)
1524         {
1525           //AliVParticle *trackMC = fMC->GetTrack(itrkMC);
1526           AliAODMCParticle *trackMC = (AliAODMCParticle *)fMC->GetTrack(itrkMC);
1527           if(!trackMC)
1528             {
1529               if(debug>0) cout<<"ERROR: Could not receive track "<<itrkMC<<" (mc loop)"<<endl;
1530               continue;
1531             }
1532
1533           float pt = trackMC->Pt();
1534           fHistPtMC->Fill(pt);
1535
1536           float phi = trackMC->Phi();
1537           float eta = trackMC->Eta();
1538           int charge = trackMC->Charge();
1539           if(charge==0) continue; // new
1540           charge /= 3; // new
1541           // int ncls = trackMC->GetMCTPCNcls();
1542           // float dedx = trackMC->GetMCTPCsignal();
1543           // float dcaxy = trackMC->DCA();
1544           // float dcaz = trackMC->ZAtDCA();
1545           bool pos = charge>0.0;
1546           bool neg = charge<0.0;
1547
1548           fHistPhiMC->Fill(phi);
1549           fHistEtaMC->Fill(eta);
1550           fHistChargeMC->Fill(charge);
1551           // fHistMCTPCnclsMC->Fill(ncls);
1552           // fHistDedxMC->Fill(dedx);
1553           // fHistDCAxyMC->Fill(dcaxy);
1554           // fHistDCAzMC->Fill(dcaz);
1555
1556           // apply kinematic cuts
1557           if(fabs(eta)>outeta) continue;
1558           if(pt<ptmin||pt>ptmax) continue;
1559
1560           //if(charge!=-3&&charge!=+3) continue;// x3 by convention
1561           if(!trackMC->IsPrimary()) continue;
1562           if(!trackMC->IsPhysicalPrimary()) continue;
1563           //if(abs(p0->GetPdgCode())==11) continue; //electrons
1564           //if(abs(p0->GetPdgCode())==13) continue; //electrons
1565
1566           ntrkMC++;
1567           if(pos) ntrkposMC++;
1568           if(neg) ntrknegMC++;
1569
1570           if(pt>ptmin&&pt<ptmax&&fabs(eta)<outeta&&fabs(eta)>excleta)
1571             {
1572               cutntrkMC++;
1573               if(pos) cutntrkposMC++;
1574               if(neg) cutntrknegMC++;
1575             }
1576
1577           if(pt>ptmin&&pt<ptmax)
1578             {
1579               // --- outer
1580               if(eta>-outeta&&eta<-ineta)
1581                 {
1582                   for(int i=1; i<hmax;i++)
1583                     {
1584                       if(charge>0)
1585                         {
1586                           MCtpcXplo[i] += cos(i*phi);
1587                           MCtpcYplo[i] += sin(i*phi);
1588                           MCtpcXpl[i] += cos(i*phi);
1589                           MCtpcYpl[i] += sin(i*phi);
1590                         }
1591                       if(charge<0)
1592                         {
1593                           MCtpcXnlo[i] += cos(i*phi);
1594                           MCtpcYnlo[i] += sin(i*phi);
1595                           MCtpcXnl[i] += cos(i*phi);
1596                           MCtpcYnl[i] += sin(i*phi);
1597                         }
1598                     }
1599                   if(charge>0)
1600                     {
1601                       MCtpcXplo[0] += 1.0;
1602                       MCtpcYplo[0] += pt;
1603                       MCtpcXpl[0] += 1.0;
1604                       MCtpcYpl[0] += pt;
1605                     }
1606                   if(charge<0)
1607                     {
1608                       MCtpcXnlo[0] += 1.0;
1609                       MCtpcYnlo[0] += pt;
1610                       MCtpcXnl[0] += 1.0;
1611                       MCtpcYnl[0] += pt;
1612                     }
1613                 } // end left half of MCTPC
1614               //----------------------------
1615               if(eta>ineta&&eta<outeta)
1616                 {
1617                   for(int i=1; i<hmax; i++)
1618                     {
1619                       if(charge>0)
1620                         {
1621                           MCtpcXpro[i] += cos(i*phi);
1622                           MCtpcYpro[i] += sin(i*phi);
1623                           MCtpcXpr[i] += cos(i*phi);
1624                           MCtpcYpr[i] += sin(i*phi);
1625                         }
1626                       if(charge<0)
1627                         {
1628                           MCtpcXnro[i] += cos(i*phi);
1629                           MCtpcYnro[i] += sin(i*phi);
1630                           MCtpcXnr[i] += cos(i*phi);
1631                           MCtpcYnr[i] += sin(i*phi);
1632                         }
1633                     }
1634                   if(charge>0)
1635                     {
1636                       MCtpcXpro[0] += 1.0;
1637                       MCtpcYpro[0] += pt;
1638                       MCtpcXpr[0] += 1.0;
1639                       MCtpcYpr[0] += pt;
1640                     }
1641                   if(charge<0)
1642                     {
1643                       MCtpcXnro[0] += 1.0;
1644                       MCtpcYnro[0] += pt;
1645                       MCtpcXnr[0] += 1.0;
1646                       MCtpcYnr[0] += pt;
1647                     }
1648                 } // end right half of MCTPC
1649               // --- inner
1650               if(eta>-ineta&&eta<-excleta)
1651                 {
1652                   for(int i=1; i<hmax;i++)
1653                     {
1654                       if(charge>0)
1655                         {
1656                           MCtpcXpli[i] += cos(i*phi);
1657                           MCtpcYpli[i] += sin(i*phi);
1658                           MCtpcXpl[i] += cos(i*phi);
1659                           MCtpcYpl[i] += sin(i*phi);
1660                         }
1661                       if(charge<0)
1662                         {
1663                           MCtpcXnli[i] += cos(i*phi);
1664                           MCtpcYnli[i] += sin(i*phi);
1665                           MCtpcXnl[i] += cos(i*phi);
1666                           MCtpcYnl[i] += sin(i*phi);
1667                         }
1668                     }
1669                   if(charge>0)
1670                     {
1671                       MCtpcXpli[0] += 1.0;
1672                       MCtpcYpli[0] += pt;
1673                       MCtpcXpl[0] += 1.0;
1674                       MCtpcYpl[0] += pt;
1675                     }
1676                   if(charge<0)
1677                     {
1678                       MCtpcXnli[0] += 1.0;
1679                       MCtpcYnli[0] += pt;
1680                       MCtpcXnl[0] += 1.0;
1681                       MCtpcYnl[0] += pt;
1682                     }
1683                 } // end left half of MCTPC
1684               //----------------------------
1685               if(eta>excleta&&eta<ineta)
1686                 {
1687                   for(int i=1; i<hmax; i++)
1688                     {
1689                       if(charge>0)
1690                         {
1691                           MCtpcXpri[i] += cos(i*phi);
1692                           MCtpcYpri[i] += sin(i*phi);
1693                           MCtpcXpr[i] += cos(i*phi);
1694                           MCtpcYpr[i] += sin(i*phi);
1695                         }
1696                       if(charge<0)
1697                         {
1698                           MCtpcXnri[i] += cos(i*phi);
1699                           MCtpcYnri[i] += sin(i*phi);
1700                           MCtpcXnr[i] += cos(i*phi);
1701                           MCtpcYnr[i] += sin(i*phi);
1702                         }
1703                     }
1704                   if(charge>0)
1705                     {
1706                       MCtpcXpri[0] += 1.0;
1707                       MCtpcYpri[0] += pt;
1708                       MCtpcXpr[0] += 1.0;
1709                       MCtpcYpr[0] += pt;
1710                     }
1711                   if(charge<0)
1712                     {
1713                       MCtpcXnri[0] += 1.0;
1714                       MCtpcYnri[0] += pt;
1715                       MCtpcXnr[0] += 1.0;
1716                       MCtpcYnr[0] += pt;
1717                     }
1718                 } // end right half of MCTPC
1719             } // end pt selection
1720
1721           h_AT_etaMC->Fill(cent,charge);
1722           h2_AT_etaMC->Fill(cent,charge);
1723
1724           MCpt3[itrkMC] = pt;
1725           MCeta3[itrkMC] = eta;
1726           //MCphi3[itrkMC] = phi;
1727           MCcharge3[itrkMC] = charge;
1728
1729         } // end of first MC track loop
1730
1731       // if(debug>0)
1732       //        {
1733       //          cout<<"number of tracks with wrong filter bit "<<badbit<<endl;
1734       //          cout<<"difference = ntrk - nbad = "<<d_ntrk-badbit<<endl;
1735       //        }
1736       float MCtpcX[9], MCtpcY[9], MCtpcQQ[9];//, qq[9];
1737       float MCtpcXp[9], MCtpcYp[9], MCtpcQQp[9];//, qqp[9]; // pos
1738       float MCtpcXn[9], MCtpcYn[9], MCtpcQQn[9];//, qqn[9]; // neg
1739
1740       float MCqasymm = float(ntrkposMC-ntrknegMC)/float(ntrkMC);
1741       if(doacuts)
1742         {
1743           MCqasymm = float(cutntrkposMC-cutntrknegMC)/float(cutntrkMC);
1744         }
1745       h_A_centMC->Fill(cent,MCqasymm);
1746       h2_A_centMC->Fill(cent,MCqasymm);
1747
1748
1749       for(int i=0; i<6; i++)
1750         {
1751           MCtpcX[i]=MCtpcXpl[i]+MCtpcXnl[i]+MCtpcXpr[i]+MCtpcXnr[i];
1752           MCtpcY[i]=MCtpcYpl[i]+MCtpcYnl[i]+MCtpcYpr[i]+MCtpcYnr[i];
1753           MCtpcQQ[i]=MCtpcX[i]*MCtpcX[i]+MCtpcY[i]*MCtpcY[i];
1754           //qq[i]=sqrt(MCtpcQQ[i]/MCtpcX[0]);
1755           // pos
1756           MCtpcXp[i]=MCtpcXpl[i]+MCtpcXpr[i];
1757           MCtpcYp[i]=MCtpcYpl[i]+MCtpcYpr[i];
1758           MCtpcQQp[i]=MCtpcXp[i]*MCtpcXp[i]+MCtpcYp[i]*MCtpcYp[i];
1759           //qqp[i]=sqrt(MCtpcQQp[i]/MCtpcXp[0]);
1760           // neg
1761           MCtpcXn[i]=MCtpcXnl[i]+MCtpcXnr[i];
1762           MCtpcYn[i]=MCtpcYnl[i]+MCtpcYnr[i];
1763           MCtpcQQn[i]=MCtpcXn[i]*MCtpcXn[i]+MCtpcYn[i]*MCtpcYn[i];
1764           //qqn[i]=sqrt(MCtpcQQn[i]/MCtpcXn[0]);
1765         }
1766
1767       float M = MCtpcX[0];
1768       float W_2 = M*(M-1);
1769       float Mp = MCtpcXp[0];
1770       float Wp_2 = Mp*(Mp-1);
1771       float Mn = MCtpcXn[0];
1772       float Wn_2 = Mn*(Mn-1);
1773
1774       float MCtpcXl2 = MCtpcXnl[2]+MCtpcXpl[2];
1775       float MCtpcYl2 = MCtpcYnl[2]+MCtpcYpl[2];
1776       float MCtpcXr2 = MCtpcXnr[2]+MCtpcXpr[2];
1777       float MCtpcYr2 = MCtpcYnr[2]+MCtpcYpr[2];
1778       float MCtpcXl0 = MCtpcXnl[0]+MCtpcXpl[0];
1779       float MCtpcXr0 = MCtpcXnr[0]+MCtpcXpr[0];
1780
1781       float MCtpcXl2o = MCtpcXnlo[2]+MCtpcXplo[2];
1782       float MCtpcYl2o = MCtpcYnlo[2]+MCtpcYplo[2];
1783       float MCtpcXr2o = MCtpcXnro[2]+MCtpcXpro[2];
1784       float MCtpcYr2o = MCtpcYnro[2]+MCtpcYpro[2];
1785       float MCtpcXl0o = MCtpcXnlo[0]+MCtpcXplo[0];
1786       float MCtpcXr0o = MCtpcXnro[0]+MCtpcXpro[0];
1787
1788       float MCq22ev = (MCtpcQQ[2]-M)/W_2;
1789       float MCq22Pev = (MCtpcQQp[2]-Mp)/Wp_2;
1790       float MCq22Nev = (MCtpcQQn[2]-Mn)/Wn_2;
1791       float MCq22gap0ev = (MCtpcXl2*MCtpcXr2+MCtpcYl2*MCtpcYr2)/(MCtpcXl0*MCtpcXr0);
1792       float MCq22gap0Pev = (MCtpcXpl[2]*MCtpcXpr[2]+MCtpcYpl[2]*MCtpcYpr[2])/(MCtpcXpl[0]*MCtpcXpr[0]);
1793       float MCq22gap0Nev = (MCtpcXnl[2]*MCtpcXnr[2]+MCtpcYnl[2]*MCtpcYnr[2])/(MCtpcXnl[0]*MCtpcXnr[0]);
1794       float MCq22gap1ev = (MCtpcXl2o*MCtpcXr2o+MCtpcYl2o*MCtpcYr2o)/(MCtpcXl0o*MCtpcXr0o);
1795       float MCq22gap1Pev = (MCtpcXplo[2]*MCtpcXpro[2]+MCtpcYplo[2]*MCtpcYpro[2])/(MCtpcXplo[0]*MCtpcXpro[0]);
1796       float MCq22gap1Nev = (MCtpcXnlo[2]*MCtpcXnro[2]+MCtpcYnlo[2]*MCtpcYnro[2])/(MCtpcXnlo[0]*MCtpcXnro[0]);
1797       h_MCq22_cent->Fill(cent,MCq22ev);
1798       h_MCq22_centP->Fill(cent,MCq22Pev);
1799       h_MCq22_centN->Fill(cent,MCq22Nev);
1800       h_MCAq22_cent->Fill(cent,MCq22ev*MCqasymm);
1801       h_MCAq22_centP->Fill(cent,MCq22Pev*MCqasymm);
1802       h_MCAq22_centN->Fill(cent,MCq22Nev*MCqasymm);
1803       h_MCq22gap0_cent->Fill(cent,MCq22gap0ev);
1804       h_MCq22gap0_centP->Fill(cent,MCq22gap0Pev);
1805       h_MCq22gap0_centN->Fill(cent,MCq22gap0Nev);
1806       h_MCq22gap1_cent->Fill(cent,MCq22gap1ev);
1807       h_MCq22gap1_centP->Fill(cent,MCq22gap1Pev);
1808       h_MCq22gap1_centN->Fill(cent,MCq22gap1Nev);
1809
1810
1811
1812       if(debug>0) cout<<"there are "<<ntrkMC<<" Monte Carlo tracks within the cuts in this event"<<endl;
1813
1814       // --- fsecond MC track loop
1815       for(int itrkMC = 0; itrkMC<d_ntrkMC; itrkMC++)
1816         {
1817           //AliVParticle *trackMC = fMC->GetTrack(itrkMC);
1818           AliAODMCParticle *trackMC = (AliAODMCParticle *)fMC->GetTrack(itrkMC);
1819           if(!trackMC)
1820             {
1821               if(debug>0) cout<<"ERROR: Could not receive track "<<itrkMC<<" (mc loop)"<<endl;
1822               continue;
1823             }
1824
1825           float pt = trackMC->Pt();
1826
1827           float phi = trackMC->Phi();
1828           float eta = trackMC->Eta();
1829           int charge = trackMC->Charge();
1830           if(charge==0) continue; // new
1831           charge /= 3; // new
1832           bool pos = charge>0.0;
1833           bool neg = charge<0.0;
1834
1835
1836           // apply kinematic cuts
1837           if(fabs(eta)>outeta) continue;
1838           if(pt<ptmin||pt>ptmax) continue;
1839
1840           //if(charge!=-3&&charge!=+3) continue;// x3 by convention
1841           if(!trackMC->IsPrimary()) continue;
1842           if(!trackMC->IsPhysicalPrimary()) continue;
1843
1844           //if(abs(p0->GetPdgCode())==11) continue; //electrons
1845           //if(abs(p0->GetPdgCode())==13) continue; //electrons
1846
1847
1848           float MCtrkXpl[hmax], MCtrkYpl[hmax], MCtrkXpr[hmax], MCtrkYpr[hmax];
1849           float MCtrkXnl[hmax], MCtrkYnl[hmax], MCtrkXnr[hmax], MCtrkYnr[hmax];
1850           for(int i=0; i<hmax;i++)
1851             {
1852               MCtrkXpl[i] = 0.0;
1853               MCtrkYpl[i] = 0.0;
1854               MCtrkXpr[i] = 0.0;
1855               MCtrkYpr[i] = 0.0;
1856               MCtrkXnl[i] = 0.0;
1857               MCtrkYnl[i] = 0.0;
1858               MCtrkXnr[i] = 0.0;
1859               MCtrkYnr[i] = 0.0;
1860             }
1861           // very similar to above, except += changed to =
1862           if(pt>ptmin&&pt<ptmax)
1863             {
1864               if(eta>-outeta&&eta<-excleta)
1865                 {
1866                   // assign values to array elements 1..8
1867                   for(int i=1; i<hmax;i++)
1868                     {
1869                       if(charge>0)
1870                         {
1871                           MCtrkXpl[i] = cos(i*phi);
1872                           MCtrkYpl[i] = sin(i*phi);
1873                         }
1874                       if(charge<0)
1875                         {
1876                           MCtrkXnl[i] = cos(i*phi);
1877                           MCtrkYnl[i] = sin(i*phi);
1878                         }
1879                     }
1880                   // assign values to array element 0
1881                   if(charge>0)
1882                     {
1883                       MCtrkXpl[0] = 1.0;
1884                       MCtrkYpl[0] = pt;
1885                     }
1886                   if(charge<0)
1887                     {
1888                       MCtrkXnl[0] = 1.0;
1889                       MCtrkYnl[0] = pt;
1890                     }
1891                 }
1892               //----------------------------
1893               if(eta>excleta&&eta<outeta) // right half of TPC
1894                 {
1895                   for(int i=1; i<hmax; i++)
1896                     {
1897                       if(charge>0)
1898                         {
1899                           MCtrkXpr[i] = cos(i*phi);
1900                           MCtrkYpr[i] = sin(i*phi);
1901                         }
1902                       if(charge<0)
1903                         {
1904                           MCtrkXnr[i] = cos(i*phi);
1905                           MCtrkYnr[i] = sin(i*phi);
1906                         }
1907                     }
1908                   if(charge>0)
1909                     {
1910                       MCtrkXpr[0] = 1.0;
1911                       MCtrkYpr[0] = pt;
1912                     }
1913                   if(charge<0)
1914                     {
1915                       MCtrkXnr[0] = 1.0;
1916                       MCtrkYnr[0] = pt;
1917                     }
1918                 } // end right half of TPC
1919             } // end pT selection
1920
1921           // COME BACK HERE
1922           float MCtrkX[9],  MCtrkY[9];
1923           float MCtrkXp[9], MCtrkYp[9]; // pos
1924           float MCtrkXn[9], MCtrkYn[9]; // neg
1925           for(int i=0; i<6; i++)
1926             {
1927               MCtrkX[i]=MCtrkXpl[i]+MCtrkXnl[i]+MCtrkXpr[i]+MCtrkXnr[i];
1928               MCtrkY[i]=MCtrkYpl[i]+MCtrkYnl[i]+MCtrkYpr[i]+MCtrkYnr[i];
1929               // pos
1930               MCtrkXp[i]=MCtrkXpl[i]+MCtrkXpr[i];
1931               MCtrkYp[i]=MCtrkYpl[i]+MCtrkYpr[i];
1932               // neg
1933               MCtrkXn[i]=MCtrkXnl[i]+MCtrkXnr[i];
1934               MCtrkYn[i]=MCtrkYnl[i]+MCtrkYnr[i];
1935             }
1936
1937           // sanity checks to prevent division by zero and other problems
1938           if(MCtrkX[0]<1) continue;
1939           if(MCtpcX[0]<2) continue;
1940           if(MCtpcXp[0]<2) continue;
1941           if(MCtpcXn[0]<2) continue;
1942
1943           // --- calculate differential cumulants from Q-vector components
1944           // --- reference flow is from whole TPC
1945           float MCdiffq22ev = ((MCtrkX[2]*MCtpcX[2])+(MCtrkY[2]*MCtpcY[2])-1)/(MCtpcX[0]-1);
1946           float MCdiffq22Pev = ((MCtrkXp[2]*MCtpcX[2])+(MCtrkYp[2]*MCtpcY[2])-1)/(MCtpcX[0]-1);
1947           float MCdiffq22Nev = ((MCtrkXn[2]*MCtpcX[2])+(MCtrkYn[2]*MCtpcY[2])-1)/(MCtpcX[0]-1);
1948           // --- third harmonic
1949           float MCdiffq32ev = ((MCtrkX[3]*MCtpcX[3])+(MCtrkY[3]*MCtpcY[3])-1)/(MCtpcX[0]-1);
1950           float MCdiffq32Pev = ((MCtrkXp[3]*MCtpcX[3])+(MCtrkYp[3]*MCtpcY[3])-1)/(MCtpcX[0]-1);
1951           float MCdiffq32Nev = ((MCtrkXn[3]*MCtpcX[3])+(MCtrkYn[3]*MCtpcY[3])-1)/(MCtpcX[0]-1);
1952           // --- fourth harmonic
1953           float MCdiffq42ev = ((MCtrkX[4]*MCtpcX[4])+(MCtrkY[4]*MCtpcY[4])-1)/(MCtpcX[0]-1);
1954           float MCdiffq42Pev = ((MCtrkXp[4]*MCtpcX[4])+(MCtrkYp[4]*MCtpcY[4])-1)/(MCtpcX[0]-1);
1955           float MCdiffq42Nev = ((MCtrkXn[4]*MCtpcX[4])+(MCtrkYn[4]*MCtpcY[4])-1)/(MCtpcX[0]-1);
1956
1957
1958           if(donested)
1959             {
1960               for(int i=0; i<d_ntrkMC; i++)
1961                 {
1962                   // ----------------------------------------------------------------------------------------
1963                   // --- want delta eta and delta pt dependence of differential cumulant weighted with charge
1964                   // --- need nested track loop to get eta1 and eta3
1965                   // ----------------------------------------------------------------------------------------
1966                   // make sure eta is in range, -99 should be autorejected by this
1967                   //if(MCeta3[i]==-99) continue;
1968                   if(MCeta3[i]<-outeta||MCeta3[i]>outeta)
1969                     {
1970                       //cout<<"MCeta3 out of range!!! "<<MCeta3[i]<<endl;
1971                       continue;
1972                     }
1973                   if(eta<-outeta||eta>outeta)
1974                     {
1975                       cout<<"eta out of range!!! "<<eta<<endl;
1976                       continue;
1977                     }
1978                   if(MCpt3[i]<ptmin||MCpt3[i]>ptmax)
1979                     {
1980                       //cout<<"MCpt3 out of range!!! "<<MCpt3[i]<<endl;
1981                       continue;
1982                     }
1983                   if(pt<ptmin||pt>ptmax)
1984                     {
1985                       cout<<"pt out of range!!!"<<endl;
1986                       continue;
1987                     }
1988                   // charge==0 should be rejected by the above cut???
1989                   if(MCcharge3[i]!=1&&MCcharge3[i]!=-1)
1990                     {
1991                       //cout<<"WTF!"<<endl;
1992                       continue;
1993                     }
1994                   // VERY IMPORTANT remove auto-correlations with 3rd particle
1995                   if(i==itrkMC) continue;
1996
1997                   float DETA = eta-MCeta3[i];
1998                   //float DPHI = phi-MCphi3[i];
1999                   //float DPT = pt-MCpt3[i];
2000
2001                   // ---
2002
2003                   hMC_AT_X_deta->Fill(DETA,MCcharge3[i]);
2004                   if(pos) hMC_AT_X_detaP->Fill(DETA,MCcharge3[i]);
2005                   if(neg) hMC_AT_X_detaN->Fill(DETA,MCcharge3[i]);
2006                   if(pos) hMC_diffq22_X_detaP->Fill(DETA,MCdiffq22Pev);
2007                   if(neg) hMC_diffq22_X_detaN->Fill(DETA,MCdiffq22Nev);
2008                   if(pos) hMC_diffq32_X_detaP->Fill(DETA,MCdiffq32Pev);
2009                   if(neg) hMC_diffq32_X_detaN->Fill(DETA,MCdiffq32Nev);
2010                   if(pos) hMC_diffq42_X_detaP->Fill(DETA,MCdiffq42Pev);
2011                   if(neg) hMC_diffq42_X_detaN->Fill(DETA,MCdiffq42Nev);
2012                   if(pos) hMC_ATdiffq22_X_detaP->Fill(DETA,MCdiffq22Pev*MCcharge3[i]);
2013                   if(neg) hMC_ATdiffq22_X_detaN->Fill(DETA,MCdiffq22Nev*MCcharge3[i]);
2014                   if(pos) hMC_ATdiffq32_X_detaP->Fill(DETA,MCdiffq32Pev*MCcharge3[i]);
2015                   if(neg) hMC_ATdiffq32_X_detaN->Fill(DETA,MCdiffq32Nev*MCcharge3[i]);
2016                   if(pos) hMC_ATdiffq42_X_detaP->Fill(DETA,MCdiffq42Pev*MCcharge3[i]);
2017                   if(neg) hMC_ATdiffq42_X_detaN->Fill(DETA,MCdiffq42Nev*MCcharge3[i]);
2018
2019                 } // nested track loop
2020
2021             } // check on donested
2022
2023         } // end of second MC track loop
2024
2025     } // check on existence of MC
2026
2027
2028
2029   // ------------------------------- //
2030   // --- Now looping over tracks --- //
2031   // ------------------------------- //
2032
2033   // --- track loop for mapping matrix
2034   TExMap *trackMap = new TExMap();
2035   for(int itrk=0; itrk<d_ntrk; itrk++)
2036     {
2037       AliAODTrack *track = fAOD->GetTrack(itrk);
2038       if(!track)
2039         {
2040           if(debug>0) cout<<"ERROR: Could not retrieve AODtrack "<<itrk<<endl;
2041           continue;
2042         }
2043       int gid = track->GetID();
2044       if(track->TestFilterBit(fbit)) trackMap->Add(gid,itrk);
2045     }
2046
2047
2048
2049   // variables for 3rd particle
2050   float pt3[d_ntrk];
2051   float eta3[d_ntrk];
2052   float phi3[d_ntrk];//new...
2053   int charge3[d_ntrk];
2054   // initial variables to out of range values
2055   // for cases when the track loop skips certain tracks
2056   for(int i=0; i<d_ntrk; i++)
2057     {
2058       pt3[i] = -99;
2059       eta3[i] = -99;
2060       phi3[i] = -99;
2061       charge3[i] = 0;
2062     }
2063
2064   int badbit = 0; // counter for number of tracks with wrong filter bit
2065   // --- main track loop
2066   for(int itrk = 0; itrk<d_ntrk; itrk++)
2067     {
2068       AliAODTrack *track = fAOD->GetTrack(itrk);
2069       if(!track)
2070         {
2071           if(debug>0) cout<<"ERROR: Could not retrieve AODtrack "<<itrk<<endl;
2072           continue;
2073         }
2074
2075       if(!track->TestFilterBit(fbit))
2076         {
2077           badbit++; // count tracks with wrong filter bit
2078           if(debug>15) cout<<"wrong filter bit for track "<<itrk<<endl;
2079           continue;
2080         }
2081
2082       int gid = track->GetID();
2083       AliAODTrack *PIDtrack;
2084       if(gid>=0) PIDtrack = track;
2085       else PIDtrack = fAOD->GetTrack(trackMap->GetValue(-1-gid));
2086
2087
2088       // if(debug>15)
2089       //        {
2090       //          cout<<"track index is "<<itrk<<endl;
2091       //          cout<<"track ID is "<<gid<<endl;
2092       //        }
2093
2094       // --- get track variables
2095       float pt = track->Pt();
2096       float phi = track->Phi();
2097       float eta = track->Eta();
2098       int charge = track->Charge();
2099       int ncls = track->GetTPCNcls();
2100       float dedx = track->GetTPCsignal();
2101
2102       float Tdcaxy = track->DCA();
2103       float Tdcaz = track->ZAtDCA();
2104
2105
2106       float dcaxy = -999;
2107       float dcax = -999;
2108       float dcay = -999;
2109       float dcaz = -999;
2110
2111       double r[3];
2112       bool dcaflag = track->GetXYZ(r);
2113       //cout<<"fbit is "<<fbit<<" GetXYZ is "<<dcaflag<<endl;
2114       //printf("fbit is %d GetXYZ is %d PropagateToDCA is %d Tdcaz is %f pt is %f \n\n",fbit,dcaflag,proptodca,Tdcaz,pt);
2115
2116       double DCA[2]; // dca
2117       double COV[3]; // covariance
2118       bool proptodca = track->PropagateToDCA(fVtx,mag,100.0,DCA,COV);
2119       if(!proptodca&&debug>0) cout<<"No DCACOV for you!"<<endl;
2120
2121       if(dcaflag)
2122         {
2123           dcaxy = r[0];
2124           dcaz = r[1];
2125           if(debug>5) cout<<"GetXYZ is true, filter bit is "<<fbit<<endl;
2126         }
2127       else
2128         {
2129           dcax = r[0] - eventX;
2130           dcay = r[1] - eventY;
2131           dcaz = r[2] - eventZ;
2132           // --- need the sign convention for dcaxy...
2133           dcaxy = sqrt(dcax*dcax+dcay*dcay);
2134           //if((float)dcaxy!=(float)fabs((float)DCA[0])&&debug>4) cout<<"hmm... "<<dcaxy<<" "<<DCA[0]<<endl;
2135           // --- set dcaxy to value from PropagateToDCA to get correct sign
2136           dcaxy = DCA[0];
2137           // --- dcaz on the other hand is unambiguous
2138           //if(dcaz!=(float)DCA[1]&&debug>4) cout<<"hmm... "<<dcaz<<" "<<DCA[1]<<endl;
2139           if(debug>5) cout<<"GetXYZ is false, filter bit is "<<fbit<<endl;
2140         }
2141
2142       if(debug>5)
2143         {
2144           cout<<"r[0] is "<<r[0]<<" ";
2145           cout<<"r[1] is "<<r[1]<<" ";
2146           cout<<"r[2] is "<<r[2]<<endl;
2147           cout<<"eventX is "<<eventX<<" ";
2148           cout<<"eventY is "<<eventY<<" ";
2149           cout<<"eventZ is "<<eventZ<<endl;
2150           cout<<"Tdcaxy is "<<Tdcaxy<<" and dcaxy is "<<dcaxy<<" and DCACOV xy is "<<DCA[0]<<endl;
2151           cout<<"Tdcaz is "<<Tdcaz<<" and dcaz is "<<dcaz<<" and DCACOV z is "<<DCA[1]<<endl;
2152         }
2153
2154       if((fbit==128||fbit==272)&&!dcaflag)
2155         {
2156           dcaxy = Tdcaxy;
2157           dcaz = Tdcaz;
2158           if(debug>4)
2159             {
2160               cout<<"GetXYZ false but should be true, flipping values"<<endl;
2161               cout<<"Tdcaxy is "<<Tdcaxy<<" and dcaxy is "<<dcaxy<<" and DCACOV xy is "<<DCA[0]<<endl;
2162               cout<<"Tdcaz is "<<Tdcaz<<" and dcaz is "<<dcaz<<" and DCACOV z is "<<DCA[1]<<endl;
2163             }
2164         }
2165
2166       // --- some diagnostic histograms for tracks
2167       fHistPt->Fill(pt);
2168       fHistPhi->Fill(phi);
2169       fHistEta->Fill(eta);
2170       fHistCharge->Fill(charge);
2171       fHistTPCncls->Fill(ncls);
2172       fHistDedx->Fill(dedx);
2173       fHistDCAxy->Fill(dcaxy);
2174       fHistDCAz->Fill(dcaz);
2175
2176       if(charge>0)
2177         {
2178           fHistPosPt->Fill(pt);
2179           fHistPosPhi->Fill(phi);
2180           fHistPosEta->Fill(eta);
2181         }
2182       if(charge<0)
2183         {
2184           fHistNegPt->Fill(pt);
2185           fHistNegPhi->Fill(phi);
2186           fHistNegEta->Fill(eta);
2187         }
2188       fProfMeanChargePt->Fill(pt,charge);
2189       fProfMeanChargeEta->Fill(eta,charge);
2190
2191       // --- some extra verbose stuff
2192       if(debug>8&&(pt<ptmin||pt>ptmax)) cout<<"pt = "<<pt<<" out of range for Q-vectors"<<endl;
2193       if(debug>9) cout<<"number of tpc clusters is "<<ncls<<endl;
2194
2195       // --- track cut on number of TPC clusters
2196       if(ncls<nclscut)
2197         {
2198           if(debug>5) cout<<"number of tpc clusters is too low "<<ncls<<endl;
2199           continue;
2200         }
2201
2202       // --- track cut on dca
2203       if(fabs(dcaxy)>dcacutxy||fabs(dcaz)>dcacutz)
2204         {
2205           if(debug>7)
2206             {
2207               cout<<"dca out of range..."<<endl;
2208               cout<<"dcaxy cut is "<<dcacutxy<<endl;
2209               cout<<"dcaz cut is "<<dcacutz<<endl;
2210               cout<<"dcaxy is "<<dcaxy<<endl;
2211               cout<<"dcaz is "<<dcaz<<endl;
2212             }
2213           if(dodcacuts) continue; // problems depending on filter bit...
2214         }
2215       // Prabhat
2216       // fHistPhi->Fill(phi);
2217       // fHistEta->Fill(eta);
2218       fHistDCAxyAfter->Fill(dcaxy);
2219       fHistDCAzAfter->Fill(dcaz);
2220
2221       float nsigmapion = 0.;
2222       float nsigmakaon = 0.;
2223       float nsigmaprot = 0.;
2224       float nsigmaelec = 0.;
2225
2226       if(!doMC){
2227         nsigmapion = fPID->NumberOfSigmasTPC(PIDtrack,(AliPID::EParticleType)AliPID::kPion);
2228         nsigmakaon = fPID->NumberOfSigmasTPC(PIDtrack,(AliPID::EParticleType)AliPID::kKaon);
2229         nsigmaprot = fPID->NumberOfSigmasTPC(PIDtrack,(AliPID::EParticleType)AliPID::kProton);
2230         nsigmaelec = fPID->NumberOfSigmasTPC(PIDtrack,(AliPID::EParticleType)AliPID::kElectron);
2231       }
2232
2233       bool isPion = fabs(nsigmapion) <= nspid;
2234       bool isKaon = fabs(nsigmakaon) <= nspid;
2235       bool isProt = fabs(nsigmaprot) <= nspid;
2236       bool isElec = fabs(nsigmaelec) <= nspid;
2237
2238       bool isPionH = isPion && (!isKaon) && (!isProt);
2239       bool isProtH = (!isPion) && (!isKaon) && isProt;
2240       bool isPionL = isPion && (!isKaon) && (!isProt) && (!isElec);
2241       bool isProtL = (!isPion) && (!isKaon) && isProt && (!isElec);
2242
2243       if(isPion) fHistPtPion->Fill(pt);
2244       if(isPionH) fHistPtPionH->Fill(pt);
2245       if(isPionL) fHistPtPionL->Fill(pt);
2246       fHistNsigmaPion->Fill(nsigmapion);
2247
2248       if(isProt) fHistPtProt->Fill(pt);
2249       if(isProtH) fHistPtProtH->Fill(pt);
2250       if(isProtL) fHistPtProtL->Fill(pt);
2251       fHistNsigmaProt->Fill(nsigmapion);
2252
2253       // --- count track numbers
2254       ntrk++;
2255       if(charge>0) ntrkpos++;
2256       if(charge<0) ntrkneg++;
2257       if(eta<0)
2258         {
2259           ntrkL++;
2260           if(charge>0) ntrkposL++;
2261           if(charge<0) ntrknegL++;
2262         }
2263       if(eta>0)
2264         {
2265           ntrkR++;
2266           if(charge>0) ntrkposR++;
2267           if(charge<0) ntrknegR++;
2268         }
2269       h_AT_eta->Fill(eta,charge);
2270       if(mag<-4.0) h_AT_etaF1->Fill(eta,charge);
2271       if(mag>4.0) h_AT_etaF3->Fill(eta,charge);
2272       h2_AT_eta->Fill(eta,charge);
2273       if(mag<-4.0) h2_AT_etaF1->Fill(eta,charge);
2274       if(mag>4.0) h2_AT_etaF3->Fill(eta,charge);
2275       //
2276       if(mag<-4.0&&charge>0) h_eta_pos_F1->Fill(eta);
2277       if(mag<-4.0&&charge<0) h_eta_neg_F1->Fill(eta);
2278       if(mag>4.0&&charge>0) h_eta_pos_F3->Fill(eta);
2279       if(mag>4.0&&charge<0) h_eta_neg_F3->Fill(eta);
2280
2281       if(doeffcorr)
2282         {
2283           if(gRandom) delete gRandom;
2284           gRandom = new TRandom3(0);
2285           gRandom->SetSeed(0);
2286           double rand = gRandom->Rndm();
2287           int ptbin = int(pt*10);
2288           if(ptbin>49) ptbin = 49;
2289           float eff = effTPC[ptbin];
2290           if(fbit==272) eff = effHYB[ptbin];
2291           if(rand>eff)
2292             {
2293               if(debug>5) cout<<"doing random throw for efficiency correction: pt is "<<pt<<" eff is "<<eff<<" rand is "<<rand<<endl;
2294               continue;
2295             }
2296         }
2297
2298       if(pt>ptmin&&pt<ptmax&&fabs(eta)<outeta&&fabs(eta)>excleta)
2299         {
2300           cutntrk++;
2301           if(charge>0) cutntrkpos++;
2302           if(charge<0) cutntrkneg++;
2303           if(eta<0)
2304             {
2305               cutntrkL++;
2306               if(charge>0) cutntrkposL++;
2307               if(charge<0) cutntrknegL++;
2308             }
2309           if(eta>0)
2310             {
2311               cutntrkR++;
2312               if(charge>0) cutntrkposR++;
2313               if(charge<0) cutntrknegR++;
2314             }
2315           h_AT_cut_eta->Fill(eta,charge);
2316           if(mag<-4.0) h_AT_cut_etaF1->Fill(eta,charge);
2317           if(mag>4.0) h_AT_cut_etaF3->Fill(eta,charge);
2318           h2_AT_cut_eta->Fill(eta,charge);
2319           if(mag<-4.0) h2_AT_cut_etaF1->Fill(eta,charge);
2320           if(mag>4.0) h2_AT_cut_etaF3->Fill(eta,charge);
2321           //
2322           if(cent>=centlo&&cent<centhi)
2323             {
2324               if(mag<-4.0&&charge>0) h_cut_eta_pos_F1->Fill(eta);
2325               if(mag<-4.0&&charge<0) h_cut_eta_neg_F1->Fill(eta);
2326               if(mag>4.0&&charge>0) h_cut_eta_pos_F3->Fill(eta);
2327               if(mag>4.0&&charge<0) h_cut_eta_neg_F3->Fill(eta);
2328             }
2329         }
2330
2331       // ------------------------------------ //
2332       // --- now fill Q-vector components --- //
2333       // ------------------------------------ //
2334
2335       if(pt>ptmin&&pt<ptmax)
2336         {
2337           // --- outer
2338           if(eta>-outeta&&eta<-ineta)
2339             {
2340               for(int i=1; i<hmax;i++)
2341                 {
2342                   if(charge>0)
2343                     {
2344                       tpcXplo[i] += cos(i*phi);
2345                       tpcYplo[i] += sin(i*phi);
2346                       tpcXpl[i] += cos(i*phi);
2347                       tpcYpl[i] += sin(i*phi);
2348                     }
2349                   if(charge<0)
2350                     {
2351                       tpcXnlo[i] += cos(i*phi);
2352                       tpcYnlo[i] += sin(i*phi);
2353                       tpcXnl[i] += cos(i*phi);
2354                       tpcYnl[i] += sin(i*phi);
2355                     }
2356                 }
2357               if(charge>0)
2358                 {
2359                   tpcXplo[0] += 1.0;
2360                   tpcYplo[0] += pt;
2361                   tpcXpl[0] += 1.0;
2362                   tpcYpl[0] += pt;
2363                 }
2364               if(charge<0)
2365                 {
2366                   tpcXnlo[0] += 1.0;
2367                   tpcYnlo[0] += pt;
2368                   tpcXnl[0] += 1.0;
2369                   tpcYnl[0] += pt;
2370                 }
2371             } // end left half of TPC
2372           //----------------------------
2373           if(eta>ineta&&eta<outeta)
2374             {
2375               for(int i=1; i<hmax; i++)
2376                 {
2377                   if(charge>0)
2378                     {
2379                       tpcXpro[i] += cos(i*phi);
2380                       tpcYpro[i] += sin(i*phi);
2381                       tpcXpr[i] += cos(i*phi);
2382                       tpcYpr[i] += sin(i*phi);
2383                     }
2384                   if(charge<0)
2385                     {
2386                       tpcXnro[i] += cos(i*phi);
2387                       tpcYnro[i] += sin(i*phi);
2388                       tpcXnr[i] += cos(i*phi);
2389                       tpcYnr[i] += sin(i*phi);
2390                     }
2391                 }
2392               if(charge>0)
2393                 {
2394                   tpcXpro[0] += 1.0;
2395                   tpcYpro[0] += pt;
2396                   tpcXpr[0] += 1.0;
2397                   tpcYpr[0] += pt;
2398                 }
2399               if(charge<0)
2400                 {
2401                   tpcXnro[0] += 1.0;
2402                   tpcYnro[0] += pt;
2403                   tpcXnr[0] += 1.0;
2404                   tpcYnr[0] += pt;
2405                 }
2406             } // end right half of TPC
2407           // --- inner
2408           if(eta>-ineta&&eta<-excleta)
2409             {
2410               for(int i=1; i<hmax;i++)
2411                 {
2412                   if(charge>0)
2413                     {
2414                       tpcXpli[i] += cos(i*phi);
2415                       tpcYpli[i] += sin(i*phi);
2416                       tpcXpl[i] += cos(i*phi);
2417                       tpcYpl[i] += sin(i*phi);
2418                     }
2419                   if(charge<0)
2420                     {
2421                       tpcXnli[i] += cos(i*phi);
2422                       tpcYnli[i] += sin(i*phi);
2423                       tpcXnl[i] += cos(i*phi);
2424                       tpcYnl[i] += sin(i*phi);
2425                     }
2426                 }
2427               if(charge>0)
2428                 {
2429                   tpcXpli[0] += 1.0;
2430                   tpcYpli[0] += pt;
2431                   tpcXpl[0] += 1.0;
2432                   tpcYpl[0] += pt;
2433                 }
2434               if(charge<0)
2435                 {
2436                   tpcXnli[0] += 1.0;
2437                   tpcYnli[0] += pt;
2438                   tpcXnl[0] += 1.0;
2439                   tpcYnl[0] += pt;
2440                 }
2441             } // end left half of TPC
2442           //----------------------------
2443           if(eta>excleta&&eta<ineta)
2444             {
2445               for(int i=1; i<hmax; i++)
2446                 {
2447                   if(charge>0)
2448                     {
2449                       tpcXpri[i] += cos(i*phi);
2450                       tpcYpri[i] += sin(i*phi);
2451                       tpcXpr[i] += cos(i*phi);
2452                       tpcYpr[i] += sin(i*phi);
2453                     }
2454                   if(charge<0)
2455                     {
2456                       tpcXnri[i] += cos(i*phi);
2457                       tpcYnri[i] += sin(i*phi);
2458                       tpcXnr[i] += cos(i*phi);
2459                       tpcYnr[i] += sin(i*phi);
2460                     }
2461                 }
2462               if(charge>0)
2463                 {
2464                   tpcXpri[0] += 1.0;
2465                   tpcYpri[0] += pt;
2466                   tpcXpr[0] += 1.0;
2467                   tpcYpr[0] += pt;
2468                 }
2469               if(charge<0)
2470                 {
2471                   tpcXnri[0] += 1.0;
2472                   tpcYnri[0] += pt;
2473                   tpcXnr[0] += 1.0;
2474                   tpcYnr[0] += pt;
2475                 }
2476             } // end right half of TPC
2477         } // end pt selection
2478
2479       // set values for 3rd particle variables
2480       pt3[itrk] = pt;
2481       eta3[itrk] = eta;
2482       phi3[itrk] = phi;
2483       charge3[itrk] = (int)charge;
2484
2485
2486       // ------------------------------------ //
2487       // --- that's it for the track loop --- //
2488       // ------------------------------------ //
2489
2490     } // end of first track loop
2491
2492   if(debug>0)
2493     {
2494       cout<<"filter bit is "<<fbit<<endl;
2495       cout<<"number of tracks with wrong filter bit "<<badbit<<endl;
2496       cout<<"difference = ntrk - nbad = "<<d_ntrk-badbit<<endl;
2497     }
2498
2499   float tpcX[9], tpcY[9], tpcQQ[9];//, qq[9];
2500   float tpcXp[9], tpcYp[9], tpcQQp[9];//, qqp[9]; // pos
2501   float tpcXn[9], tpcYn[9], tpcQQn[9];//, qqn[9]; // neg
2502
2503   for(int i=0; i<6; i++)
2504     {
2505       tpcX[i]=tpcXpl[i]+tpcXnl[i]+tpcXpr[i]+tpcXnr[i];
2506       tpcY[i]=tpcYpl[i]+tpcYnl[i]+tpcYpr[i]+tpcYnr[i];
2507       tpcQQ[i]=tpcX[i]*tpcX[i]+tpcY[i]*tpcY[i];
2508       //qq[i]=sqrt(tpcQQ[i]/tpcX[0]);
2509       // pos
2510       tpcXp[i]=tpcXpl[i]+tpcXpr[i];
2511       tpcYp[i]=tpcYpl[i]+tpcYpr[i];
2512       tpcQQp[i]=tpcXp[i]*tpcXp[i]+tpcYp[i]*tpcYp[i];
2513       //qqp[i]=sqrt(tpcQQp[i]/tpcXp[0]);
2514       // neg
2515       tpcXn[i]=tpcXnl[i]+tpcXnr[i];
2516       tpcYn[i]=tpcYnl[i]+tpcYnr[i];
2517       tpcQQn[i]=tpcXn[i]*tpcXn[i]+tpcYn[i]*tpcYn[i];
2518       //qqn[i]=sqrt(tpcQQn[i]/tpcXn[0]);
2519     }
2520
2521
2522
2523
2524
2525   float M = tpcX[0];
2526   float W_2 = M*(M-1);
2527   float Mp = tpcXp[0];
2528   float Wp_2 = Mp*(Mp-1);
2529   float Mn = tpcXn[0];
2530   float Wn_2 = Mn*(Mn-1);
2531
2532   // prevent division by zero
2533   if(M<2) return;
2534   if(Mp<2) return;
2535   if(Mn<2) return;
2536
2537   if(ntrkL<1) return;
2538   if(ntrkR<1) return;
2539
2540   float qasymm = (float)(ntrkpos-ntrkneg)/ntrk;
2541   float qasymmL = (float)(ntrkposL-ntrknegL)/ntrkL;
2542   float qasymmR = (float)(ntrkposR-ntrknegR)/ntrkR;
2543
2544   if(doacuts)
2545     {
2546       qasymm = (float)(cutntrkpos-cutntrkneg)/cutntrk;
2547       qasymmL = (float)(cutntrkposL-cutntrknegL)/cutntrkL;
2548       qasymmR = (float)(cutntrkposR-cutntrknegR)/cutntrkR;
2549     }
2550
2551   h_A_cent->Fill(cent,qasymm);
2552   if(mag<-4.0) h_A_centF1->Fill(cent,qasymm);
2553   if(mag>4.0) h_A_centF3->Fill(cent,qasymm);
2554   h2_A_cent->Fill(cent,qasymm);
2555   if(mag<-4.0) h2_A_centF1->Fill(cent,qasymm);
2556   if(mag>4.0) h2_A_centF3->Fill(cent,qasymm);
2557
2558   float Xl2 = tpcXnl[2]+tpcXpl[2];
2559   float Yl2 = tpcYnl[2]+tpcYpl[2];
2560   float Xr2 = tpcXnr[2]+tpcXpr[2];
2561   float Yr2 = tpcYnr[2]+tpcYpr[2];
2562   float Xl3 = tpcXnl[3]+tpcXpl[3];
2563   float Yl3 = tpcYnl[3]+tpcYpl[3];
2564   float Xr3 = tpcXnr[3]+tpcXpr[3];
2565   float Yr3 = tpcYnr[3]+tpcYpr[3];
2566   float Xl4 = tpcXnl[4]+tpcXpl[4];
2567   float Yl4 = tpcYnl[4]+tpcYpl[4];
2568   float Xr4 = tpcXnr[4]+tpcXpr[4];
2569   float Yr4 = tpcYnr[4]+tpcYpr[4];
2570   float Xl0 = tpcXnl[0]+tpcXpl[0];
2571   float Xr0 = tpcXnr[0]+tpcXpr[0];
2572
2573   float Xl2o = tpcXnlo[2]+tpcXplo[2];
2574   float Yl2o = tpcYnlo[2]+tpcYplo[2];
2575   float Xr2o = tpcXnro[2]+tpcXpro[2];
2576   float Yr2o = tpcYnro[2]+tpcYpro[2];
2577   float Xl3o = tpcXnlo[3]+tpcXplo[3];
2578   float Yl3o = tpcYnlo[3]+tpcYplo[3];
2579   float Xr3o = tpcXnro[3]+tpcXpro[3];
2580   float Yr3o = tpcYnro[3]+tpcYpro[3];
2581   float Xl4o = tpcXnlo[4]+tpcXplo[4];
2582   float Yl4o = tpcYnlo[4]+tpcYplo[4];
2583   float Xr4o = tpcXnro[4]+tpcXpro[4];
2584   float Yr4o = tpcYnro[4]+tpcYpro[4];
2585   float Xl0o = tpcXnlo[0]+tpcXplo[0];
2586   float Xr0o = tpcXnro[0]+tpcXpro[0];
2587
2588   // --- q vector components for recentering if needed
2589   h_X22_cent->Fill(cent,(tpcX[2]/tpcX[0]));
2590   h_X22_centP->Fill(cent,(tpcXp[2]/tpcXp[0]));
2591   h_X22_centN->Fill(cent,(tpcXn[2]/tpcXn[0]));
2592   h_Y22_cent->Fill(cent,(tpcY[2]/tpcX[0]));
2593   h_Y22_centP->Fill(cent,(tpcYp[2]/tpcXp[0]));
2594   h_Y22_centN->Fill(cent,(tpcYn[2]/tpcXn[0]));
2595   // --- outer left
2596   h_X22_lo_cent->Fill(cent,(tpcX[2]/tpcX[0]));
2597   h_X22_lo_centP->Fill(cent,(tpcXplo[2]/tpcXplo[0]));
2598   h_X22_lo_centN->Fill(cent,(tpcXnlo[2]/tpcXnlo[0]));
2599   h_Y22_lo_cent->Fill(cent,(tpcY[2]/tpcX[0]));
2600   h_Y22_lo_centP->Fill(cent,(tpcYplo[2]/tpcXplo[0]));
2601   h_Y22_lo_centN->Fill(cent,(tpcYnlo[2]/tpcXnlo[0]));
2602   // --- inner left
2603   h_X22_li_cent->Fill(cent,(tpcX[2]/tpcX[0]));
2604   h_X22_li_centP->Fill(cent,(tpcXpli[2]/tpcXpli[0]));
2605   h_X22_li_centN->Fill(cent,(tpcXnli[2]/tpcXnli[0]));
2606   h_Y22_li_cent->Fill(cent,(tpcY[2]/tpcX[0]));
2607   h_Y22_li_centP->Fill(cent,(tpcYpli[2]/tpcXpli[0]));
2608   h_Y22_li_centN->Fill(cent,(tpcYnli[2]/tpcXnli[0]));
2609   // --- inner right
2610   h_X22_ri_cent->Fill(cent,(tpcX[2]/tpcX[0]));
2611   h_X22_ri_centP->Fill(cent,(tpcXpri[2]/tpcXpri[0]));
2612   h_X22_ri_centN->Fill(cent,(tpcXnri[2]/tpcXnri[0]));
2613   h_Y22_ri_cent->Fill(cent,(tpcY[2]/tpcX[0]));
2614   h_Y22_ri_centP->Fill(cent,(tpcYpri[2]/tpcXpri[0]));
2615   h_Y22_ri_centN->Fill(cent,(tpcYnri[2]/tpcXnri[0]));
2616   // --- outer right
2617   h_X22_ro_cent->Fill(cent,(tpcX[2]/tpcX[0]));
2618   h_X22_ro_centP->Fill(cent,(tpcXpro[2]/tpcXpro[0]));
2619   h_X22_ro_centN->Fill(cent,(tpcXnro[2]/tpcXnro[0]));
2620   h_Y22_ro_cent->Fill(cent,(tpcY[2]/tpcX[0]));
2621   h_Y22_ro_centP->Fill(cent,(tpcYpro[2]/tpcXpro[0]));
2622   h_Y22_ro_centN->Fill(cent,(tpcYnro[2]/tpcXnro[0]));
2623
2624   float q22ev = (tpcQQ[2]-M)/W_2;
2625   float q22Pev = (tpcQQp[2]-Mp)/Wp_2;
2626   float q22Nev = (tpcQQn[2]-Mn)/Wn_2;
2627   float q22gap0ev = (Xl2*Xr2+Yl2*Yr2)/(Xl0*Xr0);
2628   float q22gap0Pev = (tpcXpl[2]*tpcXpr[2]+tpcYpl[2]*tpcYpr[2])/(tpcXpl[0]*tpcXpr[0]);
2629   float q22gap0Nev = (tpcXnl[2]*tpcXnr[2]+tpcYnl[2]*tpcYnr[2])/(tpcXnl[0]*tpcXnr[0]);
2630   float q22gap1ev = (Xl2o*Xr2o+Yl2o*Yr2o)/(Xl0o*Xr0o);
2631   float q22gap1Pev = (tpcXplo[2]*tpcXpro[2]+tpcYplo[2]*tpcYpro[2])/(tpcXplo[0]*tpcXpro[0]);
2632   float q22gap1Nev = (tpcXnlo[2]*tpcXnro[2]+tpcYnlo[2]*tpcYnro[2])/(tpcXnlo[0]*tpcXnro[0]);
2633   float q22SLev = (Xl2*Xr2+Yl2*Yr2)/(Xl0*Xr0);
2634   float q22SLPev = (tpcXpl[2]*Xr2+tpcYpl[2]*Yr2)/(tpcXpl[0]*Xr0);
2635   float q22SLNev = (tpcXnl[2]*Xr2+tpcYnl[2]*Yr2)/(tpcXnl[0]*Xr0);
2636   h_q22_cent->Fill(cent,q22ev);
2637   h_q22_centP->Fill(cent,q22Pev);
2638   h_q22_centN->Fill(cent,q22Nev);
2639   h_q22gap0_cent->Fill(cent,q22gap0ev);
2640   h_q22gap0_centP->Fill(cent,q22gap0Pev);
2641   h_q22gap0_centN->Fill(cent,q22gap0Nev);
2642   h_q22gap1_cent->Fill(cent,q22gap1ev);
2643   h_q22gap1_centP->Fill(cent,q22gap1Pev);
2644   h_q22gap1_centN->Fill(cent,q22gap1Nev);
2645
2646   float q32ev = (tpcQQ[3]-M)/W_2;
2647   float q32Pev = (tpcQQp[3]-Mp)/Wp_2;
2648   float q32Nev = (tpcQQn[3]-Mn)/Wn_2;
2649   float q32gap0ev = (Xl3*Xr3+Yl3*Yr3)/(Xl0*Xr0);
2650   float q32gap0Pev = (tpcXpl[3]*tpcXpr[3]+tpcYpl[3]*tpcYpr[3])/(tpcXpl[0]*tpcXpr[0]);
2651   float q32gap0Nev = (tpcXnl[3]*tpcXnr[3]+tpcYnl[3]*tpcYnr[3])/(tpcXnl[0]*tpcXnr[0]);
2652   float q32gap1ev = (Xl3o*Xr3o+Yl3o*Yr3o)/(Xl0o*Xr0o);
2653   float q32gap1Pev = (tpcXplo[3]*tpcXpro[3]+tpcYplo[3]*tpcYpro[3])/(tpcXplo[0]*tpcXpro[0]);
2654   float q32gap1Nev = (tpcXnlo[3]*tpcXnro[3]+tpcYnlo[3]*tpcYnro[3])/(tpcXnlo[0]*tpcXnro[0]);
2655   float q32SLev = (Xl3*Xr3+Yl3*Yr3)/(Xl0*Xr0);
2656   float q32SLPev = (tpcXpl[3]*Xr3+tpcYpl[3]*Yr3)/(tpcXpl[0]*Xr0);
2657   float q32SLNev = (tpcXnl[3]*Xr3+tpcYnl[3]*Yr3)/(tpcXnl[0]*Xr0);
2658
2659   float q42ev = (tpcQQ[4]-M)/W_2;
2660   float q42Pev = (tpcQQp[4]-Mp)/Wp_2;
2661   float q42Nev = (tpcQQn[4]-Mn)/Wn_2;
2662   float q42gap0ev = (Xl4*Xr4+Yl4*Yr4)/(Xl0*Xr0);
2663   float q42gap0Pev = (tpcXpl[4]*tpcXpr[4]+tpcYpl[4]*tpcYpr[4])/(tpcXpl[0]*tpcXpr[0]);
2664   float q42gap0Nev = (tpcXnl[4]*tpcXnr[4]+tpcYnl[4]*tpcYnr[4])/(tpcXnl[0]*tpcXnr[0]);
2665   float q42gap1ev = (Xl4o*Xr4o+Yl4o*Yr4o)/(Xl0o*Xr0o);
2666   float q42gap1Pev = (tpcXplo[4]*tpcXpro[4]+tpcYplo[4]*tpcYpro[4])/(tpcXplo[0]*tpcXpro[0]);
2667   float q42gap1Nev = (tpcXnlo[4]*tpcXnro[4]+tpcYnlo[4]*tpcYnro[4])/(tpcXnlo[0]*tpcXnro[0]);
2668   float q42SLev = (Xl4*Xr4+Yl4*Yr4)/(Xl0*Xr0);
2669   float q42SLPev = (tpcXpl[4]*Xr4+tpcYpl[4]*Yr4)/(tpcXpl[0]*Xr0);
2670   float q42SLNev = (tpcXnl[4]*Xr4+tpcYnl[4]*Yr4)/(tpcXnl[0]*Xr0);
2671
2672
2673   // 3p correlator
2674   float q23ev = calc3(tpcX[1],tpcY[1],tpcX[2],tpcY[2],M);
2675   float q23Pev = calc3(tpcXp[1],tpcYp[1],tpcXp[2],tpcYp[2],Mp);
2676   float q23Nev = calc3(tpcXn[1],tpcYn[1],tpcXn[2],tpcYn[2],Mn);
2677   // 4p correlator
2678   float q24ev = calc4(tpcX[2],tpcY[2],tpcX[4],tpcY[4],M);
2679   float q24Pev = calc4(tpcXp[2],tpcYp[2],tpcXp[4],tpcYp[4],Mp);
2680   float q24Nev = calc4(tpcXn[2],tpcYn[2],tpcXn[4],tpcYn[4],Mn);
2681
2682
2683
2684   h_q22_cent->Fill(cent,q22ev);
2685   h_q22_centP->Fill(cent,q22Pev);
2686   h_q22_centN->Fill(cent,q22Nev);
2687   h_q23_cent->Fill(cent,q23ev);
2688   h_q23_centP->Fill(cent,q23Pev);
2689   h_q23_centN->Fill(cent,q23Nev);
2690   h_q24_cent->Fill(cent,q24ev);
2691   h_q24_centP->Fill(cent,q24Pev);
2692   h_q24_centN->Fill(cent,q24Nev);
2693   h_q22gap0_cent->Fill(cent,q22gap0ev);
2694   h_q22gap0_centP->Fill(cent,q22gap0Pev);
2695   h_q22gap0_centN->Fill(cent,q22gap0Nev);
2696   h_q22gap1_cent->Fill(cent,q22gap1ev);
2697   h_q22gap1_centP->Fill(cent,q22gap1Pev);
2698   h_q22gap1_centN->Fill(cent,q22gap1Nev);
2699
2700   h_q32_cent->Fill(cent,q32ev);
2701   h_q32_centP->Fill(cent,q32Pev);
2702   h_q32_centN->Fill(cent,q32Nev);
2703   h_q32gap0_cent->Fill(cent,q32gap0ev);
2704   h_q32gap0_centP->Fill(cent,q32gap0Pev);
2705   h_q32gap0_centN->Fill(cent,q32gap0Nev);
2706   h_q32gap1_cent->Fill(cent,q32gap1ev);
2707   h_q32gap1_centP->Fill(cent,q32gap1Pev);
2708   h_q32gap1_centN->Fill(cent,q32gap1Nev);
2709
2710   h_q42_cent->Fill(cent,q42ev);
2711   h_q42_centP->Fill(cent,q42Pev);
2712   h_q42_centN->Fill(cent,q42Nev);
2713   h_q42gap0_cent->Fill(cent,q42gap0ev);
2714   h_q42gap0_centP->Fill(cent,q42gap0Pev);
2715   h_q42gap0_centN->Fill(cent,q42gap0Nev);
2716   h_q42gap1_cent->Fill(cent,q42gap1ev);
2717   h_q42gap1_centP->Fill(cent,q42gap1Pev);
2718   h_q42gap1_centN->Fill(cent,q42gap1Nev);
2719
2720   if(fabs(qasymm)<10.0&&fabs(qasymmL)<10.0&&fabs(qasymmR)<10.0)
2721     {
2722       h_q22qasym_cent[icent]->Fill(qasymm,q22ev);
2723       h_q22qasymP_cent[icent]->Fill(qasymm,q22Pev);
2724       h_q22qasymN_cent[icent]->Fill(qasymm,q22Nev);
2725       h_q22qasym_gap0_cent[icent]->Fill(qasymm,q22gap0ev);
2726       h_q22qasymP_gap0_cent[icent]->Fill(qasymm,q22gap0Pev);
2727       h_q22qasymN_gap0_cent[icent]->Fill(qasymm,q22gap0Nev);
2728       h_q22qasym_gap1_cent[icent]->Fill(qasymm,q22gap1ev);
2729       h_q22qasymP_gap1_cent[icent]->Fill(qasymm,q22gap1Pev);
2730       h_q22qasymN_gap1_cent[icent]->Fill(qasymm,q22gap1Nev);
2731
2732       h_q32qasym_cent[icent]->Fill(qasymm,q32ev);
2733       h_q32qasymP_cent[icent]->Fill(qasymm,q32Pev);
2734       h_q32qasymN_cent[icent]->Fill(qasymm,q32Nev);
2735       h_q32qasym_gap0_cent[icent]->Fill(qasymm,q32gap0ev);
2736       h_q32qasymP_gap0_cent[icent]->Fill(qasymm,q32gap0Pev);
2737       h_q32qasymN_gap0_cent[icent]->Fill(qasymm,q32gap0Nev);
2738       h_q32qasym_gap1_cent[icent]->Fill(qasymm,q32gap1ev);
2739       h_q32qasymP_gap1_cent[icent]->Fill(qasymm,q32gap1Pev);
2740       h_q32qasymN_gap1_cent[icent]->Fill(qasymm,q32gap1Nev);
2741
2742       h_q42qasym_cent[icent]->Fill(qasymm,q42ev);
2743       h_q42qasymP_cent[icent]->Fill(qasymm,q42Pev);
2744       h_q42qasymN_cent[icent]->Fill(qasymm,q42Nev);
2745       h_q42qasym_gap0_cent[icent]->Fill(qasymm,q42gap0ev);
2746       h_q42qasymP_gap0_cent[icent]->Fill(qasymm,q42gap0Pev);
2747       h_q42qasymN_gap0_cent[icent]->Fill(qasymm,q42gap0Nev);
2748       h_q42qasym_gap1_cent[icent]->Fill(qasymm,q42gap1ev);
2749       h_q42qasymP_gap1_cent[icent]->Fill(qasymm,q42gap1Pev);
2750       h_q42qasymN_gap1_cent[icent]->Fill(qasymm,q42gap1Nev);
2751     }
2752
2753   h_Aq22_cent->Fill(cent,q22ev*qasymm);
2754   h_Aq22_centP->Fill(cent,q22Pev*qasymm);
2755   h_Aq22_centN->Fill(cent,q22Nev*qasymm);
2756   h_Aq22_SL_cent->Fill(cent,q22SLev*qasymm);
2757   h_Aq22_SL_centP->Fill(cent,q22SLPev*qasymm);
2758   h_Aq22_SL_centN->Fill(cent,q22SLNev*qasymm);
2759   h_Aq22_SLL_cent->Fill(cent,q22SLev*qasymmL);
2760   h_Aq22_SLL_centP->Fill(cent,q22SLPev*qasymmL);
2761   h_Aq22_SLL_centN->Fill(cent,q22SLNev*qasymmL);
2762   h_Aq22_SLR_cent->Fill(cent,q22SLev*qasymmR);
2763   h_Aq22_SLR_centP->Fill(cent,q22SLPev*qasymmR);
2764   h_Aq22_SLR_centN->Fill(cent,q22SLNev*qasymmR);
2765
2766   h_Aq32_cent->Fill(cent,q32ev*qasymm);
2767   h_Aq32_centP->Fill(cent,q32Pev*qasymm);
2768   h_Aq32_centN->Fill(cent,q32Nev*qasymm);
2769   h_Aq32_SL_cent->Fill(cent,q32SLev*qasymm);
2770   h_Aq32_SL_centP->Fill(cent,q32SLPev*qasymm);
2771   h_Aq32_SL_centN->Fill(cent,q32SLNev*qasymm);
2772   h_Aq32_SLL_cent->Fill(cent,q32SLev*qasymmL);
2773   h_Aq32_SLL_centP->Fill(cent,q32SLPev*qasymmL);
2774   h_Aq32_SLL_centN->Fill(cent,q32SLNev*qasymmL);
2775   h_Aq32_SLR_cent->Fill(cent,q32SLev*qasymmR);
2776   h_Aq32_SLR_centP->Fill(cent,q32SLPev*qasymmR);
2777   h_Aq32_SLR_centN->Fill(cent,q32SLNev*qasymmR);
2778
2779   h_Aq42_cent->Fill(cent,q42ev*qasymm);
2780   h_Aq42_centP->Fill(cent,q42Pev*qasymm);
2781   h_Aq42_centN->Fill(cent,q42Nev*qasymm);
2782   h_Aq42_SL_cent->Fill(cent,q42SLev*qasymm);
2783   h_Aq42_SL_centP->Fill(cent,q42SLPev*qasymm);
2784   h_Aq42_SL_centN->Fill(cent,q42SLNev*qasymm);
2785   h_Aq42_SLL_cent->Fill(cent,q42SLev*qasymmL);
2786   h_Aq42_SLL_centP->Fill(cent,q42SLPev*qasymmL);
2787   h_Aq42_SLL_centN->Fill(cent,q42SLNev*qasymmL);
2788   h_Aq42_SLR_cent->Fill(cent,q42SLev*qasymmR);
2789   h_Aq42_SLR_centP->Fill(cent,q42SLPev*qasymmR);
2790   h_Aq42_SLR_centN->Fill(cent,q42SLNev*qasymmR);
2791
2792
2793
2794   float a1Xp[hmax], a1Yp[hmax], a2Xp[hmax], a2Yp[hmax];
2795   float a1Xn[hmax], a1Yn[hmax], a2Xn[hmax], a2Yn[hmax];
2796
2797   for(int i=0; i<hmax;i++)
2798     {
2799       a1Xp[i] = 0.0;
2800       a1Yp[i] = 0.0;
2801       a1Xn[i] = 0.0;
2802       a1Yn[i] = 0.0;
2803       a2Xn[i] = 0.0;
2804       a2Yn[i] = 0.0;
2805       a2Xp[i] = 0.0;
2806       a2Yp[i] = 0.0;
2807     }
2808
2809   // --- now second track loop
2810
2811   for(int itrk = 0; itrk<d_ntrk; itrk++)
2812     {
2813       AliAODTrack *track = fAOD->GetTrack(itrk);
2814       if(!track)
2815         {
2816           if(debug>0) cout<<"ERROR: Could not retrieve AODtrack "<<itrk<<endl;
2817           continue;
2818         }
2819
2820       if(!track->TestFilterBit(fbit))
2821         {
2822           if(debug>15) cout<<"second track loop: wrong filter bit for track "<<itrk<<endl;
2823           continue;
2824         }
2825
2826       // int gid = track->GetID();
2827       // AliAODTrack *PIDtrack;
2828       // if(gid>=0) PIDtrack = track;
2829       // else PIDtrack = fAOD->GetTrack(trackMap->GetValue(-1-gid));
2830
2831       // --- get track variables
2832       float pt = track->Pt();
2833       float phi = track->Phi();
2834       float eta = track->Eta();
2835       int charge = track->Charge();
2836       int ncls = track->GetTPCNcls();
2837       //float dedx = track->GetTPCsignal();
2838
2839       float Tdcaxy = track->DCA();
2840       float Tdcaz = track->ZAtDCA();
2841
2842
2843       float dcaxy = -999;
2844       float dcax = -999;
2845       float dcay = -999;
2846       float dcaz = -999;
2847
2848       double r[3];
2849       bool dcaflag = track->GetXYZ(r);
2850
2851       double DCA[2]; // dca
2852       double COV[3]; // covariance
2853       bool proptodca = track->PropagateToDCA(fVtx,mag,100.0,DCA,COV);
2854       if(!proptodca&&debug>0) cout<<"No DCACOV for you!"<<endl;
2855
2856       if(dcaflag)
2857         {
2858           dcaxy = r[0];
2859           dcaz = r[1];
2860           //if(debug>6) cout<<"GetXYZ is true, filter bit is "<<fbit<<endl;
2861         }
2862       else
2863         {
2864           dcax = r[0] - eventX;
2865           dcay = r[1] - eventY;
2866           dcaz = r[2] - eventZ;
2867           // --- need the sign convention for dcaxy...
2868           dcaxy = sqrt(dcax*dcax+dcay*dcay);
2869           //if((float)dcaxy!=(float)fabs((float)DCA[0])&&debug>6) cout<<"hmm... "<<dcaxy<<" "<<DCA[0]<<endl;
2870           // --- set dcaxy to value from PropagateToDCA to get correct sign
2871           dcaxy = DCA[0];
2872           // --- dcaz on the other hand is unambiguous
2873           //if(dcaz!=(float)DCA[1]&&debug>6) cout<<"hmm... "<<dcaz<<" "<<DCA[1]<<endl;
2874           //if(debug>6) cout<<"GetXYZ is false, filter bit is "<<fbit<<endl;
2875         }
2876
2877       // if(debug>6)
2878       //        {
2879       //          cout<<"r[0] is "<<r[0]<<" ";
2880       //          cout<<"r[1] is "<<r[1]<<" ";
2881       //          cout<<"r[2] is "<<r[2]<<endl;
2882       //          cout<<"eventX is "<<eventX<<" ";
2883       //          cout<<"eventY is "<<eventY<<" ";
2884       //          cout<<"eventZ is "<<eventZ<<endl;
2885       //          cout<<"Tdcaxy is "<<Tdcaxy<<" and dcaxy is "<<dcaxy<<" and DCACOV xy is "<<DCA[0]<<endl;
2886       //          cout<<"Tdcaz is "<<Tdcaz<<" and dcaz is "<<dcaz<<" and DCACOV z is "<<DCA[1]<<endl;
2887       //        }
2888
2889       if((fbit==128||fbit==272)&&!dcaflag)
2890         {
2891           dcaxy = Tdcaxy;
2892           dcaz = Tdcaz;
2893           // if(debug>4)
2894           //   {
2895           //     cout<<"GetXYZ false but should be true, flipping values"<<endl;
2896           //     cout<<"Tdcaxy is "<<Tdcaxy<<" and dcaxy is "<<dcaxy<<" and DCACOV xy is "<<DCA[0]<<endl;
2897           //     cout<<"Tdcaz is "<<Tdcaz<<" and dcaz is "<<dcaz<<" and DCACOV z is "<<DCA[1]<<endl;
2898           //   }
2899         }
2900
2901       // --- track cuts
2902       if(ncls<nclscut) continue;
2903       if(pt<ptmin||pt>ptmax) continue;
2904       if(fabs(eta)>outeta||fabs(eta)<excleta) continue;
2905       if(dodcacuts&&(fabs(dcaxy)>dcacutxy||fabs(dcaz)>dcacutz)) continue;// problems depending on filter bit
2906
2907       bool pos = charge>0.0;
2908       bool neg = charge<0.0;
2909
2910
2911       if(gRandom) delete gRandom;
2912       gRandom = new TRandom3(0);
2913       gRandom->SetSeed(0);
2914       double rand = gRandom->Rndm();
2915       bool sub1 = false;
2916       bool sub2 = false;
2917       if(rand<0.5) sub1 = true;
2918       else sub2 = true;
2919
2920       if(sub1)
2921         {
2922           ntrkA1++;
2923           if(charge>0) ntrkposA1++;
2924           if(charge<0) ntrknegA1++;
2925         }
2926       if(sub2)
2927         {
2928           ntrkA2++;
2929           if(charge>0) ntrkposA2++;
2930           if(charge<0) ntrknegA2++;
2931         }
2932
2933       if(sub1) // random selection
2934         {
2935           for(int i=1; i<hmax;i++)
2936             {
2937               if(charge>0)
2938                 {
2939                   a1Xp[i] += cos(i*phi);
2940                   a1Yp[i] += sin(i*phi);
2941                 }
2942               if(charge<0)
2943                 {
2944                   a1Xn[i] += cos(i*phi);
2945                   a1Yn[i] += sin(i*phi);
2946                 }
2947             }
2948           // assign values to array element 0
2949           if(charge>0)
2950             {
2951               a1Xp[0] += 1.0;
2952               a1Yp[0] += pt;
2953             }
2954           if(charge<0)
2955             {
2956               a1Xn[0] += 1.0;
2957               a1Yn[0] += pt;
2958             }
2959         }
2960
2961       if(sub2) // random selection
2962         {
2963           for(int i=1; i<hmax;i++)
2964             {
2965               if(charge>0)
2966                 {
2967                   a2Xp[i] += cos(i*phi);
2968                   a2Yp[i] += sin(i*phi);
2969                 }
2970               if(charge<0)
2971                 {
2972                   a2Xn[i] += cos(i*phi);
2973                   a2Yn[i] += sin(i*phi);
2974                 }
2975             }
2976           // assign values to array element 0
2977           if(charge>0)
2978             {
2979               a2Xp[0] += 1.0;
2980               a2Yp[0] += pt;
2981             }
2982           if(charge<0)
2983             {
2984               a2Xn[0] += 1.0;
2985               a2Yn[0] += pt;
2986             }
2987         }
2988
2989
2990
2991       // -------------------------------------------------------
2992       // --- now calculate differential cumulants for each track
2993       // -------------------------------------------------------
2994       float trkXpl[hmax], trkYpl[hmax], trkXpr[hmax], trkYpr[hmax];
2995       float trkXnl[hmax], trkYnl[hmax], trkXnr[hmax], trkYnr[hmax];
2996       for(int i=0; i<hmax;i++)
2997         {
2998           trkXpl[i] = 0.0;
2999           trkYpl[i] = 0.0;
3000           trkXpr[i] = 0.0;
3001           trkYpr[i] = 0.0;
3002           trkXnl[i] = 0.0;
3003           trkYnl[i] = 0.0;
3004           trkXnr[i] = 0.0;
3005           trkYnr[i] = 0.0;
3006         }
3007       // very similar to above, except += changed to =
3008       if(pt>ptmin&&pt<ptmax)
3009         {
3010           if(eta>-outeta&&eta<-excleta)
3011             {
3012               // assign values to array elements 1..8
3013               for(int i=1; i<hmax;i++)
3014                 {
3015                   if(charge>0)
3016                     {
3017                       trkXpl[i] = cos(i*phi);
3018                       trkYpl[i] = sin(i*phi);
3019                     }
3020                   if(charge<0)
3021                     {
3022                       trkXnl[i] = cos(i*phi);
3023                       trkYnl[i] = sin(i*phi);
3024                     }
3025                 }
3026               // assign values to array element 0
3027               if(charge>0)
3028                 {
3029                   trkXpl[0] = 1.0;
3030                   trkYpl[0] = pt;
3031                 }
3032               if(charge<0)
3033                 {
3034                   trkXnl[0] = 1.0;
3035                   trkYnl[0] = pt;
3036                 }
3037             }
3038           //----------------------------
3039           if(eta>excleta&&eta<outeta) // right half of TPC
3040             {
3041               for(int i=1; i<hmax; i++)
3042                 {
3043                   if(charge>0)
3044                     {
3045                       trkXpr[i] = cos(i*phi);
3046                       trkYpr[i] = sin(i*phi);
3047                     }
3048                   if(charge<0)
3049                     {
3050                       trkXnr[i] = cos(i*phi);
3051                       trkYnr[i] = sin(i*phi);
3052                     }
3053                 }
3054               if(charge>0)
3055                 {
3056                   trkXpr[0] = 1.0;
3057                   trkYpr[0] = pt;
3058                 }
3059               if(charge<0)
3060                 {
3061                   trkXnr[0] = 1.0;
3062                   trkYnr[0] = pt;
3063                 }
3064             } // end right half of TPC
3065         } // end pT selection
3066
3067       // COME BACK HERE
3068       float trkX[9],  trkY[9];
3069       float trkXp[9], trkYp[9]; // pos
3070       float trkXn[9], trkYn[9]; // neg
3071       for(int i=0; i<6; i++)
3072         {
3073           trkX[i]=trkXpl[i]+trkXnl[i]+trkXpr[i]+trkXnr[i];
3074           trkY[i]=trkYpl[i]+trkYnl[i]+trkYpr[i]+trkYnr[i];
3075           // pos
3076           trkXp[i]=trkXpl[i]+trkXpr[i];
3077           trkYp[i]=trkYpl[i]+trkYpr[i];
3078           // neg
3079           trkXn[i]=trkXnl[i]+trkXnr[i];
3080           trkYn[i]=trkYnl[i]+trkYnr[i];
3081         }
3082
3083       // sanity checks to prevent division by zero and other problems
3084       if(trkX[0]<1) continue;
3085       if(tpcX[0]<2) continue;
3086       if(tpcXp[0]<2) continue;
3087       if(tpcXn[0]<2) continue;
3088
3089       // --- calculate differential cumulants from Q-vector components
3090       // --- reference flow is from whole TPC
3091       float diffq22ev = ((trkX[2]*tpcX[2])+(trkY[2]*tpcY[2])-1)/(tpcX[0]-1);
3092       float diffq22Pev = ((trkXp[2]*tpcX[2])+(trkYp[2]*tpcY[2])-1)/(tpcX[0]-1);
3093       float diffq22Nev = ((trkXn[2]*tpcX[2])+(trkYn[2]*tpcY[2])-1)/(tpcX[0]-1);
3094       // --- third harmonic
3095       float diffq32ev = ((trkX[3]*tpcX[3])+(trkY[3]*tpcY[3])-1)/(tpcX[0]-1);
3096       float diffq32Pev = ((trkXp[3]*tpcX[3])+(trkYp[3]*tpcY[3])-1)/(tpcX[0]-1);
3097       float diffq32Nev = ((trkXn[3]*tpcX[3])+(trkYn[3]*tpcY[3])-1)/(tpcX[0]-1);
3098       // --- fourth harmonic
3099       float diffq42ev = ((trkX[4]*tpcX[4])+(trkY[4]*tpcY[4])-1)/(tpcX[0]-1);
3100       float diffq42Pev = ((trkXp[4]*tpcX[4])+(trkYp[4]*tpcY[4])-1)/(tpcX[0]-1);
3101       float diffq42Nev = ((trkXn[4]*tpcX[4])+(trkYn[4]*tpcY[4])-1)/(tpcX[0]-1);
3102
3103       // --- differential cumulants test with differing RPs
3104       // float Xdiffq22Pev = ((trkXp[2]*tpcX[2])+(trkYp[2]*tpcY[2])-1)/(tpcX[0]-1);
3105       // float Xdiffq22Nev = ((trkXn[2]*tpcX[2])+(trkYn[2]*tpcY[2])-1)/(tpcX[0]-1);
3106       // float Ydiffq22Pev = ((trkXp[2]*tpcXp[2])+(trkYp[2]*tpcYp[2])-1)/(tpcXp[0]-1);
3107       // float Ydiffq22Nev = ((trkXn[2]*tpcXn[2])+(trkYn[2]*tpcYn[2])-1)/(tpcXn[0]-1);
3108       // float Zdiffq22Pev = ((trkXp[2]*tpcXn[2])+(trkYp[2]*tpcYn[2]))/(tpcXn[0]);
3109       // float Zdiffq22Nev = ((trkXn[2]*tpcXp[2])+(trkYn[2]*tpcYp[2]))/(tpcXp[0]);
3110
3111       // --- histograms to check different reference flow
3112       // if(pos) h_Xdiffq22_ptP->Fill(pt,Xdiffq22Pev);
3113       // if(neg) h_Xdiffq22_ptN->Fill(pt,Xdiffq22Nev);
3114       // if(pos) h_Ydiffq22_ptP->Fill(pt,Ydiffq22Pev);
3115       // if(neg) h_Ydiffq22_ptN->Fill(pt,Ydiffq22Nev);
3116       // if(pos) h_Zdiffq22_ptP->Fill(pt,Zdiffq22Pev);
3117       // if(neg) h_Zdiffq22_ptN->Fill(pt,Zdiffq22Nev);
3118
3119       if(cent>=centlo&&cent<centhi)
3120         {
3121           // --- differential cumulant
3122           h_diffq22_pt->Fill(pt,diffq22ev);
3123           h_diffq22_eta->Fill(eta,diffq22ev);
3124           // --- differential cumulant weighted with qasymm
3125           h_diffAq22_pt->Fill(pt,diffq22ev*qasymm);
3126           h_diffAq22_eta->Fill(eta,diffq22ev*qasymm);
3127           if(pos)
3128             {
3129               h_diffq22_ptP->Fill(pt,diffq22Pev);
3130               h_diffq22_etaP->Fill(eta,diffq22Pev);
3131               h_diffAq22_ptP->Fill(pt,diffq22Pev*qasymm);
3132               h_diffAq22_etaP->Fill(eta,diffq22Pev*qasymm);
3133             }
3134           if(neg)
3135             {
3136               h_diffq22_ptN->Fill(pt,diffq22Nev);
3137               h_diffq22_etaN->Fill(eta,diffq22Nev);
3138               h_diffAq22_ptN->Fill(pt,diffq22Nev*qasymm);
3139               h_diffAq22_etaN->Fill(eta,diffq22Nev*qasymm);
3140             }
3141           // --- same for third harmonic
3142           h_diffq32_pt->Fill(pt,diffq32ev);
3143           h_diffq32_eta->Fill(eta,diffq32ev);
3144           h_diffAq32_pt->Fill(pt,diffq32ev*qasymm);
3145           h_diffAq32_eta->Fill(eta,diffq32ev*qasymm);
3146           if(pos)
3147             {
3148               h_diffq32_ptP->Fill(pt,diffq32Pev);
3149               h_diffq32_etaP->Fill(eta,diffq32Pev);
3150               h_diffAq32_ptP->Fill(pt,diffq32Pev*qasymm);
3151               h_diffAq32_etaP->Fill(eta,diffq32Pev*qasymm);
3152             }
3153           if(neg)
3154             {
3155               h_diffq32_ptN->Fill(pt,diffq32Nev);
3156               h_diffq32_etaN->Fill(eta,diffq32Nev);
3157               h_diffAq32_ptN->Fill(pt,diffq32Nev*qasymm);
3158               h_diffAq32_etaN->Fill(eta,diffq32Nev*qasymm);
3159             }
3160           // --- same for fourth harmonic
3161           h_diffq42_pt->Fill(pt,diffq42ev);
3162           h_diffq42_eta->Fill(eta,diffq42ev);
3163           h_diffAq42_pt->Fill(pt,diffq42ev*qasymm);
3164           h_diffAq42_eta->Fill(eta,diffq42ev*qasymm);
3165           if(pos)
3166             {
3167               h_diffq42_ptP->Fill(pt,diffq42Pev);
3168               h_diffq42_etaP->Fill(eta,diffq42Pev);
3169               h_diffAq42_ptP->Fill(pt,diffq42Pev*qasymm);
3170               h_diffAq42_etaP->Fill(eta,diffq42Pev*qasymm);
3171             }
3172           if(neg)
3173             {
3174               h_diffq42_ptN->Fill(pt,diffq42Nev);
3175               h_diffq42_etaN->Fill(eta,diffq42Nev);
3176               h_diffAq42_ptN->Fill(pt,diffq42Nev*qasymm);
3177               h_diffAq42_etaN->Fill(eta,diffq42Nev*qasymm);
3178             }
3179
3180           // selection to do nested track loop
3181           if(donested)
3182             {
3183               for(int i=0; i<d_ntrk; i++)
3184                 {
3185                   // ----------------------------------------------------------------------------------------
3186                   // --- want delta eta and delta pt dependence of differential cumulant weighted with charge
3187                   // --- need nested track loop to get eta1 and eta3
3188                   // ----------------------------------------------------------------------------------------
3189                   // make sure eta is in range, -99 should be autorejected by this
3190                   //if(eta3[i]==-99) continue;
3191                   if(eta3[i]<-outeta||eta3[i]>outeta)
3192                     {
3193                       //cout<<"eta3 out of range!!! "<<eta3[i]<<endl;
3194                       continue;
3195                     }
3196                   if(eta<-outeta||eta>outeta)
3197                     {
3198                       cout<<"eta out of range!!! "<<eta<<endl;
3199                       continue;
3200                     }
3201                   // charge==0 should be rejected by the above cut
3202                   if(charge3[i]==0) {cout<<"WTF!"<<endl;continue;}
3203                   // VERY IMPORTANT remove auto-correlations with 3rd particle
3204                   if(i==itrk) continue;
3205
3206                   float DETA = eta-eta3[i];
3207                   float DPHI = phi-phi3[i];
3208                   float DPT = pt-pt3[i];
3209
3210                   // magnetic field sign
3211                   //bSign = (eventMain->GetMagneticField() > 0) ? 1 : -1;
3212                   float bSign = 1;
3213                   if(mag<0) bSign = -1;
3214
3215                   // double deta = firstEta - secondEta[j];
3216                   // double dphi = firstPhi - secondPhi[j];
3217                   float deta = DETA;
3218                   float dphi = DPHI;
3219                   float fHBTCutValue = 0.02;
3220                   // optimization
3221                   if(dopaircut&&(fabs(deta)<fHBTCutValue*2.5*3)) //fHBTCutValue = 0.02 [default for dphicorrelations]
3222                     {
3223                       if(debug>4)
3224                         {
3225                           cout<<"inside deta cut, evaluating"<<endl;
3226                           cout<<"deta is "<<deta<<" and dphi is "<<dphi<<endl;
3227                         }
3228                       // phi in rad
3229                       // float phi1rad = firstPhi;
3230                       // float phi2rad = secondPhi[j];
3231
3232                       // check first boundaries to see if is worth to loop and find the minimum
3233                       // float dphistar1 = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, 0.8, bSign);
3234                       // float dphistar2 = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, 2.5, bSign);
3235                       float dphistar1 = GetDPhiStar(phi, pt, charge, phi3[i], pt3[i], charge3[i], 0.8, bSign);
3236                       float dphistar2 = GetDPhiStar(phi, pt, charge, phi3[i], pt3[i], charge3[i], 2.5, bSign);
3237
3238                       const float kLimit = fHBTCutValue * 3;
3239
3240                       float dphistarminabs = 1e5;
3241                       float dphistarmin = 1e5;
3242
3243                       if(fabs(dphistar1) < kLimit || fabs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0 )
3244                         {
3245                           if(debug>4)
3246                             {
3247                               cout<<"inside detadphi cut, evaluating"<<endl;
3248                               cout<<"deta is "<<deta<<" and dphi is "<<dphi<<endl;
3249                             }
3250                           for(double rad=0.8; rad<2.51; rad+=0.01)
3251                             {
3252                               //float dphistar = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, rad, bSign);
3253                               float dphistar = GetDPhiStar(phi, pt, charge, phi3[i], pt3[i], charge3[i], 0.8, bSign);
3254                               float dphistarabs = fabs(dphistar);
3255                               if(dphistarabs < dphistarminabs)
3256                                 {
3257                                   dphistarmin = dphistar;
3258                                   dphistarminabs = dphistarabs;
3259                                 }
3260                             }
3261                           if(dphistarminabs < fHBTCutValue && fabs(deta) < fHBTCutValue)
3262                             {
3263                               if(debug>4)
3264                                 printf("HBT: Removed track pair %d %d with [[%f %f]] %f %f %f | %f %f %d %f %f %d %f",
3265                                        itrk, i, deta, dphi, dphistarminabs, dphistar1, dphistar2, phi, pt, charge, phi3[i], pt3[i], charge3[i], bSign);
3266                               continue;
3267                             }
3268                         } // loop over angles
3269                     } // selection on pairs close together
3270
3271
3272
3273
3274
3275                   h_AT_X_deta->Fill(DETA,charge3[i]);
3276                   if(pos) h_AT_X_detaP->Fill(DETA,charge3[i]);
3277                   if(neg) h_AT_X_detaN->Fill(DETA,charge3[i]);
3278                   h_AT_X_dpt->Fill(DPT,charge3[i]);
3279                   if(pos) h_AT_X_dptP->Fill(DPT,charge3[i]);
3280                   if(neg) h_AT_X_dptN->Fill(DPT,charge3[i]);
3281                   // --- particles 1, 2, and 3 in whole TPC
3282                   if(pos) h_diffq22_X_detaP->Fill(DETA,diffq22Pev);
3283                   if(neg) h_diffq22_X_detaN->Fill(DETA,diffq22Nev);
3284                   if(pos) h_diffq32_X_detaP->Fill(DETA,diffq32Pev);
3285                   if(neg) h_diffq32_X_detaN->Fill(DETA,diffq32Nev);
3286                   if(pos) h_diffq42_X_detaP->Fill(DETA,diffq42Pev);
3287                   if(neg) h_diffq42_X_detaN->Fill(DETA,diffq42Nev);
3288                   if(pos) h_ATdiffq22_X_detaP->Fill(DETA,diffq22Pev*charge3[i]);
3289                   if(neg) h_ATdiffq22_X_detaN->Fill(DETA,diffq22Nev*charge3[i]);
3290                   if(pos) h_ATdiffq32_X_detaP->Fill(DETA,diffq32Pev*charge3[i]);
3291                   if(neg) h_ATdiffq32_X_detaN->Fill(DETA,diffq32Nev*charge3[i]);
3292                   if(pos) h_ATdiffq42_X_detaP->Fill(DETA,diffq42Pev*charge3[i]);
3293                   if(neg) h_ATdiffq42_X_detaN->Fill(DETA,diffq42Nev*charge3[i]);
3294                   if(pos) h_diffq22_X_dptP->Fill(DPT,diffq22Pev);
3295                   if(neg) h_diffq22_X_dptN->Fill(DPT,diffq22Nev);
3296                   if(pos) h_diffq32_X_dptP->Fill(DPT,diffq32Pev);
3297                   if(neg) h_diffq32_X_dptN->Fill(DPT,diffq32Nev);
3298                   if(pos) h_diffq42_X_dptP->Fill(DPT,diffq42Pev);
3299                   if(neg) h_diffq42_X_dptN->Fill(DPT,diffq42Nev);
3300                   if(pos) h_ATdiffq22_X_dptP->Fill(DPT,diffq22Pev*charge3[i]);
3301                   if(neg) h_ATdiffq22_X_dptN->Fill(DPT,diffq22Nev*charge3[i]);
3302                   if(pos) h_ATdiffq32_X_dptP->Fill(DPT,diffq32Pev*charge3[i]);
3303                   if(neg) h_ATdiffq32_X_dptN->Fill(DPT,diffq32Nev*charge3[i]);
3304                   if(pos) h_ATdiffq42_X_dptP->Fill(DPT,diffq42Pev*charge3[i]);
3305                   if(neg) h_ATdiffq42_X_dptN->Fill(DPT,diffq42Nev*charge3[i]);
3306
3307                   float meanchargeeta = 0;
3308                   int index = int((eta3[i]+1.0)*50);
3309                   if(index<0||index>100)
3310                     {
3311                       cout<<"out of index!!!"<<endl;
3312                       continue; // safety, should be redundant
3313                     }
3314                   //cout<<"eta is "<<eta3[i]<<" index is "<<index<<endl;
3315                   if(mag<-4.0)
3316                     {
3317                       if(fbit==1) meanchargeeta = MCEF1_fb1[index];
3318                       if(fbit==128) meanchargeeta = MCEF1_fb128[index];
3319                       if(fbit==272) meanchargeeta = MCEF1_fb272[index];
3320                       //if(fbit==1&&dodcacuts) meanchargeeta = MCEF1_dca[index];
3321                     }
3322                   if(mag>4.0)
3323                     {
3324                       if(fbit==1) meanchargeeta = MCEF3_fb1[index];
3325                       if(fbit==128) meanchargeeta = MCEF3_fb128[index];
3326                       if(fbit==272) meanchargeeta = MCEF3_fb272[index];
3327                       //if(fbit==1&&dodcacuts) meanchargeeta = MCEF3_dca[index];
3328                     }
3329                   float charge3subeta = charge3[i] - meanchargeeta;
3330
3331                   h_AT_S_deta->Fill(DETA,charge3subeta);
3332                   if(pos) h_AT_S_detaP->Fill(DETA,charge3subeta);
3333                   if(neg) h_AT_S_detaN->Fill(DETA,charge3subeta);
3334                   // --- particles 1, 2, and 3 in whole TPC
3335                   if(pos) h_diffq22_S_detaP->Fill(DETA,diffq22Pev);
3336                   if(neg) h_diffq22_S_detaN->Fill(DETA,diffq22Nev);
3337                   if(pos) h_diffq32_S_detaP->Fill(DETA,diffq32Pev);
3338                   if(neg) h_diffq32_S_detaN->Fill(DETA,diffq32Nev);
3339                   if(pos) h_diffq42_S_detaP->Fill(DETA,diffq42Pev);
3340                   if(neg) h_diffq42_S_detaN->Fill(DETA,diffq42Nev);
3341                   if(pos) h_ATdiffq22_S_detaP->Fill(DETA,diffq22Pev*charge3subeta);
3342                   if(neg) h_ATdiffq22_S_detaN->Fill(DETA,diffq22Nev*charge3subeta);
3343                   if(pos) h_ATdiffq32_S_detaP->Fill(DETA,diffq32Pev*charge3subeta);
3344                   if(neg) h_ATdiffq32_S_detaN->Fill(DETA,diffq32Nev*charge3subeta);
3345                   if(pos) h_ATdiffq42_S_detaP->Fill(DETA,diffq42Pev*charge3subeta);
3346                   if(neg) h_ATdiffq42_S_detaN->Fill(DETA,diffq42Nev*charge3subeta);
3347
3348                   // --- now field selection for delta eta stuff
3349                   if(mag<-4.0)
3350                     {
3351                       h_AT_X_deta_F1->Fill(DETA,charge3[i]);
3352                       if(pos) h_AT_X_detaP_F1->Fill(DETA,charge3[i]);
3353                       if(neg) h_AT_X_detaN_F1->Fill(DETA,charge3[i]);
3354                       if(pos) h_diffq22_X_detaP_F1->Fill(DETA,diffq22Pev);
3355                       if(neg) h_diffq22_X_detaN_F1->Fill(DETA,diffq22Nev);
3356                       if(pos) h_diffq32_X_detaP_F1->Fill(DETA,diffq32Pev);
3357                       if(neg) h_diffq32_X_detaN_F1->Fill(DETA,diffq32Nev);
3358                       if(pos) h_diffq42_X_detaP_F1->Fill(DETA,diffq42Pev);
3359                       if(neg) h_diffq42_X_detaN_F1->Fill(DETA,diffq42Nev);
3360                       if(pos) h_ATdiffq22_X_detaP_F1->Fill(DETA,diffq22Pev*charge3[i]);
3361                       if(neg) h_ATdiffq22_X_detaN_F1->Fill(DETA,diffq22Nev*charge3[i]);
3362                       if(pos) h_ATdiffq32_X_detaP_F1->Fill(DETA,diffq32Pev*charge3[i]);
3363                       if(neg) h_ATdiffq32_X_detaN_F1->Fill(DETA,diffq32Nev*charge3[i]);
3364                       if(pos) h_ATdiffq42_X_detaP_F1->Fill(DETA,diffq42Pev*charge3[i]);
3365                       if(neg) h_ATdiffq42_X_detaN_F1->Fill(DETA,diffq42Nev*charge3[i]);
3366                     }
3367                   if(mag>4.0)
3368                     {
3369                       h_AT_X_deta_F3->Fill(DETA,charge3[i]);
3370                       if(pos) h_AT_X_detaP_F3->Fill(DETA,charge3[i]);
3371                       if(neg) h_AT_X_detaN_F3->Fill(DETA,charge3[i]);
3372                       if(pos) h_diffq22_X_detaP_F3->Fill(DETA,diffq22Pev);
3373                       if(neg) h_diffq22_X_detaN_F3->Fill(DETA,diffq22Nev);
3374                       if(pos) h_diffq32_X_detaP_F3->Fill(DETA,diffq32Pev);
3375                       if(neg) h_diffq32_X_detaN_F3->Fill(DETA,diffq32Nev);
3376                       if(pos) h_diffq42_X_detaP_F3->Fill(DETA,diffq42Pev);
3377                       if(neg) h_diffq42_X_detaN_F3->Fill(DETA,diffq42Nev);
3378                       if(pos) h_ATdiffq22_X_detaP_F3->Fill(DETA,diffq22Pev*charge3[i]);
3379                       if(neg) h_ATdiffq22_X_detaN_F3->Fill(DETA,diffq22Nev*charge3[i]);
3380                       if(pos) h_ATdiffq32_X_detaP_F3->Fill(DETA,diffq32Pev*charge3[i]);
3381                       if(neg) h_ATdiffq32_X_detaN_F3->Fill(DETA,diffq32Nev*charge3[i]);
3382                       if(pos) h_ATdiffq42_X_detaP_F3->Fill(DETA,diffq42Pev*charge3[i]);
3383                       if(neg) h_ATdiffq42_X_detaN_F3->Fill(DETA,diffq42Nev*charge3[i]);
3384                     }
3385
3386                 } // end nested track loop
3387
3388             } // end if statement to check for nested track loop
3389
3390         } // centrality selection for differential cumulants
3391
3392     } // end of second track loop
3393
3394   float qasymA1 = (float)(ntrkposA1-ntrknegA1)/(ntrkA1);
3395   float qasymA2 = (float)(ntrkposA2-ntrknegA2)/(ntrkA2);
3396
3397   h_rA_cent->Fill(cent,qasymA1);
3398   h_r1A_cent->Fill(cent,qasymA1);
3399   h_r2A_cent->Fill(cent,qasymA2);
3400
3401   h2_rA_cent->Fill(cent,qasymA1);
3402   h2_r1A_cent->Fill(cent,qasymA1);
3403   h2_r2A_cent->Fill(cent,qasymA2);
3404
3405   // random subevents 1
3406   float a1XX[9],  a1YY[9],  a1QQ[9];//,  a1qq[9];
3407   float a1XXp[9], a1YYp[9], a1QQp[9];//, a1qqp[9]; // pos
3408   float a1XXn[9], a1YYn[9], a1QQn[9];//, a1qqn[9]; // neg
3409   for(int i=0; i<6; i++)
3410     {
3411       a1XX[i]=a1Xp[i]+a1Xn[i];
3412       a1YY[i]=a1Yp[i]+a1Yn[i];
3413       a1QQ[i]=a1XX[i]*a1XX[i]+a1YY[i]*a1YY[i];
3414       //a1qq[i]=sqrt(a1QQ[i]/a1XX[0]);
3415       // pos
3416       a1XXp[i]=a1Xp[i];
3417       a1YYp[i]=a1Yp[i];
3418       a1QQp[i]=a1XXp[i]*a1XXp[i]+a1YYp[i]*a1YYp[i];
3419       //a1qqp[i]=sqrt(a1QQp[i]/a1XXp[0]);
3420       // neg
3421       a1XXn[i]=a1Xn[i];
3422       a1YYn[i]=a1Yn[i];
3423       a1QQn[i]=a1XXn[i]*a1XXn[i]+a1YYn[i]*a1YYn[i];
3424       //a1qqn[i]=sqrt(a1QQn[i]/a1XXn[0]);
3425     }
3426   // random subevents 2
3427   float a2XX[9],  a2YY[9],  a2QQ[9];//,  a2qq[9];
3428   float a2XXp[9], a2YYp[9], a2QQp[9];//, a2qqp[9]; // pos
3429   float a2XXn[9], a2YYn[9], a2QQn[9];//, a2qqn[9]; // neg
3430   for(int i=0; i<6; i++)
3431     {
3432       a2XX[i]=a2Xp[i]+a2Xn[i];
3433       a2YY[i]=a2Yp[i]+a2Yn[i];
3434       a2QQ[i]=a2XX[i]*a2XX[i]+a2YY[i]*a2YY[i];
3435       //a2qq[i]=sqrt(a2QQ[i]/a2XX[0]);
3436       // pos
3437       a2XXp[i]=a2Xp[i];
3438       a2YYp[i]=a2Yp[i];
3439       a2QQp[i]=a2XXp[i]*a2XXp[i]+a2YYp[i]*a2YYp[i];
3440       //a2qqp[i]=sqrt(a2QQp[i]/a2XXp[0]);
3441       // neg
3442       a2XXn[i]=a2Xn[i];
3443       a2YYn[i]=a2Yn[i];
3444       a2QQn[i]=a2XXn[i]*a2XXn[i]+a2YYn[i]*a2YYn[i];
3445       //a2qqn[i]=sqrt(a1QQn[i]/a1XXn[0]);
3446     }
3447
3448   float a1M = a1XX[0];
3449   float a1W_2 = a1M*(a1M-(a1M/M));
3450   float a1Mp = a1XXp[0];
3451   float a1Wp_2 = a1Mp*(a1Mp-(a1Mp/Mp));
3452   float a1Mn = a1XXn[0];
3453   float a1Wn_2 = a1Mn*(a1Mn-(a1Mn/Mn));
3454
3455   float a2M = a2XX[0];
3456   float a2W_2 = a2M*(a2M-(a2M/M));
3457   float a2Mp = a2XXp[0];
3458   float a2Wp_2 = a2Mp*(a2Mp-(a2Mp/Mp));
3459   float a2Mn = a2XXn[0];
3460   float a2Wn_2 = a2Mn*(a2Mn-(a2Mn/Mn));
3461
3462   float a1q22ev = (a1QQ[2]-a1M)/a1W_2;
3463   float a1q22Pev = (a1QQp[2]-a1Mp)/a1Wp_2; // pos
3464   float a1q22Nev = (a1QQn[2]-a1Mn)/a1Wn_2; // neg
3465
3466   float a2q22ev = (a2QQ[2]-a2M)/a2W_2;
3467   float a2q22Pev = (a2QQp[2]-a2Mp)/a2Wp_2; // pos
3468   float a2q22Nev = (a2QQn[2]-a2Mn)/a2Wn_2; // neg
3469
3470   h_a1q22_cent->Fill(cent,a1q22ev);
3471   h_a1q22_centP->Fill(cent,a1q22Pev);
3472   h_a1q22_centN->Fill(cent,a1q22Nev);
3473
3474   h_a2q22_cent->Fill(cent,a2q22ev);
3475   h_a2q22_centP->Fill(cent,a2q22Pev);
3476   h_a2q22_centN->Fill(cent,a2q22Nev);
3477
3478   h_rAq22_X1_cent->Fill(cent,a1q22ev*qasymA1);
3479   h_rAq22_X1_centP->Fill(cent,a1q22Pev*qasymA1);
3480   h_rAq22_X1_centN->Fill(cent,a1q22Nev*qasymA1);
3481   h_rAq22_X2_cent->Fill(cent,a2q22ev*qasymA2);
3482   h_rAq22_X2_centP->Fill(cent,a2q22Pev*qasymA2);
3483   h_rAq22_X2_centN->Fill(cent,a2q22Nev*qasymA2);
3484   h_rAq22_X3_cent->Fill(cent,a1q22ev*qasymA2);
3485   h_rAq22_X3_centP->Fill(cent,a1q22Pev*qasymA2);
3486   h_rAq22_X3_centN->Fill(cent,a1q22Nev*qasymA2);
3487   h_rAq22_X4_cent->Fill(cent,a2q22ev*qasymA1);
3488   h_rAq22_X4_centP->Fill(cent,a2q22Pev*qasymA1);
3489   h_rAq22_X4_centN->Fill(cent,a2q22Nev*qasymA1);
3490
3491
3492
3493   // ------------------------------------ //
3494   // --- send data to the output list --- //
3495   // ------------------------------------ //
3496
3497   PostData(1,fOutputList);
3498
3499   // ------------------------- //
3500   // --- end of event loop --- //
3501   // ------------------------- //
3502
3503 }
3504
3505 void AliAnalysisTaskCMEv2A::SetParameters()
3506 {
3507
3508   cout<<"Now setting parameters!"<<endl;
3509
3510   // --- these are the TPC efficiencies
3511   effTPC[0] = 0;//pt=0.05
3512   effTPC[1] = 0;//pt=0.15
3513   effTPC[2] = 0.895961;//pt=0.25
3514   effTPC[3] = 0.92504;//pt=0.35
3515   effTPC[4] = 0.944326;//pt=0.45
3516   effTPC[5] = 0.948261;//pt=0.55
3517   effTPC[6] = 0.948522;//pt=0.65
3518   effTPC[7] = 0.946941;//pt=0.75
3519   effTPC[8] = 0.948537;//pt=0.85
3520   effTPC[9] = 0.95194;//pt=0.95
3521   effTPC[10] = 0.95228;//pt=1.05
3522   effTPC[11] = 0.954775;//pt=1.15
3523   effTPC[12] = 0.956055;//pt=1.25
3524   effTPC[13] = 0.955758;//pt=1.35
3525   effTPC[14] = 0.951268;//pt=1.45
3526   effTPC[15] = 0.952081;//pt=1.55
3527   effTPC[16] = 0.947991;//pt=1.65
3528   effTPC[17] = 0.943755;//pt=1.75
3529   effTPC[18] = 0.938106;//pt=1.85
3530   effTPC[19] = 0.932583;//pt=1.95
3531   effTPC[20] = 0.927206;//pt=2.05
3532   effTPC[21] = 0.921111;//pt=2.15
3533   effTPC[22] = 0.913127;//pt=2.25
3534   effTPC[23] = 0.910508;//pt=2.35
3535   effTPC[24] = 0.905838;//pt=2.45
3536   effTPC[25] = 0.90243;//pt=2.55
3537   effTPC[26] = 0.897893;//pt=2.65
3538   effTPC[27] = 0.895459;//pt=2.75
3539   effTPC[28] = 0.89308;//pt=2.85
3540   effTPC[29] = 0.895013;//pt=2.95
3541   effTPC[30] = 0.890872;//pt=3.05
3542   effTPC[31] = 0.892046;//pt=3.15
3543   effTPC[32] = 0.890133;//pt=3.25
3544   effTPC[33] = 0.894786;//pt=3.35
3545   effTPC[34] = 0.891886;//pt=3.45
3546   effTPC[35] = 0.893649;//pt=3.55
3547   effTPC[36] = 0.891418;//pt=3.65
3548   effTPC[37] = 0.897408;//pt=3.75
3549   effTPC[38] = 0.890273;//pt=3.85
3550   effTPC[39] = 0.896526;//pt=3.95
3551   effTPC[40] = 0.892236;//pt=4.05
3552   effTPC[41] = 0.895075;//pt=4.15
3553   effTPC[42] = 0.895159;//pt=4.25
3554   effTPC[43] = 0.889562;//pt=4.35
3555   effTPC[44] = 0.890777;//pt=4.45
3556   effTPC[45] = 0.892925;//pt=4.55
3557   effTPC[46] = 0.890358;//pt=4.65
3558   effTPC[47] = 0.897014;//pt=4.75
3559   effTPC[48] = 0.886445;//pt=4.85
3560   effTPC[49] = 0.890928;//pt=4.95
3561
3562   float max = effTPC[0];
3563   int maxbin = 0;
3564   for(int i=0; i<49; i++)
3565     {
3566       if(effTPC[i] > max)
3567         {
3568           max = effTPC[i];
3569           maxbin = i;
3570         }
3571     }
3572   cout<<"max for TPC is "<<max<<endl;
3573   cout<<"maxbin for TPC is "<<maxbin<<endl;
3574   for(int i=0; i<50; i++)
3575     {
3576       effTPC[i] /= max;
3577     }
3578
3579   effHYB[0] = 0;//pt=0.05
3580   effHYB[1] = 0;//pt=0.15
3581   effHYB[2] = 0.797902;//pt=0.25
3582   effHYB[3] = 0.842964;//pt=0.35
3583   effHYB[4] = 0.872516;//pt=0.45
3584   effHYB[5] = 0.883225;//pt=0.55
3585   effHYB[6] = 0.886472;//pt=0.65
3586   effHYB[7] = 0.886011;//pt=0.75
3587   effHYB[8] = 0.883906;//pt=0.85
3588   effHYB[9] = 0.882904;//pt=0.95
3589   effHYB[10] = 0.881162;//pt=1.05
3590   effHYB[11] = 0.883611;//pt=1.15
3591   effHYB[12] = 0.884113;//pt=1.25
3592   effHYB[13] = 0.883676;//pt=1.35
3593   effHYB[14] = 0.875165;//pt=1.45
3594   effHYB[15] = 0.87379;//pt=1.55
3595   effHYB[16] = 0.871452;//pt=1.65
3596   effHYB[17] = 0.868833;//pt=1.75
3597   effHYB[18] = 0.863982;//pt=1.85
3598   effHYB[19] = 0.858165;//pt=1.95
3599   effHYB[20] = 0.850305;//pt=2.05
3600   effHYB[21] = 0.845429;//pt=2.15
3601   effHYB[22] = 0.839418;//pt=2.25
3602   effHYB[23] = 0.836466;//pt=2.35
3603   effHYB[24] = 0.833551;//pt=2.45
3604   effHYB[25] = 0.831225;//pt=2.55
3605   effHYB[26] = 0.826823;//pt=2.65
3606   effHYB[27] = 0.822582;//pt=2.75
3607   effHYB[28] = 0.82375;//pt=2.85
3608   effHYB[29] = 0.824866;//pt=2.95
3609   effHYB[30] = 0.822234;//pt=3.05
3610   effHYB[31] = 0.825215;//pt=3.15
3611   effHYB[32] = 0.82225;//pt=3.25
3612   effHYB[33] = 0.825158;//pt=3.35
3613   effHYB[34] = 0.822447;//pt=3.45
3614   effHYB[35] = 0.825239;//pt=3.55
3615   effHYB[36] = 0.825619;//pt=3.65
3616   effHYB[37] = 0.831917;//pt=3.75
3617   effHYB[38] = 0.82425;//pt=3.85
3618   effHYB[39] = 0.830073;//pt=3.95
3619   effHYB[40] = 0.82808;//pt=4.05
3620   effHYB[41] = 0.829005;//pt=4.15
3621   effHYB[42] = 0.831138;//pt=4.25
3622   effHYB[43] = 0.827304;//pt=4.35
3623   effHYB[44] = 0.830745;//pt=4.45
3624   effHYB[45] = 0.832802;//pt=4.55
3625   effHYB[46] = 0.828073;//pt=4.65
3626   effHYB[47] = 0.832246;//pt=4.75
3627   effHYB[48] = 0.82791;//pt=4.85
3628   effHYB[49] = 0.831579;//pt=4.95
3629
3630   max = effHYB[0];
3631   maxbin = 0;
3632   for(int i=0; i<49; i++)
3633     {
3634       if(effHYB[i] > max)
3635         {
3636           max = effHYB[i];
3637           maxbin = i;
3638         }
3639     }
3640   cout<<"max for HYB is "<<max<<endl;
3641   cout<<"maxbin for HYB is "<<maxbin<<endl;
3642   for(int i=0; i<50; i++)
3643     {
3644       effHYB[i] /= max;
3645     }
3646
3647   // ------------------------------------
3648   // --- now for mean charge subtractions
3649   // ------------------------------------
3650
3651   // --- filter bit 1 (global tracks)
3652   MCEF1_fb1[0] = 0.00591299;    MCEF3_fb1[0] = 0.00383748;//eta is -0.99
3653   MCEF1_fb1[1] = 0.00588992;    MCEF3_fb1[1] = 0.003667;//eta is -0.97
3654   MCEF1_fb1[2] = 0.00570083;    MCEF3_fb1[2] = 0.00361253;//eta is -0.95
3655   MCEF1_fb1[3] = 0.00560151;    MCEF3_fb1[3] = 0.00333053;//eta is -0.93
3656   MCEF1_fb1[4] = 0.00576251;    MCEF3_fb1[4] = 0.00340385;//eta is -0.91
3657   MCEF1_fb1[5] = 0.00554351;    MCEF3_fb1[5] = 0.00292592;//eta is -0.89
3658   MCEF1_fb1[6] = 0.0056127;     MCEF3_fb1[6] = 0.00332485;//eta is -0.87
3659   MCEF1_fb1[7] = 0.00541225;    MCEF3_fb1[7] = 0.00331356;//eta is -0.85
3660   MCEF1_fb1[8] = 0.00525774;    MCEF3_fb1[8] = 0.00370341;//eta is -0.83
3661   MCEF1_fb1[9] = 0.00535911;    MCEF3_fb1[9] = 0.00380439;//eta is -0.81
3662   MCEF1_fb1[10] = 0.00532029;   MCEF3_fb1[10] = 0.00418114;//eta is -0.79
3663   MCEF1_fb1[11] = 0.00542025;   MCEF3_fb1[11] = 0.00387263;//eta is -0.77
3664   MCEF1_fb1[12] = 0.00520206;   MCEF3_fb1[12] = 0.00414251;//eta is -0.75
3665   MCEF1_fb1[13] = 0.00513343;   MCEF3_fb1[13] = 0.00400127;//eta is -0.73
3666   MCEF1_fb1[14] = 0.00543818;   MCEF3_fb1[14] = 0.00416512;//eta is -0.71
3667   MCEF1_fb1[15] = 0.00553948;   MCEF3_fb1[15] = 0.00406533;//eta is -0.69
3668   MCEF1_fb1[16] = 0.00556544;   MCEF3_fb1[16] = 0.00428908;//eta is -0.67
3669   MCEF1_fb1[17] = 0.00558886;   MCEF3_fb1[17] = 0.00435018;//eta is -0.65
3670   MCEF1_fb1[18] = 0.0057369;    MCEF3_fb1[18] = 0.00445061;//eta is -0.63
3671   MCEF1_fb1[19] = 0.00616797;   MCEF3_fb1[19] = 0.00443836;//eta is -0.61
3672   MCEF1_fb1[20] = 0.00604426;   MCEF3_fb1[20] = 0.00466612;//eta is -0.59
3673   MCEF1_fb1[21] = 0.00581836;   MCEF3_fb1[21] = 0.00473329;//eta is -0.57
3674   MCEF1_fb1[22] = 0.00573374;   MCEF3_fb1[22] = 0.00429524;//eta is -0.55
3675   MCEF1_fb1[23] = 0.00563837;   MCEF3_fb1[23] = 0.0047092;//eta is -0.53
3676   MCEF1_fb1[24] = 0.00628575;   MCEF3_fb1[24] = 0.0050631;//eta is -0.51
3677   MCEF1_fb1[25] = 0.00606812;   MCEF3_fb1[25] = 0.00480595;//eta is -0.49
3678   MCEF1_fb1[26] = 0.00614111;   MCEF3_fb1[26] = 0.00474149;//eta is -0.47
3679   MCEF1_fb1[27] = 0.0060786;    MCEF3_fb1[27] = 0.00478341;//eta is -0.45
3680   MCEF1_fb1[28] = 0.00621897;   MCEF3_fb1[28] = 0.00499422;//eta is -0.43
3681   MCEF1_fb1[29] = 0.00633416;   MCEF3_fb1[29] = 0.00498704;//eta is -0.41
3682   MCEF1_fb1[30] = 0.00638646;   MCEF3_fb1[30] = 0.00503349;//eta is -0.39
3683   MCEF1_fb1[31] = 0.00658823;   MCEF3_fb1[31] = 0.00485434;//eta is -0.37
3684   MCEF1_fb1[32] = 0.0064771;    MCEF3_fb1[32] = 0.00482841;//eta is -0.35
3685   MCEF1_fb1[33] = 0.00694801;   MCEF3_fb1[33] = 0.00471101;//eta is -0.33
3686   MCEF1_fb1[34] = 0.00636795;   MCEF3_fb1[34] = 0.00491829;//eta is -0.31
3687   MCEF1_fb1[35] = 0.00650891;   MCEF3_fb1[35] = 0.00498886;//eta is -0.29
3688   MCEF1_fb1[36] = 0.00653225;   MCEF3_fb1[36] = 0.00481843;//eta is -0.27
3689   MCEF1_fb1[37] = 0.00657161;   MCEF3_fb1[37] = 0.00511179;//eta is -0.25
3690   MCEF1_fb1[38] = 0.00640258;   MCEF3_fb1[38] = 0.00511465;//eta is -0.23
3691   MCEF1_fb1[39] = 0.00626697;   MCEF3_fb1[39] = 0.00563881;//eta is -0.21
3692   MCEF1_fb1[40] = 0.00691174;   MCEF3_fb1[40] = 0.00534971;//eta is -0.19
3693   MCEF1_fb1[41] = 0.00620562;   MCEF3_fb1[41] = 0.00557349;//eta is -0.17
3694   MCEF1_fb1[42] = 0.00649612;   MCEF3_fb1[42] = 0.00498584;//eta is -0.15
3695   MCEF1_fb1[43] = 0.0065798;    MCEF3_fb1[43] = 0.00542058;//eta is -0.13
3696   MCEF1_fb1[44] = 0.00642303;   MCEF3_fb1[44] = 0.00613666;//eta is -0.11
3697   MCEF1_fb1[45] = 0.00653218;   MCEF3_fb1[45] = 0.00546797;//eta is -0.09
3698   MCEF1_fb1[46] = 0.0064394;    MCEF3_fb1[46] = 0.00555092;//eta is -0.07
3699   MCEF1_fb1[47] = 0.00622478;   MCEF3_fb1[47] = 0.00629308;//eta is -0.05
3700   MCEF1_fb1[48] = 0.00639885;   MCEF3_fb1[48] = 0.0063706;//eta is -0.03
3701   MCEF1_fb1[49] = 0.00666821;   MCEF3_fb1[49] = 0.00631535;//eta is -0.01
3702   MCEF1_fb1[50] = 0.00653648;   MCEF3_fb1[50] = 0.00710942;//eta is 0.01
3703   MCEF1_fb1[51] = 0.00638104;   MCEF3_fb1[51] = 0.00709238;//eta is 0.03
3704   MCEF1_fb1[52] = 0.0066302;    MCEF3_fb1[52] = 0.00752929;//eta is 0.05
3705   MCEF1_fb1[53] = 0.00660631;   MCEF3_fb1[53] = 0.00756394;//eta is 0.07
3706   MCEF1_fb1[54] = 0.0068361;    MCEF3_fb1[54] = 0.00763778;//eta is 0.09
3707   MCEF1_fb1[55] = 0.00680126;   MCEF3_fb1[55] = 0.00766098;//eta is 0.11
3708   MCEF1_fb1[56] = 0.00673586;   MCEF3_fb1[56] = 0.0077562;//eta is 0.13
3709   MCEF1_fb1[57] = 0.0066831;    MCEF3_fb1[57] = 0.00749726;//eta is 0.15
3710   MCEF1_fb1[58] = 0.00632194;   MCEF3_fb1[58] = 0.00814879;//eta is 0.17
3711   MCEF1_fb1[59] = 0.00615406;   MCEF3_fb1[59] = 0.00773741;//eta is 0.19
3712   MCEF1_fb1[60] = 0.00640291;   MCEF3_fb1[60] = 0.00709483;//eta is 0.21
3713   MCEF1_fb1[61] = 0.00626011;   MCEF3_fb1[61] = 0.00748346;//eta is 0.23
3714   MCEF1_fb1[62] = 0.00588407;   MCEF3_fb1[62] = 0.00746276;//eta is 0.25
3715   MCEF1_fb1[63] = 0.00615514;   MCEF3_fb1[63] = 0.0073329;//eta is 0.27
3716   MCEF1_fb1[64] = 0.00581758;   MCEF3_fb1[64] = 0.00749753;//eta is 0.29
3717   MCEF1_fb1[65] = 0.00552913;   MCEF3_fb1[65] = 0.00758203;//eta is 0.31
3718   MCEF1_fb1[66] = 0.00618476;   MCEF3_fb1[66] = 0.00736015;//eta is 0.33
3719   MCEF1_fb1[67] = 0.00571961;   MCEF3_fb1[67] = 0.00717504;//eta is 0.35
3720   MCEF1_fb1[68] = 0.00596315;   MCEF3_fb1[68] = 0.00731385;//eta is 0.37
3721   MCEF1_fb1[69] = 0.00617447;   MCEF3_fb1[69] = 0.00715388;//eta is 0.39
3722   MCEF1_fb1[70] = 0.00548738;   MCEF3_fb1[70] = 0.00691641;//eta is 0.41
3723   MCEF1_fb1[71] = 0.00526914;   MCEF3_fb1[71] = 0.00712093;//eta is 0.43
3724   MCEF1_fb1[72] = 0.00549674;   MCEF3_fb1[72] = 0.00689862;//eta is 0.45
3725   MCEF1_fb1[73] = 0.00556518;   MCEF3_fb1[73] = 0.00701435;//eta is 0.47
3726   MCEF1_fb1[74] = 0.005665;     MCEF3_fb1[74] = 0.00676622;//eta is 0.49
3727   MCEF1_fb1[75] = 0.00552991;   MCEF3_fb1[75] = 0.00679297;//eta is 0.51
3728   MCEF1_fb1[76] = 0.00569837;   MCEF3_fb1[76] = 0.00662394;//eta is 0.53
3729   MCEF1_fb1[77] = 0.00562372;   MCEF3_fb1[77] = 0.00638628;//eta is 0.55
3730   MCEF1_fb1[78] = 0.00543001;   MCEF3_fb1[78] = 0.00657229;//eta is 0.57
3731   MCEF1_fb1[79] = 0.00524306;   MCEF3_fb1[79] = 0.00635786;//eta is 0.59
3732   MCEF1_fb1[80] = 0.00535293;   MCEF3_fb1[80] = 0.00641706;//eta is 0.61
3733   MCEF1_fb1[81] = 0.00514873;   MCEF3_fb1[81] = 0.00636306;//eta is 0.63
3734   MCEF1_fb1[82] = 0.00510695;   MCEF3_fb1[82] = 0.00638945;//eta is 0.65
3735   MCEF1_fb1[83] = 0.00484419;   MCEF3_fb1[83] = 0.0067127;//eta is 0.67
3736   MCEF1_fb1[84] = 0.00484068;   MCEF3_fb1[84] = 0.00621438;//eta is 0.69
3737   MCEF1_fb1[85] = 0.00480704;   MCEF3_fb1[85] = 0.00606007;//eta is 0.71
3738   MCEF1_fb1[86] = 0.00478109;   MCEF3_fb1[86] = 0.00596953;//eta is 0.73
3739   MCEF1_fb1[87] = 0.00488645;   MCEF3_fb1[87] = 0.00594765;//eta is 0.75
3740   MCEF1_fb1[88] = 0.00459445;   MCEF3_fb1[88] = 0.00579464;//eta is 0.77
3741   MCEF1_fb1[89] = 0.00463147;   MCEF3_fb1[89] = 0.0056488;//eta is 0.79
3742   MCEF1_fb1[90] = 0.00464974;   MCEF3_fb1[90] = 0.00567165;//eta is 0.81
3743   MCEF1_fb1[91] = 0.00458123;   MCEF3_fb1[91] = 0.00532498;//eta is 0.83
3744   MCEF1_fb1[92] = 0.0047586;    MCEF3_fb1[92] = 0.00472428;//eta is 0.85
3745   MCEF1_fb1[93] = 0.00492789;   MCEF3_fb1[93] = 0.00514747;//eta is 0.87
3746   MCEF1_fb1[94] = 0.00464254;   MCEF3_fb1[94] = 0.00544923;//eta is 0.89
3747   MCEF1_fb1[95] = 0.00453585;   MCEF3_fb1[95] = 0.00511096;//eta is 0.91
3748   MCEF1_fb1[96] = 0.00427511;   MCEF3_fb1[96] = 0.00556982;//eta is 0.93
3749   MCEF1_fb1[97] = 0.00463089;   MCEF3_fb1[97] = 0.00584684;//eta is 0.95
3750   MCEF1_fb1[98] = 0.00446822;   MCEF3_fb1[98] = 0.00615053;//eta is 0.97
3751   MCEF1_fb1[99] = 0.00433282;   MCEF3_fb1[99] = 0.0063621;//eta is 0.99
3752
3753   // --- filter bit 128 (TPC only tracks)
3754   MCEF1_fb128[0] = 0.00591299;  MCEF3_fb128[0] = 0.00383748;//eta is -0.99
3755   MCEF1_fb128[1] = 0.00588992;  MCEF3_fb128[1] = 0.003667;//eta is -0.97
3756   MCEF1_fb128[2] = 0.00570083;  MCEF3_fb128[2] = 0.00361253;//eta is -0.95
3757   MCEF1_fb128[3] = 0.00560151;  MCEF3_fb128[3] = 0.00333053;//eta is -0.93
3758   MCEF1_fb128[4] = 0.00576251;  MCEF3_fb128[4] = 0.00340385;//eta is -0.91
3759   MCEF1_fb128[5] = 0.00554351;  MCEF3_fb128[5] = 0.00292592;//eta is -0.89
3760   MCEF1_fb128[6] = 0.0056127;   MCEF3_fb128[6] = 0.00332485;//eta is -0.87
3761   MCEF1_fb128[7] = 0.00541225;  MCEF3_fb128[7] = 0.00331356;//eta is -0.85
3762   MCEF1_fb128[8] = 0.00525774;  MCEF3_fb128[8] = 0.00370341;//eta is -0.83
3763   MCEF1_fb128[9] = 0.00535911;  MCEF3_fb128[9] = 0.00380439;//eta is -0.81
3764   MCEF1_fb128[10] = 0.00532029; MCEF3_fb128[10] = 0.00418114;//eta is -0.79
3765   MCEF1_fb128[11] = 0.00542025; MCEF3_fb128[11] = 0.00387263;//eta is -0.77
3766   MCEF1_fb128[12] = 0.00520206; MCEF3_fb128[12] = 0.00414251;//eta is -0.75
3767   MCEF1_fb128[13] = 0.00513343; MCEF3_fb128[13] = 0.00400127;//eta is -0.73
3768   MCEF1_fb128[14] = 0.00543818; MCEF3_fb128[14] = 0.00416512;//eta is -0.71
3769   MCEF1_fb128[15] = 0.00553948; MCEF3_fb128[15] = 0.00406533;//eta is -0.69
3770   MCEF1_fb128[16] = 0.00556544; MCEF3_fb128[16] = 0.00428908;//eta is -0.67
3771   MCEF1_fb128[17] = 0.00558886; MCEF3_fb128[17] = 0.00435018;//eta is -0.65
3772   MCEF1_fb128[18] = 0.0057369;  MCEF3_fb128[18] = 0.00445061;//eta is -0.63
3773   MCEF1_fb128[19] = 0.00616797; MCEF3_fb128[19] = 0.00443836;//eta is -0.61
3774   MCEF1_fb128[20] = 0.00604426; MCEF3_fb128[20] = 0.00466612;//eta is -0.59
3775   MCEF1_fb128[21] = 0.00581836; MCEF3_fb128[21] = 0.00473329;//eta is -0.57
3776   MCEF1_fb128[22] = 0.00573374; MCEF3_fb128[22] = 0.00429524;//eta is -0.55
3777   MCEF1_fb128[23] = 0.00563837; MCEF3_fb128[23] = 0.0047092;//eta is -0.53
3778   MCEF1_fb128[24] = 0.00628575; MCEF3_fb128[24] = 0.0050631;//eta is -0.51
3779   MCEF1_fb128[25] = 0.00606812; MCEF3_fb128[25] = 0.00480595;//eta is -0.49
3780   MCEF1_fb128[26] = 0.00614111; MCEF3_fb128[26] = 0.00474149;//eta is -0.47
3781   MCEF1_fb128[27] = 0.0060786;  MCEF3_fb128[27] = 0.00478341;//eta is -0.45
3782   MCEF1_fb128[28] = 0.00621897; MCEF3_fb128[28] = 0.00499422;//eta is -0.43
3783   MCEF1_fb128[29] = 0.00633416; MCEF3_fb128[29] = 0.00498704;//eta is -0.41
3784   MCEF1_fb128[30] = 0.00638646; MCEF3_fb128[30] = 0.00503349;//eta is -0.39
3785   MCEF1_fb128[31] = 0.00658823; MCEF3_fb128[31] = 0.00485434;//eta is -0.37
3786   MCEF1_fb128[32] = 0.0064771;  MCEF3_fb128[32] = 0.00482841;//eta is -0.35
3787   MCEF1_fb128[33] = 0.00694801; MCEF3_fb128[33] = 0.00471101;//eta is -0.33
3788   MCEF1_fb128[34] = 0.00636795; MCEF3_fb128[34] = 0.00491829;//eta is -0.31
3789   MCEF1_fb128[35] = 0.00650891; MCEF3_fb128[35] = 0.00498886;//eta is -0.29
3790   MCEF1_fb128[36] = 0.00653225; MCEF3_fb128[36] = 0.00481843;//eta is -0.27
3791   MCEF1_fb128[37] = 0.00657161; MCEF3_fb128[37] = 0.00511179;//eta is -0.25
3792   MCEF1_fb128[38] = 0.00640258; MCEF3_fb128[38] = 0.00511465;//eta is -0.23
3793   MCEF1_fb128[39] = 0.00626697; MCEF3_fb128[39] = 0.00563881;//eta is -0.21
3794   MCEF1_fb128[40] = 0.00691174; MCEF3_fb128[40] = 0.00534971;//eta is -0.19
3795   MCEF1_fb128[41] = 0.00620562; MCEF3_fb128[41] = 0.00557349;//eta is -0.17
3796   MCEF1_fb128[42] = 0.00649612; MCEF3_fb128[42] = 0.00498584;//eta is -0.15
3797   MCEF1_fb128[43] = 0.0065798;  MCEF3_fb128[43] = 0.00542058;//eta is -0.13
3798   MCEF1_fb128[44] = 0.00642303; MCEF3_fb128[44] = 0.00613666;//eta is -0.11
3799   MCEF1_fb128[45] = 0.00653218; MCEF3_fb128[45] = 0.00546797;//eta is -0.09
3800   MCEF1_fb128[46] = 0.0064394;  MCEF3_fb128[46] = 0.00555092;//eta is -0.07
3801   MCEF1_fb128[47] = 0.00622478; MCEF3_fb128[47] = 0.00629308;//eta is -0.05
3802   MCEF1_fb128[48] = 0.00639885; MCEF3_fb128[48] = 0.0063706;//eta is -0.03
3803   MCEF1_fb128[49] = 0.00666821; MCEF3_fb128[49] = 0.00631535;//eta is -0.01
3804   MCEF1_fb128[50] = 0.00653648; MCEF3_fb128[50] = 0.00710942;//eta is 0.01
3805   MCEF1_fb128[51] = 0.00638104; MCEF3_fb128[51] = 0.00709238;//eta is 0.03
3806   MCEF1_fb128[52] = 0.0066302;  MCEF3_fb128[52] = 0.00752929;//eta is 0.05
3807   MCEF1_fb128[53] = 0.00660631; MCEF3_fb128[53] = 0.00756394;//eta is 0.07
3808   MCEF1_fb128[54] = 0.0068361;  MCEF3_fb128[54] = 0.00763778;//eta is 0.09
3809   MCEF1_fb128[55] = 0.00680126; MCEF3_fb128[55] = 0.00766098;//eta is 0.11
3810   MCEF1_fb128[56] = 0.00673586; MCEF3_fb128[56] = 0.0077562;//eta is 0.13
3811   MCEF1_fb128[57] = 0.0066831;  MCEF3_fb128[57] = 0.00749726;//eta is 0.15
3812   MCEF1_fb128[58] = 0.00632194; MCEF3_fb128[58] = 0.00814879;//eta is 0.17
3813   MCEF1_fb128[59] = 0.00615406; MCEF3_fb128[59] = 0.00773741;//eta is 0.19
3814   MCEF1_fb128[60] = 0.00640291; MCEF3_fb128[60] = 0.00709483;//eta is 0.21
3815   MCEF1_fb128[61] = 0.00626011; MCEF3_fb128[61] = 0.00748346;//eta is 0.23
3816   MCEF1_fb128[62] = 0.00588407; MCEF3_fb128[62] = 0.00746276;//eta is 0.25
3817   MCEF1_fb128[63] = 0.00615514; MCEF3_fb128[63] = 0.0073329;//eta is 0.27
3818   MCEF1_fb128[64] = 0.00581758; MCEF3_fb128[64] = 0.00749753;//eta is 0.29
3819   MCEF1_fb128[65] = 0.00552913; MCEF3_fb128[65] = 0.00758203;//eta is 0.31
3820   MCEF1_fb128[66] = 0.00618476; MCEF3_fb128[66] = 0.00736015;//eta is 0.33
3821   MCEF1_fb128[67] = 0.00571961; MCEF3_fb128[67] = 0.00717504;//eta is 0.35
3822   MCEF1_fb128[68] = 0.00596315; MCEF3_fb128[68] = 0.00731385;//eta is 0.37
3823   MCEF1_fb128[69] = 0.00617447; MCEF3_fb128[69] = 0.00715388;//eta is 0.39
3824   MCEF1_fb128[70] = 0.00548738; MCEF3_fb128[70] = 0.00691641;//eta is 0.41
3825   MCEF1_fb128[71] = 0.00526914; MCEF3_fb128[71] = 0.00712093;//eta is 0.43
3826   MCEF1_fb128[72] = 0.00549674; MCEF3_fb128[72] = 0.00689862;//eta is 0.45
3827   MCEF1_fb128[73] = 0.00556518; MCEF3_fb128[73] = 0.00701435;//eta is 0.47
3828   MCEF1_fb128[74] = 0.005665;   MCEF3_fb128[74] = 0.00676622;//eta is 0.49
3829   MCEF1_fb128[75] = 0.00552991; MCEF3_fb128[75] = 0.00679297;//eta is 0.51
3830   MCEF1_fb128[76] = 0.00569837; MCEF3_fb128[76] = 0.00662394;//eta is 0.53
3831   MCEF1_fb128[77] = 0.00562372; MCEF3_fb128[77] = 0.00638628;//eta is 0.55
3832   MCEF1_fb128[78] = 0.00543001; MCEF3_fb128[78] = 0.00657229;//eta is 0.57
3833   MCEF1_fb128[79] = 0.00524306; MCEF3_fb128[79] = 0.00635786;//eta is 0.59
3834   MCEF1_fb128[80] = 0.00535293; MCEF3_fb128[80] = 0.00641706;//eta is 0.61
3835   MCEF1_fb128[81] = 0.00514873; MCEF3_fb128[81] = 0.00636306;//eta is 0.63
3836   MCEF1_fb128[82] = 0.00510695; MCEF3_fb128[82] = 0.00638945;//eta is 0.65
3837   MCEF1_fb128[83] = 0.00484419; MCEF3_fb128[83] = 0.0067127;//eta is 0.67
3838   MCEF1_fb128[84] = 0.00484068; MCEF3_fb128[84] = 0.00621438;//eta is 0.69
3839   MCEF1_fb128[85] = 0.00480704; MCEF3_fb128[85] = 0.00606007;//eta is 0.71
3840   MCEF1_fb128[86] = 0.00478109; MCEF3_fb128[86] = 0.00596953;//eta is 0.73
3841   MCEF1_fb128[87] = 0.00488645; MCEF3_fb128[87] = 0.00594765;//eta is 0.75
3842   MCEF1_fb128[88] = 0.00459445; MCEF3_fb128[88] = 0.00579464;//eta is 0.77
3843   MCEF1_fb128[89] = 0.00463147; MCEF3_fb128[89] = 0.0056488;//eta is 0.79
3844   MCEF1_fb128[90] = 0.00464974; MCEF3_fb128[90] = 0.00567165;//eta is 0.81
3845   MCEF1_fb128[91] = 0.00458123; MCEF3_fb128[91] = 0.00532498;//eta is 0.83
3846   MCEF1_fb128[92] = 0.0047586;  MCEF3_fb128[92] = 0.00472428;//eta is 0.85
3847   MCEF1_fb128[93] = 0.00492789; MCEF3_fb128[93] = 0.00514747;//eta is 0.87
3848   MCEF1_fb128[94] = 0.00464254; MCEF3_fb128[94] = 0.00544923;//eta is 0.89
3849   MCEF1_fb128[95] = 0.00453585; MCEF3_fb128[95] = 0.00511096;//eta is 0.91
3850   MCEF1_fb128[96] = 0.00427511; MCEF3_fb128[96] = 0.00556982;//eta is 0.93
3851   MCEF1_fb128[97] = 0.00463089; MCEF3_fb128[97] = 0.00584684;//eta is 0.95
3852   MCEF1_fb128[98] = 0.00446822; MCEF3_fb128[98] = 0.00615053;//eta is 0.97
3853   MCEF1_fb128[99] = 0.00433282; MCEF3_fb128[99] = 0.0063621;//eta is 0.99
3854
3855   // --- filter bit 272 (hybrid tracks)
3856   MCEF1_fb272[0] = 0.00591299;  MCEF3_fb272[0] = 0.00383748;//eta is -0.99
3857   MCEF1_fb272[1] = 0.00588992;  MCEF3_fb272[1] = 0.003667;//eta is -0.97
3858   MCEF1_fb272[2] = 0.00570083;  MCEF3_fb272[2] = 0.00361253;//eta is -0.95
3859   MCEF1_fb272[3] = 0.00560151;  MCEF3_fb272[3] = 0.00333053;//eta is -0.93
3860   MCEF1_fb272[4] = 0.00576251;  MCEF3_fb272[4] = 0.00340385;//eta is -0.91
3861   MCEF1_fb272[5] = 0.00554351;  MCEF3_fb272[5] = 0.00292592;//eta is -0.89
3862   MCEF1_fb272[6] = 0.0056127;   MCEF3_fb272[6] = 0.00332485;//eta is -0.87
3863   MCEF1_fb272[7] = 0.00541225;  MCEF3_fb272[7] = 0.00331356;//eta is -0.85
3864   MCEF1_fb272[8] = 0.00525774;  MCEF3_fb272[8] = 0.00370341;//eta is -0.83
3865   MCEF1_fb272[9] = 0.00535911;  MCEF3_fb272[9] = 0.00380439;//eta is -0.81
3866   MCEF1_fb272[10] = 0.00532029; MCEF3_fb272[10] = 0.00418114;//eta is -0.79
3867   MCEF1_fb272[11] = 0.00542025; MCEF3_fb272[11] = 0.00387263;//eta is -0.77
3868   MCEF1_fb272[12] = 0.00520206; MCEF3_fb272[12] = 0.00414251;//eta is -0.75
3869   MCEF1_fb272[13] = 0.00513343; MCEF3_fb272[13] = 0.00400127;//eta is -0.73
3870   MCEF1_fb272[14] = 0.00543818; MCEF3_fb272[14] = 0.00416512;//eta is -0.71
3871   MCEF1_fb272[15] = 0.00553948; MCEF3_fb272[15] = 0.00406533;//eta is -0.69
3872   MCEF1_fb272[16] = 0.00556544; MCEF3_fb272[16] = 0.00428908;//eta is -0.67
3873   MCEF1_fb272[17] = 0.00558886; MCEF3_fb272[17] = 0.00435018;//eta is -0.65
3874   MCEF1_fb272[18] = 0.0057369;  MCEF3_fb272[18] = 0.00445061;//eta is -0.63
3875   MCEF1_fb272[19] = 0.00616797; MCEF3_fb272[19] = 0.00443836;//eta is -0.61
3876   MCEF1_fb272[20] = 0.00604426; MCEF3_fb272[20] = 0.00466612;//eta is -0.59
3877   MCEF1_fb272[21] = 0.00581836; MCEF3_fb272[21] = 0.00473329;//eta is -0.57
3878   MCEF1_fb272[22] = 0.00573374; MCEF3_fb272[22] = 0.00429524;//eta is -0.55
3879   MCEF1_fb272[23] = 0.00563837; MCEF3_fb272[23] = 0.0047092;//eta is -0.53
3880   MCEF1_fb272[24] = 0.00628575; MCEF3_fb272[24] = 0.0050631;//eta is -0.51
3881   MCEF1_fb272[25] = 0.00606812; MCEF3_fb272[25] = 0.00480595;//eta is -0.49
3882   MCEF1_fb272[26] = 0.00614111; MCEF3_fb272[26] = 0.00474149;//eta is -0.47
3883   MCEF1_fb272[27] = 0.0060786;  MCEF3_fb272[27] = 0.00478341;//eta is -0.45
3884   MCEF1_fb272[28] = 0.00621897; MCEF3_fb272[28] = 0.00499422;//eta is -0.43
3885   MCEF1_fb272[29] = 0.00633416; MCEF3_fb272[29] = 0.00498704;//eta is -0.41
3886   MCEF1_fb272[30] = 0.00638646; MCEF3_fb272[30] = 0.00503349;//eta is -0.39
3887   MCEF1_fb272[31] = 0.00658823; MCEF3_fb272[31] = 0.00485434;//eta is -0.37
3888   MCEF1_fb272[32] = 0.0064771;  MCEF3_fb272[32] = 0.00482841;//eta is -0.35
3889   MCEF1_fb272[33] = 0.00694801; MCEF3_fb272[33] = 0.00471101;//eta is -0.33
3890   MCEF1_fb272[34] = 0.00636795; MCEF3_fb272[34] = 0.00491829;//eta is -0.31
3891   MCEF1_fb272[35] = 0.00650891; MCEF3_fb272[35] = 0.00498886;//eta is -0.29
3892   MCEF1_fb272[36] = 0.00653225; MCEF3_fb272[36] = 0.00481843;//eta is -0.27
3893   MCEF1_fb272[37] = 0.00657161; MCEF3_fb272[37] = 0.00511179;//eta is -0.25
3894   MCEF1_fb272[38] = 0.00640258; MCEF3_fb272[38] = 0.00511465;//eta is -0.23
3895   MCEF1_fb272[39] = 0.00626697; MCEF3_fb272[39] = 0.00563881;//eta is -0.21
3896   MCEF1_fb272[40] = 0.00691174; MCEF3_fb272[40] = 0.00534971;//eta is -0.19
3897   MCEF1_fb272[41] = 0.00620562; MCEF3_fb272[41] = 0.00557349;//eta is -0.17
3898   MCEF1_fb272[42] = 0.00649612; MCEF3_fb272[42] = 0.00498584;//eta is -0.15
3899   MCEF1_fb272[43] = 0.0065798;  MCEF3_fb272[43] = 0.00542058;//eta is -0.13
3900   MCEF1_fb272[44] = 0.00642303; MCEF3_fb272[44] = 0.00613666;//eta is -0.11
3901   MCEF1_fb272[45] = 0.00653218; MCEF3_fb272[45] = 0.00546797;//eta is -0.09
3902   MCEF1_fb272[46] = 0.0064394;  MCEF3_fb272[46] = 0.00555092;//eta is -0.07
3903   MCEF1_fb272[47] = 0.00622478; MCEF3_fb272[47] = 0.00629308;//eta is -0.05
3904   MCEF1_fb272[48] = 0.00639885; MCEF3_fb272[48] = 0.0063706;//eta is -0.03
3905   MCEF1_fb272[49] = 0.00666821; MCEF3_fb272[49] = 0.00631535;//eta is -0.01
3906   MCEF1_fb272[50] = 0.00653648; MCEF3_fb272[50] = 0.00710942;//eta is 0.01
3907   MCEF1_fb272[51] = 0.00638104; MCEF3_fb272[51] = 0.00709238;//eta is 0.03
3908   MCEF1_fb272[52] = 0.0066302;  MCEF3_fb272[52] = 0.00752929;//eta is 0.05
3909   MCEF1_fb272[53] = 0.00660631; MCEF3_fb272[53] = 0.00756394;//eta is 0.07
3910   MCEF1_fb272[54] = 0.0068361;  MCEF3_fb272[54] = 0.00763778;//eta is 0.09
3911   MCEF1_fb272[55] = 0.00680126; MCEF3_fb272[55] = 0.00766098;//eta is 0.11
3912   MCEF1_fb272[56] = 0.00673586; MCEF3_fb272[56] = 0.0077562;//eta is 0.13
3913   MCEF1_fb272[57] = 0.0066831;  MCEF3_fb272[57] = 0.00749726;//eta is 0.15
3914   MCEF1_fb272[58] = 0.00632194; MCEF3_fb272[58] = 0.00814879;//eta is 0.17
3915   MCEF1_fb272[59] = 0.00615406; MCEF3_fb272[59] = 0.00773741;//eta is 0.19
3916   MCEF1_fb272[60] = 0.00640291; MCEF3_fb272[60] = 0.00709483;//eta is 0.21
3917   MCEF1_fb272[61] = 0.00626011; MCEF3_fb272[61] = 0.00748346;//eta is 0.23
3918   MCEF1_fb272[62] = 0.00588407; MCEF3_fb272[62] = 0.00746276;//eta is 0.25
3919   MCEF1_fb272[63] = 0.00615514; MCEF3_fb272[63] = 0.0073329;//eta is 0.27
3920   MCEF1_fb272[64] = 0.00581758; MCEF3_fb272[64] = 0.00749753;//eta is 0.29
3921   MCEF1_fb272[65] = 0.00552913; MCEF3_fb272[65] = 0.00758203;//eta is 0.31
3922   MCEF1_fb272[66] = 0.00618476; MCEF3_fb272[66] = 0.00736015;//eta is 0.33
3923   MCEF1_fb272[67] = 0.00571961; MCEF3_fb272[67] = 0.00717504;//eta is 0.35
3924   MCEF1_fb272[68] = 0.00596315; MCEF3_fb272[68] = 0.00731385;//eta is 0.37
3925   MCEF1_fb272[69] = 0.00617447; MCEF3_fb272[69] = 0.00715388;//eta is 0.39
3926   MCEF1_fb272[70] = 0.00548738; MCEF3_fb272[70] = 0.00691641;//eta is 0.41
3927   MCEF1_fb272[71] = 0.00526914; MCEF3_fb272[71] = 0.00712093;//eta is 0.43
3928   MCEF1_fb272[72] = 0.00549674; MCEF3_fb272[72] = 0.00689862;//eta is 0.45
3929   MCEF1_fb272[73] = 0.00556518; MCEF3_fb272[73] = 0.00701435;//eta is 0.47
3930   MCEF1_fb272[74] = 0.005665;   MCEF3_fb272[74] = 0.00676622;//eta is 0.49
3931   MCEF1_fb272[75] = 0.00552991; MCEF3_fb272[75] = 0.00679297;//eta is 0.51
3932   MCEF1_fb272[76] = 0.00569837; MCEF3_fb272[76] = 0.00662394;//eta is 0.53
3933   MCEF1_fb272[77] = 0.00562372; MCEF3_fb272[77] = 0.00638628;//eta is 0.55
3934   MCEF1_fb272[78] = 0.00543001; MCEF3_fb272[78] = 0.00657229;//eta is 0.57
3935   MCEF1_fb272[79] = 0.00524306; MCEF3_fb272[79] = 0.00635786;//eta is 0.59
3936   MCEF1_fb272[80] = 0.00535293; MCEF3_fb272[80] = 0.00641706;//eta is 0.61
3937   MCEF1_fb272[81] = 0.00514873; MCEF3_fb272[81] = 0.00636306;//eta is 0.63
3938   MCEF1_fb272[82] = 0.00510695; MCEF3_fb272[82] = 0.00638945;//eta is 0.65
3939   MCEF1_fb272[83] = 0.00484419; MCEF3_fb272[83] = 0.0067127;//eta is 0.67
3940   MCEF1_fb272[84] = 0.00484068; MCEF3_fb272[84] = 0.00621438;//eta is 0.69
3941   MCEF1_fb272[85] = 0.00480704; MCEF3_fb272[85] = 0.00606007;//eta is 0.71
3942   MCEF1_fb272[86] = 0.00478109; MCEF3_fb272[86] = 0.00596953;//eta is 0.73
3943   MCEF1_fb272[87] = 0.00488645; MCEF3_fb272[87] = 0.00594765;//eta is 0.75
3944   MCEF1_fb272[88] = 0.00459445; MCEF3_fb272[88] = 0.00579464;//eta is 0.77
3945   MCEF1_fb272[89] = 0.00463147; MCEF3_fb272[89] = 0.0056488;//eta is 0.79
3946   MCEF1_fb272[90] = 0.00464974; MCEF3_fb272[90] = 0.00567165;//eta is 0.81
3947   MCEF1_fb272[91] = 0.00458123; MCEF3_fb272[91] = 0.00532498;//eta is 0.83
3948   MCEF1_fb272[92] = 0.0047586;  MCEF3_fb272[92] = 0.00472428;//eta is 0.85
3949   MCEF1_fb272[93] = 0.00492789; MCEF3_fb272[93] = 0.00514747;//eta is 0.87
3950   MCEF1_fb272[94] = 0.00464254; MCEF3_fb272[94] = 0.00544923;//eta is 0.89
3951   MCEF1_fb272[95] = 0.00453585; MCEF3_fb272[95] = 0.00511096;//eta is 0.91
3952   MCEF1_fb272[96] = 0.00427511; MCEF3_fb272[96] = 0.00556982;//eta is 0.93
3953   MCEF1_fb272[97] = 0.00463089; MCEF3_fb272[97] = 0.00584684;//eta is 0.95
3954   MCEF1_fb272[98] = 0.00446822; MCEF3_fb272[98] = 0.00615053;//eta is 0.97
3955   MCEF1_fb272[99] = 0.00433282; MCEF3_fb272[99] = 0.0063621;//eta is 0.99
3956
3957   // --- filter bit 1 with tight DCA cut...
3958   MCEF1_dca[0] = 0.00591299;    MCEF3_dca[0] = 0.00383748;//eta is -0.99
3959   MCEF1_dca[1] = 0.00588992;    MCEF3_dca[1] = 0.003667;//eta is -0.97
3960   MCEF1_dca[2] = 0.00570083;    MCEF3_dca[2] = 0.00361253;//eta is -0.95
3961   MCEF1_dca[3] = 0.00560151;    MCEF3_dca[3] = 0.00333053;//eta is -0.93
3962   MCEF1_dca[4] = 0.00576251;    MCEF3_dca[4] = 0.00340385;//eta is -0.91
3963   MCEF1_dca[5] = 0.00554351;    MCEF3_dca[5] = 0.00292592;//eta is -0.89
3964   MCEF1_dca[6] = 0.0056127;     MCEF3_dca[6] = 0.00332485;//eta is -0.87
3965   MCEF1_dca[7] = 0.00541225;    MCEF3_dca[7] = 0.00331356;//eta is -0.85
3966   MCEF1_dca[8] = 0.00525774;    MCEF3_dca[8] = 0.00370341;//eta is -0.83
3967   MCEF1_dca[9] = 0.00535911;    MCEF3_dca[9] = 0.00380439;//eta is -0.81
3968   MCEF1_dca[10] = 0.00532029;   MCEF3_dca[10] = 0.00418114;//eta is -0.79
3969   MCEF1_dca[11] = 0.00542025;   MCEF3_dca[11] = 0.00387263;//eta is -0.77
3970   MCEF1_dca[12] = 0.00520206;   MCEF3_dca[12] = 0.00414251;//eta is -0.75
3971   MCEF1_dca[13] = 0.00513343;   MCEF3_dca[13] = 0.00400127;//eta is -0.73
3972   MCEF1_dca[14] = 0.00543818;   MCEF3_dca[14] = 0.00416512;//eta is -0.71
3973   MCEF1_dca[15] = 0.00553948;   MCEF3_dca[15] = 0.00406533;//eta is -0.69
3974   MCEF1_dca[16] = 0.00556544;   MCEF3_dca[16] = 0.00428908;//eta is -0.67
3975   MCEF1_dca[17] = 0.00558886;   MCEF3_dca[17] = 0.00435018;//eta is -0.65
3976   MCEF1_dca[18] = 0.0057369;    MCEF3_dca[18] = 0.00445061;//eta is -0.63
3977   MCEF1_dca[19] = 0.00616797;   MCEF3_dca[19] = 0.00443836;//eta is -0.61
3978   MCEF1_dca[20] = 0.00604426;   MCEF3_dca[20] = 0.00466612;//eta is -0.59
3979   MCEF1_dca[21] = 0.00581836;   MCEF3_dca[21] = 0.00473329;//eta is -0.57
3980   MCEF1_dca[22] = 0.00573374;   MCEF3_dca[22] = 0.00429524;//eta is -0.55
3981   MCEF1_dca[23] = 0.00563837;   MCEF3_dca[23] = 0.0047092;//eta is -0.53
3982   MCEF1_dca[24] = 0.00628575;   MCEF3_dca[24] = 0.0050631;//eta is -0.51
3983   MCEF1_dca[25] = 0.00606812;   MCEF3_dca[25] = 0.00480595;//eta is -0.49
3984   MCEF1_dca[26] = 0.00614111;   MCEF3_dca[26] = 0.00474149;//eta is -0.47
3985   MCEF1_dca[27] = 0.0060786;    MCEF3_dca[27] = 0.00478341;//eta is -0.45
3986   MCEF1_dca[28] = 0.00621897;   MCEF3_dca[28] = 0.00499422;//eta is -0.43
3987   MCEF1_dca[29] = 0.00633416;   MCEF3_dca[29] = 0.00498704;//eta is -0.41
3988   MCEF1_dca[30] = 0.00638646;   MCEF3_dca[30] = 0.00503349;//eta is -0.39
3989   MCEF1_dca[31] = 0.00658823;   MCEF3_dca[31] = 0.00485434;//eta is -0.37
3990   MCEF1_dca[32] = 0.0064771;    MCEF3_dca[32] = 0.00482841;//eta is -0.35
3991   MCEF1_dca[33] = 0.00694801;   MCEF3_dca[33] = 0.00471101;//eta is -0.33
3992   MCEF1_dca[34] = 0.00636795;   MCEF3_dca[34] = 0.00491829;//eta is -0.31
3993   MCEF1_dca[35] = 0.00650891;   MCEF3_dca[35] = 0.00498886;//eta is -0.29
3994   MCEF1_dca[36] = 0.00653225;   MCEF3_dca[36] = 0.00481843;//eta is -0.27
3995   MCEF1_dca[37] = 0.00657161;   MCEF3_dca[37] = 0.00511179;//eta is -0.25
3996   MCEF1_dca[38] = 0.00640258;   MCEF3_dca[38] = 0.00511465;//eta is -0.23
3997   MCEF1_dca[39] = 0.00626697;   MCEF3_dca[39] = 0.00563881;//eta is -0.21
3998   MCEF1_dca[40] = 0.00691174;   MCEF3_dca[40] = 0.00534971;//eta is -0.19
3999   MCEF1_dca[41] = 0.00620562;   MCEF3_dca[41] = 0.00557349;//eta is -0.17
4000   MCEF1_dca[42] = 0.00649612;   MCEF3_dca[42] = 0.00498584;//eta is -0.15
4001   MCEF1_dca[43] = 0.0065798;    MCEF3_dca[43] = 0.00542058;//eta is -0.13
4002   MCEF1_dca[44] = 0.00642303;   MCEF3_dca[44] = 0.00613666;//eta is -0.11
4003   MCEF1_dca[45] = 0.00653218;   MCEF3_dca[45] = 0.00546797;//eta is -0.09
4004   MCEF1_dca[46] = 0.0064394;    MCEF3_dca[46] = 0.00555092;//eta is -0.07
4005   MCEF1_dca[47] = 0.00622478;   MCEF3_dca[47] = 0.00629308;//eta is -0.05
4006   MCEF1_dca[48] = 0.00639885;   MCEF3_dca[48] = 0.0063706;//eta is -0.03
4007   MCEF1_dca[49] = 0.00666821;   MCEF3_dca[49] = 0.00631535;//eta is -0.01
4008   MCEF1_dca[50] = 0.00653648;   MCEF3_dca[50] = 0.00710942;//eta is 0.01
4009   MCEF1_dca[51] = 0.00638104;   MCEF3_dca[51] = 0.00709238;//eta is 0.03
4010   MCEF1_dca[52] = 0.0066302;    MCEF3_dca[52] = 0.00752929;//eta is 0.05
4011   MCEF1_dca[53] = 0.00660631;   MCEF3_dca[53] = 0.00756394;//eta is 0.07
4012   MCEF1_dca[54] = 0.0068361;    MCEF3_dca[54] = 0.00763778;//eta is 0.09
4013   MCEF1_dca[55] = 0.00680126;   MCEF3_dca[55] = 0.00766098;//eta is 0.11
4014   MCEF1_dca[56] = 0.00673586;   MCEF3_dca[56] = 0.0077562;//eta is 0.13
4015   MCEF1_dca[57] = 0.0066831;    MCEF3_dca[57] = 0.00749726;//eta is 0.15
4016   MCEF1_dca[58] = 0.00632194;   MCEF3_dca[58] = 0.00814879;//eta is 0.17
4017   MCEF1_dca[59] = 0.00615406;   MCEF3_dca[59] = 0.00773741;//eta is 0.19
4018   MCEF1_dca[60] = 0.00640291;   MCEF3_dca[60] = 0.00709483;//eta is 0.21
4019   MCEF1_dca[61] = 0.00626011;   MCEF3_dca[61] = 0.00748346;//eta is 0.23
4020   MCEF1_dca[62] = 0.00588407;   MCEF3_dca[62] = 0.00746276;//eta is 0.25
4021   MCEF1_dca[63] = 0.00615514;   MCEF3_dca[63] = 0.0073329;//eta is 0.27
4022   MCEF1_dca[64] = 0.00581758;   MCEF3_dca[64] = 0.00749753;//eta is 0.29
4023   MCEF1_dca[65] = 0.00552913;   MCEF3_dca[65] = 0.00758203;//eta is 0.31
4024   MCEF1_dca[66] = 0.00618476;   MCEF3_dca[66] = 0.00736015;//eta is 0.33
4025   MCEF1_dca[67] = 0.00571961;   MCEF3_dca[67] = 0.00717504;//eta is 0.35
4026   MCEF1_dca[68] = 0.00596315;   MCEF3_dca[68] = 0.00731385;//eta is 0.37
4027   MCEF1_dca[69] = 0.00617447;   MCEF3_dca[69] = 0.00715388;//eta is 0.39
4028   MCEF1_dca[70] = 0.00548738;   MCEF3_dca[70] = 0.00691641;//eta is 0.41
4029   MCEF1_dca[71] = 0.00526914;   MCEF3_dca[71] = 0.00712093;//eta is 0.43
4030   MCEF1_dca[72] = 0.00549674;   MCEF3_dca[72] = 0.00689862;//eta is 0.45
4031   MCEF1_dca[73] = 0.00556518;   MCEF3_dca[73] = 0.00701435;//eta is 0.47
4032   MCEF1_dca[74] = 0.005665;     MCEF3_dca[74] = 0.00676622;//eta is 0.49
4033   MCEF1_dca[75] = 0.00552991;   MCEF3_dca[75] = 0.00679297;//eta is 0.51
4034   MCEF1_dca[76] = 0.00569837;   MCEF3_dca[76] = 0.00662394;//eta is 0.53
4035   MCEF1_dca[77] = 0.00562372;   MCEF3_dca[77] = 0.00638628;//eta is 0.55
4036   MCEF1_dca[78] = 0.00543001;   MCEF3_dca[78] = 0.00657229;//eta is 0.57
4037   MCEF1_dca[79] = 0.00524306;   MCEF3_dca[79] = 0.00635786;//eta is 0.59
4038   MCEF1_dca[80] = 0.00535293;   MCEF3_dca[80] = 0.00641706;//eta is 0.61
4039   MCEF1_dca[81] = 0.00514873;   MCEF3_dca[81] = 0.00636306;//eta is 0.63
4040   MCEF1_dca[82] = 0.00510695;   MCEF3_dca[82] = 0.00638945;//eta is 0.65
4041   MCEF1_dca[83] = 0.00484419;   MCEF3_dca[83] = 0.0067127;//eta is 0.67
4042   MCEF1_dca[84] = 0.00484068;   MCEF3_dca[84] = 0.00621438;//eta is 0.69
4043   MCEF1_dca[85] = 0.00480704;   MCEF3_dca[85] = 0.00606007;//eta is 0.71
4044   MCEF1_dca[86] = 0.00478109;   MCEF3_dca[86] = 0.00596953;//eta is 0.73
4045   MCEF1_dca[87] = 0.00488645;   MCEF3_dca[87] = 0.00594765;//eta is 0.75
4046   MCEF1_dca[88] = 0.00459445;   MCEF3_dca[88] = 0.00579464;//eta is 0.77
4047   MCEF1_dca[89] = 0.00463147;   MCEF3_dca[89] = 0.0056488;//eta is 0.79
4048   MCEF1_dca[90] = 0.00464974;   MCEF3_dca[90] = 0.00567165;//eta is 0.81
4049   MCEF1_dca[91] = 0.00458123;   MCEF3_dca[91] = 0.00532498;//eta is 0.83
4050   MCEF1_dca[92] = 0.0047586;    MCEF3_dca[92] = 0.00472428;//eta is 0.85
4051   MCEF1_dca[93] = 0.00492789;   MCEF3_dca[93] = 0.00514747;//eta is 0.87
4052   MCEF1_dca[94] = 0.00464254;   MCEF3_dca[94] = 0.00544923;//eta is 0.89
4053   MCEF1_dca[95] = 0.00453585;   MCEF3_dca[95] = 0.00511096;//eta is 0.91
4054   MCEF1_dca[96] = 0.00427511;   MCEF3_dca[96] = 0.00556982;//eta is 0.93
4055   MCEF1_dca[97] = 0.00463089;   MCEF3_dca[97] = 0.00584684;//eta is 0.95
4056   MCEF1_dca[98] = 0.00446822;   MCEF3_dca[98] = 0.00615053;//eta is 0.97
4057   MCEF1_dca[99] = 0.00433282;   MCEF3_dca[99] = 0.0063621;//eta is 0.99
4058
4059
4060
4061 }
4062
4063
4064
4065 float AliAnalysisTaskCMEv2A::calc3(float xn, float yn, float x2n, float y2n, float M)
4066 {
4067
4068   float Qn2 = xn*xn+yn*yn;
4069   float Q2n2 = x2n*x2n+y2n*y2n;
4070   float Qn2d = xn*xn-yn*yn;
4071
4072   float W_3 = M*(M-1)*(M-2);
4073
4074   float first = ((x2n*Qn2d)-(2*Qn2))/W_3;
4075   float second = (Q2n2-(2*M))/W_3;
4076
4077   return first-second;
4078
4079 }
4080
4081
4082
4083 float AliAnalysisTaskCMEv2A::calc4(float xn, float yn, float x2n, float y2n, float M)
4084 {
4085
4086   float Qn2 = xn*xn+yn*yn;
4087   float Q2n2 = x2n*x2n+y2n*y2n;
4088   float Qn4 = Qn2*Qn2;
4089   float Qn2d = xn*xn-yn*yn;
4090
4091   float W_4 = M*(M-1)*(M-2)*(M-3);
4092
4093   float first = (Qn4+Q2n2-(2*(x2n*Qn2d)))/W_4;
4094   float second = 2*((2*(M-2)*Qn2)-(M*(M-3)))/W_4;
4095
4096   return first-second;
4097
4098 }
4099
4100
4101 // function definition
4102 float GetDPhiStar(float phi1, float pt1, float charge1, float phi2, float pt2, float charge2, float radius, float bSign)
4103 {
4104   //
4105   // calculates dphistar
4106   //
4107   float dphistar = phi1 - phi2 - charge1 * bSign * asin(0.075 * radius / pt1) + charge2 * bSign * asin(0.075 * radius / pt2);
4108
4109   //static const double kPi = pi;
4110   float kPi = pi;
4111
4112   // circularity
4113   //   if (dphistar > 2 * kPi)
4114   //     dphistar -= 2 * kPi;
4115   //   if (dphistar < -2 * kPi)
4116   //     dphistar += 2 * kPi;
4117
4118   if (dphistar > kPi)
4119     dphistar = kPi * 2 - dphistar;
4120   if (dphistar < -kPi)
4121     dphistar = -kPi * 2 - dphistar;
4122   if (dphistar > kPi) // might look funny but is needed
4123     dphistar = kPi * 2 - dphistar;
4124
4125   return dphistar;
4126
4127 }
4128
4129