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) continue;
 626     for (Int_t ipool=0; ipool<mix->GetNumberOfBins(); ++ipool){
 627       mix->MixRemaining(die, ipool);
 628     }
 629   }
 630   PostData(1, &fListHistos);
 631   PostData(2, &fListCF);
 632 }
 633
 634 //_________________________________________________________________________________
 635 void AliAnalysisTaskMultiDielectronTG::CheckGhostPairs(vector<AliDielectronSingleTG*> e1)
 636 {
 637   ////// Check whether the pairs are in ghost pairs
 638   /////  ghost means that one of the pairs is highly associated but reconstructed as a true particle
 639
 640
 641   bool reject = false;
 642   if(e1.size()>1){
 643     for(int i1=0; i1<(int)e1.size(); i1++){
 644       reject = false;
 645       for(int i2=i1+1; i2<(int)e1.size(); i2++){
 646         if( fabs(e1[i1]->Phi() - e1[i2]->Phi())<0.01 &&
 647             fabs(e1[i1]->Eta() - e1[i2]->Eta())<0.005
 648             ){
 649           reject = true;
 650           e1[i2]->SetGstFlag(0);
 651         }
 652       }
 653       if(reject==true)e1[i1]->SetGstFlag(0);
 654     }
 655   }
 656 }
 657
 658 //_________________________________________________________________________________
 659 Bool_t AliAnalysisTaskMultiDielectronTG::CheckGhost(vector<AliDielectronSingleTG*> e1, vector<AliDielectronSingleTG*> e2)
 660 {
 661   ////// To be sure whether there are no ghost pairs in h event mixing
 662
 663   if(e1.size()>0 && e2.size()>0){
 664     for(int i1=0; i1<(int)e1.size(); i1++){
 665       for(int i2=0; i2<(int)e2.size(); i2++){
 666         if( fabs(e1[i1]->Phi() - e2[i2]->Phi())<0.01 &&
 667             fabs(e1[i1]->Eta() - e2[i2]->Eta())<0.005
 668             ){
 669           return true;
 670         }
 671       }
 672     }
 673   }
 674   return false;
 675 }
 676
 677 //_________________________________________________________________________________
 678 void AliAnalysisTaskMultiDielectronTG::RandomizePool(vector<AliDielectronSingleTG*> e1, vector<AliDielectronSingleTG*> e2)
 679 {
 680
 681   ///// Randomize the pool constent to cancel the filling scheme of the single tracks
 682   ////
 683
 684
 685   int size1 = e1.size();
 686   int usedindex[1000];
 687   for(int i=0;i<1000;i++){
 688     usedindex[i] = -1;
 689   }
 690   for(int i=0;i<size1;i++){
 691     usedindex[i] = 0;
 692   }
 693
 694   for(int i=0;i<size1;i++){
 695     int j = (int)(gRandom->Uniform(0,size1));
 696     while(usedindex[j]==1){
 697       j = (int)(gRandom->Uniform(0,size1));
 698     }
 699     if( (e1[j]->GetGstFlag()==1) &&
 700         (e1[j]->GetConvFlag()==1)
 701         ){
 702       fVep.push_back(e1[j]);
 703     }
 704     usedindex[j] = 1;
 705   }
 706
 707
 708   int size2 = e2.size();
 709   for(int i=0;i<1000;i++){
 710     usedindex[i] = -1;
 711   }
 712   for(int i=0;i<size2;i++){
 713     usedindex[i] = 0;
 714   }
 715
 716   for(int i=0;i<size2;i++){
 717     int j = (int)(gRandom->Uniform(0,size2));
 718     while(usedindex[j]==1){
 719       j = (int)(gRandom->Uniform(0,size2));
 720     }
 721     if( (e2[j]->GetGstFlag()==1) &&
 722         (e2[j]->GetConvFlag()==1)
 723         ){
 724       fVem.push_back(e2[j]);
 725     }
 726     usedindex[j] = 1;
 727   }
 728 }
 729
 730
 731 //_________________________________________________________________________________
 732 void AliAnalysisTaskMultiDielectronTG::CalcPair(vector<AliDielectronSingleTG*> ve1,
 733                                                 vector<AliDielectronSingleTG*> ve2, AliDielectron *die, Int_t idie)
 734 {
 735
 736   //
 737   // main routine for the pair procesing
 738   //
 739
 740
 741   for(int i1=0; i1<(int)ve1.size(); i1++){
 742     for(int i2=0; i2<(int)ve2.size(); i2++){
 743       FillPair(ve1[i1], ve2[i2], 0, die, idie);
 744     }
 745   }
 746
 747   for(int i1=0; i1<(int)ve1.size(); i1++){
 748     for(int i2=i1+1; i2<(int)ve1.size(); i2++){
 749       FillPair(ve1[i1], ve1[i2], 1, die, idie);
 750     }
 751   }
 752
 753   for(int i1=0; i1<(int)ve2.size(); i1++){
 754     for(int i2=i1+1; i2<(int)ve2.size(); i2++){
 755       FillPair(ve2[i1], ve2[i2], 2, die, idie );
 756     }
 757   }
 758
 759
 760   double dwcent = 100.0/fgkNCENT;
 761   double dwiz = 20.0/fgkNZBIN;
 762   double dwrp = acos(-1.0)/fgkNRPBIN;
 763
 764   int icent = (int)(fgValues[AliDielectronVarManager::kCentrality]/dwcent);
 765   int izbin = (int)((fgValues[AliDielectronVarManager::kZvPrim]+10)/dwiz);
 766   int irp = (int)((fgValues[AliDielectronVarManager::kV0ACrpH2])/dwrp);
 767
 768   if(icent<0) icent=0;
 769   if(icent>=fgkNCENT) icent=fgkNCENT-1;
 770   if(izbin<0) izbin=0;
 771   if(izbin>=fgkNZBIN) izbin=fgkNZBIN-1;
 772   if(irp<0) irp=0;
 773   if(irp>=fgkNRPBIN) irp=fgkNRPBIN-1;
 774
 775   int nmixed;
 776   if(ve1.size()>0) {
 777     //
 778     // Now mixed event for +- pairs
 779     //
 780     nmixed = 0;
 781     for(int ibuf=0;(nmixed<fgkNMix);ibuf++) {
 782       int ntry = 0;
 783       while((fBGRejUnlike && CheckGhost(ve1, fvem[ibuf][idie][izbin][icent][irp])) &&  ntry<fgkMAXTRY) {
 784         ReshuffleBuffer(fvem[ibuf][idie][izbin][icent][irp],fpoolm[idie][izbin][icent][irp]);
 785         ntry++;
 786       }
 787       for(int i1=0; i1<(int)ve1.size(); i1++){
 788         for(int i2=0; i2<(int)fvem[ibuf][idie][izbin][icent][irp].size(); i2++){
 789           FillPair(ve1[i1],fvem[ibuf][idie][izbin][icent][irp][i2], 3, die, idie);
 790         }
 791       }
 792       ++nmixed;
 793     }//for(ibuf)
 794   }
 795   if(ve2.size()>0) {
 796     //
 797     // Now mixed event for -+ pairs
 798     //
 799     nmixed = 0;
 800     for(int ibuf=0;(nmixed<fgkNMix);ibuf++) {
 801       int ntry = 0;
 802       while((fBGRejUnlike && CheckGhost(ve2, fvep[ibuf][idie][izbin][icent][irp])) &&  ntry<fgkMAXTRY) {
 803         ReshuffleBuffer(fvep[ibuf][idie][izbin][icent][irp],fpoolp[idie][izbin][icent][irp]);
 804         ntry++;
 805       }
 806       for(int i1=0; i1<(int)ve2.size(); i1++){
 807         for(int i2=0; i2<(int)fvep[ibuf][idie][izbin][icent][irp].size(); i2++){
 808           FillPair(fvep[ibuf][idie][izbin][icent][irp][i2],ve2[i1],4, die, idie);
 809         }
 810       }
 811       ++nmixed;
 812     }//for(ibuf)
 813   }
 814
 815   if(ve1.size()>0) {
 816     //
 817     // Now mixed event for ++ pairs
 818     //
 819     nmixed = 0;
 820     for(int ibuf=0;(nmixed<fgkNMix);ibuf++) {
 821       int ntry = 0;
 822       while((fBGRejLike && CheckGhost(ve1, fvep[ibuf][idie][izbin][icent][irp])) &&  ntry<fgkMAXTRY) {
 823         ReshuffleBuffer(fvep[ibuf][idie][izbin][icent][irp],fpoolp[idie][izbin][icent][irp]);
 824         ntry++;
 825       }
 826       for(int i1=0; i1<(int)ve1.size(); i1++){
 827         for(int i2=0;i2<(int)fvep[ibuf][idie][izbin][icent][irp].size(); i2++){
 828           FillPair(ve1[i1],fvep[ibuf][idie][izbin][icent][irp][i2], 5, die, idie);
 829         }
 830       }
 831       ++nmixed;
 832     }//for(ibuf)
 833   }
 834
 835   if(ve2.size()>0) {
 836     //
 837     // Now mixed event for +- pairs
 838     //
 839     nmixed = 0;
 840     for(int ibuf=0;(nmixed<fgkNMix);ibuf++) {
 841       int ntry = 0;
 842       while((fBGRejLike && CheckGhost(ve2, fvem[ibuf][idie][izbin][icent][irp])) &&  ntry<fgkMAXTRY) {
 843         ReshuffleBuffer(fvem[ibuf][idie][izbin][icent][irp],fpoolm[idie][izbin][icent][irp]);
 844         ntry++;
 845       }
 846       for(int i1=0; i1<(int)ve2.size(); i1++){
 847         for(int i2=0; i2<(int)fvem[ibuf][idie][izbin][icent][irp].size(); i2++){
 848           FillPair(ve2[i1],fvem[ibuf][idie][izbin][icent][irp][i2],6, die, idie);
 849         }
 850       }
 851       ++nmixed;
 852     }//for(ibuf)
 853   }
 854
 855 }
 856
 857
 858 //_________________________________________________________________________________
 859 void AliAnalysisTaskMultiDielectronTG::FillPair(AliDielectronSingleTG *iep,
 860                                                 AliDielectronSingleTG *iem, int type, AliDielectron *die, Int_t idie)
 861 {
 862
 863   //
 864   // main routine for filling kinematics of pairs
 865   //
 866
 867
 868
 869   double dmass, dphiv, dpxpair, dpypair, dpzpair;
 870   double dptpair, depair, dphipair, detapair, dcos, dpsi;
 871
 872   CalcVars(iep, iem, dmass, dphiv, dpxpair, dpypair, dpzpair,
 873             dptpair, depair, dphipair, detapair, dcos, dpsi);
 874
 875
 876   double dopeningangle =  -9999;
 877   double dv0mass =  -9999;
 878   double dv0pxpair = -9999;
 879   double dv0pypair = -9999;
 880   double dv0pzpair = -9999;
 881   double dv0ptpair = -9999;
 882   double dv0epair = -9999;
 883   double dv0xvpair = -9999;
 884   double dv0yvpair = -9999;
 885   double dv0zvpair = -9999;
 886   double dv0phipair = -9999;
 887   double dv0etapair = -9999;
 888   double dv0rpair = -9999;
 889   double dpsipair =  -9999;
 890   double dphivpair =  -9999;
 891
 892   ////////////////////////////
 893   ///// calculate v0 ////////
 894   ///////////////////////////
 895   Bool_t v0OFF=kFALSE;
 896   /// for the moment, this doesn't work for MC
 897   if(fdv0mixing == kFALSE && (type==3 || type==4 || type==5 || type==6)){
 898     v0OFF = kTRUE;
 899   }
 900   if(die->GetHasMC()==kTRUE && (type==3 || type==4 || type==5 || type==6)){
 901     v0OFF = kTRUE;
 902   }
 903   if(type==0 || type==1 || type==2){
 904     v0OFF = kFALSE;
 905   }
 906
 907
 908   if(v0OFF==kFALSE){
 909     AliVTrack* trackob1= iep->GetTrack();
 910     AliVTrack* trackob2= iem->GetTrack();
 911
 912     if(!trackob1 || !trackob2){
 913       return;
 914     }
 915
 916     AliDielectronPair candidate;
 917     candidate.SetTracks(trackob1, (int)(11*iep->Charge()),
 918                         trackob2, (int)(11*iem->Charge()));
 919
 920     if(trackob1==trackob2){
 921       AliInfo("Same Objects!!");
 922       return;
 923     }
 924
 925     //const AliKFParticle &kfPairLeg1 = candidate.GetKFFirstDaughter();
 926     //const AliKFParticle &kfPairLeg2 = candidate.GetKFSecondDaughter();
 927
 928     const AliKFParticle &kfPair = candidate.GetKFParticle();
 929
 930     /*
 931     dv0mass = candidate.M();
 932     dv0pxpair = candidate.Px();
 933     dv0pypair = candidate.Py();
 934     dv0pzpair = candidate.Pz();
 935     dv0ptpair = candidate.Pt();
 936     dv0epair = candidate.E();
 937     dv0xvpair = candidate.Xv();
 938     dv0yvpair = candidate.Yv();
 939     dv0zvpair = candidate.Zv();
 940     dv0phipair = candidate.Phi();
 941     dv0etapair = candidate.Eta();
 942     */
 943     dv0mass = kfPair.GetMass();
 944     dv0pxpair = kfPair.GetPx();
 945     dv0pypair = kfPair.GetPy();
 946     dv0pzpair = kfPair.GetPz();
 947     dv0ptpair = kfPair.GetPt();
 948     dv0epair = kfPair.GetE();
 949     dv0xvpair = kfPair.GetX();
 950     dv0yvpair = kfPair.GetY();
 951     dv0zvpair = kfPair.GetZ();
 952     dv0phipair = kfPair.GetPhi();
 953     dv0etapair = kfPair.GetEta();
 954     dv0rpair =  kfPair.GetR();
 955
 956     dopeningangle = candidate.OpeningAngle();
 957     dpsipair = candidate.PsiPair(fbz);
 958     dphivpair = candidate.PhivPair(fbz);
 959
 960
 961   }
 962
 963   Double_t values[AliDielectronVarManager::kNMaxValues];
 964   TString  className1;
 965   TString  className2;
 966   className1.Form("MyPair_%s",kPairClassNamesTG[type]);
 967   className2.Form("MyPairV0_%s",kPairClassNamesTG[type]);
 968
 969   AliDielectronHistos *fHistos = die->GetHistoManager();
 970   Bool_t pairClass1=fHistos->GetHistogramList()->FindObject(className1.Data())!=0x0;
 971   Bool_t pairClass2=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
 972
 973   if (pairClass1 && PairTrackcut(dphiv, dcos, dmass, idie)==true){
 974     ///import pair variables to values!!
 975     values[AliDielectronVarManager::kPx] = dpxpair;
 976     values[AliDielectronVarManager::kPy] = dpypair;
 977     values[AliDielectronVarManager::kPz] = dpzpair;
 978     values[AliDielectronVarManager::kPt] = dptpair;
 979     values[AliDielectronVarManager::kXv] = dv0xvpair;
 980     values[AliDielectronVarManager::kYv] = dv0yvpair;
 981     values[AliDielectronVarManager::kZv] = dv0zvpair;
 982     values[AliDielectronVarManager::kR] = dv0rpair;
 983     values[AliDielectronVarManager::kE] = depair;
 984     values[AliDielectronVarManager::kEta] = detapair;
 985     values[AliDielectronVarManager::kM] = dmass;
 986     values[AliDielectronVarManager::kPsiPair] = dphiv;
 987     values[AliDielectronVarManager::kPhivPair] = dphiv;
 988     values[AliDielectronVarManager::kPhi]  = dphipair;
 989     values[AliDielectronVarManager::kOpeningAngle]  = dcos;
 990     values[AliDielectronVarManager::kCosPointingAngle]  = TMath::Abs(TMath::ATan2(TMath::Sin(iep->Phi()-iem->Phi()),TMath::Cos(iep->Phi()-iem->Phi())));
 991     fHistos->FillClass(className1, AliDielectronVarManager::kNMaxValues, values);
 992   }
 993
 994
 995   if (pairClass2 && PairTrackcut(dphiv, dopeningangle, dv0mass, idie)==true){
 996     values[AliDielectronVarManager::kPx] = dv0pxpair;
 997     values[AliDielectronVarManager::kPy] = dv0pypair;
 998     values[AliDielectronVarManager::kPz] = dv0pzpair;
 999     values[AliDielectronVarManager::kPt] = dv0ptpair;
1000     values[AliDielectronVarManager::kXv] = dv0xvpair;
1001     values[AliDielectronVarManager::kYv] = dv0yvpair;
1002     values[AliDielectronVarManager::kZv] = dv0zvpair;
1003     values[AliDielectronVarManager::kR] = dv0rpair;
1004     values[AliDielectronVarManager::kE] = dv0epair;
1005     values[AliDielectronVarManager::kEta] = dv0etapair;
1006     values[AliDielectronVarManager::kM] = dv0mass;
1007     values[AliDielectronVarManager::kPsiPair] = dpsipair;
1008     values[AliDielectronVarManager::kPhivPair] = dphivpair;
1009     values[AliDielectronVarManager::kPhi]  = dv0phipair;
1010     values[AliDielectronVarManager::kOpeningAngle]  = dopeningangle;
1011     values[AliDielectronVarManager::kCosPointingAngle]  = TMath::Abs(TMath::ATan2(TMath::Sin(iep->Phi()-iem->Phi()),TMath::Cos(iep->Phi()-iem->Phi())));
1012     fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
1013   }
1014
1015
1016
1017 }
1018
1019 //_________________________________________________________________________________
1020 bool AliAnalysisTaskMultiDielectronTG::PairTrackcut(double phiv, double op, double mass, int idie)
1021 {
1022
1023   //
1024   // pair-by-pair cuts
1025   //
1026
1027
1028   bool pairOK = true;
1029
1030   if(fRejectPairFlag[idie] == 1 || fRejectPairFlag[idie] == 2 ||
1031      fRejectPairFlag[idie] == 3 || fRejectPairFlag[idie] == 4 ){
1032     if(fRejectPairFlag[idie] == 2 || fRejectPairFlag[idie] == 4 ){
1033       if(fbz>0 && (phiv>fdconvphiv && mass < fdconvMee) ){
1034         pairOK = false;
1035       }else if(fbz<0 && phiv<acos(-1.0)-fdconvphiv  && mass < fdconvMee){
1036         pairOK = false;
1037       }
1038     }else if(fRejectPairFlag[idie] == 1 || fRejectPairFlag[idie] == 3){
1039       if(op<fdop){
1040         pairOK = false;
1041       }
1042     }
1043   }
1044
1045   return pairOK;
1046
1047 }
1048
1049
1050 //_________________________________________________________________________________
1051 void AliAnalysisTaskMultiDielectronTG::CalcVars(AliDielectronSingleTG *iep, AliDielectronSingleTG *iem,
1052                                                 double &mass, double &phiv, double &px, double &py, double&pz,
1053                                                 double &pt, double &e, double &phi,
1054                                                 double &eta, double &cos, double &psi)
1055 {
1056
1057
1058   //
1059   // standalone calculator for the pair variables
1060   //
1061
1062   px = iep->Px()+iem->Px();
1063   py = iep->Py()+iem->Py();
1064   pz = iep->Pz()+iem->Pz();
1065   pt = sqrt(px*px+py*py);
1066   double dppair = sqrt(pt*pt+pz*pz);
1067   static const double me=0.0005109989;
1068   e = sqrt(me*me+iep->Px()*iep->Px()+iep->Py()*iep->Py()+iep->Pz()*iep->Pz())
1069     + sqrt(me*me+iem->Px()*iem->Px()+iem->Py()*iem->Py()+iem->Pz()*iem->Pz());
1070
1071   mass =  e*e-px*px-py*py-pz*pz;
1072   if(mass>=0){
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           Double_t mee = GetMass(e1[i1], e2[i2]);
1184           if(fbz>0 && ( phiv>fdconvphiv && mee < fdconvMee) ){
1185             e1[i1]->SetConvFlag(0);
1186             e2[i2]->SetConvFlag(0);
1187           }else if(fbz<0 && phiv<acos(-1.0)-fdconvphiv && mee < fdconvMee){
1188             e1[i1]->SetConvFlag(0);
1189             e2[i2]->SetConvFlag(0);
1190           }
1191         }
1192       }
1193     }
1194   }
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           Double_t mee = GetMass(e1[i1], e1[i2]);
1207           if(fbz>0 && phiv>fdconvphiv && mee < fdconvMee){
1208             e1[i1]->SetConvFlag(0);
1209             e1[i2]->SetConvFlag(0);
1210           }else if(fbz<0 && phiv<acos(-1.0)-fdconvphiv && mee < fdconvMee){
1211             e1[i1]->SetConvFlag(0);
1212             e1[i2]->SetConvFlag(0);
1213           }
1214         }
1215       }
1216     }
1217   }
1218
1219   if(e2.size()>0){
1220     for(int i1=0; i1<(int)e2.size(); i1++){
1221       for(int i2=i1+1; i2<(int)e2.size(); i2++){
1222         if(fRejectPairFlag[idie]==1){
1223           Double_t cosine = GetOpeningAngle(e2[i1], e2[i2]);
1224           if(cosine<fdop){
1225             e2[i1]->SetConvFlag(0);
1226             e2[i2]->SetConvFlag(0);
1227           }
1228         }else if(fRejectPairFlag[idie]==2){
1229           Double_t phiv = GetPhiv(e2[i1], e2[i2]);
1230           Double_t mee = GetMass(e2[i1], e2[i2]);
1231           if(fbz>0 && phiv>fdconvphiv && mee < fdconvMee){
1232             e2[i1]->SetConvFlag(0);
1233             e2[i2]->SetConvFlag(0);
1234           }else if(fbz<0 && phiv<acos(-1.0)-fdconvphiv && mee < fdconvMee){
1235             e2[i1]->SetConvFlag(0);
1236             e2[i2]->SetConvFlag(0);
1237           }
1238         }
1239       }
1240     }
1241   }
1242 }
1243
1244
1245 //_________________________________________________________________________________
1246 Double_t AliAnalysisTaskMultiDielectronTG::GetOpeningAngle(AliDielectronSingleTG* e1, AliDielectronSingleTG* e2){
1247
1248   //////////////////////
1249   //////// calculate pairs and get opening angle
1250   //////////////////////
1251
1252   double dmass, dphiv, dpxpair, dpypair, dpzpair;
1253   double dptpair, depair, dphipair, detapair, dcos, dpsi;
1254
1255   CalcVars(e1, e2, dmass, dphiv, dpxpair, dpypair, dpzpair,
1256             dptpair, depair, dphipair, detapair, dcos, dpsi);
1257
1258   return dcos;
1259 }
1260
1261 //_________________________________________________________________________________
1262 Double_t AliAnalysisTaskMultiDielectronTG::GetPhiv(AliDielectronSingleTG* e1, AliDielectronSingleTG* e2){
1263
1264   //////////////////////
1265   //////// calculate pairs and get phiv
1266   //////////////////////
1267
1268   double dmass, dphiv, dpxpair, dpypair, dpzpair;
1269   double dptpair, depair, dphipair, detapair, dcos, dpsi;
1270
1271   CalcVars(e1, e2, dmass, dphiv, dpxpair, dpypair, dpzpair,
1272             dptpair, depair, dphipair, detapair, dcos, dpsi);
1273
1274   return dphiv;
1275 }
1276
1277 //_________________________________________________________________________________
1278 Double_t AliAnalysisTaskMultiDielectronTG::GetMass(AliDielectronSingleTG* e1, AliDielectronSingleTG* e2){
1279
1280   //////////////////////
1281   //////// calculate pairs and get mass
1282   //////////////////////
1283
1284   double dmass, dphiv, dpxpair, dpypair, dpzpair;
1285   double dptpair, depair, dphipair, detapair, dcos, dpsi;
1286
1287   CalcVars(e1, e2, dmass, dphiv, dpxpair, dpypair, dpzpair,
1288             dptpair, depair, dphipair, detapair, dcos, dpsi);
1289
1290   return dmass;
1291 }