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