PWGDQ/dielectron/AliDielectron.cxx

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