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