PWGDQ/dielectron/AliDielectron.cxx

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