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