PWGDQ/dielectron/AliDielectron.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 ///////////////////////////////////////////////////////////////////////////
  17 //                Dielectron Analysis Main class                         //
  18 //                                                                       //
  19 /*
  20 Framework to perform event selectoin, single track selection and track pair
  21 selection.
  22
  23 Convention for the signs of the pair in fPairCandidates:
  24 The names are available via the function PairClassName(Int_t i)
  25
  26 0: ev1+ ev1+  (same event like sign +)
  27 1: ev1+ ev1-  (same event unlike sign)
  28 2: ev1- ev1-  (same event like sign -)
  29
  30 3: ev1+ ev2+  (mixed event like sign +)
  31 4: ev1- ev2+  (mixed event unlike sign -+)
  32 6: ev1+ ev2-  (mixed event unlike sign +-)
  33 7: ev1- ev2-  (mixed event like sign -)
  34
  35 5: ev2+ ev2+  (same event like sign +)
  36 8: ev2+ ev2-  (same event unlike sign)
  37 9: ev2- ev2-  (same event like sign -)
  38
  39 10: ev1+ ev1- (same event track rotation)
  40
  41 */
  42 //                                                                       //
  43 ///////////////////////////////////////////////////////////////////////////
  44
  45 #include <TString.h>
  46 #include <TList.h>
  47 #include <TMath.h>
  48 #include <TObject.h>
  49
  50 #include <AliESDEvent.h>
  51 #include <AliESDtrack.h>
  52 #include <AliKFParticle.h>
  53
  54 #include <AliEventplane.h>
  55 #include <AliVEvent.h>
  56 #include <AliVParticle.h>
  57 #include <AliVTrack.h>
  58 #include "AliDielectronPair.h"
  59 #include "AliDielectronHistos.h"
  60 #include "AliDielectronCF.h"
  61 #include "AliDielectronMC.h"
  62 #include "AliDielectronVarManager.h"
  63 #include "AliDielectronTrackRotator.h"
  64 #include "AliDielectronDebugTree.h"
  65 #include "AliDielectronSignalMC.h"
  66 #include "AliDielectronMixingHandler.h"
  67
  68 #include "AliDielectron.h"
  69
  70 ClassImp(AliDielectron)
  71
  72 const char* AliDielectron::fgkTrackClassNames[4] = {
  73   "ev1+",
  74   "ev1-",
  75   "ev2+",
  76   "ev2-"
  77 };
  78
  79 const char* AliDielectron::fgkPairClassNames[11] = {
  80   "ev1+_ev1+",
  81   "ev1+_ev1-",
  82   "ev1-_ev1-",
  83   "ev1+_ev2+",
  84   "ev1-_ev2+",
  85   "ev2+_ev2+",
  86   "ev1+_ev2-",
  87   "ev1-_ev2-",
  88   "ev2+_ev2-",
  89   "ev2-_ev2-",
  90   "ev1+_ev1-_TR"
  91 };
  92
  93 //________________________________________________________________
  94 AliDielectron::AliDielectron() :
  95   TNamed("AliDielectron","AliDielectron"),
  96   fEventFilter("EventFilter"),
  97   fTrackFilter("TrackFilter"),
  98   fPairPreFilter("PairPreFilter"),
  99   fPairPreFilterLegs("PairPreFilterLegs"),
 100   fPairFilter("PairFilter"),
 101   fEventPlanePreFilter("EventPlanePreFilter"),
 102   fEventPlanePOIPreFilter("EventPlanePOIPreFilter"),
 103   fPdgMother(443),
 104   fPdgLeg1(11),
 105   fPdgLeg2(11),
 106   fSignalsMC(0x0),
 107   fNoPairing(kFALSE),
 108   fHistos(0x0),
 109   fPairCandidates(new TObjArray(11)),
 110   fCfManagerPair(0x0),
 111   fTrackRotator(0x0),
 112   fDebugTree(0x0),
 113   fMixing(0x0),
 114   fPreFilterEventPlane(kFALSE),
 115   fLikeSignSubEvents(kFALSE),
 116   fPreFilterUnlikeOnly(kFALSE),
 117   fPreFilterAllSigns(kFALSE),
 118   fHasMC(kFALSE),
 119   fStoreRotatedPairs(kFALSE),
 120   fDontClearArrays(kFALSE),
 121   fEstimatorFilename(""),
 122   fTRDpidCorrectionFilename(""),
 123   fVZEROCalibrationFilename(""),
 124   fVZERORecenteringFilename("")
 125 {
 126   //
 127   // Default constructor
 128   //
 129
 130 }
 131
 132 //________________________________________________________________
 133 AliDielectron::AliDielectron(const char* name, const char* title) :
 134   TNamed(name,title),
 135   fEventFilter("EventFilter"),
 136   fTrackFilter("TrackFilter"),
 137   fPairPreFilter("PairPreFilter"),
 138   fPairPreFilterLegs("PairPreFilterLegs"),
 139   fPairFilter("PairFilter"),
 140   fEventPlanePreFilter("EventPlanePreFilter"),
 141   fEventPlanePOIPreFilter("EventPlanePOIPreFilter"),
 142   fPdgMother(443),
 143   fPdgLeg1(11),
 144   fPdgLeg2(11),
 145   fSignalsMC(0x0),
 146   fNoPairing(kFALSE),
 147   fHistos(0x0),
 148   fPairCandidates(new TObjArray(11)),
 149   fCfManagerPair(0x0),
 150   fTrackRotator(0x0),
 151   fDebugTree(0x0),
 152   fMixing(0x0),
 153   fPreFilterEventPlane(kFALSE),
 154   fLikeSignSubEvents(kFALSE),
 155   fPreFilterUnlikeOnly(kFALSE),
 156   fPreFilterAllSigns(kFALSE),
 157   fHasMC(kFALSE),
 158   fStoreRotatedPairs(kFALSE),
 159   fDontClearArrays(kFALSE),
 160   fEstimatorFilename(""),
 161   fTRDpidCorrectionFilename(""),
 162   fVZEROCalibrationFilename(""),
 163   fVZERORecenteringFilename("")
 164 {
 165   //
 166   // Named constructor
 167   //
 168
 169 }
 170
 171 //________________________________________________________________
 172 AliDielectron::~AliDielectron()
 173 {
 174   //
 175   // Default destructor
 176   //
 177   if (fHistos) delete fHistos;
 178   if (fPairCandidates) delete fPairCandidates;
 179   if (fDebugTree) delete fDebugTree;
 180   if (fMixing) delete fMixing;
 181   if (fSignalsMC) delete fSignalsMC;
 182   if (fCfManagerPair) delete fCfManagerPair;
 183 }
 184
 185 //________________________________________________________________
 186 void AliDielectron::Init()
 187 {
 188   //
 189   // Initialise objects
 190   //
 191
 192   if(GetHasMC()) AliDielectronMC::Instance()->SetHasMC(GetHasMC());
 193
 194   InitPairCandidateArrays();
 195
 196   if (fCfManagerPair) {
 197     fCfManagerPair->SetSignalsMC(fSignalsMC);
 198     fCfManagerPair->InitialiseContainer(fPairFilter);
 199   }
 200   if (fTrackRotator)  {
 201     fTrackRotator->SetTrackArrays(&fTracks[0],&fTracks[1]);
 202     fTrackRotator->SetPdgLegs(fPdgLeg1,fPdgLeg2);
 203   }
 204   if (fDebugTree) fDebugTree->SetDielectron(this);
 205   if(fEstimatorFilename.Contains(".root")) AliDielectronVarManager::InitEstimatorAvg(fEstimatorFilename.Data());
 206   if(fTRDpidCorrectionFilename.Contains(".root")) AliDielectronVarManager::InitTRDpidEffHistograms(fTRDpidCorrectionFilename.Data());
 207   if(fVZEROCalibrationFilename.Contains(".root")) AliDielectronVarManager::SetVZEROCalibrationFile(fVZEROCalibrationFilename.Data());
 208   if(fVZERORecenteringFilename.Contains(".root")) AliDielectronVarManager::SetVZERORecenteringFile(fVZERORecenteringFilename.Data());
 209
 210   if (fMixing) fMixing->Init(this);
 211 }
 212
 213 //________________________________________________________________
 214 void AliDielectron::Process(AliVEvent *ev1, AliVEvent *ev2)
 215 {
 216   //
 217   // Process the events
 218   //
 219
 220   //at least first event is needed!
 221   if (!ev1){
 222     AliError("At least first event must be set!");
 223     return;
 224   }
 225   AliDielectronVarManager::SetEvent(ev1);
 226   if (fMixing){
 227     //set mixing bin to event data
 228     Int_t bin=fMixing->FindBin(AliDielectronVarManager::GetData());
 229     AliDielectronVarManager::SetValue(AliDielectronVarManager::kMixingBin,bin);
 230   }
 231
 232   //in case we have MC load the MC event and process the MC particles
 233   if (AliDielectronMC::Instance()->ConnectMCEvent()){
 234     ProcessMC(ev1);
 235   }
 236
 237   //if candidate array doesn't exist, create it
 238   if (!fPairCandidates->UncheckedAt(0)) {
 239     InitPairCandidateArrays();
 240   } else {
 241     ClearArrays();
 242   }
 243
 244   //mask used to require that all cuts are fulfilled
 245   UInt_t selectedMask=(1<<fEventFilter.GetCuts()->GetEntries())-1;
 246
 247   //apply event cuts
 248     if ((ev1&&fEventFilter.IsSelected(ev1)!=selectedMask) ||
 249         (ev2&&fEventFilter.IsSelected(ev2)!=selectedMask)) return;
 250
 251 //   AliDielectronVarManager::SetEvent(ev1); // why a second time???
 252
 253   //fill track arrays for the first event
 254   if (ev1){
 255     FillTrackArrays(ev1);
 256     if (((fPreFilterAllSigns)||(fPreFilterUnlikeOnly)) && ( fPairPreFilter.GetCuts()->GetEntries()>0 )) PairPreFilter(0, 1, fTracks[0], fTracks[1]);
 257   }
 258
 259
 260   //fill track arrays for the second event
 261   if (ev2) {
 262     FillTrackArrays(ev2,1);
 263     if (((fPreFilterAllSigns)||(fPreFilterUnlikeOnly)) && ( fPairPreFilter.GetCuts()->GetEntries()>0 )) PairPreFilter(2, 3, fTracks[2], fTracks[3]);
 264   }
 265
 266   // TPC event plane correction
 267   AliEventplane *cevplane = new AliEventplane();
 268   if (ev1 && fPreFilterEventPlane && ( fEventPlanePreFilter.GetCuts()->GetEntries()>0 || fEventPlanePOIPreFilter.GetCuts()->GetEntries()>0))
 269     EventPlanePreFilter(0, 1, fTracks[0], fTracks[1], ev1, cevplane);
 270
 271   if (!fNoPairing){
 272     // create pairs and fill pair candidate arrays
 273     for (Int_t itrackArr1=0; itrackArr1<4; ++itrackArr1){
 274       for (Int_t itrackArr2=itrackArr1; itrackArr2<4; ++itrackArr2){
 275         FillPairArrays(itrackArr1, itrackArr2);
 276       }
 277     }
 278
 279     //track rotation
 280     if (fTrackRotator) {
 281       fTrackRotator->SetEvent(ev1);
 282       FillPairArrayTR();
 283     }
 284   }
 285
 286   //fill debug tree if a manager is attached
 287   if (fDebugTree) FillDebugTree();
 288
 289   //process event mixing
 290   if (fMixing) {
 291     fMixing->Fill(ev1,this);
 292 //     FillHistograms(0x0,kTRUE);
 293   }
 294
 295   //in case there is a histogram manager, fill the QA histograms
 296   if (fHistos) FillHistograms(ev1);
 297
 298   // clear arrays
 299   if (!fDontClearArrays) ClearArrays();
 300   AliDielectronVarManager::SetTPCEventPlane(0x0);
 301   delete cevplane;
 302 }
 303
 304 //________________________________________________________________
 305 void AliDielectron::ProcessMC(AliVEvent *ev1)
 306 {
 307   //
 308   // Process the MC data
 309   //
 310
 311   AliDielectronMC *dieMC=AliDielectronMC::Instance();
 312
 313   if (fHistos) FillHistogramsMC(dieMC->GetMCEvent(), ev1);
 314
 315   if(!fSignalsMC) return;
 316   //loop over all MC data and Fill the CF container if it exist
 317   if (!fCfManagerPair) return;
 318   fCfManagerPair->SetPdgMother(fPdgMother);
 319   if(!fCfManagerPair->GetStepForMCtruth()) return;
 320
 321   // signals to be studied
 322   Int_t nSignals = fSignalsMC->GetEntries();
 323
 324   // initialize 2D arrays of labels for particles from each MC signal
 325   Int_t** labels1;      // labels for particles satisfying branch 1
 326   Int_t** labels2;      // labels for particles satisfying branch 2
 327   Int_t** labels12;     // labels for particles satisfying both branches
 328   labels1 = new Int_t*[nSignals];
 329   labels2 = new Int_t*[nSignals];
 330   labels12 = new Int_t*[nSignals];
 331   Int_t* indexes1=new Int_t[nSignals];
 332   Int_t* indexes2=new Int_t[nSignals];
 333   Int_t* indexes12=new Int_t[nSignals];
 334   for(Int_t isig=0;isig<nSignals;++isig) {
 335     *(labels1+isig) = new Int_t[dieMC->GetNMCTracks()];
 336     *(labels2+isig) = new Int_t[dieMC->GetNMCTracks()];
 337     *(labels12+isig) = new Int_t[dieMC->GetNMCTracks()];
 338     for(Int_t ip=0; ip<dieMC->GetNMCTracks();++ip) {
 339       labels1[isig][ip] = -1;
 340       labels2[isig][ip] = -1;
 341       labels12[isig][ip] = -1;
 342     }
 343     indexes1[isig]=0;
 344     indexes2[isig]=0;
 345     indexes12[isig]=0;
 346   }
 347
 348   Bool_t truth1=kFALSE;
 349   Bool_t truth2=kFALSE;
 350   // loop over the MC tracks
 351   for(Int_t ipart=0; ipart<dieMC->GetNMCTracks(); ++ipart) {
 352     for(Int_t isig=0; isig<nSignals; ++isig) {       // loop over signals
 353       // Proceed only if this signal is required in the pure MC step
 354       // NOTE: Some signals can be satisfied by many particles and this leads to high
 355       //       computation times (e.g. secondary electrons from the GEANT transport). Be aware of this!!
 356       if(!((AliDielectronSignalMC*)fSignalsMC->At(isig))->GetFillPureMCStep()) continue;
 357
 358       truth1 = dieMC->IsMCTruth(ipart, (AliDielectronSignalMC*)fSignalsMC->At(isig), 1);
 359       truth2 = dieMC->IsMCTruth(ipart, (AliDielectronSignalMC*)fSignalsMC->At(isig), 2);
 360
 361       // particles satisfying both branches are treated separately to avoid double counting during pairing
 362       if(truth1 && truth2) {
 363         labels12[isig][indexes12[isig]] = ipart;
 364         ++indexes12[isig];
 365       }
 366       else {
 367         if(truth1) {
 368           labels1[isig][indexes1[isig]] = ipart;
 369           ++indexes1[isig];
 370         }
 371         if(truth2) {
 372           labels2[isig][indexes2[isig]] = ipart;
 373           ++indexes2[isig];
 374         }
 375       }
 376     }
 377   }  // end loop over MC particles
 378
 379   // Do the pairing and fill the CF container with pure MC info
 380   for(Int_t isig=0; isig<nSignals; ++isig) {
 381     // mix the particles which satisfy only one of the signal branches
 382     for(Int_t i1=0;i1<indexes1[isig];++i1) {
 383       for(Int_t i2=0;i2<indexes2[isig];++i2) {
 384         fCfManagerPair->FillMC(labels1[isig][i1], labels2[isig][i2], isig);
 385       }
 386     }
 387     // mix the particles which satisfy both branches
 388     for(Int_t i1=0;i1<indexes12[isig];++i1) {
 389       for(Int_t i2=0; i2<i1; ++i2) {
 390         fCfManagerPair->FillMC(labels12[isig][i1], labels12[isig][i2], isig);
 391       }
 392     }
 393   }    // end loop over signals
 394
 395   // release the memory
 396   for(Int_t isig=0;isig<nSignals;++isig) {
 397     delete [] *(labels1+isig);
 398     delete [] *(labels2+isig);
 399     delete [] *(labels12+isig);
 400   }
 401   delete [] labels1;
 402   delete [] labels2;
 403   delete [] labels12;
 404   delete [] indexes1;
 405   delete [] indexes2;
 406   delete [] indexes12;
 407 }
 408
 409 //________________________________________________________________
 410 void AliDielectron::FillHistogramsTracks(TObjArray **tracks)
 411 {
 412   //
 413   // Fill Histogram information for tracks after prefilter
 414   // ignore mixed events - for prefilter, only single tracks +/- are relevant
 415   //
 416
 417   TString  className,className2;
 418   Double_t values[AliDielectronVarManager::kNMaxValues];
 419
 420   //Fill track information, separately for the track array candidates
 421   for (Int_t i=0; i<2; ++i){
 422     className.Form("Pre_%s",fgkTrackClassNames[i]);
 423     if (!fHistos->GetHistogramList()->FindObject(className.Data())) continue;
 424     Int_t ntracks=tracks[i]->GetEntriesFast();
 425     for (Int_t itrack=0; itrack<ntracks; ++itrack){
 426       AliDielectronVarManager::Fill(tracks[i]->UncheckedAt(itrack), values);
 427       fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
 428     }
 429   }
 430 }
 431
 432
 433 //________________________________________________________________
 434 void AliDielectron::FillHistogramsMC(const AliMCEvent *ev, AliVEvent *ev1)
 435 {
 436   //
 437   // Fill Histogram information for MCEvents
 438   //
 439
 440   Double_t values[AliDielectronVarManager::kNMaxValues]={0.};
 441   // Fill event information
 442   AliDielectronVarManager::Fill(ev1, values);    // ESD/AOD information
 443   AliDielectronVarManager::Fill(ev, values);     // MC truth info
 444   if (fHistos->GetHistogramList()->FindObject("MCEvent"))
 445     fHistos->FillClass("MCEvent", AliDielectronVarManager::kNMaxValues, values);
 446 }
 447
 448
 449 //________________________________________________________________
 450 void AliDielectron::FillHistograms(const AliVEvent *ev, Bool_t pairInfoOnly)
 451 {
 452   //
 453   // Fill Histogram information for tracks and pairs
 454   //
 455
 456   TString  className,className2;
 457   Double_t values[AliDielectronVarManager::kNMaxValues]={0.};
 458   //Fill event information
 459   if (ev){  //TODO: Why not use GetData() ??? See below event plane stuff!!!
 460     AliDielectronVarManager::Fill(ev, values); //data should already be stored in AliDielectronVarManager from SetEvent, does EV plane correction rely on this???
 461       if (fMixing){
 462     //set mixing bin to event data
 463     Int_t bin=fMixing->FindBin(values);
 464     values[AliDielectronVarManager::kMixingBin]=bin;
 465   }
 466
 467     if (fHistos->GetHistogramList()->FindObject("Event"))
 468 //       fHistos->FillClass("Event", AliDielectronVarManager::kNMaxValues, AliDielectronVarManager::GetData());
 469       fHistos->FillClass("Event", AliDielectronVarManager::kNMaxValues, values); //check event plane stuff and replace with above...
 470   }
 471
 472   //Fill track information, separately for the track array candidates
 473   if (!pairInfoOnly){
 474     for (Int_t i=0; i<4; ++i){
 475       className.Form("Track_%s",fgkTrackClassNames[i]);
 476       if (!fHistos->GetHistogramList()->FindObject(className.Data())) continue;
 477       Int_t ntracks=fTracks[i].GetEntriesFast();
 478       for (Int_t itrack=0; itrack<ntracks; ++itrack){
 479         AliDielectronVarManager::Fill(fTracks[i].UncheckedAt(itrack), values);
 480         fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
 481       }
 482     }
 483   }
 484
 485   //Fill Pair information, separately for all pair candidate arrays and the legs
 486   TObjArray arrLegs(100);
 487   for (Int_t i=0; i<10; ++i){
 488     className.Form("Pair_%s",fgkPairClassNames[i]);
 489     className2.Form("Track_Legs_%s",fgkPairClassNames[i]);
 490     Bool_t pairClass=fHistos->GetHistogramList()->FindObject(className.Data())!=0x0;
 491     Bool_t legClass=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
 492     if (!pairClass&&!legClass) continue;
 493     Int_t ntracks=PairArray(i)->GetEntriesFast();
 494     for (Int_t ipair=0; ipair<ntracks; ++ipair){
 495       AliDielectronPair *pair=static_cast<AliDielectronPair*>(PairArray(i)->UncheckedAt(ipair));
 496
 497       //fill pair information
 498       if (pairClass){
 499         AliDielectronVarManager::Fill(pair, values);
 500         fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
 501       }
 502
 503       //fill leg information, don't fill the information twice
 504       if (legClass){
 505         AliVParticle *d1=pair->GetFirstDaughter();
 506         AliVParticle *d2=pair->GetSecondDaughter();
 507         if (!arrLegs.FindObject(d1)){
 508           AliDielectronVarManager::Fill(d1, values);
 509           fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
 510           arrLegs.Add(d1);
 511         }
 512         if (!arrLegs.FindObject(d2)){
 513           AliDielectronVarManager::Fill(d2, values);
 514           fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
 515           arrLegs.Add(d2);
 516         }
 517       }
 518     }
 519     if (legClass) arrLegs.Clear();
 520   }
 521
 522 }
 523 //________________________________________________________________
 524 void AliDielectron::FillHistogramsPair(AliDielectronPair *pair,Bool_t fromPreFilter/*=kFALSE*/)
 525 {
 526   //
 527   // Fill Histogram information for pairs and the track in the pair
 528   // NOTE: in this funtion the leg information may be filled multiple
 529   //       times. This funtion is used in the track rotation pairing
 530   //       and those legs are not saved!
 531   //
 532   TString  className,className2;
 533   Double_t values[AliDielectronVarManager::kNMaxValues];
 534
 535   //Fill Pair information, separately for all pair candidate arrays and the legs
 536   TObjArray arrLegs(100);
 537   const Int_t type=pair->GetType();
 538   if (fromPreFilter) {
 539     className.Form("RejPair_%s",fgkPairClassNames[type]);
 540     className2.Form("RejTrack_%s",fgkPairClassNames[type]);
 541   } else {
 542     className.Form("Pair_%s",fgkPairClassNames[type]);
 543     className2.Form("Track_Legs_%s",fgkPairClassNames[type]);
 544   }
 545
 546   Bool_t pairClass=fHistos->GetHistogramList()->FindObject(className.Data())!=0x0;
 547   Bool_t legClass=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
 548
 549   //fill pair information
 550   if (pairClass){
 551     AliDielectronVarManager::Fill(pair, values);
 552     fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
 553   }
 554
 555   if (legClass){
 556     AliVParticle *d1=pair->GetFirstDaughter();
 557     AliDielectronVarManager::Fill(d1, values);
 558     fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
 559
 560     AliVParticle *d2=pair->GetSecondDaughter();
 561     AliDielectronVarManager::Fill(d2, values);
 562     fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
 563   }
 564 }
 565
 566 //________________________________________________________________
 567 void AliDielectron::FillTrackArrays(AliVEvent * const ev, Int_t eventNr)
 568 {
 569   //
 570   // select tracks and fill track candidate arrays
 571   // eventNr = 0: First  event, use track arrays 0 and 1
 572   // eventNr = 1: Second event, use track arrays 2 and 3
 573   //
 574
 575   Int_t ntracks=ev->GetNumberOfTracks();
 576   UInt_t selectedMask=(1<<fTrackFilter.GetCuts()->GetEntries())-1;
 577   for (Int_t itrack=0; itrack<ntracks; ++itrack){
 578     //get particle
 579     AliVParticle *particle=ev->GetTrack(itrack);
 580     //TODO: temporary solution, perhaps think about a better implementation
 581     //      This is needed to use AliESDpidCuts, which relies on the ESD event
 582     //      is set as a AliESDtrack attribute... somehow ugly!
 583     if (ev->IsA()==AliESDEvent::Class()){
 584       AliESDtrack *track=static_cast<AliESDtrack*>(particle);
 585       track->SetESDEvent(static_cast<AliESDEvent*>(ev)); //only in trunk...
 586     }
 587
 588     //apply track cuts
 589     if (fTrackFilter.IsSelected(particle)!=selectedMask) continue;
 590
 591     //fill selected particle into the corresponding track arrays
 592     Short_t charge=particle->Charge();
 593     if (charge>0)      fTracks[eventNr*2].Add(particle);
 594     else if (charge<0) fTracks[eventNr*2+1].Add(particle);
 595   }
 596 }
 597
 598 //________________________________________________________________
 599 void AliDielectron::EventPlanePreFilter(Int_t arr1, Int_t arr2, TObjArray arrTracks1, TObjArray arrTracks2, const AliVEvent *ev, AliEventplane *cevplane)
 600 {
 601   //
 602   // Prefilter tracks and tracks from pairs
 603   // Needed for rejection in the Q-Vector of the event plane
 604   // remove contribution of all tracks to the Q-vector that are in invariant mass window
 605   //
 606   AliEventplane *evplane = const_cast<AliVEvent *>(ev)->GetEventplane();
 607 //   AliEventplane *evplane = ev->GetEventplane();
 608   if(!evplane) return;
 609
 610   // do not change these vectors qref
 611   TVector2 * const qref  = evplane->GetQVector();
 612   if(!qref) return;
 613   // random subevents
 614   TVector2 *qrsub1 = evplane->GetQsub1();
 615   TVector2 *qrsub2 = evplane->GetQsub2();
 616
 617   // copy references
 618   TVector2 *qcorr  = new TVector2(*qref);
 619   TVector2 *qcsub1 = 0x0;
 620   TVector2 *qcsub2 = 0x0;
 621   //  printf("qrsub1 %p %f \n",qrsub1,qrsub1->X());
 622
 623
 624   // eta gap ?
 625   Bool_t etagap = kFALSE;
 626   for (Int_t iCut=0; iCut<fEventPlanePreFilter.GetCuts()->GetEntries();++iCut) {
 627     TString cutName=fEventPlanePreFilter.GetCuts()->At(iCut)->GetName();
 628     if(cutName.Contains("eta") || cutName.Contains("Eta"))  etagap=kTRUE;
 629   }
 630
 631   // subevent separation
 632   if(fLikeSignSubEvents || etagap) {
 633     qcsub1 = new TVector2(*qcorr);
 634     qcsub2 = new TVector2(*qcorr);
 635
 636     Int_t ntracks=ev->GetNumberOfTracks();
 637
 638     // track removals
 639     for (Int_t itrack=0; itrack<ntracks; ++itrack){
 640       AliVParticle *particle=ev->GetTrack(itrack);
 641       AliVTrack *track= static_cast<AliVTrack*>(particle);
 642       if (!track) continue;
 643
 644       Double_t cQX     = evplane->GetQContributionX(track);
 645       Double_t cQY     = evplane->GetQContributionY(track);
 646
 647       // by charge sub1+ sub2-
 648       if(fLikeSignSubEvents) {
 649         Short_t charge=track->Charge();
 650         if (charge<0) qcsub1->Set(qcsub1->X()-cQX, qcsub1->Y()-cQY);
 651         if (charge>0) qcsub2->Set(qcsub2->X()-cQX, qcsub2->Y()-cQY);
 652       }
 653       // by eta sub1+ sub2-
 654       if(etagap) {
 655         Double_t eta=track->Eta();
 656         if (eta<0.0) qcsub1->Set(qcsub1->X()-cQX, qcsub1->Y()-cQY);
 657         if (eta>0.0) qcsub2->Set(qcsub2->X()-cQX, qcsub2->Y()-cQY);
 658       }
 659     }
 660   }
 661   else {
 662     // by a random
 663     qcsub1 = new TVector2(*qrsub1);
 664     qcsub2 = new TVector2(*qrsub2);
 665   }
 666
 667   // apply cuts, e.g. etagap
 668   if(fEventPlanePreFilter.GetCuts()->GetEntries()) {
 669     UInt_t selectedMask=(1<<fEventPlanePreFilter.GetCuts()->GetEntries())-1;
 670     Int_t ntracks=ev->GetNumberOfTracks();
 671     for (Int_t itrack=0; itrack<ntracks; ++itrack){
 672       AliVParticle *particle=ev->GetTrack(itrack);
 673       AliVTrack *track= static_cast<AliVTrack*>(particle);
 674       if (!track) continue;
 675
 676       //event plane cuts
 677       UInt_t cutMask=fEventPlanePreFilter.IsSelected(track);
 678       //apply cut
 679       if (cutMask==selectedMask) continue;
 680
 681       Double_t cQX     = 0.0;
 682       Double_t cQY     = 0.0;
 683       if(!etagap) {
 684         cQX = evplane->GetQContributionX(track);
 685         cQY = evplane->GetQContributionY(track);
 686       }
 687       Double_t cQXsub1 = evplane->GetQContributionXsub1(track);
 688       Double_t cQYsub1 = evplane->GetQContributionYsub1(track);
 689       Double_t cQXsub2 = evplane->GetQContributionXsub2(track);
 690       Double_t cQYsub2 = evplane->GetQContributionYsub2(track);
 691
 692       // update Q vectors
 693       qcorr->Set(qcorr->X()-cQX, qcorr->Y()-cQY);
 694       qcsub1->Set(qcsub1->X()-cQXsub1, qcsub1->Y()-cQYsub1);
 695       qcsub2->Set(qcsub2->X()-cQXsub2, qcsub2->Y()-cQYsub2);
 696     }
 697   }
 698
 699   // POI (particle of interest) rejection
 700   Int_t pairIndex=GetPairIndex(arr1,arr2);
 701
 702   Int_t ntrack1=arrTracks1.GetEntriesFast();
 703   Int_t ntrack2=arrTracks2.GetEntriesFast();
 704   AliDielectronPair candidate;
 705
 706   UInt_t selectedMask=(1<<fEventPlanePOIPreFilter.GetCuts()->GetEntries())-1;
 707   for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1){
 708     Int_t end=ntrack2;
 709     if (arr1==arr2) end=itrack1;
 710     Bool_t accepted=kFALSE;
 711     for (Int_t itrack2=0; itrack2<end; ++itrack2){
 712       TObject *track1=arrTracks1.UncheckedAt(itrack1);
 713       TObject *track2=arrTracks2.UncheckedAt(itrack2);
 714       if (!track1 || !track2) continue;
 715       //create the pair
 716       candidate.SetTracks(static_cast<AliVTrack*>(track1), fPdgLeg1,
 717                           static_cast<AliVTrack*>(track2), fPdgLeg2);
 718
 719       candidate.SetType(pairIndex);
 720       candidate.SetLabel(AliDielectronMC::Instance()->GetLabelMotherWithPdg(&candidate,fPdgMother));
 721
 722       //event plane cuts
 723       UInt_t cutMask=fEventPlanePOIPreFilter.IsSelected(&candidate);
 724       //apply cut
 725       if (cutMask==selectedMask) continue;
 726
 727       accepted=kTRUE;
 728       //remove the tracks from the Track arrays
 729       arrTracks2.AddAt(0x0,itrack2);
 730     }
 731       if ( accepted ) arrTracks1.AddAt(0x0,itrack1);
 732   }
 733   //compress the track arrays
 734   arrTracks1.Compress();
 735   arrTracks2.Compress();
 736
 737
 738   //Modify the components: subtract the tracks
 739   ntrack1=arrTracks1.GetEntriesFast();
 740   ntrack2=arrTracks2.GetEntriesFast();
 741
 742   // remove leg1 contribution
 743   for (Int_t itrack=0; itrack<ntrack1; ++itrack){
 744     AliVTrack *track= static_cast<AliVTrack*>(arrTracks1.UncheckedAt(itrack));
 745     if (!track) continue;
 746
 747     Double_t cQX     = evplane->GetQContributionX(track);
 748     Double_t cQY     = evplane->GetQContributionY(track);
 749     Double_t cQXsub1 = evplane->GetQContributionXsub1(track);
 750     Double_t cQYsub1 = evplane->GetQContributionYsub1(track);
 751     Double_t cQXsub2 = evplane->GetQContributionXsub2(track);
 752     Double_t cQYsub2 = evplane->GetQContributionYsub2(track);
 753
 754     // update Q vectors
 755     qcorr->Set(qcorr->X()-cQX, qcorr->Y()-cQY);
 756     qcsub1->Set(qcsub1->X()-cQXsub1, qcsub1->Y()-cQYsub1);
 757     qcsub2->Set(qcsub2->X()-cQXsub2, qcsub2->Y()-cQYsub2);
 758   }
 759   // remove leg2 contribution
 760   for (Int_t itrack=0; itrack<ntrack2; ++itrack){
 761     AliVTrack *track= static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack));
 762     if (!track) continue;
 763
 764     Double_t cQX     = evplane->GetQContributionX(track);
 765     Double_t cQY     = evplane->GetQContributionY(track);
 766     Double_t cQXsub1 = evplane->GetQContributionXsub1(track);
 767     Double_t cQYsub1 = evplane->GetQContributionYsub1(track);
 768     Double_t cQXsub2 = evplane->GetQContributionXsub2(track);
 769     Double_t cQYsub2 = evplane->GetQContributionYsub2(track);
 770
 771     // update Q vectors
 772     qcorr->Set(qcorr->X()-cQX, qcorr->Y()-cQY);
 773     qcsub1->Set(qcsub1->X()-cQXsub1, qcsub1->Y()-cQYsub1);
 774     qcsub2->Set(qcsub2->X()-cQXsub2, qcsub2->Y()-cQYsub2);
 775   }
 776
 777   //  printf("qrsub1 %p %f \t qcsub1 %p %f \n",qrsub1,qrsub1->X(),qcsub1,qcsub1->X());
 778   // set AliEventplane with corrected values
 779   cevplane->SetQVector(qcorr);
 780   cevplane->SetQsub(qcsub1, qcsub2);
 781   AliDielectronVarManager::SetTPCEventPlane(cevplane);
 782 }
 783
 784 //________________________________________________________________
 785 void AliDielectron::PairPreFilter(Int_t arr1, Int_t arr2, TObjArray &arrTracks1, TObjArray &arrTracks2)
 786 {
 787   //
 788   // Prefilter tracks from pairs
 789   // Needed for datlitz rejections
 790   // remove all tracks from the Single track arrays that pass the cuts in this filter
 791   //
 792
 793   Int_t ntrack1=arrTracks1.GetEntriesFast();
 794   Int_t ntrack2=arrTracks2.GetEntriesFast();
 795   AliDielectronPair candidate;
 796
 797   // flag arrays for track removal
 798   Bool_t *bTracks1 = new Bool_t[ntrack1];
 799   for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1) bTracks1[itrack1]=kFALSE;
 800   Bool_t *bTracks2 = new Bool_t[ntrack2];
 801   for (Int_t itrack2=0; itrack2<ntrack2; ++itrack2) bTracks2[itrack2]=kFALSE;
 802
 803   UInt_t selectedMask=(1<<fPairPreFilter.GetCuts()->GetEntries())-1;
 804   UInt_t selectedMaskPair=(1<<fPairFilter.GetCuts()->GetEntries())-1;
 805
 806   Int_t nRejPasses = 1; //for fPreFilterUnlikeOnly and no set flag
 807   if (fPreFilterAllSigns) nRejPasses = 3;
 808
 809   for (Int_t iRP=0; iRP < nRejPasses; ++iRP) {
 810         Int_t arr1RP=arr1, arr2RP=arr2;
 811         TObjArray *arrTracks1RP=&arrTracks1;
 812         TObjArray *arrTracks2RP=&arrTracks2;
 813         Bool_t *bTracks1RP = bTracks1;
 814         Bool_t *bTracks2RP = bTracks2;
 815         switch (iRP) {
 816                 case 1: arr1RP=arr1;arr2RP=arr1;
 817                                 arrTracks1RP=&arrTracks1;
 818                                 arrTracks2RP=&arrTracks1;
 819                                 bTracks1RP = bTracks1;
 820                                 bTracks2RP = bTracks1;
 821                                 break;
 822                 case 2: arr1RP=arr2;arr2RP=arr2;
 823                                 arrTracks1RP=&arrTracks2;
 824                                 arrTracks2RP=&arrTracks2;
 825                                 bTracks1RP = bTracks2;
 826                                 bTracks2RP = bTracks2;
 827                                 break;
 828                 default: ;//nothing to do
 829         }
 830         Int_t ntrack1RP=(*arrTracks1RP).GetEntriesFast();
 831         Int_t ntrack2RP=(*arrTracks2RP).GetEntriesFast();
 832
 833         Int_t pairIndex=GetPairIndex(arr1RP,arr2RP);
 834
 835         for (Int_t itrack1=0; itrack1<ntrack1RP; ++itrack1){
 836           Int_t end=ntrack2RP;
 837           if (arr1RP==arr2RP) end=itrack1;
 838           for (Int_t itrack2=0; itrack2<end; ++itrack2){
 839                 TObject *track1=(*arrTracks1RP).UncheckedAt(itrack1);
 840                 TObject *track2=(*arrTracks2RP).UncheckedAt(itrack2);
 841                 if (!track1 || !track2) continue;
 842                 //create the pair
 843                 candidate.SetTracks(static_cast<AliVTrack*>(track1), fPdgLeg1,
 844                         static_cast<AliVTrack*>(track2), fPdgLeg2);
 845
 846                 candidate.SetType(pairIndex);
 847                 candidate.SetLabel(AliDielectronMC::Instance()->GetLabelMotherWithPdg(&candidate,fPdgMother));
 848                 //relate to the production vertex
 849                 //       if (AliDielectronVarManager::GetKFVertex()) candidate.SetProductionVertex(*AliDielectronVarManager::GetKFVertex());
 850
 851                 //pair cuts
 852                 UInt_t cutMask=fPairPreFilter.IsSelected(&candidate);
 853
 854                 //apply cut
 855                 if (cutMask!=selectedMask) continue;
 856                 if (fCfManagerPair) fCfManagerPair->Fill(selectedMaskPair+1 ,&candidate);
 857                 if (fHistos) FillHistogramsPair(&candidate,kTRUE);
 858                 //set flags for track removal
 859                 bTracks1RP[itrack1]=kTRUE;
 860                 bTracks2RP[itrack2]=kTRUE;
 861           }
 862         }
 863   }
 864
 865   //remove the tracks from the Track arrays
 866   for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1){
 867     if(bTracks1[itrack1]) arrTracks1.AddAt(0x0, itrack1);
 868   }
 869   for (Int_t itrack2=0; itrack2<ntrack2; ++itrack2){
 870     if(bTracks2[itrack2]) arrTracks2.AddAt(0x0, itrack2);
 871   }
 872
 873   // clean up
 874   delete [] bTracks1;
 875   delete [] bTracks2;
 876
 877   //compress the track arrays
 878   arrTracks1.Compress();
 879   arrTracks2.Compress();
 880
 881   //apply leg cuts after the pre filter
 882   if ( fPairPreFilterLegs.GetCuts()->GetEntries()>0 ) {
 883     selectedMask=(1<<fPairPreFilterLegs.GetCuts()->GetEntries())-1;
 884     //loop over tracks from array 1
 885     for (Int_t itrack=0; itrack<arrTracks1.GetEntriesFast();++itrack){
 886       //test cuts
 887       UInt_t cutMask=fPairPreFilterLegs.IsSelected(arrTracks1.UncheckedAt(itrack));
 888
 889       //apply cut
 890       if (cutMask!=selectedMask) arrTracks1.AddAt(0x0,itrack);;
 891     }
 892     arrTracks1.Compress();
 893
 894     //in case of like sign don't loop over second array
 895     if (arr1==arr2) {
 896       arrTracks2=arrTracks1;
 897     } else {
 898
 899       //loop over tracks from array 2
 900       for (Int_t itrack=0; itrack<arrTracks2.GetEntriesFast();++itrack){
 901       //test cuts
 902         UInt_t cutMask=fPairPreFilterLegs.IsSelected(arrTracks2.UncheckedAt(itrack));
 903       //apply cut
 904         if (cutMask!=selectedMask) arrTracks2.AddAt(0x0,itrack);
 905       }
 906       arrTracks2.Compress();
 907
 908     }
 909   }
 910   //For unlike-sign monitor track-cuts:
 911   if (arr1!=arr2&&fHistos) {
 912     TObjArray *unlikesignArray[2] = {&arrTracks1,&arrTracks2};
 913     FillHistogramsTracks(unlikesignArray);
 914   }
 915 }
 916
 917 //________________________________________________________________
 918 void AliDielectron::FillPairArrays(Int_t arr1, Int_t arr2)
 919 {
 920   //
 921   // select pairs and fill pair candidate arrays
 922   //
 923
 924   TObjArray arrTracks1=fTracks[arr1];
 925   TObjArray arrTracks2=fTracks[arr2];
 926
 927   //process pre filter if set
 928   if ((!fPreFilterAllSigns) && (!fPreFilterUnlikeOnly) && ( fPairPreFilter.GetCuts()->GetEntries()>0 ))  PairPreFilter(arr1, arr2, arrTracks1, arrTracks2);
 929
 930   Int_t pairIndex=GetPairIndex(arr1,arr2);
 931
 932   Int_t ntrack1=arrTracks1.GetEntriesFast();
 933   Int_t ntrack2=arrTracks2.GetEntriesFast();
 934
 935   AliDielectronPair *candidate=new AliDielectronPair;
 936
 937   UInt_t selectedMask=(1<<fPairFilter.GetCuts()->GetEntries())-1;
 938
 939   for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1){
 940     Int_t end=ntrack2;
 941     if (arr1==arr2) end=itrack1;
 942     for (Int_t itrack2=0; itrack2<end; ++itrack2){
 943       //create the pair
 944       candidate->SetTracks(static_cast<AliVTrack*>(arrTracks1.UncheckedAt(itrack1)), fPdgLeg1,
 945                              static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack2)), fPdgLeg2);
 946       candidate->SetType(pairIndex);
 947       Int_t label=AliDielectronMC::Instance()->GetLabelMotherWithPdg(candidate,fPdgMother);
 948       candidate->SetLabel(label);
 949       if (label>-1) candidate->SetPdgCode(fPdgMother);
 950
 951       //pair cuts
 952       UInt_t cutMask=fPairFilter.IsSelected(candidate);
 953
 954       //CF manager for the pair
 955       if (fCfManagerPair) fCfManagerPair->Fill(cutMask,candidate);
 956
 957       //apply cut
 958       if (cutMask!=selectedMask) continue;
 959
 960       //add the candidate to the candidate array
 961       PairArray(pairIndex)->Add(candidate);
 962       //get a new candidate
 963       candidate=new AliDielectronPair;
 964     }
 965   }
 966   //delete the surplus candidate
 967   delete candidate;
 968 }
 969
 970 //________________________________________________________________
 971 void AliDielectron::FillPairArrayTR()
 972 {
 973   //
 974   // select pairs and fill pair candidate arrays
 975   //
 976   UInt_t selectedMask=(1<<fPairFilter.GetCuts()->GetEntries())-1;
 977
 978   while ( fTrackRotator->NextCombination() ){
 979     AliDielectronPair candidate;
 980     candidate.SetTracks(&fTrackRotator->GetKFTrackP(), &fTrackRotator->GetKFTrackN(),
 981                         fTrackRotator->GetVTrackP(),fTrackRotator->GetVTrackN());
 982     candidate.SetType(kEv1PMRot);
 983
 984     //pair cuts
 985     UInt_t cutMask=fPairFilter.IsSelected(&candidate);
 986
 987     //CF manager for the pair
 988     if (fCfManagerPair) fCfManagerPair->Fill(cutMask,&candidate);
 989
 990     //apply cut
 991     if (cutMask==selectedMask) {
 992      if(fHistos) FillHistogramsPair(&candidate);
 993      if(fStoreRotatedPairs) PairArray(kEv1PMRot)->Add(new AliDielectronPair(candidate));
 994     }
 995   }
 996 }
 997
 998 //________________________________________________________________
 999 void AliDielectron::FillDebugTree()
1000 {
1001   //
1002   // Fill Histogram information for tracks and pairs
1003   //
1004
1005   //Fill Debug tree
1006   for (Int_t i=0; i<10; ++i){
1007     Int_t ntracks=PairArray(i)->GetEntriesFast();
1008     for (Int_t ipair=0; ipair<ntracks; ++ipair){
1009       fDebugTree->Fill(static_cast<AliDielectronPair*>(PairArray(i)->UncheckedAt(ipair)));
1010     }
1011   }
1012 }
1013
1014 //________________________________________________________________
1015 void AliDielectron::SaveDebugTree()
1016 {
1017   //
1018   // delete the debug tree, this will also write the tree
1019   //
1020   if (fDebugTree) fDebugTree->DeleteStreamer();
1021 }
1022
1023
1024 //__________________________________________________________________
1025 void AliDielectron::AddSignalMC(AliDielectronSignalMC* signal) {
1026   //
1027   //  Add an MC signal to the signals list
1028   //
1029   if(!fSignalsMC) {
1030     fSignalsMC = new TObjArray();
1031     fSignalsMC->SetOwner();
1032   }
1033   fSignalsMC->Add(signal);
1034 }