PWGDQ/dielectron/AliDielectron.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 //                Dielectron Analysis Main class                         //
  18 //                                                                       //
  19 /*
  20 Framework to perform event selectoin, single track selection and track pair
  21 selection.
  22
  23 Convention for the signs of the pair in fPairCandidates:
  24 The names are available via the function PairClassName(Int_t i)
  25
  26 0: ev1+ ev1+  (same event like sign +)
  27 1: ev1+ ev1-  (same event unlike sign)
  28 2: ev1- ev1-  (same event like sign -)
  29
  30 3: ev1+ ev2+  (mixed event like sign +)
  31 4: ev1- ev2+  (mixed event unlike sign -+)
  32 6: ev1+ ev2-  (mixed event unlike sign +-)
  33 7: ev1- ev2-  (mixed event like sign -)
  34
  35 5: ev2+ ev2+  (same event like sign +)
  36 8: ev2+ ev2-  (same event unlike sign)
  37 9: ev2- ev2-  (same event like sign -)
  38
  39 10: ev1+ ev1- (same event track rotation)
  40
  41 */
  42 //                                                                       //
  43 ///////////////////////////////////////////////////////////////////////////
  44
  45 #include <TString.h>
  46 #include <TList.h>
  47 #include <TMath.h>
  48 #include <TObject.h>
  49 #include <TGrid.h>
  50
  51 #include <AliKFParticle.h>
  52
  53 #include <AliESDInputHandler.h>
  54 #include <AliAnalysisManager.h>
  55 #include <AliEPSelectionTask.h>
  56 #include <AliEventplane.h>
  57 #include <AliVEvent.h>
  58 #include <AliVParticle.h>
  59 #include <AliVTrack.h>
  60 #include "AliDielectronPair.h"
  61 #include "AliDielectronHistos.h"
  62 #include "AliDielectronCF.h"
  63 #include "AliDielectronMC.h"
  64 #include "AliDielectronVarManager.h"
  65 #include "AliDielectronTrackRotator.h"
  66 #include "AliDielectronDebugTree.h"
  67 #include "AliDielectronSignalMC.h"
  68 #include "AliDielectronMixingHandler.h"
  69 #include "AliDielectronV0Cuts.h"
  70 #include "AliDielectronPID.h"
  71
  72 #include "AliDielectron.h"
  73
  74 ClassImp(AliDielectron)
  75
  76 const char* AliDielectron::fgkTrackClassNames[4] = {
  77   "ev1+",
  78   "ev1-",
  79   "ev2+",
  80   "ev2-"
  81 };
  82
  83 const char* AliDielectron::fgkPairClassNames[11] = {
  84   "ev1+_ev1+",
  85   "ev1+_ev1-",
  86   "ev1-_ev1-",
  87   "ev1+_ev2+",
  88   "ev1-_ev2+",
  89   "ev2+_ev2+",
  90   "ev1+_ev2-",
  91   "ev1-_ev2-",
  92   "ev2+_ev2-",
  93   "ev2-_ev2-",
  94   "ev1+_ev1-_TR"
  95 };
  96
  97 //________________________________________________________________
  98 AliDielectron::AliDielectron() :
  99   TNamed("AliDielectron","AliDielectron"),
 100   fCutQA(kFALSE),
 101   fQAmonitor(0x0),
 102   fPostPIDCntrdCorr(0x0),
 103   fPostPIDWdthCorr(0x0),
 104   fLegEffMap(0x0),
 105   fPairEffMap(0x0),
 106   fEventFilter("EventFilter"),
 107   fTrackFilter("TrackFilter"),
 108   fPairPreFilter("PairPreFilter"),
 109   fPairPreFilterLegs("PairPreFilterLegs"),
 110   fPairFilter("PairFilter"),
 111   fEventPlanePreFilter("EventPlanePreFilter"),
 112   fEventPlanePOIPreFilter("EventPlanePOIPreFilter"),
 113   fPdgMother(443),
 114   fPdgLeg1(11),
 115   fPdgLeg2(11),
 116   fSignalsMC(0x0),
 117   fNoPairing(kFALSE),
 118   fProcessLS(kTRUE),
 119   fUseKF(kTRUE),
 120   fHistoArray(0x0),
 121   fHistos(0x0),
 122   fUsedVars(new TBits(AliDielectronVarManager::kNMaxValues)),
 123   fPairCandidates(new TObjArray(11)),
 124   fCfManagerPair(0x0),
 125   fTrackRotator(0x0),
 126   fDebugTree(0x0),
 127   fMixing(0x0),
 128   fPreFilterEventPlane(kFALSE),
 129   fLikeSignSubEvents(kFALSE),
 130   fPreFilterUnlikeOnly(kFALSE),
 131   fPreFilterAllSigns(kFALSE),
 132   fHasMC(kFALSE),
 133   fStoreRotatedPairs(kFALSE),
 134   fDontClearArrays(kFALSE),
 135   fEstimatorFilename(""),
 136   fTRDpidCorrectionFilename(""),
 137   fVZEROCalibrationFilename(""),
 138   fVZERORecenteringFilename(""),
 139   fZDCRecenteringFilename("")
 140
 141 {
 142   //
 143   // Default constructor
 144   //
 145
 146 }
 147
 148 //________________________________________________________________
 149 AliDielectron::AliDielectron(const char* name, const char* title) :
 150   TNamed(name,title),
 151   fCutQA(kFALSE),
 152   fQAmonitor(0x0),
 153   fPostPIDCntrdCorr(0x0),
 154   fPostPIDWdthCorr(0x0),
 155   fLegEffMap(0x0),
 156   fPairEffMap(0x0),
 157   fEventFilter("EventFilter"),
 158   fTrackFilter("TrackFilter"),
 159   fPairPreFilter("PairPreFilter"),
 160   fPairPreFilterLegs("PairPreFilterLegs"),
 161   fPairFilter("PairFilter"),
 162   fEventPlanePreFilter("EventPlanePreFilter"),
 163   fEventPlanePOIPreFilter("EventPlanePOIPreFilter"),
 164   fPdgMother(443),
 165   fPdgLeg1(11),
 166   fPdgLeg2(11),
 167   fSignalsMC(0x0),
 168   fNoPairing(kFALSE),
 169   fProcessLS(kTRUE),
 170   fUseKF(kTRUE),
 171   fHistoArray(0x0),
 172   fHistos(0x0),
 173   fUsedVars(new TBits(AliDielectronVarManager::kNMaxValues)),
 174   fPairCandidates(new TObjArray(11)),
 175   fCfManagerPair(0x0),
 176   fTrackRotator(0x0),
 177   fDebugTree(0x0),
 178   fMixing(0x0),
 179   fPreFilterEventPlane(kFALSE),
 180   fLikeSignSubEvents(kFALSE),
 181   fPreFilterUnlikeOnly(kFALSE),
 182   fPreFilterAllSigns(kFALSE),
 183   fHasMC(kFALSE),
 184   fStoreRotatedPairs(kFALSE),
 185   fDontClearArrays(kFALSE),
 186   fEstimatorFilename(""),
 187   fTRDpidCorrectionFilename(""),
 188   fVZEROCalibrationFilename(""),
 189   fVZERORecenteringFilename(""),
 190   fZDCRecenteringFilename("")
 191 {
 192   //
 193   // Named constructor
 194   //
 195
 196 }
 197
 198 //________________________________________________________________
 199 AliDielectron::~AliDielectron()
 200 {
 201   //
 202   // Default destructor
 203   //
 204   if (fQAmonitor) delete fQAmonitor;
 205   if (fPostPIDCntrdCorr) delete fPostPIDCntrdCorr;
 206   if (fPostPIDWdthCorr) delete fPostPIDWdthCorr;
 207   if (fLegEffMap) delete fLegEffMap;
 208   if (fPairEffMap) delete fPairEffMap;
 209   if (fHistos) delete fHistos;
 210   if (fUsedVars) delete fUsedVars;
 211   if (fPairCandidates) delete fPairCandidates;
 212   if (fDebugTree) delete fDebugTree;
 213   if (fMixing) delete fMixing;
 214   if (fSignalsMC) delete fSignalsMC;
 215   if (fCfManagerPair) delete fCfManagerPair;
 216   if (fHistoArray) delete fHistoArray;
 217 }
 218
 219 //________________________________________________________________
 220 void AliDielectron::Init()
 221 {
 222   //
 223   // Initialise objects
 224   //
 225
 226   if(GetHasMC()) AliDielectronMC::Instance()->SetHasMC(GetHasMC());
 227
 228   InitPairCandidateArrays();
 229
 230   if (fCfManagerPair) {
 231     fCfManagerPair->SetSignalsMC(fSignalsMC);
 232     fCfManagerPair->InitialiseContainer(fPairFilter);
 233   }
 234   if (fTrackRotator)  {
 235     fTrackRotator->SetTrackArrays(&fTracks[0],&fTracks[1]);
 236     fTrackRotator->SetPdgLegs(fPdgLeg1,fPdgLeg2);
 237   }
 238   if (fDebugTree) fDebugTree->SetDielectron(this);
 239
 240   if(fEstimatorFilename.Contains(".root"))        AliDielectronVarManager::InitEstimatorAvg(fEstimatorFilename.Data());
 241   if(fTRDpidCorrectionFilename.Contains(".root")) AliDielectronVarManager::InitTRDpidEffHistograms(fTRDpidCorrectionFilename.Data());
 242   if(fVZEROCalibrationFilename.Contains(".root")) AliDielectronVarManager::SetVZEROCalibrationFile(fVZEROCalibrationFilename.Data());
 243   if(fVZERORecenteringFilename.Contains(".root")) AliDielectronVarManager::SetVZERORecenteringFile(fVZERORecenteringFilename.Data());
 244   if(fZDCRecenteringFilename.Contains(".root")) AliDielectronVarManager::SetZDCRecenteringFile(fZDCRecenteringFilename.Data());
 245
 246   if(fLegEffMap)  AliDielectronVarManager::SetLegEffMap(fLegEffMap);
 247   if(fPairEffMap) AliDielectronVarManager::SetPairEffMap(fPairEffMap);
 248
 249
 250   if (fMixing) fMixing->Init(this);
 251   if (fHistoArray) {
 252     fHistoArray->SetSignalsMC(fSignalsMC);
 253     fHistoArray->Init();
 254   }
 255
 256   if(fPostPIDCntrdCorr) AliDielectronPID::SetCentroidCorrFunction(fPostPIDCntrdCorr);
 257   if(fPostPIDWdthCorr)  AliDielectronPID::SetWidthCorrFunction(fPostPIDWdthCorr);
 258
 259   if (fCutQA) {
 260     fQAmonitor = new AliDielectronCutQA(Form("QAcuts_%s",GetName()),"QAcuts");
 261     fQAmonitor->AddTrackFilter(&fTrackFilter);
 262     fQAmonitor->AddPairFilter(&fPairFilter);
 263     fQAmonitor->AddEventFilter(&fEventFilter);
 264     fQAmonitor->Init();
 265   }
 266
 267   if(fHistos) {
 268     (*fUsedVars)|= (*fHistos->GetUsedVars());
 269   }
 270
 271 }
 272
 273 //________________________________________________________________
 274 void AliDielectron::Process(AliVEvent *ev1, AliVEvent *ev2)
 275 {
 276   //
 277   // Process the events
 278   //
 279
 280   //at least first event is needed!
 281   if (!ev1){
 282     AliError("At least first event must be set!");
 283     return;
 284   }
 285
 286   // modify event numbers in MC so that we can identify new events
 287   // in AliDielectronV0Cuts (not neeeded for collision data)
 288   if(GetHasMC()) {
 289     ev1->SetBunchCrossNumber(1);
 290     ev1->SetOrbitNumber(1);
 291     ev1->SetPeriodNumber(1);
 292   }
 293
 294   // set event
 295   AliDielectronVarManager::SetEvent(ev1);
 296   if (fMixing){
 297     //set mixing bin to event data
 298     Int_t bin=fMixing->FindBin(AliDielectronVarManager::GetData());
 299     AliDielectronVarManager::SetValue(AliDielectronVarManager::kMixingBin,bin);
 300   }
 301
 302   //in case we have MC load the MC event and process the MC particles
 303   // why do not apply the event cuts first ????
 304   if (AliDielectronMC::Instance()->ConnectMCEvent()){
 305     ProcessMC(ev1);
 306   }
 307
 308   //if candidate array doesn't exist, create it
 309   if (!fPairCandidates->UncheckedAt(0)) {
 310     InitPairCandidateArrays();
 311   } else {
 312     ClearArrays();
 313   }
 314
 315   //mask used to require that all cuts are fulfilled
 316   UInt_t selectedMask=(1<<fEventFilter.GetCuts()->GetEntries())-1;
 317
 318   //apply event cuts
 319   UInt_t cutmask = fEventFilter.IsSelected(ev1);
 320   if(fCutQA) fQAmonitor->FillAll(ev1);
 321   if(fCutQA) fQAmonitor->Fill(cutmask,ev1);
 322   if ((ev1&&cutmask!=selectedMask) ||
 323       (ev2&&fEventFilter.IsSelected(ev2)!=selectedMask)) return;
 324
 325   //fill track arrays for the first event
 326   if (ev1){
 327     FillTrackArrays(ev1);
 328     if (((fPreFilterAllSigns)||(fPreFilterUnlikeOnly)) && ( fPairPreFilter.GetCuts()->GetEntries()>0 )) PairPreFilter(0, 1, fTracks[0], fTracks[1]);
 329   }
 330
 331
 332   //fill track arrays for the second event
 333   if (ev2) {
 334     FillTrackArrays(ev2,1);
 335     if (((fPreFilterAllSigns)||(fPreFilterUnlikeOnly)) && ( fPairPreFilter.GetCuts()->GetEntries()>0 )) PairPreFilter(2, 3, fTracks[2], fTracks[3]);
 336   }
 337
 338   // TPC event plane correction
 339   if (ev1 && fPreFilterEventPlane && ( fEventPlanePreFilter.GetCuts()->GetEntries()>0 || fEventPlanePOIPreFilter.GetCuts()->GetEntries()>0))
 340     EventPlanePreFilter(0, 1, fTracks[0], fTracks[1], ev1);
 341
 342   if (!fNoPairing){
 343     // create pairs and fill pair candidate arrays
 344     for (Int_t itrackArr1=0; itrackArr1<4; ++itrackArr1){
 345       for (Int_t itrackArr2=itrackArr1; itrackArr2<4; ++itrackArr2){
 346         if(!fProcessLS && GetPairIndex(itrackArr1,itrackArr2)!=kEv1PM) continue;
 347         FillPairArrays(itrackArr1, itrackArr2);
 348       }
 349     }
 350
 351     //track rotation
 352     if (fTrackRotator) {
 353       fTrackRotator->SetEvent(ev1);
 354       FillPairArrayTR();
 355     }
 356   }
 357
 358   //fill debug tree if a manager is attached
 359   if (fDebugTree) FillDebugTree();
 360
 361   //process event mixing
 362   if (fMixing) {
 363     fMixing->Fill(ev1,this);
 364     //     FillHistograms(0x0,kTRUE);
 365   }
 366
 367   // fill candidate variables
 368   Double_t ntracks = fTracks[0].GetEntriesFast() + fTracks[1].GetEntriesFast();
 369   Double_t npairs  = PairArray(AliDielectron::kEv1PM)->GetEntriesFast();
 370   AliDielectronVarManager::SetValue(AliDielectronVarManager::kTracks, ntracks);
 371   AliDielectronVarManager::SetValue(AliDielectronVarManager::kPairs,  npairs);
 372
 373   //in case there is a histogram manager, fill the QA histograms
 374   if (fHistos && fSignalsMC) FillMCHistograms(ev1);
 375   if (fHistos) FillHistograms(ev1);
 376
 377
 378   // clear arrays
 379   if (!fDontClearArrays) ClearArrays();
 380
 381   // reset TPC EP and unique identifiers for v0 cut class
 382   AliDielectronVarManager::SetTPCEventPlane(0x0);
 383   if(GetHasMC()) { // only for MC needed
 384     for (Int_t iCut=0; iCut<fTrackFilter.GetCuts()->GetEntries();++iCut) {
 385       if ( fTrackFilter.GetCuts()->At(iCut)->IsA() == AliDielectronV0Cuts::Class() )
 386         ((AliDielectronV0Cuts*)fTrackFilter.GetCuts()->At(iCut))->ResetUniqueEventNumbers();
 387     }
 388   }
 389
 390 }
 391
 392 //________________________________________________________________
 393 void AliDielectron::ProcessMC(AliVEvent *ev1)
 394 {
 395   //
 396   // Process the MC data
 397   //
 398
 399   AliDielectronMC *dieMC=AliDielectronMC::Instance();
 400
 401   if (fHistos) FillHistogramsMC(dieMC->GetMCEvent(), ev1);
 402
 403   // mc tracks
 404   if(!dieMC->GetNMCTracks()) return;
 405
 406   // signals to be studied
 407   if(!fSignalsMC) return;
 408   Int_t nSignals = fSignalsMC->GetEntries();
 409   if(!nSignals) return;
 410
 411   //loop over all MC data and Fill the HF, CF containers and histograms if they exist
 412   if(fCfManagerPair) fCfManagerPair->SetPdgMother(fPdgMother);
 413
 414   Bool_t bFillCF   = (fCfManagerPair ? fCfManagerPair->GetStepForMCtruth()  : kFALSE);
 415   Bool_t bFillHF   = (fHistoArray    ? fHistoArray->GetStepForMCGenerated() : kFALSE);
 416   Bool_t bFillHist = kFALSE;
 417   if(fHistos) {
 418     const THashList *histlist =  fHistos->GetHistogramList();
 419     for(Int_t isig=0;isig<nSignals;isig++) {
 420       TString sigName = fSignalsMC->At(isig)->GetName();
 421       bFillHist |= histlist->FindObject(Form("Pair_%s_MCtruth",sigName.Data()))!=0x0;
 422       bFillHist |= histlist->FindObject(Form("Track_Leg_%s_MCtruth",sigName.Data()))!=0x0;
 423       bFillHist |= histlist->FindObject(Form("Track_%s_%s_MCtruth",fgkPairClassNames[1],sigName.Data()))!=0x0;
 424       if(bFillHist) break;
 425     }
 426   }
 427   // check if there is anything to fill
 428   if(!bFillCF && !bFillHF && !bFillHist) return;
 429
 430
 431   // initialize 2D arrays of labels for particles from each MC signal
 432   Int_t** labels1;      // labels for particles satisfying branch 1
 433   Int_t** labels2;      // labels for particles satisfying branch 2
 434   Int_t** labels12;     // labels for particles satisfying both branches
 435   labels1 = new Int_t*[nSignals];
 436   labels2 = new Int_t*[nSignals];
 437   labels12 = new Int_t*[nSignals];
 438   Int_t* indexes1=new Int_t[nSignals];
 439   Int_t* indexes2=new Int_t[nSignals];
 440   Int_t* indexes12=new Int_t[nSignals];
 441   for(Int_t isig=0;isig<nSignals;++isig) {
 442     *(labels1+isig) = new Int_t[dieMC->GetNMCTracks()];
 443     *(labels2+isig) = new Int_t[dieMC->GetNMCTracks()];
 444     *(labels12+isig) = new Int_t[dieMC->GetNMCTracks()];
 445     for(Int_t ip=0; ip<dieMC->GetNMCTracks();++ip) {
 446       labels1[isig][ip] = -1;
 447       labels2[isig][ip] = -1;
 448       labels12[isig][ip] = -1;
 449     }
 450     indexes1[isig]=0;
 451     indexes2[isig]=0;
 452     indexes12[isig]=0;
 453   }
 454
 455   Bool_t truth1=kFALSE;
 456   Bool_t truth2=kFALSE;
 457   // loop over the MC tracks
 458   for(Int_t ipart=0; ipart<dieMC->GetNMCTracks(); ++ipart) {
 459     for(Int_t isig=0; isig<nSignals; ++isig) {       // loop over signals
 460       // Proceed only if this signal is required in the pure MC step
 461       // NOTE: Some signals can be satisfied by many particles and this leads to high
 462       //       computation times (e.g. secondary electrons from the GEANT transport). Be aware of this!!
 463       if(!((AliDielectronSignalMC*)fSignalsMC->At(isig))->GetFillPureMCStep()) continue;
 464
 465       truth1 = dieMC->IsMCTruth(ipart, (AliDielectronSignalMC*)fSignalsMC->At(isig), 1);
 466       truth2 = dieMC->IsMCTruth(ipart, (AliDielectronSignalMC*)fSignalsMC->At(isig), 2);
 467
 468       // particles satisfying both branches are treated separately to avoid double counting during pairing
 469       if(truth1 && truth2) {
 470         labels12[isig][indexes12[isig]] = ipart;
 471         ++indexes12[isig];
 472       }
 473       else {
 474         if(truth1) {
 475           labels1[isig][indexes1[isig]] = ipart;
 476           ++indexes1[isig];
 477         }
 478         if(truth2) {
 479           labels2[isig][indexes2[isig]] = ipart;
 480           ++indexes2[isig];
 481         }
 482       }
 483     }
 484   }  // end loop over MC particles
 485
 486   // Do the pairing and fill the CF container with pure MC info
 487   for(Int_t isig=0; isig<nSignals; ++isig) {
 488     //    printf("INDEXES: %d-%d both%d\n",indexes1[isig],indexes2[isig],indexes12[isig]);
 489     // mix the particles which satisfy only one of the signal branches
 490     for(Int_t i1=0;i1<indexes1[isig];++i1) {
 491       if(!indexes2[isig]) FillMCHistograms(labels1[isig][i1], -1, isig); // (e.g. single electrons only, no pairs)
 492       for(Int_t i2=0;i2<indexes2[isig];++i2) {
 493         if(bFillCF) fCfManagerPair->FillMC(labels1[isig][i1], labels2[isig][i2], isig);
 494         if(bFillHF) fHistoArray->Fill(labels1[isig][i1], labels2[isig][i2], isig);
 495         FillMCHistograms(labels1[isig][i1], labels2[isig][i2], isig);
 496       }
 497     }
 498     // mix the particles which satisfy both branches
 499     for(Int_t i1=0;i1<indexes12[isig];++i1) {
 500       for(Int_t i2=0; i2<i1; ++i2) {
 501         if(bFillCF) fCfManagerPair->FillMC(labels12[isig][i1], labels12[isig][i2], isig);
 502         if(bFillHF) fHistoArray->Fill(labels12[isig][i1], labels12[isig][i2], isig);
 503         FillMCHistograms(labels12[isig][i1], labels12[isig][i2], isig);
 504       }
 505     }
 506   }    // end loop over signals
 507
 508   // release the memory
 509   for(Int_t isig=0;isig<nSignals;++isig) {
 510     delete [] *(labels1+isig);
 511     delete [] *(labels2+isig);
 512     delete [] *(labels12+isig);
 513   }
 514   delete [] labels1;
 515   delete [] labels2;
 516   delete [] labels12;
 517   delete [] indexes1;
 518   delete [] indexes2;
 519   delete [] indexes12;
 520 }
 521
 522 //________________________________________________________________
 523 void AliDielectron::FillHistogramsTracks(TObjArray **tracks)
 524 {
 525   //
 526   // Fill Histogram information for tracks after prefilter
 527   // ignore mixed events - for prefilter, only single tracks +/- are relevant
 528   //
 529
 530   TString  className,className2;
 531   Double_t values[AliDielectronVarManager::kNMaxValues];
 532   AliDielectronVarManager::SetFillMap(fUsedVars);
 533
 534   //Fill track information, separately for the track array candidates
 535   for (Int_t i=0; i<2; ++i){
 536     className.Form("Pre_%s",fgkTrackClassNames[i]);
 537     if (!fHistos->GetHistogramList()->FindObject(className.Data())) continue;
 538     Int_t ntracks=tracks[i]->GetEntriesFast();
 539     for (Int_t itrack=0; itrack<ntracks; ++itrack){
 540       AliDielectronVarManager::Fill(tracks[i]->UncheckedAt(itrack), values);
 541       fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
 542     }
 543   }
 544 }
 545
 546
 547 //________________________________________________________________
 548 void AliDielectron::FillHistogramsMC(const AliMCEvent *ev, AliVEvent *ev1)
 549 {
 550   //
 551   // Fill Histogram information for MCEvents
 552   //
 553
 554   Double_t values[AliDielectronVarManager::kNMaxValues]={0.};
 555   AliDielectronVarManager::SetFillMap(fUsedVars);
 556
 557   // Fill event information
 558   AliDielectronVarManager::Fill(ev1, values);    // ESD/AOD information
 559   AliDielectronVarManager::Fill(ev, values);     // MC truth info
 560   if (fHistos->GetHistogramList()->FindObject("MCEvent"))
 561     fHistos->FillClass("MCEvent", AliDielectronVarManager::kNMaxValues, values);
 562 }
 563
 564
 565 //________________________________________________________________
 566 void AliDielectron::FillHistograms(const AliVEvent *ev, Bool_t pairInfoOnly)
 567 {
 568   //
 569   // Fill Histogram information for tracks and pairs
 570   //
 571
 572   TString  className,className2;
 573   Double_t values[AliDielectronVarManager::kNMaxValues]={0.};
 574   AliDielectronVarManager::SetFillMap(fUsedVars);
 575
 576   //Fill event information
 577   if (ev){
 578     if (fHistos->GetHistogramList()->FindObject("Event")) {
 579       fHistos->FillClass("Event", AliDielectronVarManager::kNMaxValues, AliDielectronVarManager::GetData());
 580     }
 581   }
 582
 583   //Fill track information, separately for the track array candidates
 584   if (!pairInfoOnly){
 585     className2.Form("Track_%s",fgkPairClassNames[1]);  // unlike sign, SE only
 586     for (Int_t i=0; i<4; ++i){
 587       className.Form("Track_%s",fgkTrackClassNames[i]);
 588       Bool_t mergedtrkClass=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
 589       Bool_t trkClass=fHistos->GetHistogramList()->FindObject(className.Data())!=0x0;
 590       if (!trkClass && !mergedtrkClass) continue;
 591       Int_t ntracks=fTracks[i].GetEntriesFast();
 592       for (Int_t itrack=0; itrack<ntracks; ++itrack){
 593         AliDielectronVarManager::Fill(fTracks[i].UncheckedAt(itrack), values);
 594         if(trkClass)
 595           fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
 596         if(mergedtrkClass && i<2)
 597           fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values); //only ev1
 598       }
 599     }
 600   }
 601
 602   //Fill Pair information, separately for all pair candidate arrays and the legs
 603   TObjArray arrLegs(100);
 604   for (Int_t i=0; i<10; ++i){
 605     className.Form("Pair_%s",fgkPairClassNames[i]);
 606     className2.Form("Track_Legs_%s",fgkPairClassNames[i]);
 607     Bool_t pairClass=fHistos->GetHistogramList()->FindObject(className.Data())!=0x0;
 608     Bool_t legClass=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
 609     if (!pairClass&&!legClass) continue;
 610     Int_t ntracks=PairArray(i)->GetEntriesFast();
 611     for (Int_t ipair=0; ipair<ntracks; ++ipair){
 612       AliDielectronPair *pair=static_cast<AliDielectronPair*>(PairArray(i)->UncheckedAt(ipair));
 613
 614       //fill pair information
 615       if (pairClass){
 616         AliDielectronVarManager::Fill(pair, values);
 617         fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
 618       }
 619
 620       //fill leg information, don't fill the information twice
 621       if (legClass){
 622         AliVParticle *d1=pair->GetFirstDaughter();
 623         AliVParticle *d2=pair->GetSecondDaughter();
 624         if (!arrLegs.FindObject(d1)){
 625           AliDielectronVarManager::Fill(d1, values);
 626           fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
 627           arrLegs.Add(d1);
 628         }
 629         if (!arrLegs.FindObject(d2)){
 630           AliDielectronVarManager::Fill(d2, values);
 631           fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
 632           arrLegs.Add(d2);
 633         }
 634       }
 635     }
 636     if (legClass) arrLegs.Clear();
 637   }
 638
 639 }
 640 //________________________________________________________________
 641 void AliDielectron::FillHistogramsPair(AliDielectronPair *pair,Bool_t fromPreFilter/*=kFALSE*/)
 642 {
 643   //
 644   // Fill Histogram information for pairs and the track in the pair
 645   // NOTE: in this funtion the leg information may be filled multiple
 646   //       times. This funtion is used in the track rotation pairing
 647   //       and those legs are not saved!
 648   //
 649   TString  className,className2;
 650   Double_t values[AliDielectronVarManager::kNMaxValues];
 651   AliDielectronVarManager::SetFillMap(fUsedVars);
 652
 653   //Fill Pair information, separately for all pair candidate arrays and the legs
 654   TObjArray arrLegs(100);
 655   const Int_t type=pair->GetType();
 656   if (fromPreFilter) {
 657     className.Form("RejPair_%s",fgkPairClassNames[type]);
 658     className2.Form("RejTrack_%s",fgkPairClassNames[type]);
 659   } else {
 660     className.Form("Pair_%s",fgkPairClassNames[type]);
 661     className2.Form("Track_Legs_%s",fgkPairClassNames[type]);
 662   }
 663
 664   Bool_t pairClass=fHistos->GetHistogramList()->FindObject(className.Data())!=0x0;
 665   Bool_t legClass=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
 666
 667   //fill pair information
 668   if (pairClass){
 669     AliDielectronVarManager::Fill(pair, values);
 670     fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
 671   }
 672
 673   if (legClass){
 674     AliVParticle *d1=pair->GetFirstDaughter();
 675     AliDielectronVarManager::Fill(d1, values);
 676     fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
 677
 678     AliVParticle *d2=pair->GetSecondDaughter();
 679     AliDielectronVarManager::Fill(d2, values);
 680     fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
 681   }
 682 }
 683
 684 //________________________________________________________________
 685 void AliDielectron::FillTrackArrays(AliVEvent * const ev, Int_t eventNr)
 686 {
 687   //
 688   // select tracks and fill track candidate arrays
 689   // eventNr = 0: First  event, use track arrays 0 and 1
 690   // eventNr = 1: Second event, use track arrays 2 and 3
 691   //
 692
 693   Int_t ntracks=ev->GetNumberOfTracks();
 694
 695   UInt_t selectedMask=(1<<fTrackFilter.GetCuts()->GetEntries())-1;
 696   for (Int_t itrack=0; itrack<ntracks; ++itrack){
 697     //get particle
 698     AliVParticle *particle=ev->GetTrack(itrack);
 699
 700     //apply track cuts
 701     UInt_t cutmask=fTrackFilter.IsSelected(particle);
 702     //fill cut QA
 703     if(fCutQA) fQAmonitor->FillAll(particle);
 704     if(fCutQA) fQAmonitor->Fill(cutmask,particle);
 705
 706     if (cutmask!=selectedMask) continue;
 707
 708     //fill selected particle into the corresponding track arrays
 709     Short_t charge=particle->Charge();
 710     if (charge>0)      fTracks[eventNr*2].Add(particle);
 711     else if (charge<0) fTracks[eventNr*2+1].Add(particle);
 712
 713
 714   }
 715 }
 716
 717 //________________________________________________________________
 718 void AliDielectron::EventPlanePreFilter(Int_t arr1, Int_t arr2, TObjArray arrTracks1, TObjArray arrTracks2, const AliVEvent *ev)
 719 {
 720   //
 721   // Prefilter tracks and tracks from pairs
 722   // Needed for rejection in the Q-Vector of the event plane
 723   // remove contribution of all tracks to the Q-vector that are in invariant mass window
 724   //
 725
 726   AliEventplane *evplane = const_cast<AliVEvent *>(ev)->GetEventplane();
 727   if(!evplane) { // nanoAODs , here we do NOT have sub event reaction planes
 728     //  if(1) {
 729     // get the EPselectionTask for recalculation of weighting factors
 730     AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 731     AliEPSelectionTask *eptask = dynamic_cast<AliEPSelectionTask *>(man->GetTask("EventplaneSelection"));
 732     if(!eptask) return;
 733
 734     // track mapping
 735     TMap mapRemovedTracks;
 736
 737
 738     Double_t cQX=0., cQY=0.;
 739     // apply cuts to the tracks, e.g. etagap
 740     if(fEventPlanePreFilter.GetCuts()->GetEntries()) {
 741       UInt_t selectedMask=(1<<fEventPlanePreFilter.GetCuts()->GetEntries())-1;
 742       Int_t ntracks=ev->GetNumberOfTracks();
 743       for (Int_t itrack=0; itrack<ntracks; ++itrack){
 744         AliVParticle *particle=ev->GetTrack(itrack);
 745         AliVTrack *track= static_cast<AliVTrack*>(particle);
 746         if (!track) continue;
 747         //event plane cuts
 748         UInt_t cutMask=fEventPlanePreFilter.IsSelected(track);
 749         //apply cut
 750         if (cutMask==selectedMask) continue;
 751
 752         mapRemovedTracks.Add(track,track);
 753         cQX += (eptask->GetWeight(track) * TMath::Cos(2*track->Phi()));
 754         cQY += (eptask->GetWeight(track) * TMath::Sin(2*track->Phi()));
 755       }
 756     }
 757
 758     // POI (particle of interest) rejection
 759     Int_t pairIndex=GetPairIndex(arr1,arr2);
 760
 761     Int_t ntrack1=arrTracks1.GetEntriesFast();
 762     Int_t ntrack2=arrTracks2.GetEntriesFast();
 763     AliDielectronPair candidate;
 764     candidate.SetKFUsage(fUseKF);
 765
 766     UInt_t selectedMask=(1<<fEventPlanePOIPreFilter.GetCuts()->GetEntries())-1;
 767     for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1){
 768       Int_t end=ntrack2;
 769       if (arr1==arr2) end=itrack1;
 770       Bool_t accepted=kFALSE;
 771       for (Int_t itrack2=0; itrack2<end; ++itrack2){
 772         TObject *track1=arrTracks1.UncheckedAt(itrack1);
 773         TObject *track2=arrTracks2.UncheckedAt(itrack2);
 774         if (!track1 || !track2) continue;
 775         //create the pair
 776         candidate.SetTracks(static_cast<AliVTrack*>(track1), fPdgLeg1,
 777                             static_cast<AliVTrack*>(track2), fPdgLeg2);
 778         candidate.SetType(pairIndex);
 779         candidate.SetLabel(AliDielectronMC::Instance()->GetLabelMotherWithPdg(&candidate,fPdgMother));
 780
 781         //event plane pair cuts
 782         UInt_t cutMask=fEventPlanePOIPreFilter.IsSelected(&candidate);
 783         //apply cut
 784         if (cutMask==selectedMask) continue;
 785
 786         accepted=kTRUE;
 787         //remove the tracks from the Track arrays
 788         arrTracks2.AddAt(0x0,itrack2);
 789       }
 790       if ( accepted ) arrTracks1.AddAt(0x0,itrack1);
 791     }
 792     //compress the track arrays
 793     arrTracks1.Compress();
 794     arrTracks2.Compress();
 795
 796     //Modify the components: subtract the tracks
 797     ntrack1=arrTracks1.GetEntriesFast();
 798     ntrack2=arrTracks2.GetEntriesFast();
 799     // remove leg1 contribution
 800     for (Int_t itrack=0; itrack<ntrack1; ++itrack){
 801       AliVTrack *track= static_cast<AliVTrack*>(arrTracks1.UncheckedAt(itrack));
 802       if (!track) continue;
 803       // track contribution was already removed
 804       if (mapRemovedTracks.FindObject(track)) continue;
 805       else mapRemovedTracks.Add(track,track);
 806
 807       cQX += (eptask->GetWeight(track) * TMath::Cos(2*track->Phi()));
 808       cQY += (eptask->GetWeight(track) * TMath::Sin(2*track->Phi()));
 809     }
 810     // remove leg2 contribution
 811     for (Int_t itrack=0; itrack<ntrack2; ++itrack){
 812       AliVTrack *track= static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack));
 813       if (!track) continue;
 814       // track contribution was already removed
 815       if (mapRemovedTracks.FindObject(track)) continue;
 816       else mapRemovedTracks.Add(track,track);
 817
 818       cQX += (eptask->GetWeight(track) * TMath::Cos(2*track->Phi()));
 819       cQY += (eptask->GetWeight(track) * TMath::Sin(2*track->Phi()));
 820     }
 821
 822     // build a corrected alieventplane using the values from the var manager
 823     // these uncorrected values are filled using the stored magnitude and angle  in the header
 824     TVector2 qcorr;
 825     qcorr.Set(AliDielectronVarManager::GetValue(AliDielectronVarManager::kTPCxH2uc)-cQX,
 826               AliDielectronVarManager::GetValue(AliDielectronVarManager::kTPCyH2uc)-cQY);
 827     // fill alieventplane
 828     AliEventplane cevplane;
 829     cevplane.SetQVector(&qcorr);
 830     AliDielectronVarManager::SetTPCEventPlane(&cevplane);
 831     cevplane.SetQVector(0);
 832     return;
 833   } //end: nanoAODs
 834   else
 835     {
 836     // this is done in case of ESDs or AODs
 837     Bool_t isESD=(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->IsA()==AliESDInputHandler::Class());
 838     // copy event plane object
 839     AliEventplane cevplane(*evplane);
 840     //    Int_t nMaxID=cevplane->GetQContributionXArray()->GetSize();
 841
 842     TVector2 *qcorr  = cevplane.GetQVector();
 843     if(!qcorr) return;
 844     TVector2 *qcsub1 = 0x0;
 845     TVector2 *qcsub2 = 0x0;
 846
 847     // eta gap ?
 848     Bool_t etagap = kFALSE;
 849     for (Int_t iCut=0; iCut<fEventPlanePreFilter.GetCuts()->GetEntries();++iCut) {
 850       TString cutName=fEventPlanePreFilter.GetCuts()->At(iCut)->GetName();
 851       if(cutName.Contains("eta") || cutName.Contains("Eta"))  etagap=kTRUE;
 852     }
 853
 854     // subevent configuration for eta gap or LS (default is rndm)
 855     if(fLikeSignSubEvents && etagap) {
 856       // start with the full Qvector/event in both sub events
 857       qcsub1 = new TVector2(*qcorr);
 858       qcsub2 = new TVector2(*qcorr);
 859       cevplane.SetQsub(qcsub1,qcsub2);
 860
 861       Int_t ntracks=ev->GetNumberOfTracks();
 862       // track removals
 863       for (Int_t itrack=0; itrack<ntracks; ++itrack){
 864         AliVParticle *particle=ev->GetTrack(itrack);
 865         AliVTrack *track= static_cast<AliVTrack*>(particle);
 866         if (!track) continue;
 867         if (track->GetID()>=0 && !isESD) continue;
 868         Int_t tmpID = isESD ? track->GetID() : track->GetID()*-1 - 1;
 869
 870         // set contributions to zero
 871         // charge sub1+ sub2-
 872         if(fLikeSignSubEvents) {
 873           Short_t charge=track->Charge();
 874           if (charge<0) {
 875             cevplane.GetQContributionXArraysub1()->SetAt(0.0, tmpID);
 876             cevplane.GetQContributionYArraysub1()->SetAt(0.0, tmpID);
 877           }
 878           if (charge>0) {
 879             cevplane.GetQContributionXArraysub2()->SetAt(0.0, tmpID);
 880             cevplane.GetQContributionYArraysub2()->SetAt(0.0, tmpID);
 881           }
 882         }
 883         // eta sub1+ sub2-
 884         if(etagap) {
 885           Double_t eta=track->Eta();
 886           if (eta<0.0) {
 887             cevplane.GetQContributionXArraysub1()->SetAt(0.0, tmpID);
 888             cevplane.GetQContributionYArraysub1()->SetAt(0.0, tmpID);
 889           }
 890           if (eta>0.0) {
 891             cevplane.GetQContributionXArraysub2()->SetAt(0.0, tmpID);
 892             cevplane.GetQContributionYArraysub2()->SetAt(0.0, tmpID);
 893           }
 894         }
 895       } // end: loop over tracks
 896     } // end: sub event configuration
 897
 898     // apply cuts, e.g. etagap
 899     if(fEventPlanePreFilter.GetCuts()->GetEntries()) {
 900       UInt_t selectedMask=(1<<fEventPlanePreFilter.GetCuts()->GetEntries())-1;
 901       Int_t ntracks=ev->GetNumberOfTracks();
 902       for (Int_t itrack=0; itrack<ntracks; ++itrack){
 903         AliVParticle *particle=ev->GetTrack(itrack);
 904         AliVTrack *track= static_cast<AliVTrack*>(particle);
 905         if (!track) continue;
 906         if (track->GetID()>=0 && !isESD) continue;
 907         Int_t tmpID = isESD ? track->GetID() : track->GetID()*-1 - 1;
 908
 909         //event plane cuts
 910         UInt_t cutMask=fEventPlanePreFilter.IsSelected(track);
 911         //apply cut
 912         if (cutMask==selectedMask) continue;
 913
 914         // set contributions to zero
 915         cevplane.GetQContributionXArray()->SetAt(0.0, tmpID);
 916         cevplane.GetQContributionYArray()->SetAt(0.0, tmpID);
 917         cevplane.GetQContributionXArraysub1()->SetAt(0.0, tmpID);
 918         cevplane.GetQContributionYArraysub1()->SetAt(0.0, tmpID);
 919         cevplane.GetQContributionXArraysub2()->SetAt(0.0, tmpID);
 920         cevplane.GetQContributionYArraysub2()->SetAt(0.0, tmpID);
 921       }
 922     } // end: track cuts
 923
 924     // POI (particle of interest) rejection
 925     Int_t pairIndex=GetPairIndex(arr1,arr2);
 926     Int_t ntrack1=arrTracks1.GetEntriesFast();
 927     Int_t ntrack2=arrTracks2.GetEntriesFast();
 928     AliDielectronPair candidate;
 929     candidate.SetKFUsage(fUseKF);
 930
 931     UInt_t selectedMask=(1<<fEventPlanePOIPreFilter.GetCuts()->GetEntries())-1;
 932     for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1){
 933       Int_t end=ntrack2;
 934       if (arr1==arr2) end=itrack1;
 935       Bool_t accepted=kFALSE;
 936       for (Int_t itrack2=0; itrack2<end; ++itrack2){
 937         TObject *track1=arrTracks1.UncheckedAt(itrack1);
 938         TObject *track2=arrTracks2.UncheckedAt(itrack2);
 939         if (!track1 || !track2) continue;
 940         //create the pair
 941         candidate.SetTracks(static_cast<AliVTrack*>(track1), fPdgLeg1,
 942                             static_cast<AliVTrack*>(track2), fPdgLeg2);
 943
 944         candidate.SetType(pairIndex);
 945         candidate.SetLabel(AliDielectronMC::Instance()->GetLabelMotherWithPdg(&candidate,fPdgMother));
 946
 947         //event plane cuts
 948         UInt_t cutMask=fEventPlanePOIPreFilter.IsSelected(&candidate);
 949         //apply cut
 950         if (cutMask==selectedMask) continue;
 951
 952         accepted=kTRUE;
 953         //remove the tracks from the Track arrays
 954         arrTracks2.AddAt(0x0,itrack2);
 955       }
 956       if ( accepted ) arrTracks1.AddAt(0x0,itrack1);
 957     }
 958     //compress the track arrays
 959     arrTracks1.Compress();
 960     arrTracks2.Compress();
 961
 962     //Modify the components: subtract the tracks
 963     ntrack1=arrTracks1.GetEntriesFast();
 964     ntrack2=arrTracks2.GetEntriesFast();
 965     // remove leg1 contribution
 966     for (Int_t itrack=0; itrack<ntrack1; ++itrack){
 967       AliVTrack *track= static_cast<AliVTrack*>(arrTracks1.UncheckedAt(itrack));
 968       if (!track) continue;
 969       if (track->GetID()>=0 && !isESD) continue;
 970       Int_t tmpID = isESD ? track->GetID() : track->GetID()*-1 - 1;
 971       // set contributions to zero
 972       cevplane.GetQContributionXArray()->SetAt(0.0, tmpID);
 973       cevplane.GetQContributionYArray()->SetAt(0.0, tmpID);
 974       cevplane.GetQContributionXArraysub1()->SetAt(0.0, tmpID);
 975       cevplane.GetQContributionYArraysub1()->SetAt(0.0, tmpID);
 976       cevplane.GetQContributionXArraysub2()->SetAt(0.0, tmpID);
 977       cevplane.GetQContributionYArraysub2()->SetAt(0.0, tmpID);
 978     }
 979     // remove leg2 contribution
 980     for (Int_t itrack=0; itrack<ntrack2; ++itrack){
 981       AliVTrack *track= static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack));
 982       if (!track) continue;
 983       if (track->GetID()>=0 && !isESD) continue;
 984       Int_t tmpID = isESD ? track->GetID() : track->GetID()*-1 - 1;
 985       // set contributions to zero
 986       cevplane.GetQContributionXArray()->SetAt(0.0, tmpID);
 987       cevplane.GetQContributionYArray()->SetAt(0.0, tmpID);
 988       cevplane.GetQContributionXArraysub1()->SetAt(0.0, tmpID);
 989       cevplane.GetQContributionYArraysub1()->SetAt(0.0, tmpID);
 990       cevplane.GetQContributionXArraysub2()->SetAt(0.0, tmpID);
 991       cevplane.GetQContributionYArraysub2()->SetAt(0.0, tmpID);
 992     }
 993
 994     // set corrected AliEventplane and fill variables with corrected values
 995     AliDielectronVarManager::SetTPCEventPlane(&cevplane);
 996     delete qcsub1;
 997     delete qcsub2;
 998   } // end: ESD or AOD case
 999
1000 }
1001 //________________________________________________________________
1002 void AliDielectron::PairPreFilter(Int_t arr1, Int_t arr2, TObjArray &arrTracks1, TObjArray &arrTracks2)
1003 {
1004   //
1005   // Prefilter tracks from pairs
1006   // Needed for datlitz rejections
1007   // remove all tracks from the Single track arrays that pass the cuts in this filter
1008   //
1009
1010   Int_t ntrack1=arrTracks1.GetEntriesFast();
1011   Int_t ntrack2=arrTracks2.GetEntriesFast();
1012   AliDielectronPair candidate;
1013   candidate.SetKFUsage(fUseKF);
1014   // flag arrays for track removal
1015   Bool_t *bTracks1 = new Bool_t[ntrack1];
1016   for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1) bTracks1[itrack1]=kFALSE;
1017   Bool_t *bTracks2 = new Bool_t[ntrack2];
1018   for (Int_t itrack2=0; itrack2<ntrack2; ++itrack2) bTracks2[itrack2]=kFALSE;
1019
1020   UInt_t selectedMask=(1<<fPairPreFilter.GetCuts()->GetEntries())-1;
1021   UInt_t selectedMaskPair=(1<<fPairFilter.GetCuts()->GetEntries())-1;
1022
1023   Int_t nRejPasses = 1; //for fPreFilterUnlikeOnly and no set flag
1024   if (fPreFilterAllSigns) nRejPasses = 3;
1025
1026   for (Int_t iRP=0; iRP < nRejPasses; ++iRP) {
1027         Int_t arr1RP=arr1, arr2RP=arr2;
1028         TObjArray *arrTracks1RP=&arrTracks1;
1029         TObjArray *arrTracks2RP=&arrTracks2;
1030         Bool_t *bTracks1RP = bTracks1;
1031         Bool_t *bTracks2RP = bTracks2;
1032         switch (iRP) {
1033                 case 1: arr1RP=arr1;arr2RP=arr1;
1034                                 arrTracks1RP=&arrTracks1;
1035                                 arrTracks2RP=&arrTracks1;
1036                                 bTracks1RP = bTracks1;
1037                                 bTracks2RP = bTracks1;
1038                                 break;
1039                 case 2: arr1RP=arr2;arr2RP=arr2;
1040                                 arrTracks1RP=&arrTracks2;
1041                                 arrTracks2RP=&arrTracks2;
1042                                 bTracks1RP = bTracks2;
1043                                 bTracks2RP = bTracks2;
1044                                 break;
1045                 default: ;//nothing to do
1046         }
1047         Int_t ntrack1RP=(*arrTracks1RP).GetEntriesFast();
1048         Int_t ntrack2RP=(*arrTracks2RP).GetEntriesFast();
1049
1050         Int_t pairIndex=GetPairIndex(arr1RP,arr2RP);
1051
1052         for (Int_t itrack1=0; itrack1<ntrack1RP; ++itrack1){
1053           Int_t end=ntrack2RP;
1054           if (arr1RP==arr2RP) end=itrack1;
1055           for (Int_t itrack2=0; itrack2<end; ++itrack2){
1056                 TObject *track1=(*arrTracks1RP).UncheckedAt(itrack1);
1057                 TObject *track2=(*arrTracks2RP).UncheckedAt(itrack2);
1058                 if (!track1 || !track2) continue;
1059                 //create the pair
1060                 candidate.SetTracks(static_cast<AliVTrack*>(track1), fPdgLeg1,
1061                         static_cast<AliVTrack*>(track2), fPdgLeg2);
1062
1063                 candidate.SetType(pairIndex);
1064                 candidate.SetLabel(AliDielectronMC::Instance()->GetLabelMotherWithPdg(&candidate,fPdgMother));
1065                 //relate to the production vertex
1066                 //       if (AliDielectronVarManager::GetKFVertex()) candidate.SetProductionVertex(*AliDielectronVarManager::GetKFVertex());
1067
1068                 //pair cuts
1069                 UInt_t cutMask=fPairPreFilter.IsSelected(&candidate);
1070
1071                 //apply cut
1072                 if (cutMask!=selectedMask) continue;
1073                 if (fCfManagerPair) fCfManagerPair->Fill(selectedMaskPair+1 ,&candidate);
1074                 if (fHistos) FillHistogramsPair(&candidate,kTRUE);
1075                 //set flags for track removal
1076                 bTracks1RP[itrack1]=kTRUE;
1077                 bTracks2RP[itrack2]=kTRUE;
1078           }
1079         }
1080   }
1081
1082   //remove the tracks from the Track arrays
1083   for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1){
1084     if(bTracks1[itrack1]) arrTracks1.AddAt(0x0, itrack1);
1085   }
1086   for (Int_t itrack2=0; itrack2<ntrack2; ++itrack2){
1087     if(bTracks2[itrack2]) arrTracks2.AddAt(0x0, itrack2);
1088   }
1089
1090   // clean up
1091   delete [] bTracks1;
1092   delete [] bTracks2;
1093
1094   //compress the track arrays
1095   arrTracks1.Compress();
1096   arrTracks2.Compress();
1097
1098   //apply leg cuts after the pre filter
1099   if ( fPairPreFilterLegs.GetCuts()->GetEntries()>0 ) {
1100     selectedMask=(1<<fPairPreFilterLegs.GetCuts()->GetEntries())-1;
1101     //loop over tracks from array 1
1102     for (Int_t itrack=0; itrack<arrTracks1.GetEntriesFast();++itrack){
1103       //test cuts
1104       UInt_t cutMask=fPairPreFilterLegs.IsSelected(arrTracks1.UncheckedAt(itrack));
1105
1106       //apply cut
1107       if (cutMask!=selectedMask) arrTracks1.AddAt(0x0,itrack);
1108     }
1109     arrTracks1.Compress();
1110
1111     //in case of like sign don't loop over second array
1112     if (arr1==arr2) {
1113       arrTracks2=arrTracks1;
1114     } else {
1115
1116       //loop over tracks from array 2
1117       for (Int_t itrack=0; itrack<arrTracks2.GetEntriesFast();++itrack){
1118       //test cuts
1119         UInt_t cutMask=fPairPreFilterLegs.IsSelected(arrTracks2.UncheckedAt(itrack));
1120       //apply cut
1121         if (cutMask!=selectedMask) arrTracks2.AddAt(0x0,itrack);
1122       }
1123       arrTracks2.Compress();
1124
1125     }
1126   }
1127   //For unlike-sign monitor track-cuts:
1128   if (arr1!=arr2&&fHistos) {
1129     TObjArray *unlikesignArray[2] = {&arrTracks1,&arrTracks2};
1130     FillHistogramsTracks(unlikesignArray);
1131   }
1132 }
1133
1134 //________________________________________________________________
1135 void AliDielectron::FillPairArrays(Int_t arr1, Int_t arr2)
1136 {
1137   //
1138   // select pairs and fill pair candidate arrays
1139   //
1140
1141   TObjArray arrTracks1=fTracks[arr1];
1142   TObjArray arrTracks2=fTracks[arr2];
1143
1144   //process pre filter if set
1145   if ((!fPreFilterAllSigns) && (!fPreFilterUnlikeOnly) && ( fPairPreFilter.GetCuts()->GetEntries()>0 ))  PairPreFilter(arr1, arr2, arrTracks1, arrTracks2);
1146
1147   Int_t pairIndex=GetPairIndex(arr1,arr2);
1148
1149   Int_t ntrack1=arrTracks1.GetEntriesFast();
1150   Int_t ntrack2=arrTracks2.GetEntriesFast();
1151
1152   AliDielectronPair *candidate=new AliDielectronPair;
1153   candidate->SetKFUsage(fUseKF);
1154
1155   UInt_t selectedMask=(1<<fPairFilter.GetCuts()->GetEntries())-1;
1156
1157   for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1){
1158     Int_t end=ntrack2;
1159     if (arr1==arr2) end=itrack1;
1160     for (Int_t itrack2=0; itrack2<end; ++itrack2){
1161       //create the pair
1162       candidate->SetTracks(static_cast<AliVTrack*>(arrTracks1.UncheckedAt(itrack1)), fPdgLeg1,
1163                            static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack2)), fPdgLeg2);
1164       candidate->SetType(pairIndex);
1165       Int_t label=AliDielectronMC::Instance()->GetLabelMotherWithPdg(candidate,fPdgMother);
1166       candidate->SetLabel(label);
1167       if (label>-1) candidate->SetPdgCode(fPdgMother);
1168       else candidate->SetPdgCode(0);
1169
1170       // check for gamma kf particle
1171       label=AliDielectronMC::Instance()->GetLabelMotherWithPdg(candidate,22);
1172       if (label>-1) {
1173         candidate->SetGammaTracks(static_cast<AliVTrack*>(arrTracks1.UncheckedAt(itrack1)), fPdgLeg1,
1174                                   static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack2)), fPdgLeg2);
1175       // should we set the pdgmothercode and the label
1176       }
1177
1178       //pair cuts
1179       UInt_t cutMask=fPairFilter.IsSelected(candidate);
1180
1181       //CF manager for the pair
1182       if (fCfManagerPair) fCfManagerPair->Fill(cutMask,candidate);
1183
1184       // cut qa
1185       if(pairIndex==kEv1PM && fCutQA) {
1186         fQAmonitor->FillAll(candidate);
1187         fQAmonitor->Fill(cutMask,candidate);
1188       }
1189
1190       //apply cut
1191       if (cutMask!=selectedMask) continue;
1192
1193       //histogram array for the pair
1194       if (fHistoArray) fHistoArray->Fill(pairIndex,candidate);
1195
1196       //add the candidate to the candidate array
1197       PairArray(pairIndex)->Add(candidate);
1198       //get a new candidate
1199       candidate=new AliDielectronPair;
1200           candidate->SetKFUsage(fUseKF);
1201     }
1202   }
1203   //delete the surplus candidate
1204   delete candidate;
1205 }
1206
1207 //________________________________________________________________
1208 void AliDielectron::FillPairArrayTR()
1209 {
1210   //
1211   // select pairs and fill pair candidate arrays
1212   //
1213   UInt_t selectedMask=(1<<fPairFilter.GetCuts()->GetEntries())-1;
1214
1215   while ( fTrackRotator->NextCombination() ){
1216     AliDielectronPair candidate;
1217     candidate.SetKFUsage(fUseKF);
1218     candidate.SetTracks(&fTrackRotator->GetKFTrackP(), &fTrackRotator->GetKFTrackN(),
1219                         fTrackRotator->GetVTrackP(),fTrackRotator->GetVTrackN());
1220     candidate.SetType(kEv1PMRot);
1221
1222     //pair cuts
1223     UInt_t cutMask=fPairFilter.IsSelected(&candidate);
1224
1225     //CF manager for the pair
1226     if (fCfManagerPair) fCfManagerPair->Fill(cutMask,&candidate);
1227
1228     //apply cut
1229     if (cutMask==selectedMask) {
1230
1231       //histogram array for the pair
1232       if (fHistoArray) fHistoArray->Fill((Int_t)kEv1PMRot,&candidate);
1233
1234       if(fHistos) FillHistogramsPair(&candidate);
1235       if(fStoreRotatedPairs) PairArray(kEv1PMRot)->Add(new AliDielectronPair(candidate));
1236     }
1237   }
1238 }
1239
1240 //________________________________________________________________
1241 void AliDielectron::FillDebugTree()
1242 {
1243   //
1244   // Fill Histogram information for tracks and pairs
1245   //
1246
1247   //Fill Debug tree
1248   for (Int_t i=0; i<10; ++i){
1249     Int_t ntracks=PairArray(i)->GetEntriesFast();
1250     for (Int_t ipair=0; ipair<ntracks; ++ipair){
1251       fDebugTree->Fill(static_cast<AliDielectronPair*>(PairArray(i)->UncheckedAt(ipair)));
1252     }
1253   }
1254 }
1255
1256 //________________________________________________________________
1257 void AliDielectron::SaveDebugTree()
1258 {
1259   //
1260   // delete the debug tree, this will also write the tree
1261   //
1262   if (fDebugTree) fDebugTree->DeleteStreamer();
1263 }
1264
1265
1266 //__________________________________________________________________
1267 void AliDielectron::AddSignalMC(AliDielectronSignalMC* signal) {
1268   //
1269   //  Add an MC signal to the signals list
1270   //
1271   if(!fSignalsMC) {
1272     fSignalsMC = new TObjArray();
1273     fSignalsMC->SetOwner();
1274   }
1275   fSignalsMC->Add(signal);
1276 }
1277
1278 //________________________________________________________________
1279 void AliDielectron::FillMCHistograms(Int_t label1, Int_t label2, Int_t nSignal) {
1280   //
1281   // fill QA MC TRUTH histograms for pairs and legs of all added mc signals
1282   //
1283
1284   TString className,className2,className3;
1285   className.Form("Pair_%s_MCtruth",fSignalsMC->At(nSignal)->GetName());
1286   className2.Form("Track_Legs_%s_MCtruth",fSignalsMC->At(nSignal)->GetName());
1287   className3.Form("Track_%s_%s_MCtruth",fgkPairClassNames[1],fSignalsMC->At(nSignal)->GetName());
1288   Bool_t pairClass=fHistos->GetHistogramList()->FindObject(className.Data())!=0x0;
1289   Bool_t legClass=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
1290   Bool_t trkClass=fHistos->GetHistogramList()->FindObject(className3.Data())!=0x0;
1291   //  printf("fill signal %d: pair %d legs %d trk %d \n",nSignal,pairClass,legClass,trkClass);
1292   if(!pairClass && !legClass && !trkClass) return;
1293
1294   //  printf("leg labels: %d-%d \n",label1,label2);
1295   AliVParticle* part1 = AliDielectronMC::Instance()->GetMCTrackFromMCEvent(label1);
1296   AliVParticle* part2 = AliDielectronMC::Instance()->GetMCTrackFromMCEvent(label2);
1297   if(!part1 && !part2) return;
1298   if(part1&&part2) {
1299     // fill only unlike sign (and only SE)
1300     if(part1->Charge()*part2->Charge()>=0) return;
1301   }
1302
1303
1304   AliDielectronMC* dieMC = AliDielectronMC::Instance();
1305
1306   Int_t mLabel1 = dieMC->GetMothersLabel(label1);    // should work for both ESD and AOD
1307   Int_t mLabel2 = dieMC->GetMothersLabel(label2);
1308
1309   // check the same mother option
1310   AliDielectronSignalMC* sigMC = (AliDielectronSignalMC*)fSignalsMC->At(nSignal);
1311   if(sigMC->GetMothersRelation()==AliDielectronSignalMC::kSame && mLabel1!=mLabel2) return;
1312   if(sigMC->GetMothersRelation()==AliDielectronSignalMC::kDifferent && mLabel1==mLabel2) return;
1313
1314   // fill event values
1315   Double_t values[AliDielectronVarManager::kNMaxValues];
1316   AliDielectronVarManager::SetFillMap(fUsedVars);
1317   AliDielectronVarManager::Fill(dieMC->GetMCEvent(), values); // get event informations
1318
1319   // fill the leg variables
1320   //  printf("leg:%d trk:%d part1:%p part2:%p \n",legClass,trkClass,part1,part2);
1321   if (legClass || trkClass) {
1322     if(part1) AliDielectronVarManager::Fill(part1,values);
1323     if(part1 && trkClass)          fHistos->FillClass(className3, AliDielectronVarManager::kNMaxValues, values);
1324     if(part1 && part2 && legClass) fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
1325     if(part2) AliDielectronVarManager::Fill(part2,values);
1326     if(part2 && trkClass)          fHistos->FillClass(className3, AliDielectronVarManager::kNMaxValues, values);
1327     if(part1 && part2 && legClass) fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
1328   }
1329
1330   //fill pair information
1331   if (pairClass && part1 && part2) {
1332     AliDielectronVarManager::FillVarMCParticle2(part1,part2,values);
1333     fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
1334   }
1335
1336 }
1337
1338 //________________________________________________________________
1339 void AliDielectron::FillMCHistograms(const AliVEvent *ev) {
1340   //
1341   // fill QA MC histograms for pairs and legs of all added mc signals
1342   //
1343
1344   if (!fSignalsMC) return;
1345   TString className,className2,className3;
1346   Double_t values[AliDielectronVarManager::kNMaxValues]={0.};
1347   AliDielectronVarManager::SetFillMap(fUsedVars);
1348   AliDielectronVarManager::Fill(ev, values); // get event informations
1349   //loop over all added mc signals
1350   for(Int_t isig=0; isig<fSignalsMC->GetEntries(); isig++) {
1351
1352     //check if and what to fill
1353     className.Form("Pair_%s",fSignalsMC->At(isig)->GetName());
1354     className2.Form("Track_Legs_%s",fSignalsMC->At(isig)->GetName());
1355     className3.Form("Track_%s_%s",fgkPairClassNames[1],fSignalsMC->At(isig)->GetName());  // unlike sign, SE only
1356     Bool_t pairClass=fHistos->GetHistogramList()->FindObject(className.Data())!=0x0;
1357     Bool_t legClass=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
1358     Bool_t mergedtrkClass=fHistos->GetHistogramList()->FindObject(className3.Data())!=0x0;
1359     if(!pairClass && !legClass && !mergedtrkClass) continue;
1360
1361     // fill pair and/or their leg variables
1362     if(pairClass || legClass) {
1363       Int_t npairs=PairArray(AliDielectron::kEv1PM)->GetEntriesFast(); // only SE +-
1364       for (Int_t ipair=0; ipair<npairs; ++ipair){
1365         AliDielectronPair *pair=static_cast<AliDielectronPair*>(PairArray(AliDielectron::kEv1PM)->UncheckedAt(ipair));
1366
1367         Bool_t isMCtruth = AliDielectronMC::Instance()->IsMCTruth(pair, (AliDielectronSignalMC*)fSignalsMC->At(isig));
1368         if(isMCtruth) {
1369           //fill pair information
1370           if (pairClass){
1371             AliDielectronVarManager::Fill(pair, values);
1372             fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
1373           }
1374           //fill leg information, both + and - in the same histo
1375           if (legClass){
1376             AliDielectronVarManager::Fill(pair->GetFirstDaughter(),values);
1377             fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
1378             AliDielectronVarManager::Fill(pair->GetSecondDaughter(),values);
1379             fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
1380           }
1381         } //is signal
1382       } //loop: pairs
1383     }
1384
1385     // fill single tracks of signals
1386     if(!mergedtrkClass) continue;
1387     // loop over SE track arrays
1388     for (Int_t i=0; i<2; ++i){
1389       Int_t ntracks=fTracks[i].GetEntriesFast();
1390       for (Int_t itrack=0; itrack<ntracks; ++itrack){
1391         Int_t label=((AliVParticle*)fTracks[i].UncheckedAt(itrack))->GetLabel();
1392         Bool_t isMCtruth1 = AliDielectronMC::Instance()->IsMCTruth(label, (AliDielectronSignalMC*)fSignalsMC->At(isig), 1);
1393         Bool_t isMCtruth2 = AliDielectronMC::Instance()->IsMCTruth(label, (AliDielectronSignalMC*)fSignalsMC->At(isig), 2);
1394         // skip if track does not correspond to the signal
1395         if(!isMCtruth1 && !isMCtruth2) continue;
1396         AliDielectronVarManager::Fill(fTracks[i].UncheckedAt(itrack), values);
1397         fHistos->FillClass(className3, AliDielectronVarManager::kNMaxValues, values);
1398       } //loop: tracks
1399     } //loop: arrays
1400
1401   } //loop: MCsignals
1402
1403 }
1404
1405 //______________________________________________
1406 void AliDielectron::SetCentroidCorrFunction(TF1 *fun, UInt_t varx, UInt_t vary, UInt_t varz)
1407 {
1408   fun->GetHistogram()->GetXaxis()->SetUniqueID(varx);
1409   fun->GetHistogram()->GetYaxis()->SetUniqueID(vary);
1410   fun->GetHistogram()->GetZaxis()->SetUniqueID(varz);
1411   // clone temporare histogram since otherwise it will not be streamed to file!
1412   TString key = Form("cntrd%d%d%d",varx,vary,varz);
1413   fPostPIDCntrdCorr = (TH1*)fun->GetHistogram()->Clone(key.Data());
1414   fPostPIDCntrdCorr->GetListOfFunctions()->AddAt(fun,0);
1415 }
1416 //______________________________________________
1417 void AliDielectron::SetWidthCorrFunction(TF1 *fun, UInt_t varx, UInt_t vary, UInt_t varz)
1418 {
1419   fun->GetHistogram()->GetXaxis()->SetUniqueID(varx);
1420   fun->GetHistogram()->GetYaxis()->SetUniqueID(vary);
1421   fun->GetHistogram()->GetZaxis()->SetUniqueID(varz);
1422   // clone temporare histogram since otherwise it will not be streamed to file!
1423   TString key = Form("wdth%d%d%d",varx,vary,varz);
1424   fPostPIDWdthCorr = (TH1*)fun->GetHistogram()->Clone(key.Data());
1425   fPostPIDWdthCorr->GetListOfFunctions()->AddAt(fun,0);
1426 }
1427
1428 //______________________________________________
1429 THnBase* AliDielectron::InitEffMap(TString filename)
1430 {
1431   // init an efficiency object for on-the-fly correction calculations
1432   if(filename.Contains("alien://") && !gGrid) TGrid::Connect("alien://",0,0,"t");
1433
1434   TFile* file=TFile::Open(filename.Data());
1435   if(!file) return 0x0;
1436   THnBase *hGen = (THnBase*) file->Get("hGenerated");
1437   THnBase *hFnd = (THnBase*) file->Get("hFound");
1438   if(!hFnd || !hGen) return 0x0;
1439
1440   hFnd->Divide(hGen);
1441   printf("[I] AliDielectron::InitEffMap efficiency maps %s with %d dimensions loaded! \n",filename.Data(),hFnd->GetNdimensions());
1442   return ((THnBase*) hFnd->Clone("effMap"));
1443 }