PWGDQ/dielectron/AliAnalysisTaskMultiDielectronTG.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 ///////////////////////////////////////////////////////////////////////////
  17 //                                                                       //
  18 //                        Basic Analysis Task                            //
  19 //                                                                       //
  20 ///////////////////////////////////////////////////////////////////////////
  21
  22 #include <TChain.h>
  23 #include <TH1D.h>
  24
  25 #include <AliCFContainer.h>
  26 #include <AliInputEventHandler.h>
  27 #include <AliESDInputHandler.h>
  28 #include <AliAODInputHandler.h>
  29 #include <AliAnalysisManager.h>
  30 #include <AliVEvent.h>
  31 #include <AliTriggerAnalysis.h>
  32 #include <AliPIDResponse.h>
  33 #include <AliTPCPIDResponse.h>
  34
  35 #include "AliDielectron.h"
  36 #include "AliDielectronHistos.h"
  37 #include "AliDielectronCF.h"
  38 #include "AliDielectronMC.h"
  39 #include "AliDielectronMixingHandler.h"
  40 #include "AliAnalysisTaskMultiDielectronTG.h"
  41
  42 #include "AliESDtrack.h"
  43 #include "AliESDtrackCuts.h"
  44
  45 #include "AliLog.h"
  46
  47 #include <vector>
  48 #include <deque>
  49 #include <cstdlib>
  50 #include <TRandom.h>
  51
  52 using std::vector;
  53 using std::deque;
  54
  55 const char* kPairClassNamesTG[7] = {
  56   "unlike",
  57   "like_pp",
  58   "like_ee",
  59   "mixunlike_pe",
  60   "mixunlike_ep",
  61   "mixlike_pp",
  62   "mixlike_ee"
  63 };
  64
  65
  66 ClassImp(AliDielectronSingleTG)
  67 ClassImp(AliAnalysisTaskMultiDielectronTG)
  68
  69 //_________________________________________________________________________________
  70 AliAnalysisTaskMultiDielectronTG::AliAnalysisTaskMultiDielectronTG() :
  71   AliAnalysisTaskSE(),
  72   fListDielectron(),
  73   fListHistos(),
  74   fListCF(),
  75   fQAElectron(),
  76   fSelectPhysics(kFALSE),
  77   fTriggerMask(AliVEvent::kMB),
  78   fExcludeTriggerMask(0),
  79   fTriggerOnV0AND(kFALSE),
  80   fRejectPileup(kFALSE),
  81   fTriggerLogic(kAny),
  82   fTriggerAnalysis(0x0),
  83   fEventFilter(0x0),
  84   fCutsMother(0x0),
  85   fEventStat(0x0),
  86   fEvent(0x0),
  87   fdEdXvsPt(0x0),
  88   fdEdXnSigmaElecvsPt(0x0),
  89   fdEdXvsPtTOF(0x0),
  90   fdEdXnSigmaElecvsPtTOF(0x0),
  91   fTOFbetavsPt(0x0),
  92   fTOFnSigmaElecvsPt(0x0),
  93   fNCrossedRowsTPC(0x0),
  94   fChi2ClusTPC(0x0),
  95   fRatioCrossClusTPC(0x0),
  96   fVem(0x0),
  97   fVep(0x0),
  98   fVemtmp(0x0),
  99   fVeptmp(0x0),
 100   fdconvphiv(acos(-1.0)),
 101   fdop(0),
 102   fbz(0),
 103   fdv0mixing(kTRUE),
 104   fBGRejUnlike(kFALSE),
 105   fBGRejLike(kTRUE)
 106 {
 107   //
 108   // Constructor
 109   //
 110 }
 111
 112 //_________________________________________________________________________________
 113 AliAnalysisTaskMultiDielectronTG::AliAnalysisTaskMultiDielectronTG(const char *name) :
 114   AliAnalysisTaskSE(name),
 115   fListDielectron(),
 116   fListHistos(),
 117   fListCF(),
 118   fQAElectron(),
 119   fSelectPhysics(kFALSE),
 120   fTriggerMask(AliVEvent::kMB),
 121   fExcludeTriggerMask(0),
 122   fTriggerOnV0AND(kFALSE),
 123   fRejectPileup(kFALSE),
 124   fTriggerLogic(kAny),
 125   fTriggerAnalysis(0x0),
 126   fEventFilter(0x0),
 127   fCutsMother(0x0),
 128   fEventStat(0x0),
 129   fEvent(0x0),
 130   fdEdXvsPt(0x0),
 131   fdEdXnSigmaElecvsPt(0x0),
 132   fdEdXvsPtTOF(0x0),
 133   fdEdXnSigmaElecvsPtTOF(0x0),
 134   fTOFbetavsPt(0x0),
 135   fTOFnSigmaElecvsPt(0x0),
 136   fNCrossedRowsTPC(0x0),
 137   fChi2ClusTPC(0x0),
 138   fRatioCrossClusTPC(0x0),
 139   fVem(0x0),
 140   fVep(0x0),
 141   fVemtmp(0x0),
 142   fVeptmp(0x0),
 143   fdconvphiv(acos(-1.0)),
 144   fdop(0),
 145   fbz(0),
 146   fdv0mixing(kTRUE),
 147   fBGRejUnlike(kFALSE),
 148   fBGRejLike(kTRUE)
 149 {
 150   //
 151   // Constructor
 152   //
 153   DefineInput(0,TChain::Class());
 154   DefineOutput(1, TList::Class());
 155   DefineOutput(2, TList::Class());
 156   DefineOutput(3, TH1D::Class());
 157   fListHistos.SetName("Dielectron_Histos_Multi");
 158   fListCF.SetName("Dielectron_CF_Multi");
 159   fListDielectron.SetOwner();
 160   fListHistos.SetOwner();
 161   fListCF.SetOwner();
 162
 163   ///////////////
 164   for(int i=0;i<fgkNDIE; i++){
 165     for(int j=0;j<fgkNZBIN;j++){
 166       for(int k=0;k<fgkNCENT;k++){
 167         for(int l=0; l<fgkNRPBIN; l++){
 168           fibuf[i][j][k][l] = 0;
 169           fpoolp[i][j][k][l].clear();
 170           fpoolm[i][j][k][l].clear();
 171           for(int ib=0;ib<fgkNBUF; ib++){
 172             fvep[ib][i][j][k][l].clear();
 173             fvem[ib][i][j][k][l].clear();
 174           }
 175         }
 176       }
 177     }
 178   }
 179
 180
 181
 182 }
 183
 184 //_________________________________________________________________________________
 185 AliAnalysisTaskMultiDielectronTG::~AliAnalysisTaskMultiDielectronTG()
 186 {
 187   //
 188   // Destructor
 189   //
 190
 191   //histograms and CF are owned by the dielectron framework.
 192   //however they are streamed to file, so in the first place the
 193   //lists need to be owner...
 194   fListHistos.SetOwner(kFALSE);
 195   fListCF.SetOwner(kFALSE);
 196
 197   for(int i=0;i<fgkNDIE; i++){
 198     for(int j=0;j<fgkNZBIN;j++){
 199       for(int k=0;k<fgkNCENT;k++){
 200         for(int l=0; l<fgkNRPBIN; l++){
 201           fibuf[i][j][k][l] = 0;
 202           fpoolp[i][j][k][l].clear();
 203           fpoolm[i][j][k][l].clear();
 204           for(int ib=0;ib<fgkNBUF; ib++){
 205             fvep[ib][i][j][k][l].clear();
 206             fvem[ib][i][j][k][l].clear();
 207           }
 208         }
 209       }
 210     }
 211   }
 212 }
 213 //_________________________________________________________________________________
 214 void AliAnalysisTaskMultiDielectronTG::UserCreateOutputObjects()
 215 {
 216   //
 217   // Add all histogram manager histogram lists to the output TList
 218   //
 219
 220   if (!fListHistos.IsEmpty()||!fListCF.IsEmpty()) return; //already initialised
 221
 222 //   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 223 //   Bool_t isESD=man->GetInputEventHandler()->IsA()==AliESDInputHandler::Class();
 224 //   Bool_t isAOD=man->GetInputEventHandler()->IsA()==AliAODInputHandler::Class();
 225
 226   TIter nextDie(&fListDielectron);
 227   AliDielectron *die=0;
 228   while ( (die=static_cast<AliDielectron*>(nextDie())) ){
 229     die->Init();
 230     if (die->GetHistogramList()) fListHistos.Add(const_cast<THashList*>(die->GetHistogramList()));
 231     if (die->GetCFManagerPair()) fListCF.Add(const_cast<AliCFContainer*>(die->GetCFManagerPair()->GetContainer()));
 232   }
 233
 234   Int_t cuts=fListDielectron.GetEntries();
 235   Int_t nbins=kNbinsEvent+2*cuts;
 236   if (!fEventStat){
 237     fEventStat=new TH1D("hEventStat","Event statistics",nbins,0,nbins);
 238     fEventStat->GetXaxis()->SetBinLabel(1,"Before Phys. Sel.");
 239     fEventStat->GetXaxis()->SetBinLabel(2,"After Phys. Sel.");
 240
 241     //default names
 242     fEventStat->GetXaxis()->SetBinLabel(3,"Bin3 not used");
 243     fEventStat->GetXaxis()->SetBinLabel(4,"Bin4 not used");
 244     fEventStat->GetXaxis()->SetBinLabel(5,"Bin5 not used");
 245
 246     if(fTriggerOnV0AND) fEventStat->GetXaxis()->SetBinLabel(3,"V0and triggers");
 247     if (fEventFilter) fEventStat->GetXaxis()->SetBinLabel(4,"After Event Filter");
 248     if (fRejectPileup) fEventStat->GetXaxis()->SetBinLabel(5,"After Pileup rejection");
 249
 250     for (Int_t i=0; i<cuts; ++i){
 251       fEventStat->GetXaxis()->SetBinLabel((kNbinsEvent+1)+2*i,Form("#splitline{1 candidate}{%s}",fListDielectron.At(i)->GetName()));
 252       fEventStat->GetXaxis()->SetBinLabel((kNbinsEvent+2)+2*i,Form("#splitline{With >1 candidate}{%s}",fListDielectron.At(i)->GetName()));
 253     }
 254   }
 255
 256   if (!fTriggerAnalysis) fTriggerAnalysis=new AliTriggerAnalysis;
 257   fTriggerAnalysis->EnableHistograms();
 258   fTriggerAnalysis->SetAnalyzeMC(AliDielectronMC::Instance()->HasMC());
 259
 260
 261
 262   int nbinx=400;
 263   float maxx=20;
 264   float minx=0.2;
 265   float binw = (TMath::Log(maxx)-TMath::Log(minx))/nbinx;
 266   double xbin[401];
 267   for(int ii=0;ii<nbinx+1;ii++){
 268     xbin[ii] = TMath::Exp(TMath::Log(minx) + 0.5*binw+binw*ii);
 269   }
 270
 271
 272   fQAElectron = new TList();
 273   fQAElectron->SetName("QAElectron");
 274   fQAElectron->SetOwner();
 275
 276
 277   fEvent = new TH1D("Event","centrality",   100,0,100);
 278   fQAElectron->Add(fEvent);
 279   fdEdXvsPt = new TH2D("dEdXvsPt","dE/dX vs. PT of TPC", nbinx, xbin, 2000,0,200);
 280   fQAElectron->Add(fdEdXvsPt);
 281   fdEdXnSigmaElecvsPt = new TH2D("fdEdXnSigmaElecvsPt"," dE/dX normalized to electron vs. pT of TPC",
 282                                  nbinx, xbin, 2000, -10, 10);
 283   fQAElectron->Add(fdEdXnSigmaElecvsPt);
 284
 285   fdEdXvsPtTOF = new TH2D("dEdXvsPtTOF","dE/dX vs. PT of TPC", nbinx, xbin, 2000,0,200);
 286   fQAElectron->Add(fdEdXvsPtTOF);
 287   fdEdXnSigmaElecvsPtTOF = new TH2D("fdEdXnSigmaElecvsPtTOF"," dE/dX normalized to electron vs. pT of TPC",
 288                                  nbinx, xbin, 2000, -10, 10);
 289   fQAElectron->Add(fdEdXnSigmaElecvsPtTOF);
 290
 291
 292
 293   fTOFbetavsPt = new TH2D("fTOFbetavsPt","TOF beta vs. p", 400, 0, 20, 1200, 0, 1.2);
 294   fQAElectron->Add(fTOFbetavsPt);
 295   fTOFnSigmaElecvsPt = new TH2D("fTOFnSigmaElecvsPt","TOF nsigma for electron", 400, 0, 20, 2000, -10, 10);
 296   fQAElectron->Add(fTOFnSigmaElecvsPt);
 297
 298   fNCrossedRowsTPC = new TH2F("fNCrossedRowsTPC", "TPC nCrossed Rows vs. pT", 200, 0, 20, 200, 0, 200);
 299   fQAElectron->Add(fNCrossedRowsTPC);
 300   fChi2ClusTPC = new TH2F("fChi2ClusTPC", "hChi2ClusTPC vs. pT", 200, 0, 20, 200, 0, 10);
 301   fQAElectron->Add(fChi2ClusTPC);
 302
 303   fRatioCrossClusTPC = new TH2F("fRatioCrossClusTPC", "hRatioCrossClusTPC vs. pT", 200, 0, 20, 200, 0, 10);
 304   fQAElectron->Add(fRatioCrossClusTPC);
 305
 306   fListHistos.Add(fQAElectron);
 307
 308   fListHistos.SetOwner();
 309
 310   PostData(1, &fListHistos);
 311   PostData(2, &fListCF);
 312   PostData(3, fEventStat);
 313
 314   fCutsMother = new AliESDtrackCuts;
 315   fCutsMother->SetDCAToVertex2D(kTRUE);
 316   fCutsMother->SetMaxDCAToVertexZ(3.0);
 317   fCutsMother->SetMaxDCAToVertexXY(1.0);
 318   fCutsMother->SetPtRange(  0.05 , 200.0);
 319   fCutsMother->SetEtaRange( -0.84 , 0.84 );
 320   fCutsMother->SetAcceptKinkDaughters(kFALSE);
 321   fCutsMother->SetRequireITSRefit(kTRUE);
 322   fCutsMother->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kFirst);
 323   fCutsMother->SetMinNClustersITS(3);
 324   fCutsMother->SetRequireTPCRefit(kTRUE);
 325
 326
 327   AliInfo(Form("PairCutType = %d %d %d %d %d",
 328                fRejectPairFlag[0],
 329                fRejectPairFlag[1],
 330                fRejectPairFlag[2],
 331                fRejectPairFlag[3],
 332                fRejectPairFlag[4]));
 333
 334 }
 335
 336 //_________________________________________________________________________________
 337 void AliAnalysisTaskMultiDielectronTG::UserExec(Option_t *)
 338 {
 339   //
 340   // Main loop. Called for every event
 341   //
 342
 343   if (fListHistos.IsEmpty()&&fListCF.IsEmpty()) return;
 344
 345   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 346   Bool_t isESD=man->GetInputEventHandler()->IsA()==AliESDInputHandler::Class();
 347   Bool_t isAOD=man->GetInputEventHandler()->IsA()==AliAODInputHandler::Class();
 348
 349   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 350   if (!inputHandler) return;
 351
 352 //   AliPIDResponse *pidRes=inputHandler->GetPIDResponse();
 353   if ( inputHandler->GetPIDResponse() ){
 354     // for the 2.76 pass2 MC private train. Together with a sigma shift of -0.169
 355 //    pidRes->GetTPCResponse().SetSigma(4.637e-3,2.41332105409873257e+04);
 356     AliDielectronVarManager::SetPIDResponse( inputHandler->GetPIDResponse() );
 357   } else {
 358     AliFatal("This task needs the PID response attached to the input event handler!");
 359   }
 360
 361   // Was event selected ?
 362   ULong64_t isSelected = AliVEvent::kAny;
 363   Bool_t isRejected = kFALSE;
 364   if( fSelectPhysics && inputHandler){
 365     if((isESD && inputHandler->GetEventSelection()) || isAOD){
 366       isSelected = inputHandler->IsEventSelected();
 367       if (fExcludeTriggerMask && (isSelected&fExcludeTriggerMask)) isRejected=kTRUE;
 368       if (fTriggerLogic==kAny) isSelected&=fTriggerMask;
 369       else if (fTriggerLogic==kExact) isSelected=((isSelected&fTriggerMask)==fTriggerMask);
 370     }
 371    }
 372
 373
 374   //Before physics selection
 375   fEventStat->Fill(kAllEvents);
 376   if (isSelected==0||isRejected) {
 377     PostData(3,fEventStat);
 378     return;
 379   }
 380   //after physics selection
 381   fEventStat->Fill(kSelectedEvents);
 382
 383   //V0and
 384   if(fTriggerOnV0AND){
 385     if(isESD){if (!fTriggerAnalysis->IsOfflineTriggerFired(static_cast<AliESDEvent*>(InputEvent()), AliTriggerAnalysis::kV0AND))
 386             return;}
 387     if(isAOD){if(!((static_cast<AliAODEvent*>(InputEvent()))->GetVZEROData()->GetV0ADecision() == AliVVZERO::kV0BB &&
 388                    (static_cast<AliAODEvent*>(InputEvent()))->GetVZEROData()->GetV0CDecision() == AliVVZERO::kV0BB) )
 389         return;}
 390   }
 391
 392
 393   fEventStat->Fill(kV0andEvents);
 394
 395   //Fill Event histograms before the event filter
 396   TIter nextDie(&fListDielectron);
 397   AliDielectron *die=0;
 398   Double_t values[AliDielectronVarManager::kNMaxValues]={0};
 399   Double_t valuesMC[AliDielectronVarManager::kNMaxValues]={0};
 400   AliDielectronVarManager::SetEvent(InputEvent());
 401   AliDielectronVarManager::Fill(InputEvent(),values);
 402   AliDielectronVarManager::Fill(InputEvent(),valuesMC);
 403   Bool_t hasMC=AliDielectronMC::Instance()->HasMC();
 404   if (hasMC) {
 405     if (AliDielectronMC::Instance()->ConnectMCEvent())
 406       AliDielectronVarManager::Fill(AliDielectronMC::Instance()->GetMCEvent(),valuesMC);
 407   }
 408
 409   while ( (die=static_cast<AliDielectron*>(nextDie())) ){
 410     AliDielectronHistos *h=die->GetHistoManager();
 411     if (h){
 412       if (h->GetHistogramList()->FindObject("Event_noCuts"))
 413         h->FillClass("Event_noCuts",AliDielectronVarManager::kNMaxValues,values);
 414       if (hasMC && h->GetHistogramList()->FindObject("MCEvent_noCuts"))
 415         h->FillClass("Event_noCuts",AliDielectronVarManager::kNMaxValues,valuesMC);
 416     }
 417   }
 418   nextDie.Reset();
 419
 420   //event filter
 421   if (fEventFilter) {
 422     if (!fEventFilter->IsSelected(InputEvent())) return;
 423   }
 424   fEventStat->Fill(kFilteredEvents);
 425
 426   //pileup
 427   if (fRejectPileup){
 428     if (InputEvent()->IsPileupFromSPD(3,0.8,3.,2.,5.)) return;
 429   }
 430   fEventStat->Fill(kPileupEvents);
 431
 432   //bz for AliKF
 433   fbz = InputEvent()->GetMagneticField();
 434   AliKFParticle::SetField( fbz );
 435   AliDielectronPID::SetCorrVal((Double_t)InputEvent()->GetRunNumber());
 436
 437   //Process event in all AliDielectron instances
 438   //   TIter nextDie(&fListDielectron);
 439   //   AliDielectron *die=0;
 440   Int_t idie=0;
 441   while ( (die=static_cast<AliDielectron*>(nextDie())) ){
 442     //AliInfo(Form(" **************** die->Process(InputEvent()) : %d **************************", idie));
 443     die->SetDontClearArrays(kTRUE);
 444     die->Process(InputEvent());
 445     if (die->HasCandidates()){
 446       Int_t ncandidates=die->GetPairArray(1)->GetEntriesFast();
 447       if (ncandidates==1) fEventStat->Fill((kNbinsEvent)+2*idie);
 448       else if (ncandidates>1) fEventStat->Fill((kNbinsEvent+1)+2*idie);
 449     }
 450
 451
 452     AliDielectronVarManager::Fill(InputEvent(), fgValues);
 453     for (Int_t ii=0; ii<2; ++ii){
 454       TObjArray *obj = (TObjArray*)die->GetTrackArray(ii);
 455       Int_t ntracks=obj->GetEntriesFast();
 456       //AliInfo(Form(" ************** # of tracks = %d", ntracks));
 457       for (Int_t itrack=0; itrack<ntracks; ++itrack){
 458
 459         ////////////////////////////////////////////////////////////////////
 460         AliDielectronVarManager::Fill(obj->UncheckedAt(itrack), fgValues);
 461         ////////////////////////////////////////////////////////////////////
 462
 463         if(fgValues[AliDielectronVarManager::kCharge]>0){
 464           fVeptmp.push_back(new  AliDielectronSingleTG(1,
 465                                                        fgValues[AliDielectronVarManager::kCentrality],
 466                                                        fgValues[AliDielectronVarManager::kXv],
 467                                                        fgValues[AliDielectronVarManager::kYv],
 468                                                        fgValues[AliDielectronVarManager::kZv],
 469                                                        fgValues[AliDielectronVarManager::kPx],
 470                                                        fgValues[AliDielectronVarManager::kPy],
 471                                                        fgValues[AliDielectronVarManager::kPz],
 472                                                        fgValues[AliDielectronVarManager::kPt],
 473                                                        fgValues[AliDielectronVarManager::kEta],
 474                                                        fgValues[AliDielectronVarManager::kPhi],
 475                                                        fgValues[AliDielectronVarManager::kTheta],
 476                                                        1, 1, static_cast<AliVTrack*>(obj->UncheckedAt(itrack)))
 477                             );
 478         }else if(fgValues[AliDielectronVarManager::kCharge]<0){
 479           fVemtmp.push_back(new  AliDielectronSingleTG(-1,
 480                                                        fgValues[AliDielectronVarManager::kCentrality],
 481                                                        fgValues[AliDielectronVarManager::kXv],
 482                                                        fgValues[AliDielectronVarManager::kYv],
 483                                                        fgValues[AliDielectronVarManager::kZv],
 484                                                        fgValues[AliDielectronVarManager::kPx],
 485                                                        fgValues[AliDielectronVarManager::kPy],
 486                                                        fgValues[AliDielectronVarManager::kPz],
 487                                                        fgValues[AliDielectronVarManager::kPt],
 488                                                        fgValues[AliDielectronVarManager::kEta],
 489                                                        fgValues[AliDielectronVarManager::kPhi],
 490                                                        fgValues[AliDielectronVarManager::kTheta],
 491                                                        1, 1, static_cast<AliVTrack*>(obj->UncheckedAt(itrack)))
 492                             );
 493         }
 494       }
 495     }
 496     //AliInfo(Form("size of e and p = %d %d", (int)fVeptmp.size(), (int)fVemtmp.size()));
 497
 498
 499     CheckGhostPairs(fVeptmp);
 500     CheckGhostPairs(fVemtmp);
 501     if(fRejectPairFlag[idie]==1 || fRejectPairFlag[idie]==2){
 502       RejectPairs(fVeptmp, fVemtmp, idie);
 503     }
 504     RandomizePool(fVeptmp, fVemtmp);
 505     CalcPair(fVep, fVem, die, idie);
 506
 507     //    AliInfo(Form("size of e and p (after) = %d %d", (int)fVep.size(), (int)fVem.size()));
 508
 509     double dwcent = 100.0/fgkNCENT;
 510     double dwiz = 20.0/fgkNZBIN;
 511     double dwrp = acos(-1.0)/fgkNRPBIN;
 512
 513     int icent = (int)(fgValues[AliDielectronVarManager::kCentrality]/dwcent);
 514     int izbin = (int)((fgValues[AliDielectronVarManager::kZvPrim]+10)/dwiz);
 515     int irp = (int)((fgValues[AliDielectronVarManager::kV0ACrpH2])/dwrp);
 516
 517     if(icent<0) icent=0;
 518     if(icent>=fgkNCENT) icent=fgkNCENT-1;
 519     if(izbin<0) izbin=0;
 520     if(izbin>=fgkNZBIN) izbin=fgkNZBIN-1;
 521     if(irp<0) irp=0;
 522     if(irp>=fgkNRPBIN) irp=fgkNRPBIN-1;
 523
 524     fvep[fibuf[idie][izbin][icent][irp]][idie][izbin][icent][irp].clear();
 525     for(int iep = 0; iep<(int)fVep.size();iep++) {
 526       fvep[fibuf[idie][izbin][icent][irp]][idie][izbin][icent][irp].push_back(fVep[iep]);
 527       fpoolp[idie][izbin][icent][irp].push_back(fVep[iep]);
 528       if(fpoolp[idie][izbin][icent][irp].size()>fgkMAXPOOL) {
 529         fpoolp[idie][izbin][icent][irp].pop_front();
 530       }
 531     }
 532     fvem[fibuf[idie][izbin][icent][irp]][idie][izbin][icent][irp].clear();
 533     for(int iem = 0; iem<(int)fVem.size();iem++) {
 534       fvem[fibuf[idie][izbin][icent][irp]][idie][izbin][icent][irp].push_back(fVem[iem]);
 535       fpoolm[idie][izbin][icent][irp].push_back(fVem[iem]);
 536       if(fpoolm[idie][izbin][icent][irp].size()>fgkMAXPOOL) {
 537         fpoolm[idie][izbin][icent][irp].pop_front();
 538       }
 539     }
 540
 541
 542     fibuf[idie][izbin][icent][irp]++;
 543     if(fibuf[idie][izbin][icent][irp]>= fgkNBUF) fibuf[idie][izbin][icent][irp]=0;
 544
 545
 546     fVeptmp.clear();
 547     fVemtmp.clear();
 548     fVep.clear();
 549     fVem.clear();
 550
 551
 552     ++idie;
 553   }
 554
 555
 556   AliESDEvent *fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 557   fEvent->Fill(values[AliDielectronVarManager::kCentrality]);
 558   for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
 559     AliESDtrack* track = fESD->GetTrack(iTracks);
 560     if (!track) {
 561       Printf("ERROR: Could not receive track %d", iTracks);
 562       continue;
 563     }
 564     if(!fCutsMother->AcceptTrack(track)) continue;
 565     fdEdXvsPt->Fill(track->GetTPCmomentum(), track->GetTPCsignal());
 566     fdEdXnSigmaElecvsPt->Fill(track->GetTPCmomentum(),
 567                               AliDielectronVarManager::GetPIDResponse()->NumberOfSigmasTPC(track,
 568                                                                                       AliPID::kElectron)
 569                               -AliDielectronPID::GetCorrVal());
 570     /// for beta caliculaton
 571     Double_t l = track->GetIntegratedLength();  // cm
 572     Double_t t = track->GetTOFsignal();
 573     Double_t t0 = AliDielectronVarManager::GetPIDResponse()->GetTOFResponse().GetTimeZero(); // ps
 574     Double_t beta = 0;
 575     if( (l < 360. || l > 800.) || (t <= 0.) || (t0 >999990.0) ) {
 576       beta=-9999;
 577     }
 578     else {
 579       t -= t0; // subtract the T0
 580       l *= 0.01;  // cm ->m
 581       t *= 1e-12; //ps -> s
 582
 583       Double_t v = l / t;
 584       beta = v / TMath::C();
 585     }
 586
 587     fTOFbetavsPt->Fill(track->GetTPCmomentum(), beta);
 588     fTOFnSigmaElecvsPt->Fill(track->GetTPCmomentum(),
 589                              AliDielectronVarManager::GetPIDResponse()->NumberOfSigmasTOF(track,
 590                                                                                      AliPID::kElectron));
 591     ////rough PID is required
 592     if( fabs(AliDielectronVarManager::GetPIDResponse()->NumberOfSigmasTOF(track, AliPID::kElectron))<3){
 593       fdEdXvsPtTOF->Fill(track->GetTPCmomentum(), track->GetTPCsignal());
 594       fdEdXnSigmaElecvsPtTOF->Fill(track->GetTPCmomentum(),
 595                                    AliDielectronVarManager::GetPIDResponse()->NumberOfSigmasTPC(track,
 596                                                                                                 AliPID::kElectron)
 597                                    -AliDielectronPID::GetCorrVal());
 598
 599
 600       if(track->GetTPCsignal()>70 && track->GetTPCsignal()<90){
 601         fNCrossedRowsTPC->Fill(track->GetTPCmomentum(), track->GetTPCCrossedRows());
 602         //fChi2ClusTPC->Fill(track->GetTPCmomentum(), track->GetTPCchi2()/track->GetTPCNcls());
 603         fChi2ClusTPC->Fill(track->GetTPCmomentum(), track->GetTPCchi2());
 604         fRatioCrossClusTPC->Fill(track->GetTPCmomentum(), track->GetTPCCrossedRows()/track->GetTPCNclsF());
 605       }
 606     }
 607   }
 608
 609   PostData(1, &fListHistos);
 610   PostData(2, &fListCF);
 611   PostData(3,fEventStat);
 612 }
 613
 614 //_________________________________________________________________________________
 615 void AliAnalysisTaskMultiDielectronTG::FinishTaskOutput()
 616 {
 617   //
 618   // Write debug tree
 619   //
 620   TIter nextDie(&fListDielectron);
 621   AliDielectron *die=0;
 622   while ( (die=static_cast<AliDielectron*>(nextDie())) ){
 623     die->SaveDebugTree();
 624     AliDielectronMixingHandler *mix=die->GetMixingHandler();
 625 //    printf("\n\n\n===============\ncall mix in Terminate: %p (%p)\n=================\n\n",mix,die);
 626     if (mix) mix->MixRemaining(die);
 627   }
 628   PostData(1, &fListHistos);
 629   PostData(2, &fListCF);
 630 }
 631
 632 //_________________________________________________________________________________
 633 void AliAnalysisTaskMultiDielectronTG::CheckGhostPairs(vector<AliDielectronSingleTG*> e1)
 634 {
 635   ////// Check whether the pairs are in ghost pairs
 636   /////  ghost means that one of the pairs is highly associated but reconstructed as a true particle
 637
 638
 639   bool reject = false;
 640   if(e1.size()>1){
 641     for(int i1=0; i1<(int)e1.size(); i1++){
 642       reject = false;
 643       for(int i2=i1+1; i2<(int)e1.size(); i2++){
 644         if( fabs(e1[i1]->Phi() - e1[i2]->Phi())<0.01 &&
 645             fabs(e1[i1]->Eta() - e1[i2]->Eta())<0.005
 646             ){
 647           reject = true;
 648           e1[i2]->SetGstFlag(0);
 649         }
 650       }
 651       if(reject==true)e1[i1]->SetGstFlag(0);
 652     }
 653   }
 654 }
 655
 656 //_________________________________________________________________________________
 657 Bool_t AliAnalysisTaskMultiDielectronTG::CheckGhost(vector<AliDielectronSingleTG*> e1, vector<AliDielectronSingleTG*> e2)
 658 {
 659   ////// To be sure whether there are no ghost pairs in h event mixing
 660
 661   if(e1.size()>0 && e2.size()>0){
 662     for(int i1=0; i1<(int)e1.size(); i1++){
 663       for(int i2=0; i2<(int)e2.size(); i2++){
 664         if( fabs(e1[i1]->Phi() - e2[i2]->Phi())<0.01 &&
 665             fabs(e1[i1]->Eta() - e2[i2]->Eta())<0.005
 666             ){
 667           return true;
 668         }
 669       }
 670     }
 671   }
 672   return false;
 673 }
 674
 675 //_________________________________________________________________________________
 676 void AliAnalysisTaskMultiDielectronTG::RandomizePool(vector<AliDielectronSingleTG*> e1, vector<AliDielectronSingleTG*> e2)
 677 {
 678
 679   ///// Randomize the pool constent to cancel the filling scheme of the single tracks
 680   ////
 681
 682
 683   int size1 = e1.size();
 684   int usedindex[1000];
 685   for(int i=0;i<1000;i++){
 686     usedindex[i] = -1;
 687   }
 688   for(int i=0;i<size1;i++){
 689     usedindex[i] = 0;
 690   }
 691
 692   for(int i=0;i<size1;i++){
 693     int j = (int)(gRandom->Uniform(0,size1));
 694     while(usedindex[j]==1){
 695       j = (int)(gRandom->Uniform(0,size1));
 696     }
 697     if( (e1[j]->GetGstFlag()==1) &&
 698         (e1[j]->GetConvFlag()==1)
 699         ){
 700       fVep.push_back(e1[j]);
 701     }
 702     usedindex[j] = 1;
 703   }
 704
 705
 706   int size2 = e2.size();
 707   for(int i=0;i<1000;i++){
 708     usedindex[i] = -1;
 709   }
 710   for(int i=0;i<size2;i++){
 711     usedindex[i] = 0;
 712   }
 713
 714   for(int i=0;i<size2;i++){
 715     int j = (int)(gRandom->Uniform(0,size2));
 716     while(usedindex[j]==1){
 717       j = (int)(gRandom->Uniform(0,size2));
 718     }
 719     if( (e2[j]->GetGstFlag()==1) &&
 720         (e2[j]->GetConvFlag()==1)
 721         ){
 722       fVem.push_back(e2[j]);
 723     }
 724     usedindex[j] = 1;
 725   }
 726 }
 727
 728
 729 //_________________________________________________________________________________
 730 void AliAnalysisTaskMultiDielectronTG::CalcPair(vector<AliDielectronSingleTG*> ve1,
 731                                                 vector<AliDielectronSingleTG*> ve2, AliDielectron *die, Int_t idie)
 732 {
 733
 734   //
 735   // main routine for the pair procesing
 736   //
 737
 738
 739   for(int i1=0; i1<(int)ve1.size(); i1++){
 740     for(int i2=0; i2<(int)ve2.size(); i2++){
 741       FillPair(ve1[i1], ve2[i2], 0, die, idie);
 742     }
 743   }
 744
 745   for(int i1=0; i1<(int)ve1.size(); i1++){
 746     for(int i2=i1+1; i2<(int)ve1.size(); i2++){
 747       FillPair(ve1[i1], ve1[i2], 1, die, idie);
 748     }
 749   }
 750
 751   for(int i1=0; i1<(int)ve2.size(); i1++){
 752     for(int i2=i1+1; i2<(int)ve2.size(); i2++){
 753       FillPair(ve2[i1], ve2[i2], 2, die, idie );
 754     }
 755   }
 756
 757
 758   double dwcent = 100.0/fgkNCENT;
 759   double dwiz = 20.0/fgkNZBIN;
 760   double dwrp = acos(-1.0)/fgkNRPBIN;
 761
 762   int icent = (int)(fgValues[AliDielectronVarManager::kCentrality]/dwcent);
 763   int izbin = (int)((fgValues[AliDielectronVarManager::kZvPrim]+10)/dwiz);
 764   int irp = (int)((fgValues[AliDielectronVarManager::kV0ACrpH2])/dwrp);
 765
 766   if(icent<0) icent=0;
 767   if(icent>=fgkNCENT) icent=fgkNCENT-1;
 768   if(izbin<0) izbin=0;
 769   if(izbin>=fgkNZBIN) izbin=fgkNZBIN-1;
 770   if(irp<0) irp=0;
 771   if(irp>=fgkNRPBIN) irp=fgkNRPBIN-1;
 772
 773   int nmixed;
 774   if(ve1.size()>0) {
 775     //
 776     // Now mixed event for +- pairs
 777     //
 778     nmixed = 0;
 779     for(int ibuf=0;(nmixed<fgkNMix);ibuf++) {
 780       int ntry = 0;
 781       while((fBGRejUnlike && CheckGhost(ve1, fvem[ibuf][idie][izbin][icent][irp])) &&  ntry<fgkMAXTRY) {
 782         ReshuffleBuffer(fvem[ibuf][idie][izbin][icent][irp],fpoolm[idie][izbin][icent][irp]);
 783         ntry++;
 784       }
 785       for(int i1=0; i1<(int)ve1.size(); i1++){
 786         for(int i2=0; i2<(int)fvem[ibuf][idie][izbin][icent][irp].size(); i2++){
 787           FillPair(ve1[i1],fvem[ibuf][idie][izbin][icent][irp][i2], 3, die, idie);
 788         }
 789       }
 790       ++nmixed;
 791     }//for(ibuf)
 792   }
 793   if(ve2.size()>0) {
 794     //
 795     // Now mixed event for -+ pairs
 796     //
 797     nmixed = 0;
 798     for(int ibuf=0;(nmixed<fgkNMix);ibuf++) {
 799       int ntry = 0;
 800       while((fBGRejUnlike && CheckGhost(ve2, fvep[ibuf][idie][izbin][icent][irp])) &&  ntry<fgkMAXTRY) {
 801         ReshuffleBuffer(fvep[ibuf][idie][izbin][icent][irp],fpoolp[idie][izbin][icent][irp]);
 802         ntry++;
 803       }
 804       for(int i1=0; i1<(int)ve2.size(); i1++){
 805         for(int i2=0; i2<(int)fvep[ibuf][idie][izbin][icent][irp].size(); i2++){
 806           FillPair(fvep[ibuf][idie][izbin][icent][irp][i2],ve2[i1],4, die, idie);
 807         }
 808       }
 809       ++nmixed;
 810     }//for(ibuf)
 811   }
 812
 813   if(ve1.size()>0) {
 814     //
 815     // Now mixed event for ++ pairs
 816     //
 817     nmixed = 0;
 818     for(int ibuf=0;(nmixed<fgkNMix);ibuf++) {
 819       int ntry = 0;
 820       while((fBGRejLike && CheckGhost(ve1, fvep[ibuf][idie][izbin][icent][irp])) &&  ntry<fgkMAXTRY) {
 821         ReshuffleBuffer(fvep[ibuf][idie][izbin][icent][irp],fpoolp[idie][izbin][icent][irp]);
 822         ntry++;
 823       }
 824       for(int i1=0; i1<(int)ve1.size(); i1++){
 825         for(int i2=0;i2<(int)fvep[ibuf][idie][izbin][icent][irp].size(); i2++){
 826           FillPair(ve1[i1],fvep[ibuf][idie][izbin][icent][irp][i2], 5, die, idie);
 827         }
 828       }
 829       ++nmixed;
 830     }//for(ibuf)
 831   }
 832
 833   if(ve2.size()>0) {
 834     //
 835     // Now mixed event for +- pairs
 836     //
 837     nmixed = 0;
 838     for(int ibuf=0;(nmixed<fgkNMix);ibuf++) {
 839       int ntry = 0;
 840       while((fBGRejLike && CheckGhost(ve2, fvem[ibuf][idie][izbin][icent][irp])) &&  ntry<fgkMAXTRY) {
 841         ReshuffleBuffer(fvem[ibuf][idie][izbin][icent][irp],fpoolm[idie][izbin][icent][irp]);
 842         ntry++;
 843       }
 844       for(int i1=0; i1<(int)ve2.size(); i1++){
 845         for(int i2=0; i2<(int)fvem[ibuf][idie][izbin][icent][irp].size(); i2++){
 846           FillPair(ve2[i1],fvem[ibuf][idie][izbin][icent][irp][i2],6, die, idie);
 847         }
 848       }
 849       ++nmixed;
 850     }//for(ibuf)
 851   }
 852
 853 }
 854
 855
 856 //_________________________________________________________________________________
 857 void AliAnalysisTaskMultiDielectronTG::FillPair(AliDielectronSingleTG *iep,
 858                                                 AliDielectronSingleTG *iem, int type, AliDielectron *die, Int_t idie)
 859 {
 860
 861   //
 862   // main routine for filling kinematics of pairs
 863   //
 864
 865
 866
 867   double dmass, dphiv, dpxpair, dpypair, dpzpair;
 868   double dptpair, depair, dphipair, detapair, dcos, dpsi;
 869
 870   CalcVars(iep, iem, dmass, dphiv, dpxpair, dpypair, dpzpair,
 871             dptpair, depair, dphipair, detapair, dcos, dpsi);
 872
 873
 874   double dopeningangle =  -9999;
 875   double dv0mass =  -9999;
 876   double dv0pxpair = -9999;
 877   double dv0pypair = -9999;
 878   double dv0pzpair = -9999;
 879   double dv0ptpair = -9999;
 880   double dv0epair = -9999;
 881   double dv0xvpair = -9999;
 882   double dv0yvpair = -9999;
 883   double dv0zvpair = -9999;
 884   double dv0phipair = -9999;
 885   double dv0etapair = -9999;
 886   double dv0rpair = -9999;
 887   double dpsipair =  -9999;
 888   double dphivpair =  -9999;
 889
 890   ////////////////////////////
 891   ///// calculate v0 ////////
 892   ///////////////////////////
 893   Bool_t v0OFF=kFALSE;
 894   /// for the moment, this doesn't work for MC
 895   if(fdv0mixing == kFALSE && (type==3 || type==4 || type==5 || type==6)){
 896     v0OFF = kTRUE;
 897   }
 898   if(die->GetHasMC()==kTRUE && (type==3 || type==4 || type==5 || type==6)){
 899     v0OFF = kTRUE;
 900   }
 901   if(type==0 || type==1 || type==2){
 902     v0OFF = kFALSE;
 903   }
 904
 905
 906   if(v0OFF==kFALSE){
 907     AliVTrack* trackob1= iep->GetTrack();
 908     AliVTrack* trackob2= iem->GetTrack();
 909
 910     if(!trackob1 || !trackob2){
 911       return;
 912     }
 913
 914     AliDielectronPair candidate;
 915     candidate.SetTracks(trackob1, (int)(11*iep->Charge()),
 916                         trackob2, (int)(11*iem->Charge()));
 917
 918     if(trackob1==trackob2){
 919       AliInfo("Same Objects!!");
 920       return;
 921     }
 922
 923     //const AliKFParticle &kfPairLeg1 = candidate.GetKFFirstDaughter();
 924     //const AliKFParticle &kfPairLeg2 = candidate.GetKFSecondDaughter();
 925
 926     const AliKFParticle &kfPair = candidate.GetKFParticle();
 927
 928     /*
 929     dv0mass = candidate.M();
 930     dv0pxpair = candidate.Px();
 931     dv0pypair = candidate.Py();
 932     dv0pzpair = candidate.Pz();
 933     dv0ptpair = candidate.Pt();
 934     dv0epair = candidate.E();
 935     dv0xvpair = candidate.Xv();
 936     dv0yvpair = candidate.Yv();
 937     dv0zvpair = candidate.Zv();
 938     dv0phipair = candidate.Phi();
 939     dv0etapair = candidate.Eta();
 940     */
 941     dv0mass = kfPair.GetMass();
 942     dv0pxpair = kfPair.GetPx();
 943     dv0pypair = kfPair.GetPy();
 944     dv0pzpair = kfPair.GetPz();
 945     dv0ptpair = kfPair.GetPt();
 946     dv0epair = kfPair.GetE();
 947     dv0xvpair = kfPair.GetX();
 948     dv0yvpair = kfPair.GetY();
 949     dv0zvpair = kfPair.GetZ();
 950     dv0phipair = kfPair.GetPhi();
 951     dv0etapair = kfPair.GetEta();
 952     dv0rpair =  kfPair.GetR();
 953
 954     dopeningangle = candidate.OpeningAngle();
 955     dpsipair = candidate.PsiPair(fbz);
 956     dphivpair = candidate.PhivPair(fbz);
 957
 958
 959   }
 960
 961   Double_t values[AliDielectronVarManager::kNMaxValues];
 962   TString  className1;
 963   TString  className2;
 964   className1.Form("MyPair_%s",kPairClassNamesTG[type]);
 965   className2.Form("MyPairV0_%s",kPairClassNamesTG[type]);
 966
 967   AliDielectronHistos *fHistos = die->GetHistoManager();
 968   Bool_t pairClass1=fHistos->GetHistogramList()->FindObject(className1.Data())!=0x0;
 969   Bool_t pairClass2=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
 970
 971   if (pairClass1 && PairTrackcut(dphiv, dcos, idie)==true){
 972     ///import pair variables to values!!
 973     values[AliDielectronVarManager::kPx] = dpxpair;
 974     values[AliDielectronVarManager::kPy] = dpypair;
 975     values[AliDielectronVarManager::kPz] = dpzpair;
 976     values[AliDielectronVarManager::kPt] = dptpair;
 977     values[AliDielectronVarManager::kXv] = dv0xvpair;
 978     values[AliDielectronVarManager::kYv] = dv0yvpair;
 979     values[AliDielectronVarManager::kZv] = dv0zvpair;
 980     values[AliDielectronVarManager::kR] = dv0rpair;
 981     values[AliDielectronVarManager::kE] = depair;
 982     values[AliDielectronVarManager::kEta] = detapair;
 983     values[AliDielectronVarManager::kM] = dmass;
 984     values[AliDielectronVarManager::kPsiPair] = dphiv;
 985     values[AliDielectronVarManager::kPhivPair] = dphiv;
 986     values[AliDielectronVarManager::kPhi]  = dphipair;
 987     values[AliDielectronVarManager::kOpeningAngle]  = dcos;
 988     values[AliDielectronVarManager::kCosPointingAngle]  = TMath::Abs(TMath::ATan2(TMath::Sin(iep->Phi()-iem->Phi()),TMath::Cos(iep->Phi()-iem->Phi())));
 989     fHistos->FillClass(className1, AliDielectronVarManager::kNMaxValues, values);
 990   }
 991
 992
 993   if (pairClass2 && PairTrackcut(dphiv, dopeningangle, idie)==true){
 994     values[AliDielectronVarManager::kPx] = dv0pxpair;
 995     values[AliDielectronVarManager::kPy] = dv0pypair;
 996     values[AliDielectronVarManager::kPz] = dv0pzpair;
 997     values[AliDielectronVarManager::kPt] = dv0ptpair;
 998     values[AliDielectronVarManager::kXv] = dv0xvpair;
 999     values[AliDielectronVarManager::kYv] = dv0yvpair;
1000     values[AliDielectronVarManager::kZv] = dv0zvpair;
1001     values[AliDielectronVarManager::kR] = dv0rpair;
1002     values[AliDielectronVarManager::kE] = dv0epair;
1003     values[AliDielectronVarManager::kEta] = dv0etapair;
1004     values[AliDielectronVarManager::kM] = dv0mass;
1005     values[AliDielectronVarManager::kPsiPair] = dpsipair;
1006     values[AliDielectronVarManager::kPhivPair] = dphivpair;
1007     values[AliDielectronVarManager::kPhi]  = dv0phipair;
1008     values[AliDielectronVarManager::kOpeningAngle]  = dopeningangle;
1009     values[AliDielectronVarManager::kCosPointingAngle]  = TMath::Abs(TMath::ATan2(TMath::Sin(iep->Phi()-iem->Phi()),TMath::Cos(iep->Phi()-iem->Phi())));
1010     fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
1011   }
1012
1013
1014
1015 }
1016
1017 //_________________________________________________________________________________
1018 bool AliAnalysisTaskMultiDielectronTG::PairTrackcut(double phiv, double op, int idie)
1019 {
1020
1021   //
1022   // pair-by-pair cuts
1023   //
1024
1025
1026   bool pairOK = true;
1027
1028   if(fRejectPairFlag[idie] == 1 || fRejectPairFlag[idie] == 2 ||
1029      fRejectPairFlag[idie] == 3 || fRejectPairFlag[idie] == 4 ){
1030     if(fRejectPairFlag[idie] == 2 || fRejectPairFlag[idie] == 4 ){
1031       if(fbz>0 && phiv>fdconvphiv){
1032         pairOK = false;
1033       }else if(fbz<0 && phiv<acos(-1.0)-fdconvphiv){
1034         pairOK = false;
1035       }
1036     }else if(fRejectPairFlag[idie] == 1 || fRejectPairFlag[idie] == 3){
1037       if(op<fdop){
1038         pairOK = false;
1039       }
1040     }
1041   }
1042
1043   return pairOK;
1044
1045 }
1046
1047
1048 //_________________________________________________________________________________
1049 void AliAnalysisTaskMultiDielectronTG::CalcVars(AliDielectronSingleTG *iep, AliDielectronSingleTG *iem,
1050                                                 double &mass, double &phiv, double &px, double &py, double&pz,
1051                                                 double &pt, double &e, double &phi,
1052                                                 double &eta, double &cos, double &psi)
1053 {
1054
1055
1056   //
1057   // standalone calculator for the pair variables
1058   //
1059
1060   px = iep->Px()+iem->Px();
1061   py = iep->Py()+iem->Py();
1062   pz = iep->Pz()+iem->Pz();
1063   pt = sqrt(px*px+py*py);
1064   double dppair = sqrt(pt*pt+pz*pz);
1065   static const double me=0.0005109989;
1066   e = sqrt(me*me+iep->Px()*iep->Px()+iep->Py()*iep->Py()+iep->Pz()*iep->Pz())
1067     + sqrt(me*me+iem->Px()*iem->Px()+iem->Py()*iem->Py()+iem->Pz()*iem->Pz());
1068
1069   mass =  e*e-px*px-py*py-pz*pz;
1070   if(mass<0){
1071     mass = mass;
1072   }else{
1073     mass = sqrt(mass);
1074   }
1075
1076
1077   phi = atan2(py, px);
1078   eta = -0.5*TMath::Log((dppair+pz)/(dppair-pz));
1079   double p1 = sqrt(pow(iep->Px(),2)+pow(iep->Py(),2)+pow(iep->Pz(),2));
1080   double p2 = sqrt(pow(iem->Px(),2)+pow(iem->Py(),2)+pow(iem->Pz(),2));
1081   cos = acos((iep->Px()*iem->Px()+iep->Py()*iem->Py()+iep->Pz()*iem->Pz())/(p1*p2));
1082
1083   double dtheta = iep->Theta()-iem->Theta();
1084   psi = asin(dtheta/cos);
1085
1086
1087   //unit vector of (pep+pem)
1088   double pl = dppair;
1089   double ux = px/pl;
1090   double uy = py/pl;
1091   double uz = pz/pl;
1092   double ax = uy/sqrt(ux*ux+uy*uy);
1093   double ay = -ux/sqrt(ux*ux+uy*uy);
1094
1095   //momentum of e+ and e- in (ax,ay,az) axis. Note that az=0 by
1096   //definition.
1097   //double ptep = iep->Px()*ax + iep->Py()*ay;
1098   //double ptem = iem->Px()*ax + iem->Py()*ay;
1099
1100   double pxep = iep->Px();
1101   double pyep = iep->Py();
1102   double pzep = iep->Pz();
1103   double pxem = iem->Px();
1104   double pyem = iem->Py();
1105   double pzem = iem->Pz();
1106
1107
1108   //vector product of pep X pem
1109   double vpx = pyep*pzem - pzep*pyem;
1110   double vpy = pzep*pxem - pxep*pzem;
1111   double vpz = pxep*pyem - pyep*pxem;
1112   double vp = sqrt(vpx*vpx+vpy*vpy+vpz*vpz);
1113   //double thev = acos(vpz/vp);
1114
1115   //unit vector of pep X pem
1116   double vx = vpx/vp;
1117   double vy = vpy/vp;
1118   double vz = vpz/vp;
1119
1120   //The third axis defined by vector product (ux,uy,uz)X(vx,vy,vz)
1121   double wx = uy*vz - uz*vy;
1122   double wy = uz*vx - ux*vz;
1123   double wz = ux*vy - uy*vx;
1124   double wl = sqrt(wx*wx+wy*wy+wz*wz);
1125   // by construction, (wx,wy,wz) must be a unit vector.
1126   if(fabs(wl - 1.0) > 0.00001) std::cout << "Calculation error in W vector"<<std::endl;
1127   // measure angle between (wx,wy,wz) and (ax,ay,0). The angle between them
1128   // should be small if the pair is conversion
1129   //
1130   double cosPhiV = wx*ax + wy*ay;
1131   phiv = acos(cosPhiV);
1132
1133 }
1134
1135 //_________________________________________________________________________________
1136 void AliAnalysisTaskMultiDielectronTG::ReshuffleBuffer(vector<AliDielectronSingleTG*> ve, deque<AliDielectronSingleTG*> pool)
1137 {
1138   //If there is not enough electron in the pool, give up
1139   //
1140   // ReshuffleBuffer for th event mixing
1141   //
1142
1143   unsigned int ne = ve.size();
1144   unsigned int poolsize = pool.size();
1145   int used[fgkMAXPOOL];
1146   for(int i=0;i<(int)fgkMAXPOOL;i++){
1147     used[i]=0;
1148   }
1149
1150   if(poolsize < ne) {
1151     std::cout <<" pool size="<<poolsize<<" ne"<<ne<<std::endl;
1152     return;
1153   }
1154   for(unsigned int ie=0; ie < ne; ie++) {
1155     int j = rand()%poolsize;
1156     while(used[j]==1){
1157       j = rand()%poolsize;
1158     }
1159     ve[ie] = pool[j];
1160     used[j]=1;
1161   }
1162
1163 }
1164
1165 //_________________________________________________________________________________
1166 void AliAnalysisTaskMultiDielectronTG::RejectPairs(vector<AliDielectronSingleTG*> e1, vector<AliDielectronSingleTG*> e2, Int_t idie){
1167
1168   ////////////////////////////////////
1169   ///// to reject pairs from track arrays
1170   ///// conversions, ghost ..
1171   ///////////////////////////////////
1172   if(e1.size()>0 && e2.size()>0){
1173     for(int i1=0; i1<(int)e1.size(); i1++){
1174       for(int i2=0; i2<(int)e2.size(); i2++){
1175         if(fRejectPairFlag[idie]==1){
1176           Double_t cosine = GetOpeningAngle(e1[i1], e2[i2]);
1177           if(cosine<fdop){
1178             e1[i1]->SetConvFlag(0);
1179             e2[i2]->SetConvFlag(0);
1180           }
1181         }else if(fRejectPairFlag[idie]==2){
1182           Double_t phiv = GetPhiv(e1[i1], e2[i2]);
1183           if(fbz>0 && phiv>fdconvphiv){
1184             e1[i1]->SetConvFlag(0);
1185             e2[i2]->SetConvFlag(0);
1186           }else if(fbz<0 && phiv<acos(-1.0)-fdconvphiv){
1187             e1[i1]->SetConvFlag(0);
1188             e2[i2]->SetConvFlag(0);
1189           }
1190         }
1191       }
1192     }
1193   }
1194   /////// this is really necessary??????
1195   if(e1.size()>0){
1196     for(int i1=0; i1<(int)e1.size(); i1++){
1197       for(int i2=i1+1; i2<(int)e1.size(); i2++){
1198         if(fRejectPairFlag[idie]==1){
1199           Double_t cosine = GetOpeningAngle(e1[i1], e1[i2]);
1200           if(cosine<fdop){
1201             e1[i1]->SetConvFlag(0);
1202             e1[i2]->SetConvFlag(0);
1203           }
1204         }else if(fRejectPairFlag[idie]==2){
1205           Double_t phiv = GetPhiv(e1[i1], e1[i2]);
1206           if(fbz>0 && phiv>fdconvphiv){
1207             e1[i1]->SetConvFlag(0);
1208             e1[i2]->SetConvFlag(0);
1209           }else if(fbz<0 && phiv<acos(-1.0)-fdconvphiv){
1210             e1[i1]->SetConvFlag(0);
1211             e1[i2]->SetConvFlag(0);
1212           }
1213         }
1214       }
1215     }
1216   }
1217
1218   if(e2.size()>0){
1219     for(int i1=0; i1<(int)e2.size(); i1++){
1220       for(int i2=i1+1; i2<(int)e2.size(); i2++){
1221         if(fRejectPairFlag[idie]==1){
1222           Double_t cosine = GetOpeningAngle(e2[i1], e2[i2]);
1223           if(cosine<fdop){
1224             e2[i1]->SetConvFlag(0);
1225             e2[i2]->SetConvFlag(0);
1226           }
1227         }else if(fRejectPairFlag[idie]==2){
1228           Double_t phiv = GetPhiv(e2[i1], e2[i2]);
1229           if(fbz>0 && phiv>fdconvphiv){
1230             e2[i1]->SetConvFlag(0);
1231             e2[i2]->SetConvFlag(0);
1232           }else if(fbz<0 && phiv<acos(-1.0)-fdconvphiv){
1233             e2[i1]->SetConvFlag(0);
1234             e2[i2]->SetConvFlag(0);
1235           }
1236         }
1237       }
1238     }
1239   }
1240 }
1241
1242
1243 //_________________________________________________________________________________
1244 Double_t AliAnalysisTaskMultiDielectronTG::GetOpeningAngle(AliDielectronSingleTG* e1, AliDielectronSingleTG* e2){
1245
1246   //////////////////////
1247   //////// calculate pairs and get opening angle
1248   //////////////////////
1249
1250   double dmass, dphiv, dpxpair, dpypair, dpzpair;
1251   double dptpair, depair, dphipair, detapair, dcos, dpsi;
1252
1253   CalcVars(e1, e2, dmass, dphiv, dpxpair, dpypair, dpzpair,
1254             dptpair, depair, dphipair, detapair, dcos, dpsi);
1255
1256   return dcos;
1257 }
1258
1259 //_________________________________________________________________________________
1260 Double_t AliAnalysisTaskMultiDielectronTG::GetPhiv(AliDielectronSingleTG* e1, AliDielectronSingleTG* e2){
1261
1262   //////////////////////
1263   //////// calculate pairs and get opening angle
1264   //////////////////////
1265
1266   double dmass, dphiv, dpxpair, dpypair, dpzpair;
1267   double dptpair, depair, dphipair, detapair, dcos, dpsi;
1268
1269   CalcVars(e1, e2, dmass, dphiv, dpxpair, dpypair, dpzpair,
1270             dptpair, depair, dphipair, detapair, dcos, dpsi);
1271
1272   return dphiv;
1273 }