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