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 HF, CF containers and histograms if they exist
 362   if(fCfManagerPair) fCfManagerPair->SetPdgMother(fPdgMother);
 363
 364   // signals to be studied
 365   Int_t nSignals = fSignalsMC->GetEntries();
 366
 367   Bool_t bFillCF   = (fCfManagerPair ? fCfManagerPair->GetStepForMCtruth()  : kFALSE);
 368   Bool_t bFillHF   = (fHistoArray    ? fHistoArray->GetStepForMCGenerated() : kFALSE);
 369   Bool_t bFillHist = kFALSE;
 370   for(Int_t isig=0;isig<nSignals;isig++) {
 371     TString sigName = fSignalsMC->At(isig)->GetName();
 372     const THashList *histlist =  fHistos->GetHistogramList();
 373     bFillHist |= histlist->FindObject(Form("Pair_%s_MCtruth",sigName.Data()))!=0x0;
 374     bFillHist |= histlist->FindObject(Form("Track_Leg_%s_MCtruth",sigName.Data()))!=0x0;
 375     bFillHist |= histlist->FindObject(Form("Track_%s_%s_MCtruth",fgkPairClassNames[1],sigName.Data()))!=0x0;
 376     if(bFillHist) break;
 377   }
 378   // check if there is anything to fill
 379   if(!bFillCF && !bFillHF && !bFillHist) return;
 380
 381
 382   // initialize 2D arrays of labels for particles from each MC signal
 383   Int_t** labels1;      // labels for particles satisfying branch 1
 384   Int_t** labels2;      // labels for particles satisfying branch 2
 385   Int_t** labels12;     // labels for particles satisfying both branches
 386   labels1 = new Int_t*[nSignals];
 387   labels2 = new Int_t*[nSignals];
 388   labels12 = new Int_t*[nSignals];
 389   Int_t* indexes1=new Int_t[nSignals];
 390   Int_t* indexes2=new Int_t[nSignals];
 391   Int_t* indexes12=new Int_t[nSignals];
 392   for(Int_t isig=0;isig<nSignals;++isig) {
 393     *(labels1+isig) = new Int_t[dieMC->GetNMCTracks()];
 394     *(labels2+isig) = new Int_t[dieMC->GetNMCTracks()];
 395     *(labels12+isig) = new Int_t[dieMC->GetNMCTracks()];
 396     for(Int_t ip=0; ip<dieMC->GetNMCTracks();++ip) {
 397       labels1[isig][ip] = -1;
 398       labels2[isig][ip] = -1;
 399       labels12[isig][ip] = -1;
 400     }
 401     indexes1[isig]=0;
 402     indexes2[isig]=0;
 403     indexes12[isig]=0;
 404   }
 405
 406   Bool_t truth1=kFALSE;
 407   Bool_t truth2=kFALSE;
 408   // loop over the MC tracks
 409   for(Int_t ipart=0; ipart<dieMC->GetNMCTracks(); ++ipart) {
 410     for(Int_t isig=0; isig<nSignals; ++isig) {       // loop over signals
 411       // Proceed only if this signal is required in the pure MC step
 412       // NOTE: Some signals can be satisfied by many particles and this leads to high
 413       //       computation times (e.g. secondary electrons from the GEANT transport). Be aware of this!!
 414       if(!((AliDielectronSignalMC*)fSignalsMC->At(isig))->GetFillPureMCStep()) continue;
 415
 416       truth1 = dieMC->IsMCTruth(ipart, (AliDielectronSignalMC*)fSignalsMC->At(isig), 1);
 417       truth2 = dieMC->IsMCTruth(ipart, (AliDielectronSignalMC*)fSignalsMC->At(isig), 2);
 418
 419       // particles satisfying both branches are treated separately to avoid double counting during pairing
 420       if(truth1 && truth2) {
 421         labels12[isig][indexes12[isig]] = ipart;
 422         ++indexes12[isig];
 423       }
 424       else {
 425         if(truth1) {
 426           labels1[isig][indexes1[isig]] = ipart;
 427           ++indexes1[isig];
 428         }
 429         if(truth2) {
 430           labels2[isig][indexes2[isig]] = ipart;
 431           ++indexes2[isig];
 432         }
 433       }
 434     }
 435   }  // end loop over MC particles
 436
 437   // Do the pairing and fill the CF container with pure MC info
 438   for(Int_t isig=0; isig<nSignals; ++isig) {
 439     //    printf("INDEXES: %d-%d both%d\n",indexes1[isig],indexes2[isig],indexes12[isig]);
 440     // mix the particles which satisfy only one of the signal branches
 441     for(Int_t i1=0;i1<indexes1[isig];++i1) {
 442       if(!indexes2[isig]) FillMCHistograms(labels1[isig][i1], -1, isig); // (e.g. single electrons only, no pairs)
 443       for(Int_t i2=0;i2<indexes2[isig];++i2) {
 444         if(bFillCF) fCfManagerPair->FillMC(labels1[isig][i1], labels2[isig][i2], isig);
 445         if(bFillHF) fHistoArray->Fill(labels1[isig][i1], labels2[isig][i2], isig);
 446         FillMCHistograms(labels1[isig][i1], labels2[isig][i2], isig);
 447       }
 448     }
 449     // mix the particles which satisfy both branches
 450     for(Int_t i1=0;i1<indexes12[isig];++i1) {
 451       for(Int_t i2=0; i2<i1; ++i2) {
 452         if(bFillCF) fCfManagerPair->FillMC(labels12[isig][i1], labels12[isig][i2], isig);
 453         if(bFillHF) fHistoArray->Fill(labels12[isig][i1], labels12[isig][i2], isig);
 454         FillMCHistograms(labels12[isig][i1], labels12[isig][i2], isig);
 455       }
 456     }
 457   }    // end loop over signals
 458
 459   // release the memory
 460   for(Int_t isig=0;isig<nSignals;++isig) {
 461     delete [] *(labels1+isig);
 462     delete [] *(labels2+isig);
 463     delete [] *(labels12+isig);
 464   }
 465   delete [] labels1;
 466   delete [] labels2;
 467   delete [] labels12;
 468   delete [] indexes1;
 469   delete [] indexes2;
 470   delete [] indexes12;
 471 }
 472
 473 //________________________________________________________________
 474 void AliDielectron::FillHistogramsTracks(TObjArray **tracks)
 475 {
 476   //
 477   // Fill Histogram information for tracks after prefilter
 478   // ignore mixed events - for prefilter, only single tracks +/- are relevant
 479   //
 480
 481   TString  className,className2;
 482   Double_t values[AliDielectronVarManager::kNMaxValues];
 483
 484   //Fill track information, separately for the track array candidates
 485   for (Int_t i=0; i<2; ++i){
 486     className.Form("Pre_%s",fgkTrackClassNames[i]);
 487     if (!fHistos->GetHistogramList()->FindObject(className.Data())) continue;
 488     Int_t ntracks=tracks[i]->GetEntriesFast();
 489     for (Int_t itrack=0; itrack<ntracks; ++itrack){
 490       AliDielectronVarManager::Fill(tracks[i]->UncheckedAt(itrack), values);
 491       fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
 492     }
 493   }
 494 }
 495
 496
 497 //________________________________________________________________
 498 void AliDielectron::FillHistogramsMC(const AliMCEvent *ev, AliVEvent *ev1)
 499 {
 500   //
 501   // Fill Histogram information for MCEvents
 502   //
 503
 504   Double_t values[AliDielectronVarManager::kNMaxValues]={0.};
 505   // Fill event information
 506   AliDielectronVarManager::Fill(ev1, values);    // ESD/AOD information
 507   AliDielectronVarManager::Fill(ev, values);     // MC truth info
 508   if (fHistos->GetHistogramList()->FindObject("MCEvent"))
 509     fHistos->FillClass("MCEvent", AliDielectronVarManager::kNMaxValues, values);
 510 }
 511
 512
 513 //________________________________________________________________
 514 void AliDielectron::FillHistograms(const AliVEvent *ev, Bool_t pairInfoOnly)
 515 {
 516   //
 517   // Fill Histogram information for tracks and pairs
 518   //
 519
 520   TString  className,className2;
 521   Double_t values[AliDielectronVarManager::kNMaxValues]={0.};
 522
 523   //Fill event information
 524   if (ev){
 525     if (fHistos->GetHistogramList()->FindObject("Event")) {
 526       fHistos->FillClass("Event", AliDielectronVarManager::kNMaxValues, AliDielectronVarManager::GetData());
 527     }
 528   }
 529
 530   //Fill track information, separately for the track array candidates
 531   if (!pairInfoOnly){
 532     className2.Form("Track_%s",fgkPairClassNames[1]);  // unlike sign, SE only
 533     for (Int_t i=0; i<4; ++i){
 534       className.Form("Track_%s",fgkTrackClassNames[i]);
 535       Bool_t mergedtrkClass=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
 536       Bool_t trkClass=fHistos->GetHistogramList()->FindObject(className.Data())!=0x0;
 537       if (!trkClass && !mergedtrkClass) continue;
 538       Int_t ntracks=fTracks[i].GetEntriesFast();
 539       for (Int_t itrack=0; itrack<ntracks; ++itrack){
 540         AliDielectronVarManager::Fill(fTracks[i].UncheckedAt(itrack), values);
 541         if(trkClass)
 542           fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
 543         if(mergedtrkClass && i<2)
 544           fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values); //only ev1
 545       }
 546     }
 547   }
 548
 549   //Fill Pair information, separately for all pair candidate arrays and the legs
 550   TObjArray arrLegs(100);
 551   for (Int_t i=0; i<10; ++i){
 552     className.Form("Pair_%s",fgkPairClassNames[i]);
 553     className2.Form("Track_Legs_%s",fgkPairClassNames[i]);
 554     Bool_t pairClass=fHistos->GetHistogramList()->FindObject(className.Data())!=0x0;
 555     Bool_t legClass=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
 556     if (!pairClass&&!legClass) continue;
 557     Int_t ntracks=PairArray(i)->GetEntriesFast();
 558     for (Int_t ipair=0; ipair<ntracks; ++ipair){
 559       AliDielectronPair *pair=static_cast<AliDielectronPair*>(PairArray(i)->UncheckedAt(ipair));
 560
 561       //fill pair information
 562       if (pairClass){
 563         AliDielectronVarManager::Fill(pair, values);
 564         fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
 565       }
 566
 567       //fill leg information, don't fill the information twice
 568       if (legClass){
 569         AliVParticle *d1=pair->GetFirstDaughter();
 570         AliVParticle *d2=pair->GetSecondDaughter();
 571         if (!arrLegs.FindObject(d1)){
 572           AliDielectronVarManager::Fill(d1, values);
 573           fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
 574           arrLegs.Add(d1);
 575         }
 576         if (!arrLegs.FindObject(d2)){
 577           AliDielectronVarManager::Fill(d2, values);
 578           fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
 579           arrLegs.Add(d2);
 580         }
 581       }
 582     }
 583     if (legClass) arrLegs.Clear();
 584   }
 585
 586 }
 587 //________________________________________________________________
 588 void AliDielectron::FillHistogramsPair(AliDielectronPair *pair,Bool_t fromPreFilter/*=kFALSE*/)
 589 {
 590   //
 591   // Fill Histogram information for pairs and the track in the pair
 592   // NOTE: in this funtion the leg information may be filled multiple
 593   //       times. This funtion is used in the track rotation pairing
 594   //       and those legs are not saved!
 595   //
 596   TString  className,className2;
 597   Double_t values[AliDielectronVarManager::kNMaxValues];
 598
 599   //Fill Pair information, separately for all pair candidate arrays and the legs
 600   TObjArray arrLegs(100);
 601   const Int_t type=pair->GetType();
 602   if (fromPreFilter) {
 603     className.Form("RejPair_%s",fgkPairClassNames[type]);
 604     className2.Form("RejTrack_%s",fgkPairClassNames[type]);
 605   } else {
 606     className.Form("Pair_%s",fgkPairClassNames[type]);
 607     className2.Form("Track_Legs_%s",fgkPairClassNames[type]);
 608   }
 609
 610   Bool_t pairClass=fHistos->GetHistogramList()->FindObject(className.Data())!=0x0;
 611   Bool_t legClass=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
 612
 613   //fill pair information
 614   if (pairClass){
 615     AliDielectronVarManager::Fill(pair, values);
 616     fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
 617   }
 618
 619   if (legClass){
 620     AliVParticle *d1=pair->GetFirstDaughter();
 621     AliDielectronVarManager::Fill(d1, values);
 622     fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
 623
 624     AliVParticle *d2=pair->GetSecondDaughter();
 625     AliDielectronVarManager::Fill(d2, values);
 626     fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
 627   }
 628 }
 629
 630 //________________________________________________________________
 631 void AliDielectron::FillTrackArrays(AliVEvent * const ev, Int_t eventNr)
 632 {
 633   //
 634   // select tracks and fill track candidate arrays
 635   // eventNr = 0: First  event, use track arrays 0 and 1
 636   // eventNr = 1: Second event, use track arrays 2 and 3
 637   //
 638
 639   Int_t ntracks=ev->GetNumberOfTracks();
 640
 641   UInt_t selectedMask=(1<<fTrackFilter.GetCuts()->GetEntries())-1;
 642   for (Int_t itrack=0; itrack<ntracks; ++itrack){
 643     //get particle
 644     AliVParticle *particle=ev->GetTrack(itrack);
 645
 646     //apply track cuts
 647     UInt_t cutmask=fTrackFilter.IsSelected(particle);
 648     //fill cut QA
 649     if(fCutQA) fQAmonitor->FillAll(particle);
 650     if(fCutQA) fQAmonitor->Fill(cutmask,particle);
 651
 652     if (cutmask!=selectedMask) continue;
 653
 654     //fill selected particle into the corresponding track arrays
 655     Short_t charge=particle->Charge();
 656     if (charge>0)      fTracks[eventNr*2].Add(particle);
 657     else if (charge<0) fTracks[eventNr*2+1].Add(particle);
 658
 659
 660   }
 661 }
 662
 663 //________________________________________________________________
 664 void AliDielectron::EventPlanePreFilter(Int_t arr1, Int_t arr2, TObjArray arrTracks1, TObjArray arrTracks2, const AliVEvent *ev)
 665 {
 666   //
 667   // Prefilter tracks and tracks from pairs
 668   // Needed for rejection in the Q-Vector of the event plane
 669   // remove contribution of all tracks to the Q-vector that are in invariant mass window
 670   //
 671
 672   AliEventplane *evplane = const_cast<AliVEvent *>(ev)->GetEventplane();
 673   if(!evplane) { // nanoAODs , here we do NOT have sub event reaction planes
 674     //  if(1) {
 675     // get the EPselectionTask for recalculation of weighting factors
 676     AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 677     AliEPSelectionTask *eptask = dynamic_cast<AliEPSelectionTask *>(man->GetTask("EventplaneSelection"));
 678     if(!eptask) return;
 679
 680     // track mapping
 681     TMap mapRemovedTracks;
 682
 683     // eta gap ?
 684     Bool_t etagap = kFALSE;
 685     for (Int_t iCut=0; iCut<fEventPlanePreFilter.GetCuts()->GetEntries();++iCut) {
 686       TString cutName=fEventPlanePreFilter.GetCuts()->At(iCut)->GetName();
 687       if(cutName.Contains("eta") || cutName.Contains("Eta"))  etagap=kTRUE;
 688     }
 689
 690     Double_t cQX=0., cQY=0.;
 691     // apply cuts to the tracks, e.g. etagap
 692     if(fEventPlanePreFilter.GetCuts()->GetEntries()) {
 693       UInt_t selectedMask=(1<<fEventPlanePreFilter.GetCuts()->GetEntries())-1;
 694       Int_t ntracks=ev->GetNumberOfTracks();
 695       for (Int_t itrack=0; itrack<ntracks; ++itrack){
 696         AliVParticle *particle=ev->GetTrack(itrack);
 697         AliVTrack *track= static_cast<AliVTrack*>(particle);
 698         if (!track) continue;
 699         //event plane cuts
 700         UInt_t cutMask=fEventPlanePreFilter.IsSelected(track);
 701         //apply cut
 702         if (cutMask==selectedMask) continue;
 703
 704         mapRemovedTracks.Add(track,track);
 705         cQX += (eptask->GetWeight(track) * TMath::Cos(2*track->Phi()));
 706         cQY += (eptask->GetWeight(track) * TMath::Sin(2*track->Phi()));
 707       }
 708     }
 709
 710     // POI (particle of interest) rejection
 711     Int_t pairIndex=GetPairIndex(arr1,arr2);
 712
 713     Int_t ntrack1=arrTracks1.GetEntriesFast();
 714     Int_t ntrack2=arrTracks2.GetEntriesFast();
 715     AliDielectronPair candidate;
 716     candidate.SetKFUsage(fUseKF);
 717
 718     UInt_t selectedMask=(1<<fEventPlanePOIPreFilter.GetCuts()->GetEntries())-1;
 719     for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1){
 720       Int_t end=ntrack2;
 721       if (arr1==arr2) end=itrack1;
 722       Bool_t accepted=kFALSE;
 723       for (Int_t itrack2=0; itrack2<end; ++itrack2){
 724         TObject *track1=arrTracks1.UncheckedAt(itrack1);
 725         TObject *track2=arrTracks2.UncheckedAt(itrack2);
 726         if (!track1 || !track2) continue;
 727         //create the pair
 728         candidate.SetTracks(static_cast<AliVTrack*>(track1), fPdgLeg1,
 729                             static_cast<AliVTrack*>(track2), fPdgLeg2);
 730         candidate.SetType(pairIndex);
 731         candidate.SetLabel(AliDielectronMC::Instance()->GetLabelMotherWithPdg(&candidate,fPdgMother));
 732
 733         //event plane pair cuts
 734         UInt_t cutMask=fEventPlanePOIPreFilter.IsSelected(&candidate);
 735         //apply cut
 736         if (cutMask==selectedMask) continue;
 737
 738         accepted=kTRUE;
 739         //remove the tracks from the Track arrays
 740         arrTracks2.AddAt(0x0,itrack2);
 741       }
 742       if ( accepted ) arrTracks1.AddAt(0x0,itrack1);
 743     }
 744     //compress the track arrays
 745     arrTracks1.Compress();
 746     arrTracks2.Compress();
 747
 748     //Modify the components: subtract the tracks
 749     ntrack1=arrTracks1.GetEntriesFast();
 750     ntrack2=arrTracks2.GetEntriesFast();
 751     // remove leg1 contribution
 752     for (Int_t itrack=0; itrack<ntrack1; ++itrack){
 753       AliVTrack *track= static_cast<AliVTrack*>(arrTracks1.UncheckedAt(itrack));
 754       if (!track) continue;
 755       // track contribution was already removed
 756       if (mapRemovedTracks.FindObject(track)) continue;
 757       else mapRemovedTracks.Add(track,track);
 758
 759       cQX += (eptask->GetWeight(track) * TMath::Cos(2*track->Phi()));
 760       cQY += (eptask->GetWeight(track) * TMath::Sin(2*track->Phi()));
 761     }
 762     // remove leg2 contribution
 763     for (Int_t itrack=0; itrack<ntrack2; ++itrack){
 764       AliVTrack *track= static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack));
 765       if (!track) continue;
 766       // track contribution was already removed
 767       if (mapRemovedTracks.FindObject(track)) continue;
 768       else mapRemovedTracks.Add(track,track);
 769
 770       cQX += (eptask->GetWeight(track) * TMath::Cos(2*track->Phi()));
 771       cQY += (eptask->GetWeight(track) * TMath::Sin(2*track->Phi()));
 772     }
 773
 774     // build a corrected alieventplane using the values from the var manager
 775     // these uncorrected values are filled using the stored magnitude and angle  in the header
 776     TVector2 qcorr;
 777     qcorr.Set(AliDielectronVarManager::GetValue(AliDielectronVarManager::kTPCxH2uc)-cQX,
 778               AliDielectronVarManager::GetValue(AliDielectronVarManager::kTPCyH2uc)-cQY);
 779     // fill alieventplane
 780     AliEventplane cevplane;
 781     cevplane.SetQVector(&qcorr);
 782     AliDielectronVarManager::SetTPCEventPlane(&cevplane);
 783     cevplane.SetQVector(0);
 784     return;
 785   } //end: nanoAODs
 786   else
 787     {
 788     // this is done in case of ESDs or AODs
 789     Bool_t isESD=(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->IsA()==AliESDInputHandler::Class());
 790     // copy event plane object
 791     AliEventplane cevplane(*evplane);
 792     //    Int_t nMaxID=cevplane->GetQContributionXArray()->GetSize();
 793
 794     TVector2 *qcorr  = cevplane.GetQVector();
 795     if(!qcorr) return;
 796     TVector2 *qcsub1 = 0x0;
 797     TVector2 *qcsub2 = 0x0;
 798
 799     // eta gap ?
 800     Bool_t etagap = kFALSE;
 801     for (Int_t iCut=0; iCut<fEventPlanePreFilter.GetCuts()->GetEntries();++iCut) {
 802       TString cutName=fEventPlanePreFilter.GetCuts()->At(iCut)->GetName();
 803       if(cutName.Contains("eta") || cutName.Contains("Eta"))  etagap=kTRUE;
 804     }
 805
 806     // subevent configuration for eta gap or LS (default is rndm)
 807     if(fLikeSignSubEvents && etagap) {
 808       // start with the full Qvector/event in both sub events
 809       qcsub1 = new TVector2(*qcorr);
 810       qcsub2 = new TVector2(*qcorr);
 811       cevplane.SetQsub(qcsub1,qcsub2);
 812
 813       Int_t ntracks=ev->GetNumberOfTracks();
 814       // track removals
 815       for (Int_t itrack=0; itrack<ntracks; ++itrack){
 816         AliVParticle *particle=ev->GetTrack(itrack);
 817         AliVTrack *track= static_cast<AliVTrack*>(particle);
 818         if (!track) continue;
 819         if (track->GetID()>=0 && !isESD) continue;
 820         Int_t tmpID = isESD ? track->GetID() : track->GetID()*-1 - 1;
 821
 822         // set contributions to zero
 823         // charge sub1+ sub2-
 824         if(fLikeSignSubEvents) {
 825           Short_t charge=track->Charge();
 826           if (charge<0) {
 827             cevplane.GetQContributionXArraysub1()->SetAt(0.0, tmpID);
 828             cevplane.GetQContributionYArraysub1()->SetAt(0.0, tmpID);
 829           }
 830           if (charge>0) {
 831             cevplane.GetQContributionXArraysub2()->SetAt(0.0, tmpID);
 832             cevplane.GetQContributionYArraysub2()->SetAt(0.0, tmpID);
 833           }
 834         }
 835         // eta sub1+ sub2-
 836         if(etagap) {
 837           Double_t eta=track->Eta();
 838           if (eta<0.0) {
 839             cevplane.GetQContributionXArraysub1()->SetAt(0.0, tmpID);
 840             cevplane.GetQContributionYArraysub1()->SetAt(0.0, tmpID);
 841           }
 842           if (eta>0.0) {
 843             cevplane.GetQContributionXArraysub2()->SetAt(0.0, tmpID);
 844             cevplane.GetQContributionYArraysub2()->SetAt(0.0, tmpID);
 845           }
 846         }
 847       } // end: loop over tracks
 848     } // end: sub event configuration
 849
 850     // apply cuts, e.g. etagap
 851     if(fEventPlanePreFilter.GetCuts()->GetEntries()) {
 852       UInt_t selectedMask=(1<<fEventPlanePreFilter.GetCuts()->GetEntries())-1;
 853       Int_t ntracks=ev->GetNumberOfTracks();
 854       for (Int_t itrack=0; itrack<ntracks; ++itrack){
 855         AliVParticle *particle=ev->GetTrack(itrack);
 856         AliVTrack *track= static_cast<AliVTrack*>(particle);
 857         if (!track) continue;
 858         if (track->GetID()>=0 && !isESD) continue;
 859         Int_t tmpID = isESD ? track->GetID() : track->GetID()*-1 - 1;
 860
 861         //event plane cuts
 862         UInt_t cutMask=fEventPlanePreFilter.IsSelected(track);
 863         //apply cut
 864         if (cutMask==selectedMask) continue;
 865
 866         // set contributions to zero
 867         cevplane.GetQContributionXArray()->SetAt(0.0, tmpID);
 868         cevplane.GetQContributionYArray()->SetAt(0.0, tmpID);
 869         cevplane.GetQContributionXArraysub1()->SetAt(0.0, tmpID);
 870         cevplane.GetQContributionYArraysub1()->SetAt(0.0, tmpID);
 871         cevplane.GetQContributionXArraysub2()->SetAt(0.0, tmpID);
 872         cevplane.GetQContributionYArraysub2()->SetAt(0.0, tmpID);
 873       }
 874     } // end: track cuts
 875
 876     // POI (particle of interest) rejection
 877     Int_t pairIndex=GetPairIndex(arr1,arr2);
 878     Int_t ntrack1=arrTracks1.GetEntriesFast();
 879     Int_t ntrack2=arrTracks2.GetEntriesFast();
 880     AliDielectronPair candidate;
 881     candidate.SetKFUsage(fUseKF);
 882
 883     UInt_t selectedMask=(1<<fEventPlanePOIPreFilter.GetCuts()->GetEntries())-1;
 884     for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1){
 885       Int_t end=ntrack2;
 886       if (arr1==arr2) end=itrack1;
 887       Bool_t accepted=kFALSE;
 888       for (Int_t itrack2=0; itrack2<end; ++itrack2){
 889         TObject *track1=arrTracks1.UncheckedAt(itrack1);
 890         TObject *track2=arrTracks2.UncheckedAt(itrack2);
 891         if (!track1 || !track2) continue;
 892         //create the pair
 893         candidate.SetTracks(static_cast<AliVTrack*>(track1), fPdgLeg1,
 894                             static_cast<AliVTrack*>(track2), fPdgLeg2);
 895
 896         candidate.SetType(pairIndex);
 897         candidate.SetLabel(AliDielectronMC::Instance()->GetLabelMotherWithPdg(&candidate,fPdgMother));
 898
 899         //event plane cuts
 900         UInt_t cutMask=fEventPlanePOIPreFilter.IsSelected(&candidate);
 901         //apply cut
 902         if (cutMask==selectedMask) continue;
 903
 904         accepted=kTRUE;
 905         //remove the tracks from the Track arrays
 906         arrTracks2.AddAt(0x0,itrack2);
 907       }
 908       if ( accepted ) arrTracks1.AddAt(0x0,itrack1);
 909     }
 910     //compress the track arrays
 911     arrTracks1.Compress();
 912     arrTracks2.Compress();
 913
 914     //Modify the components: subtract the tracks
 915     ntrack1=arrTracks1.GetEntriesFast();
 916     ntrack2=arrTracks2.GetEntriesFast();
 917     // remove leg1 contribution
 918     for (Int_t itrack=0; itrack<ntrack1; ++itrack){
 919       AliVTrack *track= static_cast<AliVTrack*>(arrTracks1.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     // remove leg2 contribution
 932     for (Int_t itrack=0; itrack<ntrack2; ++itrack){
 933       AliVTrack *track= static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack));
 934       if (!track) continue;
 935       if (track->GetID()>=0 && !isESD) continue;
 936       Int_t tmpID = isESD ? track->GetID() : track->GetID()*-1 - 1;
 937       // set contributions to zero
 938       cevplane.GetQContributionXArray()->SetAt(0.0, tmpID);
 939       cevplane.GetQContributionYArray()->SetAt(0.0, tmpID);
 940       cevplane.GetQContributionXArraysub1()->SetAt(0.0, tmpID);
 941       cevplane.GetQContributionYArraysub1()->SetAt(0.0, tmpID);
 942       cevplane.GetQContributionXArraysub2()->SetAt(0.0, tmpID);
 943       cevplane.GetQContributionYArraysub2()->SetAt(0.0, tmpID);
 944     }
 945
 946     // set corrected AliEventplane and fill variables with corrected values
 947     AliDielectronVarManager::SetTPCEventPlane(&cevplane);
 948     delete qcsub1;
 949     delete qcsub2;
 950   } // end: ESD or AOD case
 951
 952 }
 953 //________________________________________________________________
 954 void AliDielectron::PairPreFilter(Int_t arr1, Int_t arr2, TObjArray &arrTracks1, TObjArray &arrTracks2)
 955 {
 956   //
 957   // Prefilter tracks from pairs
 958   // Needed for datlitz rejections
 959   // remove all tracks from the Single track arrays that pass the cuts in this filter
 960   //
 961
 962   Int_t ntrack1=arrTracks1.GetEntriesFast();
 963   Int_t ntrack2=arrTracks2.GetEntriesFast();
 964   AliDielectronPair candidate;
 965   candidate.SetKFUsage(fUseKF);
 966   // flag arrays for track removal
 967   Bool_t *bTracks1 = new Bool_t[ntrack1];
 968   for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1) bTracks1[itrack1]=kFALSE;
 969   Bool_t *bTracks2 = new Bool_t[ntrack2];
 970   for (Int_t itrack2=0; itrack2<ntrack2; ++itrack2) bTracks2[itrack2]=kFALSE;
 971
 972   UInt_t selectedMask=(1<<fPairPreFilter.GetCuts()->GetEntries())-1;
 973   UInt_t selectedMaskPair=(1<<fPairFilter.GetCuts()->GetEntries())-1;
 974
 975   Int_t nRejPasses = 1; //for fPreFilterUnlikeOnly and no set flag
 976   if (fPreFilterAllSigns) nRejPasses = 3;
 977
 978   for (Int_t iRP=0; iRP < nRejPasses; ++iRP) {
 979         Int_t arr1RP=arr1, arr2RP=arr2;
 980         TObjArray *arrTracks1RP=&arrTracks1;
 981         TObjArray *arrTracks2RP=&arrTracks2;
 982         Bool_t *bTracks1RP = bTracks1;
 983         Bool_t *bTracks2RP = bTracks2;
 984         switch (iRP) {
 985                 case 1: arr1RP=arr1;arr2RP=arr1;
 986                                 arrTracks1RP=&arrTracks1;
 987                                 arrTracks2RP=&arrTracks1;
 988                                 bTracks1RP = bTracks1;
 989                                 bTracks2RP = bTracks1;
 990                                 break;
 991                 case 2: arr1RP=arr2;arr2RP=arr2;
 992                                 arrTracks1RP=&arrTracks2;
 993                                 arrTracks2RP=&arrTracks2;
 994                                 bTracks1RP = bTracks2;
 995                                 bTracks2RP = bTracks2;
 996                                 break;
 997                 default: ;//nothing to do
 998         }
 999         Int_t ntrack1RP=(*arrTracks1RP).GetEntriesFast();
1000         Int_t ntrack2RP=(*arrTracks2RP).GetEntriesFast();
1001
1002         Int_t pairIndex=GetPairIndex(arr1RP,arr2RP);
1003
1004         for (Int_t itrack1=0; itrack1<ntrack1RP; ++itrack1){
1005           Int_t end=ntrack2RP;
1006           if (arr1RP==arr2RP) end=itrack1;
1007           for (Int_t itrack2=0; itrack2<end; ++itrack2){
1008                 TObject *track1=(*arrTracks1RP).UncheckedAt(itrack1);
1009                 TObject *track2=(*arrTracks2RP).UncheckedAt(itrack2);
1010                 if (!track1 || !track2) continue;
1011                 //create the pair
1012                 candidate.SetTracks(static_cast<AliVTrack*>(track1), fPdgLeg1,
1013                         static_cast<AliVTrack*>(track2), fPdgLeg2);
1014
1015                 candidate.SetType(pairIndex);
1016                 candidate.SetLabel(AliDielectronMC::Instance()->GetLabelMotherWithPdg(&candidate,fPdgMother));
1017                 //relate to the production vertex
1018                 //       if (AliDielectronVarManager::GetKFVertex()) candidate.SetProductionVertex(*AliDielectronVarManager::GetKFVertex());
1019
1020                 //pair cuts
1021                 UInt_t cutMask=fPairPreFilter.IsSelected(&candidate);
1022
1023                 //apply cut
1024                 if (cutMask!=selectedMask) continue;
1025                 if (fCfManagerPair) fCfManagerPair->Fill(selectedMaskPair+1 ,&candidate);
1026                 if (fHistos) FillHistogramsPair(&candidate,kTRUE);
1027                 //set flags for track removal
1028                 bTracks1RP[itrack1]=kTRUE;
1029                 bTracks2RP[itrack2]=kTRUE;
1030           }
1031         }
1032   }
1033
1034   //remove the tracks from the Track arrays
1035   for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1){
1036     if(bTracks1[itrack1]) arrTracks1.AddAt(0x0, itrack1);
1037   }
1038   for (Int_t itrack2=0; itrack2<ntrack2; ++itrack2){
1039     if(bTracks2[itrack2]) arrTracks2.AddAt(0x0, itrack2);
1040   }
1041
1042   // clean up
1043   delete [] bTracks1;
1044   delete [] bTracks2;
1045
1046   //compress the track arrays
1047   arrTracks1.Compress();
1048   arrTracks2.Compress();
1049
1050   //apply leg cuts after the pre filter
1051   if ( fPairPreFilterLegs.GetCuts()->GetEntries()>0 ) {
1052     selectedMask=(1<<fPairPreFilterLegs.GetCuts()->GetEntries())-1;
1053     //loop over tracks from array 1
1054     for (Int_t itrack=0; itrack<arrTracks1.GetEntriesFast();++itrack){
1055       //test cuts
1056       UInt_t cutMask=fPairPreFilterLegs.IsSelected(arrTracks1.UncheckedAt(itrack));
1057
1058       //apply cut
1059       if (cutMask!=selectedMask) arrTracks1.AddAt(0x0,itrack);;
1060     }
1061     arrTracks1.Compress();
1062
1063     //in case of like sign don't loop over second array
1064     if (arr1==arr2) {
1065       arrTracks2=arrTracks1;
1066     } else {
1067
1068       //loop over tracks from array 2
1069       for (Int_t itrack=0; itrack<arrTracks2.GetEntriesFast();++itrack){
1070       //test cuts
1071         UInt_t cutMask=fPairPreFilterLegs.IsSelected(arrTracks2.UncheckedAt(itrack));
1072       //apply cut
1073         if (cutMask!=selectedMask) arrTracks2.AddAt(0x0,itrack);
1074       }
1075       arrTracks2.Compress();
1076
1077     }
1078   }
1079   //For unlike-sign monitor track-cuts:
1080   if (arr1!=arr2&&fHistos) {
1081     TObjArray *unlikesignArray[2] = {&arrTracks1,&arrTracks2};
1082     FillHistogramsTracks(unlikesignArray);
1083   }
1084 }
1085
1086 //________________________________________________________________
1087 void AliDielectron::FillPairArrays(Int_t arr1, Int_t arr2)
1088 {
1089   //
1090   // select pairs and fill pair candidate arrays
1091   //
1092
1093   TObjArray arrTracks1=fTracks[arr1];
1094   TObjArray arrTracks2=fTracks[arr2];
1095
1096   //process pre filter if set
1097   if ((!fPreFilterAllSigns) && (!fPreFilterUnlikeOnly) && ( fPairPreFilter.GetCuts()->GetEntries()>0 ))  PairPreFilter(arr1, arr2, arrTracks1, arrTracks2);
1098
1099   Int_t pairIndex=GetPairIndex(arr1,arr2);
1100
1101   Int_t ntrack1=arrTracks1.GetEntriesFast();
1102   Int_t ntrack2=arrTracks2.GetEntriesFast();
1103
1104   AliDielectronPair *candidate=new AliDielectronPair;
1105   candidate->SetKFUsage(fUseKF);
1106
1107   UInt_t selectedMask=(1<<fPairFilter.GetCuts()->GetEntries())-1;
1108
1109   for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1){
1110     Int_t end=ntrack2;
1111     if (arr1==arr2) end=itrack1;
1112     for (Int_t itrack2=0; itrack2<end; ++itrack2){
1113       //create the pair
1114       candidate->SetTracks(static_cast<AliVTrack*>(arrTracks1.UncheckedAt(itrack1)), fPdgLeg1,
1115                              static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack2)), fPdgLeg2);
1116       candidate->SetType(pairIndex);
1117       Int_t label=AliDielectronMC::Instance()->GetLabelMotherWithPdg(candidate,fPdgMother);
1118       candidate->SetLabel(label);
1119       if (label>-1) candidate->SetPdgCode(fPdgMother);
1120
1121       // check for gamma kf particle
1122       label=AliDielectronMC::Instance()->GetLabelMotherWithPdg(candidate,22);
1123       if (label>-1) {
1124         candidate->SetGammaTracks(static_cast<AliVTrack*>(arrTracks1.UncheckedAt(itrack1)), fPdgLeg1,
1125                                   static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack2)), fPdgLeg2);
1126       // should we set the pdgmothercode and the label
1127       }
1128
1129       //pair cuts
1130       UInt_t cutMask=fPairFilter.IsSelected(candidate);
1131
1132       //CF manager for the pair
1133       if (fCfManagerPair) fCfManagerPair->Fill(cutMask,candidate);
1134       //histogram array for the pair
1135       if (fHistoArray) fHistoArray->Fill(pairIndex,candidate);
1136       // cut qa
1137       if(pairIndex==kEv1PM && fCutQA) {
1138         fQAmonitor->FillAll(candidate);
1139         fQAmonitor->Fill(cutMask,candidate);
1140       }
1141
1142       //apply cut
1143       if (cutMask!=selectedMask) continue;
1144
1145       //add the candidate to the candidate array
1146       PairArray(pairIndex)->Add(candidate);
1147       //get a new candidate
1148       candidate=new AliDielectronPair;
1149           candidate->SetKFUsage(fUseKF);
1150     }
1151   }
1152   //delete the surplus candidate
1153   delete candidate;
1154 }
1155
1156 //________________________________________________________________
1157 void AliDielectron::FillPairArrayTR()
1158 {
1159   //
1160   // select pairs and fill pair candidate arrays
1161   //
1162   UInt_t selectedMask=(1<<fPairFilter.GetCuts()->GetEntries())-1;
1163
1164   while ( fTrackRotator->NextCombination() ){
1165     AliDielectronPair candidate;
1166     candidate.SetKFUsage(fUseKF);
1167     candidate.SetTracks(&fTrackRotator->GetKFTrackP(), &fTrackRotator->GetKFTrackN(),
1168                         fTrackRotator->GetVTrackP(),fTrackRotator->GetVTrackN());
1169     candidate.SetType(kEv1PMRot);
1170
1171     //pair cuts
1172     UInt_t cutMask=fPairFilter.IsSelected(&candidate);
1173
1174     //CF manager for the pair
1175     if (fCfManagerPair) fCfManagerPair->Fill(cutMask,&candidate);
1176     //histogram array for the pair
1177     if (fHistoArray) fHistoArray->Fill((Int_t)kEv1PMRot,&candidate);
1178
1179     //apply cut
1180     if (cutMask==selectedMask) {
1181      if(fHistos) FillHistogramsPair(&candidate);
1182      if(fStoreRotatedPairs) PairArray(kEv1PMRot)->Add(new AliDielectronPair(candidate));
1183     }
1184   }
1185 }
1186
1187 //________________________________________________________________
1188 void AliDielectron::FillDebugTree()
1189 {
1190   //
1191   // Fill Histogram information for tracks and pairs
1192   //
1193
1194   //Fill Debug tree
1195   for (Int_t i=0; i<10; ++i){
1196     Int_t ntracks=PairArray(i)->GetEntriesFast();
1197     for (Int_t ipair=0; ipair<ntracks; ++ipair){
1198       fDebugTree->Fill(static_cast<AliDielectronPair*>(PairArray(i)->UncheckedAt(ipair)));
1199     }
1200   }
1201 }
1202
1203 //________________________________________________________________
1204 void AliDielectron::SaveDebugTree()
1205 {
1206   //
1207   // delete the debug tree, this will also write the tree
1208   //
1209   if (fDebugTree) fDebugTree->DeleteStreamer();
1210 }
1211
1212
1213 //__________________________________________________________________
1214 void AliDielectron::AddSignalMC(AliDielectronSignalMC* signal) {
1215   //
1216   //  Add an MC signal to the signals list
1217   //
1218   if(!fSignalsMC) {
1219     fSignalsMC = new TObjArray();
1220     fSignalsMC->SetOwner();
1221   }
1222   fSignalsMC->Add(signal);
1223 }
1224
1225 //________________________________________________________________
1226 void AliDielectron::FillMCHistograms(Int_t label1, Int_t label2, Int_t nSignal) {
1227   //
1228   // fill QA MC TRUTH histograms for pairs and legs of all added mc signals
1229   //
1230
1231   TString className,className2,className3;
1232   className.Form("Pair_%s_MCtruth",fSignalsMC->At(nSignal)->GetName());
1233   className2.Form("Track_Legs_%s_MCtruth",fSignalsMC->At(nSignal)->GetName());
1234   className3.Form("Track_%s_%s_MCtruth",fgkPairClassNames[1],fSignalsMC->At(nSignal)->GetName());
1235   Bool_t pairClass=fHistos->GetHistogramList()->FindObject(className.Data())!=0x0;
1236   Bool_t legClass=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
1237   Bool_t trkClass=fHistos->GetHistogramList()->FindObject(className3.Data())!=0x0;
1238   //  printf("fill signal %d: pair %d legs %d trk %d \n",nSignal,pairClass,legClass,trkClass);
1239   if(!pairClass && !legClass && !trkClass) return;
1240
1241   //  printf("leg labels: %d-%d \n",label1,label2);
1242   AliVParticle* part1 = AliDielectronMC::Instance()->GetMCTrackFromMCEvent(label1);
1243   AliVParticle* part2 = AliDielectronMC::Instance()->GetMCTrackFromMCEvent(label2);
1244   if(!part1 && !part2) return;
1245   if(part1&&part2) {
1246     // fill only unlike sign (and only SE)
1247     if(part1->Charge()*part2->Charge()>=0) return;
1248   }
1249
1250
1251   AliDielectronMC* dieMC = AliDielectronMC::Instance();
1252
1253   Int_t mLabel1 = dieMC->GetMothersLabel(label1);    // should work for both ESD and AOD
1254   Int_t mLabel2 = dieMC->GetMothersLabel(label2);
1255
1256   // check the same mother option
1257   AliDielectronSignalMC* sigMC = (AliDielectronSignalMC*)fSignalsMC->At(nSignal);
1258   if(sigMC->GetMothersRelation()==AliDielectronSignalMC::kSame && mLabel1!=mLabel2) return;
1259   if(sigMC->GetMothersRelation()==AliDielectronSignalMC::kDifferent && mLabel1==mLabel2) return;
1260
1261   // fill event values
1262   Double_t values[AliDielectronVarManager::kNMaxValues];
1263   AliDielectronVarManager::Fill(dieMC->GetMCEvent(), values); // get event informations
1264
1265   // fill the leg variables
1266   //  printf("leg:%d trk:%d part1:%p part2:%p \n",legClass,trkClass,part1,part2);
1267   if (legClass || trkClass) {
1268     if(part1) AliDielectronVarManager::Fill(part1,values);
1269     if(part1 && trkClass)          fHistos->FillClass(className3, AliDielectronVarManager::kNMaxValues, values);
1270     if(part1 && part2 && legClass) fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
1271     if(part2) AliDielectronVarManager::Fill(part2,values);
1272     if(part2 && trkClass)          fHistos->FillClass(className3, AliDielectronVarManager::kNMaxValues, values);
1273     if(part1 && part2 && legClass) fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
1274   }
1275
1276   //fill pair information
1277   if (pairClass && part1 && part2) {
1278     AliDielectronVarManager::FillVarMCParticle2(part1,part2,values);
1279     fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
1280   }
1281
1282 }
1283
1284 //________________________________________________________________
1285 void AliDielectron::FillMCHistograms(const AliVEvent *ev) {
1286   //
1287   // fill QA MC histograms for pairs and legs of all added mc signals
1288   //
1289
1290   if (!fSignalsMC) return;
1291   TString className,className2,className3;
1292   Double_t values[AliDielectronVarManager::kNMaxValues]={0.};
1293   AliDielectronVarManager::Fill(ev, values); // get event informations
1294   //loop over all added mc signals
1295   for(Int_t isig=0; isig<fSignalsMC->GetEntries(); isig++) {
1296
1297     //check if and what to fill
1298     className.Form("Pair_%s",fSignalsMC->At(isig)->GetName());
1299     className2.Form("Track_Legs_%s",fSignalsMC->At(isig)->GetName());
1300     className3.Form("Track_%s_%s",fgkPairClassNames[1],fSignalsMC->At(isig)->GetName());  // unlike sign, SE only
1301     Bool_t pairClass=fHistos->GetHistogramList()->FindObject(className.Data())!=0x0;
1302     Bool_t legClass=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
1303     Bool_t mergedtrkClass=fHistos->GetHistogramList()->FindObject(className3.Data())!=0x0;
1304     if(!pairClass && !legClass && !mergedtrkClass) continue;
1305
1306     // fill pair and/or their leg variables
1307     if(pairClass || legClass) {
1308       Int_t npairs=PairArray(AliDielectron::kEv1PM)->GetEntriesFast(); // only SE +-
1309       for (Int_t ipair=0; ipair<npairs; ++ipair){
1310         AliDielectronPair *pair=static_cast<AliDielectronPair*>(PairArray(AliDielectron::kEv1PM)->UncheckedAt(ipair));
1311
1312         Bool_t isMCtruth = AliDielectronMC::Instance()->IsMCTruth(pair, (AliDielectronSignalMC*)fSignalsMC->At(isig));
1313         if(isMCtruth) {
1314           //fill pair information
1315           if (pairClass){
1316             AliDielectronVarManager::Fill(pair, values);
1317             fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
1318           }
1319           //fill leg information, both + and - in the same histo
1320           if (legClass){
1321             AliDielectronVarManager::Fill(pair->GetFirstDaughter(),values);
1322             fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
1323             AliDielectronVarManager::Fill(pair->GetSecondDaughter(),values);
1324             fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
1325           }
1326         } //is signal
1327       } //loop: pairs
1328     }
1329
1330     // fill single tracks of signals
1331     if(!mergedtrkClass) continue;
1332     // loop over SE track arrays
1333     for (Int_t i=0; i<2; ++i){
1334       Int_t ntracks=fTracks[i].GetEntriesFast();
1335       for (Int_t itrack=0; itrack<ntracks; ++itrack){
1336         Int_t label=((AliVParticle*)fTracks[i].UncheckedAt(itrack))->GetLabel();
1337         Bool_t isMCtruth1 = AliDielectronMC::Instance()->IsMCTruth(label, (AliDielectronSignalMC*)fSignalsMC->At(isig), 1);
1338         Bool_t isMCtruth2 = AliDielectronMC::Instance()->IsMCTruth(label, (AliDielectronSignalMC*)fSignalsMC->At(isig), 2);
1339         // skip if track does not correspond to the signal
1340         if(!isMCtruth1 && !isMCtruth2) continue;
1341         AliDielectronVarManager::Fill(fTracks[i].UncheckedAt(itrack), values);
1342         fHistos->FillClass(className3, AliDielectronVarManager::kNMaxValues, values);
1343       } //loop: tracks
1344     } //loop: arrays
1345
1346   } //loop: MCsignals
1347
1348 }