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