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