PWGDQ/dielectron/AliAnalysisTaskMultiDielectronTG.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 ///////////////////////////////////////////////////////////////////////////
  17 //                                                                       //
  18 //                        Basic Analysis Task                            //
  19 //                                                                       //
  20 ///////////////////////////////////////////////////////////////////////////
  21
  22 #include <TChain.h>
  23 #include <TH1D.h>
  24
  25 #include <AliCFContainer.h>
  26 #include <AliInputEventHandler.h>
  27 #include <AliESDInputHandler.h>
  28 #include <AliAODInputHandler.h>
  29 #include <AliAnalysisManager.h>
  30 #include <AliVEvent.h>
  31 #include <AliTriggerAnalysis.h>
  32 #include <AliPIDResponse.h>
  33 #include <AliTPCPIDResponse.h>
  34
  35 #include "AliDielectron.h"
  36 #include "AliDielectronHistos.h"
  37 #include "AliDielectronCF.h"
  38 #include "AliDielectronMC.h"
  39 #include "AliDielectronMixingHandler.h"
  40 #include "AliAnalysisTaskMultiDielectronTG.h"
  41
  42 #include "AliESDtrack.h"
  43 #include "AliESDtrackCuts.h"
  44
  45 #include "AliLog.h"
  46
  47 #include <vector>
  48 #include <deque>
  49 #include <cstdlib>
  50 #include <TRandom.h>
  51
  52 using std::vector;
  53 using std::deque;
  54
  55 const char* kPairClassNamesTG[7] = {
  56   "unlike",
  57   "like_pp",
  58   "like_ee",
  59   "mixunlike_pe",
  60   "mixunlike_ep",
  61   "mixlike_pp",
  62   "mixlike_ee"
  63 };
  64
  65
  66 ClassImp(AliDielectronSingleTG)
  67 ClassImp(AliAnalysisTaskMultiDielectronTG)
  68
  69 //_________________________________________________________________________________
  70 AliAnalysisTaskMultiDielectronTG::AliAnalysisTaskMultiDielectronTG() :
  71   AliAnalysisTaskSE(),
  72   fListDielectron(),
  73   fListHistos(),
  74   fListCF(),
  75   fQAElectron(),
  76   fSelectPhysics(kFALSE),
  77   fTriggerMask(AliVEvent::kMB),
  78   fExcludeTriggerMask(0),
  79   fTriggerOnV0AND(kFALSE),
  80   fRejectPileup(kFALSE),
  81   fTriggerLogic(kAny),
  82   fTriggerAnalysis(0x0),
  83   fEventFilter(0x0),
  84   fCutsMother(0x0),
  85   fEventStat(0x0),
  86   fEvent(0x0),
  87   fdEdXvsPt(0x0),
  88   fdEdXnSigmaElecvsPt(0x0),
  89   fdEdXvsPtTOF(0x0),
  90   fdEdXnSigmaElecvsPtTOF(0x0),
  91   fTOFbetavsPt(0x0),
  92   fTOFnSigmaElecvsPt(0x0),
  93   fNCrossedRowsTPC(0x0),
  94   fChi2ClusTPC(0x0),
  95   fRatioCrossClusTPC(0x0),
  96   fVem(0x0),
  97   fVep(0x0),
  98   fVemtmp(0x0),
  99   fVeptmp(0x0),
 100   fdconvphiv(acos(-1.0)),
 101   fdconvMee(100),
 102   fdop(0),
 103   fbz(0),
 104   fdv0mixing(kTRUE),
 105   fBGRejUnlike(kFALSE),
 106   fBGRejLike(kTRUE)
 107 {
 108   //
 109   // Constructor
 110   //
 111 }
 112
 113 //_________________________________________________________________________________
 114 AliAnalysisTaskMultiDielectronTG::AliAnalysisTaskMultiDielectronTG(const char *name) :
 115   AliAnalysisTaskSE(name),
 116   fListDielectron(),
 117   fListHistos(),
 118   fListCF(),
 119   fQAElectron(),
 120   fSelectPhysics(kFALSE),
 121   fTriggerMask(AliVEvent::kMB),
 122   fExcludeTriggerMask(0),
 123   fTriggerOnV0AND(kFALSE),
 124   fRejectPileup(kFALSE),
 125   fTriggerLogic(kAny),
 126   fTriggerAnalysis(0x0),
 127   fEventFilter(0x0),
 128   fCutsMother(0x0),
 129   fEventStat(0x0),
 130   fEvent(0x0),
 131   fdEdXvsPt(0x0),
 132   fdEdXnSigmaElecvsPt(0x0),
 133   fdEdXvsPtTOF(0x0),
 134   fdEdXnSigmaElecvsPtTOF(0x0),
 135   fTOFbetavsPt(0x0),
 136   fTOFnSigmaElecvsPt(0x0),
 137   fNCrossedRowsTPC(0x0),
 138   fChi2ClusTPC(0x0),
 139   fRatioCrossClusTPC(0x0),
 140   fVem(0x0),
 141   fVep(0x0),
 142   fVemtmp(0x0),
 143   fVeptmp(0x0),
 144   fdconvphiv(acos(-1.0)),
 145   fdconvMee(100),
 146   fdop(0),
 147   fbz(0),
 148   fdv0mixing(kTRUE),
 149   fBGRejUnlike(kFALSE),
 150   fBGRejLike(kTRUE)
 151 {
 152   //
 153   // Constructor
 154   //
 155   DefineInput(0,TChain::Class());
 156   DefineOutput(1, TList::Class());
 157   DefineOutput(2, TList::Class());
 158   DefineOutput(3, TH1D::Class());
 159   fListHistos.SetName("Dielectron_Histos_Multi");
 160   fListCF.SetName("Dielectron_CF_Multi");
 161   fListDielectron.SetOwner();
 162   fListHistos.SetOwner();
 163   fListCF.SetOwner();
 164
 165   ///////////////
 166   for(int i=0;i<fgkNDIE; i++){
 167     for(int j=0;j<fgkNZBIN;j++){
 168       for(int k=0;k<fgkNCENT;k++){
 169         for(int l=0; l<fgkNRPBIN; l++){
 170           fibuf[i][j][k][l] = 0;
 171           fpoolp[i][j][k][l].clear();
 172           fpoolm[i][j][k][l].clear();
 173           for(int ib=0;ib<fgkNBUF; ib++){
 174             fvep[ib][i][j][k][l].clear();
 175             fvem[ib][i][j][k][l].clear();
 176           }
 177         }
 178       }
 179     }
 180   }
 181
 182
 183
 184 }
 185
 186 //_________________________________________________________________________________
 187 AliAnalysisTaskMultiDielectronTG::~AliAnalysisTaskMultiDielectronTG()
 188 {
 189   //
 190   // Destructor
 191   //
 192
 193   //histograms and CF are owned by the dielectron framework.
 194   //however they are streamed to file, so in the first place the
 195   //lists need to be owner...
 196   fListHistos.SetOwner(kFALSE);
 197   fListCF.SetOwner(kFALSE);
 198
 199   for(int i=0;i<fgkNDIE; i++){
 200     for(int j=0;j<fgkNZBIN;j++){
 201       for(int k=0;k<fgkNCENT;k++){
 202         for(int l=0; l<fgkNRPBIN; l++){
 203           fibuf[i][j][k][l] = 0;
 204           fpoolp[i][j][k][l].clear();
 205           fpoolm[i][j][k][l].clear();
 206           for(int ib=0;ib<fgkNBUF; ib++){
 207             fvep[ib][i][j][k][l].clear();
 208             fvem[ib][i][j][k][l].clear();
 209           }
 210         }
 211       }
 212     }
 213   }
 214 }
 215 //_________________________________________________________________________________
 216 void AliAnalysisTaskMultiDielectronTG::UserCreateOutputObjects()
 217 {
 218   //
 219   // Add all histogram manager histogram lists to the output TList
 220   //
 221
 222   if (!fListHistos.IsEmpty()||!fListCF.IsEmpty()) return; //already initialised
 223
 224 //   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 225 //   Bool_t isESD=man->GetInputEventHandler()->IsA()==AliESDInputHandler::Class();
 226 //   Bool_t isAOD=man->GetInputEventHandler()->IsA()==AliAODInputHandler::Class();
 227
 228   TIter nextDie(&fListDielectron);
 229   AliDielectron *die=0;
 230   while ( (die=static_cast<AliDielectron*>(nextDie())) ){
 231     die->Init();
 232     if (die->GetHistogramList()) fListHistos.Add(const_cast<THashList*>(die->GetHistogramList()));
 233     if (die->GetCFManagerPair()) fListCF.Add(const_cast<AliCFContainer*>(die->GetCFManagerPair()->GetContainer()));
 234   }
 235
 236   Int_t cuts=fListDielectron.GetEntries();
 237   Int_t nbins=kNbinsEvent+2*cuts;
 238   if (!fEventStat){
 239     fEventStat=new TH1D("hEventStat","Event statistics",nbins,0,nbins);
 240     fEventStat->GetXaxis()->SetBinLabel(1,"Before Phys. Sel.");
 241     fEventStat->GetXaxis()->SetBinLabel(2,"After Phys. Sel.");
 242
 243     //default names
 244     fEventStat->GetXaxis()->SetBinLabel(3,"Bin3 not used");
 245     fEventStat->GetXaxis()->SetBinLabel(4,"Bin4 not used");
 246     fEventStat->GetXaxis()->SetBinLabel(5,"Bin5 not used");
 247
 248     if(fTriggerOnV0AND) fEventStat->GetXaxis()->SetBinLabel(3,"V0and triggers");
 249     if (fEventFilter) fEventStat->GetXaxis()->SetBinLabel(4,"After Event Filter");
 250     if (fRejectPileup) fEventStat->GetXaxis()->SetBinLabel(5,"After Pileup rejection");
 251
 252     for (Int_t i=0; i<cuts; ++i){
 253       fEventStat->GetXaxis()->SetBinLabel((kNbinsEvent+1)+2*i,Form("#splitline{1 candidate}{%s}",fListDielectron.At(i)->GetName()));
 254       fEventStat->GetXaxis()->SetBinLabel((kNbinsEvent+2)+2*i,Form("#splitline{With >1 candidate}{%s}",fListDielectron.At(i)->GetName()));
 255     }
 256   }
 257
 258   if (!fTriggerAnalysis) fTriggerAnalysis=new AliTriggerAnalysis;
 259   fTriggerAnalysis->EnableHistograms();
 260   fTriggerAnalysis->SetAnalyzeMC(AliDielectronMC::Instance()->HasMC());
 261
 262
 263
 264   int nbinx=400;
 265   float maxx=20;
 266   float minx=0.2;
 267   float binw = (TMath::Log(maxx)-TMath::Log(minx))/nbinx;
 268   double xbin[401];
 269   for(int ii=0;ii<nbinx+1;ii++){
 270     xbin[ii] = TMath::Exp(TMath::Log(minx) + 0.5*binw+binw*ii);
 271   }
 272
 273
 274   fQAElectron = new TList();
 275   fQAElectron->SetName("QAElectron");
 276   fQAElectron->SetOwner();
 277
 278
 279   fEvent = new TH1D("Event","centrality",   100,0,100);
 280   fQAElectron->Add(fEvent);
 281   fdEdXvsPt = new TH2D("dEdXvsPt","dE/dX vs. PT of TPC", nbinx, xbin, 2000,0,200);
 282   fQAElectron->Add(fdEdXvsPt);
 283   fdEdXnSigmaElecvsPt = new TH2D("fdEdXnSigmaElecvsPt"," dE/dX normalized to electron vs. pT of TPC",
 284                                  nbinx, xbin, 2000, -10, 10);
 285   fQAElectron->Add(fdEdXnSigmaElecvsPt);
 286
 287   fdEdXvsPtTOF = new TH2D("dEdXvsPtTOF","dE/dX vs. PT of TPC", nbinx, xbin, 2000,0,200);
 288   fQAElectron->Add(fdEdXvsPtTOF);
 289   fdEdXnSigmaElecvsPtTOF = new TH2D("fdEdXnSigmaElecvsPtTOF"," dE/dX normalized to electron vs. pT of TPC",
 290                                  nbinx, xbin, 2000, -10, 10);
 291   fQAElectron->Add(fdEdXnSigmaElecvsPtTOF);
 292
 293
 294
 295   fTOFbetavsPt = new TH2D("fTOFbetavsPt","TOF beta vs. p", 400, 0, 20, 1200, 0, 1.2);
 296   fQAElectron->Add(fTOFbetavsPt);
 297   fTOFnSigmaElecvsPt = new TH2D("fTOFnSigmaElecvsPt","TOF nsigma for electron", 400, 0, 20, 2000, -10, 10);
 298   fQAElectron->Add(fTOFnSigmaElecvsPt);
 299
 300   fNCrossedRowsTPC = new TH2F("fNCrossedRowsTPC", "TPC nCrossed Rows vs. pT", 200, 0, 20, 200, 0, 200);
 301   fQAElectron->Add(fNCrossedRowsTPC);
 302   fChi2ClusTPC = new TH2F("fChi2ClusTPC", "hChi2ClusTPC vs. pT", 200, 0, 20, 200, 0, 10);
 303   fQAElectron->Add(fChi2ClusTPC);
 304
 305   fRatioCrossClusTPC = new TH2F("fRatioCrossClusTPC", "hRatioCrossClusTPC vs. pT", 200, 0, 20, 200, 0, 10);
 306   fQAElectron->Add(fRatioCrossClusTPC);
 307
 308   fListHistos.Add(fQAElectron);
 309
 310   fListHistos.SetOwner();
 311
 312   PostData(1, &fListHistos);
 313   PostData(2, &fListCF);
 314   PostData(3, fEventStat);
 315
 316   fCutsMother = new AliESDtrackCuts;
 317   fCutsMother->SetDCAToVertex2D(kTRUE);
 318   fCutsMother->SetMaxDCAToVertexZ(3.0);
 319   fCutsMother->SetMaxDCAToVertexXY(1.0);
 320   fCutsMother->SetPtRange(  0.05 , 200.0);
 321   fCutsMother->SetEtaRange( -0.84 , 0.84 );
 322   fCutsMother->SetAcceptKinkDaughters(kFALSE);
 323   fCutsMother->SetRequireITSRefit(kTRUE);
 324   fCutsMother->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kFirst);
 325   fCutsMother->SetMinNClustersITS(3);
 326   fCutsMother->SetRequireTPCRefit(kTRUE);
 327
 328
 329   AliInfo(Form("PairCutType = %d %d %d %d %d",
 330                fRejectPairFlag[0],
 331                fRejectPairFlag[1],
 332                fRejectPairFlag[2],
 333                fRejectPairFlag[3],
 334                fRejectPairFlag[4]));
 335
 336 }
 337
 338 //_________________________________________________________________________________
 339 void AliAnalysisTaskMultiDielectronTG::UserExec(Option_t *)
 340 {
 341   //
 342   // Main loop. Called for every event
 343   //
 344
 345   if (fListHistos.IsEmpty()&&fListCF.IsEmpty()) return;
 346
 347   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 348   Bool_t isESD=man->GetInputEventHandler()->IsA()==AliESDInputHandler::Class();
 349   Bool_t isAOD=man->GetInputEventHandler()->IsA()==AliAODInputHandler::Class();
 350
 351   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 352   if (!inputHandler) return;
 353
 354 //   AliPIDResponse *pidRes=inputHandler->GetPIDResponse();
 355   if ( inputHandler->GetPIDResponse() ){
 356     // for the 2.76 pass2 MC private train. Together with a sigma shift of -0.169
 357 //    pidRes->GetTPCResponse().SetSigma(4.637e-3,2.41332105409873257e+04);
 358     AliDielectronVarManager::SetPIDResponse( inputHandler->GetPIDResponse() );
 359   } else {
 360     AliFatal("This task needs the PID response attached to the input event handler!");
 361   }
 362
 363   // Was event selected ?
 364   ULong64_t isSelected = AliVEvent::kAny;
 365   Bool_t isRejected = kFALSE;
 366   if( fSelectPhysics && inputHandler){
 367     if((isESD && inputHandler->GetEventSelection()) || isAOD){
 368       isSelected = inputHandler->IsEventSelected();
 369       if (fExcludeTriggerMask && (isSelected&fExcludeTriggerMask)) isRejected=kTRUE;
 370       if (fTriggerLogic==kAny) isSelected&=fTriggerMask;
 371       else if (fTriggerLogic==kExact) isSelected=((isSelected&fTriggerMask)==fTriggerMask);
 372     }
 373    }
 374
 375
 376   //Before physics selection
 377   fEventStat->Fill(kAllEvents);
 378   if (isSelected==0||isRejected) {
 379     PostData(3,fEventStat);
 380     return;
 381   }
 382   //after physics selection
 383   fEventStat->Fill(kSelectedEvents);
 384
 385   //V0and
 386   if(fTriggerOnV0AND){
 387     if(isESD){if (!fTriggerAnalysis->IsOfflineTriggerFired(static_cast<AliESDEvent*>(InputEvent()), AliTriggerAnalysis::kV0AND))
 388             return;}
 389     if(isAOD){if(!((static_cast<AliAODEvent*>(InputEvent()))->GetVZEROData()->GetV0ADecision() == AliVVZERO::kV0BB &&
 390                    (static_cast<AliAODEvent*>(InputEvent()))->GetVZEROData()->GetV0CDecision() == AliVVZERO::kV0BB) )
 391         return;}
 392   }
 393
 394
 395   fEventStat->Fill(kV0andEvents);
 396
 397   //Fill Event histograms before the event filter
 398   TIter nextDie(&fListDielectron);
 399   AliDielectron *die=0;
 400   Double_t values[AliDielectronVarManager::kNMaxValues]={0};
 401   Double_t valuesMC[AliDielectronVarManager::kNMaxValues]={0};
 402   AliDielectronVarManager::SetEvent(InputEvent());
 403   AliDielectronVarManager::Fill(InputEvent(),values);
 404   AliDielectronVarManager::Fill(InputEvent(),valuesMC);
 405   Bool_t hasMC=AliDielectronMC::Instance()->HasMC();
 406   if (hasMC) {
 407     if (AliDielectronMC::Instance()->ConnectMCEvent())
 408       AliDielectronVarManager::Fill(AliDielectronMC::Instance()->GetMCEvent(),valuesMC);
 409   }
 410
 411   while ( (die=static_cast<AliDielectron*>(nextDie())) ){
 412     AliDielectronHistos *h=die->GetHistoManager();
 413     if (h){
 414       if (h->GetHistogramList()->FindObject("Event_noCuts"))
 415         h->FillClass("Event_noCuts",AliDielectronVarManager::kNMaxValues,values);
 416       if (hasMC && h->GetHistogramList()->FindObject("MCEvent_noCuts"))
 417         h->FillClass("Event_noCuts",AliDielectronVarManager::kNMaxValues,valuesMC);
 418     }
 419   }
 420   nextDie.Reset();
 421
 422   //event filter
 423   if (fEventFilter) {
 424     if (!fEventFilter->IsSelected(InputEvent())) return;
 425   }
 426   fEventStat->Fill(kFilteredEvents);
 427
 428   //pileup
 429   if (fRejectPileup){
 430     if (InputEvent()->IsPileupFromSPD(3,0.8,3.,2.,5.)) return;
 431   }
 432   fEventStat->Fill(kPileupEvents);
 433
 434   //bz for AliKF
 435   fbz = InputEvent()->GetMagneticField();
 436   AliKFParticle::SetField( fbz );
 437   AliDielectronPID::SetCorrVal((Double_t)InputEvent()->GetRunNumber());
 438
 439   //Process event in all AliDielectron instances
 440   //   TIter nextDie(&fListDielectron);
 441   //   AliDielectron *die=0;
 442   Int_t idie=0;
 443   while ( (die=static_cast<AliDielectron*>(nextDie())) ){
 444     //AliInfo(Form(" **************** die->Process(InputEvent()) : %d **************************", idie));
 445     die->SetDontClearArrays(kTRUE);
 446     die->Process(InputEvent());
 447     if (die->HasCandidates()){
 448       Int_t ncandidates=die->GetPairArray(1)->GetEntriesFast();
 449       if (ncandidates==1) fEventStat->Fill((kNbinsEvent)+2*idie);
 450       else if (ncandidates>1) fEventStat->Fill((kNbinsEvent+1)+2*idie);
 451     }
 452
 453
 454     AliDielectronVarManager::Fill(InputEvent(), fgValues);
 455     for (Int_t ii=0; ii<2; ++ii){
 456       TObjArray *obj = (TObjArray*)die->GetTrackArray(ii);
 457       Int_t ntracks=obj->GetEntriesFast();
 458       //AliInfo(Form(" ************** # of tracks = %d", ntracks));
 459       for (Int_t itrack=0; itrack<ntracks; ++itrack){
 460
 461         ////////////////////////////////////////////////////////////////////
 462         AliDielectronVarManager::Fill(obj->UncheckedAt(itrack), fgValues);
 463         ////////////////////////////////////////////////////////////////////
 464
 465         if(fgValues[AliDielectronVarManager::kCharge]>0){
 466           fVeptmp.push_back(new  AliDielectronSingleTG(1,
 467                                                        fgValues[AliDielectronVarManager::kCentrality],
 468                                                        fgValues[AliDielectronVarManager::kXv],
 469                                                        fgValues[AliDielectronVarManager::kYv],
 470                                                        fgValues[AliDielectronVarManager::kZv],
 471                                                        fgValues[AliDielectronVarManager::kPx],
 472                                                        fgValues[AliDielectronVarManager::kPy],
 473                                                        fgValues[AliDielectronVarManager::kPz],
 474                                                        fgValues[AliDielectronVarManager::kPt],
 475                                                        fgValues[AliDielectronVarManager::kEta],
 476                                                        fgValues[AliDielectronVarManager::kPhi],
 477                                                        fgValues[AliDielectronVarManager::kTheta],
 478                                                        1, 1, static_cast<AliVTrack*>(obj->UncheckedAt(itrack)))
 479                             );
 480         }else if(fgValues[AliDielectronVarManager::kCharge]<0){
 481           fVemtmp.push_back(new  AliDielectronSingleTG(-1,
 482                                                        fgValues[AliDielectronVarManager::kCentrality],
 483                                                        fgValues[AliDielectronVarManager::kXv],
 484                                                        fgValues[AliDielectronVarManager::kYv],
 485                                                        fgValues[AliDielectronVarManager::kZv],
 486                                                        fgValues[AliDielectronVarManager::kPx],
 487                                                        fgValues[AliDielectronVarManager::kPy],
 488                                                        fgValues[AliDielectronVarManager::kPz],
 489                                                        fgValues[AliDielectronVarManager::kPt],
 490                                                        fgValues[AliDielectronVarManager::kEta],
 491                                                        fgValues[AliDielectronVarManager::kPhi],
 492                                                        fgValues[AliDielectronVarManager::kTheta],
 493                                                        1, 1, static_cast<AliVTrack*>(obj->UncheckedAt(itrack)))
 494                             );
 495         }
 496       }
 497     }
 498     //AliInfo(Form("size of e and p = %d %d", (int)fVeptmp.size(), (int)fVemtmp.size()));
 499
 500
 501     CheckGhostPairs(fVeptmp);
 502     CheckGhostPairs(fVemtmp);
 503     if(fRejectPairFlag[idie]==1 || fRejectPairFlag[idie]==2){
 504       RejectPairs(fVeptmp, fVemtmp, idie);
 505     }
 506     RandomizePool(fVeptmp, fVemtmp);
 507     CalcPair(fVep, fVem, die, idie);
 508
 509     //    AliInfo(Form("size of e and p (after) = %d %d", (int)fVep.size(), (int)fVem.size()));
 510
 511     double dwcent = 100.0/fgkNCENT;
 512     double dwiz = 20.0/fgkNZBIN;
 513     double dwrp = acos(-1.0)/fgkNRPBIN;
 514
 515     int icent = (int)(fgValues[AliDielectronVarManager::kCentrality]/dwcent);
 516     int izbin = (int)((fgValues[AliDielectronVarManager::kZvPrim]+10)/dwiz);
 517     int irp = (int)((fgValues[AliDielectronVarManager::kV0ACrpH2])/dwrp);
 518
 519     if(icent<0) icent=0;
 520     if(icent>=fgkNCENT) icent=fgkNCENT-1;
 521     if(izbin<0) izbin=0;
 522     if(izbin>=fgkNZBIN) izbin=fgkNZBIN-1;
 523     if(irp<0) irp=0;
 524     if(irp>=fgkNRPBIN) irp=fgkNRPBIN-1;
 525
 526     fvep[fibuf[idie][izbin][icent][irp]][idie][izbin][icent][irp].clear();
 527     for(int iep = 0; iep<(int)fVep.size();iep++) {
 528       fvep[fibuf[idie][izbin][icent][irp]][idie][izbin][icent][irp].push_back(fVep[iep]);
 529       fpoolp[idie][izbin][icent][irp].push_back(fVep[iep]);
 530       if(fpoolp[idie][izbin][icent][irp].size()>fgkMAXPOOL) {
 531         fpoolp[idie][izbin][icent][irp].pop_front();
 532       }
 533     }
 534     fvem[fibuf[idie][izbin][icent][irp]][idie][izbin][icent][irp].clear();
 535     for(int iem = 0; iem<(int)fVem.size();iem++) {
 536       fvem[fibuf[idie][izbin][icent][irp]][idie][izbin][icent][irp].push_back(fVem[iem]);
 537       fpoolm[idie][izbin][icent][irp].push_back(fVem[iem]);
 538       if(fpoolm[idie][izbin][icent][irp].size()>fgkMAXPOOL) {
 539         fpoolm[idie][izbin][icent][irp].pop_front();
 540       }
 541     }
 542
 543
 544     fibuf[idie][izbin][icent][irp]++;
 545     if(fibuf[idie][izbin][icent][irp]>= fgkNBUF) fibuf[idie][izbin][icent][irp]=0;
 546
 547
 548     fVeptmp.clear();
 549     fVemtmp.clear();
 550     fVep.clear();
 551     fVem.clear();
 552
 553
 554     ++idie;
 555   }
 556
 557
 558   AliESDEvent *fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 559   fEvent->Fill(values[AliDielectronVarManager::kCentrality]);
 560   for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
 561     AliESDtrack* track = fESD->GetTrack(iTracks);
 562     if (!track) {
 563       Printf("ERROR: Could not receive track %d", iTracks);
 564       continue;
 565     }
 566     if(!fCutsMother->AcceptTrack(track)) continue;
 567     fdEdXvsPt->Fill(track->GetTPCmomentum(), track->GetTPCsignal());
 568     fdEdXnSigmaElecvsPt->Fill(track->GetTPCmomentum(),
 569                               AliDielectronVarManager::GetPIDResponse()->NumberOfSigmasTPC(track,
 570                                                                                       AliPID::kElectron)
 571                               -AliDielectronPID::GetCorrVal());
 572     /// for beta caliculaton
 573     Double_t l = track->GetIntegratedLength();  // cm
 574     Double_t t = track->GetTOFsignal();
 575     Double_t t0 = AliDielectronVarManager::GetPIDResponse()->GetTOFResponse().GetTimeZero(); // ps
 576     Double_t beta = 0;
 577     if( (l < 360. || l > 800.) || (t <= 0.) || (t0 >999990.0) ) {
 578       beta=-9999;
 579     }
 580     else {
 581       t -= t0; // subtract the T0
 582       l *= 0.01;  // cm ->m
 583       t *= 1e-12; //ps -> s
 584
 585       Double_t v = l / t;
 586       beta = v / TMath::C();
 587     }
 588
 589     fTOFbetavsPt->Fill(track->GetTPCmomentum(), beta);
 590     fTOFnSigmaElecvsPt->Fill(track->GetTPCmomentum(),
 591                              AliDielectronVarManager::GetPIDResponse()->NumberOfSigmasTOF(track,
 592                                                                                      AliPID::kElectron));
 593     ////rough PID is required
 594     if( fabs(AliDielectronVarManager::GetPIDResponse()->NumberOfSigmasTOF(track, AliPID::kElectron))<3){
 595       fdEdXvsPtTOF->Fill(track->GetTPCmomentum(), track->GetTPCsignal());
 596       fdEdXnSigmaElecvsPtTOF->Fill(track->GetTPCmomentum(),
 597                                    AliDielectronVarManager::GetPIDResponse()->NumberOfSigmasTPC(track,
 598                                                                                                 AliPID::kElectron)
 599                                    -AliDielectronPID::GetCorrVal());
 600
 601
 602       if(track->GetTPCsignal()>70 && track->GetTPCsignal()<90){
 603         fNCrossedRowsTPC->Fill(track->GetTPCmomentum(), track->GetTPCCrossedRows());
 604         //fChi2ClusTPC->Fill(track->GetTPCmomentum(), track->GetTPCchi2()/track->GetTPCNcls());
 605         fChi2ClusTPC->Fill(track->GetTPCmomentum(), track->GetTPCchi2());
 606         fRatioCrossClusTPC->Fill(track->GetTPCmomentum(), track->GetTPCCrossedRows()/track->GetTPCNclsF());
 607       }
 608     }
 609   }
 610
 611   PostData(1, &fListHistos);
 612   PostData(2, &fListCF);
 613   PostData(3,fEventStat);
 614 }
 615
 616 //_________________________________________________________________________________
 617 void AliAnalysisTaskMultiDielectronTG::FinishTaskOutput()
 618 {
 619   //
 620   // Write debug tree
 621   //
 622   TIter nextDie(&fListDielectron);
 623   AliDielectron *die=0;
 624   while ( (die=static_cast<AliDielectron*>(nextDie())) ){
 625     die->SaveDebugTree();
 626     AliDielectronMixingHandler *mix=die->GetMixingHandler();
 627 //    printf("\n\n\n===============\ncall mix in Terminate: %p (%p)\n=================\n\n",mix,die);
 628     if (mix) mix->MixRemaining(die);
 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 = mass;
1074   }else{
1075     mass = sqrt(mass);
1076   }
1077
1078
1079   phi = atan2(py, px);
1080   eta = -0.5*TMath::Log((dppair+pz)/(dppair-pz));
1081   double p1 = sqrt(pow(iep->Px(),2)+pow(iep->Py(),2)+pow(iep->Pz(),2));
1082   double p2 = sqrt(pow(iem->Px(),2)+pow(iem->Py(),2)+pow(iem->Pz(),2));
1083   cos = acos((iep->Px()*iem->Px()+iep->Py()*iem->Py()+iep->Pz()*iem->Pz())/(p1*p2));
1084
1085   double dtheta = iep->Theta()-iem->Theta();
1086   psi = asin(dtheta/cos);
1087
1088
1089   //unit vector of (pep+pem)
1090   double pl = dppair;
1091   double ux = px/pl;
1092   double uy = py/pl;
1093   double uz = pz/pl;
1094   double ax = uy/sqrt(ux*ux+uy*uy);
1095   double ay = -ux/sqrt(ux*ux+uy*uy);
1096
1097   //momentum of e+ and e- in (ax,ay,az) axis. Note that az=0 by
1098   //definition.
1099   //double ptep = iep->Px()*ax + iep->Py()*ay;
1100   //double ptem = iem->Px()*ax + iem->Py()*ay;
1101
1102   double pxep = iep->Px();
1103   double pyep = iep->Py();
1104   double pzep = iep->Pz();
1105   double pxem = iem->Px();
1106   double pyem = iem->Py();
1107   double pzem = iem->Pz();
1108
1109
1110   //vector product of pep X pem
1111   double vpx = pyep*pzem - pzep*pyem;
1112   double vpy = pzep*pxem - pxep*pzem;
1113   double vpz = pxep*pyem - pyep*pxem;
1114   double vp = sqrt(vpx*vpx+vpy*vpy+vpz*vpz);
1115   //double thev = acos(vpz/vp);
1116
1117   //unit vector of pep X pem
1118   double vx = vpx/vp;
1119   double vy = vpy/vp;
1120   double vz = vpz/vp;
1121
1122   //The third axis defined by vector product (ux,uy,uz)X(vx,vy,vz)
1123   double wx = uy*vz - uz*vy;
1124   double wy = uz*vx - ux*vz;
1125   double wz = ux*vy - uy*vx;
1126   double wl = sqrt(wx*wx+wy*wy+wz*wz);
1127   // by construction, (wx,wy,wz) must be a unit vector.
1128   if(fabs(wl - 1.0) > 0.00001) std::cout << "Calculation error in W vector"<<std::endl;
1129   // measure angle between (wx,wy,wz) and (ax,ay,0). The angle between them
1130   // should be small if the pair is conversion
1131   //
1132   double cosPhiV = wx*ax + wy*ay;
1133   phiv = acos(cosPhiV);
1134
1135 }
1136
1137 //_________________________________________________________________________________
1138 void AliAnalysisTaskMultiDielectronTG::ReshuffleBuffer(vector<AliDielectronSingleTG*> ve, deque<AliDielectronSingleTG*> pool)
1139 {
1140   //If there is not enough electron in the pool, give up
1141   //
1142   // ReshuffleBuffer for th event mixing
1143   //
1144
1145   unsigned int ne = ve.size();
1146   unsigned int poolsize = pool.size();
1147   int used[fgkMAXPOOL];
1148   for(int i=0;i<(int)fgkMAXPOOL;i++){
1149     used[i]=0;
1150   }
1151
1152   if(poolsize < ne) {
1153     std::cout <<" pool size="<<poolsize<<" ne"<<ne<<std::endl;
1154     return;
1155   }
1156   for(unsigned int ie=0; ie < ne; ie++) {
1157     int j = rand()%poolsize;
1158     while(used[j]==1){
1159       j = rand()%poolsize;
1160     }
1161     ve[ie] = pool[j];
1162     used[j]=1;
1163   }
1164
1165 }
1166
1167 //_________________________________________________________________________________
1168 void AliAnalysisTaskMultiDielectronTG::RejectPairs(vector<AliDielectronSingleTG*> e1, vector<AliDielectronSingleTG*> e2, Int_t idie){
1169
1170   ////////////////////////////////////
1171   ///// to reject pairs from track arrays
1172   ///// conversions, ghost ..
1173   ///////////////////////////////////
1174   if(e1.size()>0 && e2.size()>0){
1175     for(int i1=0; i1<(int)e1.size(); i1++){
1176       for(int i2=0; i2<(int)e2.size(); i2++){
1177         if(fRejectPairFlag[idie]==1){
1178           Double_t cosine = GetOpeningAngle(e1[i1], e2[i2]);
1179           if(cosine<fdop){
1180             e1[i1]->SetConvFlag(0);
1181             e2[i2]->SetConvFlag(0);
1182           }
1183         }else if(fRejectPairFlag[idie]==2){
1184           Double_t phiv = GetPhiv(e1[i1], e2[i2]);
1185           Double_t mee = GetMass(e1[i1], e2[i2]);
1186           if(fbz>0 && ( phiv>fdconvphiv && mee < fdconvMee) ){
1187             e1[i1]->SetConvFlag(0);
1188             e2[i2]->SetConvFlag(0);
1189           }else if(fbz<0 && phiv<acos(-1.0)-fdconvphiv && mee < fdconvMee){
1190             e1[i1]->SetConvFlag(0);
1191             e2[i2]->SetConvFlag(0);
1192           }
1193         }
1194       }
1195     }
1196   }
1197   if(e1.size()>0){
1198     for(int i1=0; i1<(int)e1.size(); i1++){
1199       for(int i2=i1+1; i2<(int)e1.size(); i2++){
1200         if(fRejectPairFlag[idie]==1){
1201           Double_t cosine = GetOpeningAngle(e1[i1], e1[i2]);
1202           if(cosine<fdop){
1203             e1[i1]->SetConvFlag(0);
1204             e1[i2]->SetConvFlag(0);
1205           }
1206         }else if(fRejectPairFlag[idie]==2){
1207           Double_t phiv = GetPhiv(e1[i1], e1[i2]);
1208           Double_t mee = GetMass(e1[i1], e1[i2]);
1209           if(fbz>0 && phiv>fdconvphiv && mee < fdconvMee){
1210             e1[i1]->SetConvFlag(0);
1211             e1[i2]->SetConvFlag(0);
1212           }else if(fbz<0 && phiv<acos(-1.0)-fdconvphiv && mee < fdconvMee){
1213             e1[i1]->SetConvFlag(0);
1214             e1[i2]->SetConvFlag(0);
1215           }
1216         }
1217       }
1218     }
1219   }
1220
1221   if(e2.size()>0){
1222     for(int i1=0; i1<(int)e2.size(); i1++){
1223       for(int i2=i1+1; i2<(int)e2.size(); i2++){
1224         if(fRejectPairFlag[idie]==1){
1225           Double_t cosine = GetOpeningAngle(e2[i1], e2[i2]);
1226           if(cosine<fdop){
1227             e2[i1]->SetConvFlag(0);
1228             e2[i2]->SetConvFlag(0);
1229           }
1230         }else if(fRejectPairFlag[idie]==2){
1231           Double_t phiv = GetPhiv(e2[i1], e2[i2]);
1232           Double_t mee = GetMass(e2[i1], e2[i2]);
1233           if(fbz>0 && phiv>fdconvphiv && mee < fdconvMee){
1234             e2[i1]->SetConvFlag(0);
1235             e2[i2]->SetConvFlag(0);
1236           }else if(fbz<0 && phiv<acos(-1.0)-fdconvphiv && mee < fdconvMee){
1237             e2[i1]->SetConvFlag(0);
1238             e2[i2]->SetConvFlag(0);
1239           }
1240         }
1241       }
1242     }
1243   }
1244 }
1245
1246
1247 //_________________________________________________________________________________
1248 Double_t AliAnalysisTaskMultiDielectronTG::GetOpeningAngle(AliDielectronSingleTG* e1, AliDielectronSingleTG* e2){
1249
1250   //////////////////////
1251   //////// calculate pairs and get opening angle
1252   //////////////////////
1253
1254   double dmass, dphiv, dpxpair, dpypair, dpzpair;
1255   double dptpair, depair, dphipair, detapair, dcos, dpsi;
1256
1257   CalcVars(e1, e2, dmass, dphiv, dpxpair, dpypair, dpzpair,
1258             dptpair, depair, dphipair, detapair, dcos, dpsi);
1259
1260   return dcos;
1261 }
1262
1263 //_________________________________________________________________________________
1264 Double_t AliAnalysisTaskMultiDielectronTG::GetPhiv(AliDielectronSingleTG* e1, AliDielectronSingleTG* e2){
1265
1266   //////////////////////
1267   //////// calculate pairs and get phiv
1268   //////////////////////
1269
1270   double dmass, dphiv, dpxpair, dpypair, dpzpair;
1271   double dptpair, depair, dphipair, detapair, dcos, dpsi;
1272
1273   CalcVars(e1, e2, dmass, dphiv, dpxpair, dpypair, dpzpair,
1274             dptpair, depair, dphipair, detapair, dcos, dpsi);
1275
1276   return dphiv;
1277 }
1278
1279 //_________________________________________________________________________________
1280 Double_t AliAnalysisTaskMultiDielectronTG::GetMass(AliDielectronSingleTG* e1, AliDielectronSingleTG* e2){
1281
1282   //////////////////////
1283   //////// calculate pairs and get mass
1284   //////////////////////
1285
1286   double dmass, dphiv, dpxpair, dpypair, dpzpair;
1287   double dptpair, depair, dphipair, detapair, dcos, dpsi;
1288
1289   CalcVars(e1, e2, dmass, dphiv, dpxpair, dpypair, dpzpair,
1290             dptpair, depair, dphipair, detapair, dcos, dpsi);
1291
1292   return dmass;
1293 }