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