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