PWG3/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
  49 #include <AliESDEvent.h>
  50 #include <AliESDtrack.h>
  51
  52 #include <AliVEvent.h>
  53 #include <AliVParticle.h>
  54 #include <AliVTrack.h>
  55 #include "AliDielectronPair.h"
  56 #include "AliDielectronHistos.h"
  57 #include "AliDielectronCF.h"
  58 #include "AliDielectronMC.h"
  59 #include "AliDielectronVarManager.h"
  60 #include "AliDielectronTrackRotator.h"
  61 #include "AliDielectronDebugTree.h"
  62
  63 #include "AliDielectron.h"
  64
  65 ClassImp(AliDielectron)
  66
  67 const char* AliDielectron::fgkTrackClassNames[4] = {
  68   "ev1+",
  69   "ev1-",
  70   "ev2+",
  71   "ev2-"
  72 };
  73
  74 const char* AliDielectron::fgkPairClassNames[11] = {
  75   "ev1+_ev1+",
  76   "ev1+_ev1-",
  77   "ev1-_ev1-",
  78   "ev1+_ev2+",
  79   "ev1-_ev2+",
  80   "ev2+_ev2+",
  81   "ev1+_ev2-",
  82   "ev1-_ev2-",
  83   "ev2+_ev2-",
  84   "ev2-_ev2-",
  85   "ev1+_ev1-_TR"
  86 };
  87
  88 //________________________________________________________________
  89 AliDielectron::AliDielectron() :
  90   TNamed("AliDielectron","AliDielectron"),
  91   fEventFilter("EventFilter"),
  92   fTrackFilter("TrackFilter"),
  93   fPairPreFilter("PairPreFilter"),
  94   fPairPreFilterLegs("PairPreFilterLegs"),
  95   fPairFilter("PairFilter"),
  96   fPdgMother(443),
  97   fPdgLeg1(11),
  98   fPdgLeg2(11),
  99   fHistos(0x0),
 100   fPairCandidates(new TObjArray(10)),
 101   fCfManagerPair(0x0),
 102   fTrackRotator(0x0),
 103   fDebugTree(0x0)
 104 {
 105   //
 106   // Default constructor
 107   //
 108
 109 }
 110
 111 //________________________________________________________________
 112 AliDielectron::AliDielectron(const char* name, const char* title) :
 113   TNamed(name,title),
 114   fEventFilter("EventFilter"),
 115   fTrackFilter("TrackFilter"),
 116   fPairPreFilter("PairPreFilter"),
 117   fPairPreFilterLegs("PairPreFilterLegs"),
 118   fPairFilter("PairFilter"),
 119   fPdgMother(443),
 120   fPdgLeg1(11),
 121   fPdgLeg2(11),
 122   fHistos(0x0),
 123   fPairCandidates(new TObjArray(10)),
 124   fCfManagerPair(0x0),
 125   fTrackRotator(0x0),
 126   fDebugTree(0x0)
 127 {
 128   //
 129   // Named constructor
 130   //
 131
 132 }
 133
 134 //________________________________________________________________
 135 AliDielectron::~AliDielectron()
 136 {
 137   //
 138   // Default destructor
 139   //
 140   if (fHistos) delete fHistos;
 141   if (fPairCandidates) delete fPairCandidates;
 142   if (fDebugTree) delete fDebugTree;
 143 }
 144
 145 //________________________________________________________________
 146 void AliDielectron::Init()
 147 {
 148   //
 149   // Initialise objects
 150   //
 151   if (fCfManagerPair) fCfManagerPair->InitialiseContainer(fPairFilter);
 152   if (fTrackRotator)  fTrackRotator->SetTrackArrays(&fTracks[0],&fTracks[1]);
 153   if (fDebugTree) fDebugTree->SetDielectron(this);
 154 }
 155
 156 //________________________________________________________________
 157 void AliDielectron::Process(AliVEvent *ev1, AliVEvent *ev2)
 158 {
 159   //
 160   // Process the events
 161   //
 162
 163   AliDielectronVarManager::SetEvent(ev1);
 164
 165   //in case we have MC load the MC event and process the MC particles
 166   if (AliDielectronMC::Instance()->HasMC()) {
 167     if (!AliDielectronMC::Instance()->ConnectMCEvent()){
 168       AliError("Could not properly connect the MC event, skipping this event!");
 169       return;
 170     }
 171     ProcessMC();
 172   }
 173
 174   //if candidate array doesn't exist, create it
 175   if (!fPairCandidates->UncheckedAt(0)) {
 176     InitPairCandidateArrays();
 177   } else {
 178     ClearArrays();
 179   }
 180
 181   //mask used to require that all cuts are fulfilled
 182   UInt_t selectedMask=(1<<fEventFilter.GetCuts()->GetEntries())-1;
 183
 184   //apply event cuts
 185     if ((ev1&&fEventFilter.IsSelected(ev1)!=selectedMask) ||
 186         (ev2&&fEventFilter.IsSelected(ev2)!=selectedMask)) return;
 187
 188   AliDielectronVarManager::SetEvent(ev1);
 189
 190   //fill track arrays for the first event
 191   if (ev1) FillTrackArrays(ev1);
 192
 193   //fill track arrays for the second event
 194   if (ev2) FillTrackArrays(ev2,1);
 195
 196   // create pairs and fill pair candidate arrays
 197   for (Int_t itrackArr1=0; itrackArr1<4; ++itrackArr1){
 198     for (Int_t itrackArr2=itrackArr1; itrackArr2<4; ++itrackArr2){
 199       FillPairArrays(itrackArr1, itrackArr2);
 200     }
 201   }
 202
 203   //track rotation
 204   if (fTrackRotator) FillPairArrayTR();
 205
 206   //in case there is a histogram manager, fill the QA histograms
 207   if (fHistos) FillHistograms(ev1);
 208
 209   //fill debug tree if a manager is attached
 210   if (fDebugTree) FillDebugTree();
 211 }
 212
 213 //________________________________________________________________
 214 void AliDielectron::ProcessMC()
 215 {
 216   //
 217   // Process the MC data
 218   //
 219
 220   //loop over all MC data and Fill the CF container if it exist
 221   if (!fCfManagerPair) return;
 222   fCfManagerPair->SetPdgMother(fPdgMother);
 223   AliDielectronMC *dieMC=AliDielectronMC::Instance();
 224   for (Int_t ipart=0; ipart<dieMC->GetNMCTracks();++ipart){
 225     //TODO: MC truth cut properly!!!
 226     AliVParticle *mcPart=dieMC->GetMCTrackFromMCEvent(ipart);
 227     if (!dieMC->IsMCMotherToEE(mcPart, fPdgMother)) continue;
 228     fCfManagerPair->FillMC(mcPart);
 229   }
 230 }
 231
 232 //________________________________________________________________
 233 void AliDielectron::FillHistograms(const AliVEvent *ev)
 234 {
 235   //
 236   // Fill Histogram information for tracks and pairs
 237   //
 238
 239   TString  className,className2;
 240   Double_t values[AliDielectronVarManager::kNMaxValues];
 241   //Fill event information
 242   AliDielectronVarManager::Fill(ev, values);
 243   if (fHistos->GetHistogramList()->FindObject("Event"))
 244     fHistos->FillClass("Event", AliDielectronVarManager::kNMaxValues, values);
 245
 246   //Fill track information, separately for the track array candidates
 247   for (Int_t i=0; i<4; ++i){
 248     className.Form("Track_%s",fgkTrackClassNames[i]);
 249     if (!fHistos->GetHistogramList()->FindObject(className.Data())) continue;
 250     Int_t ntracks=fTracks[i].GetEntriesFast();
 251     for (Int_t itrack=0; itrack<ntracks; ++itrack){
 252       AliDielectronVarManager::Fill(fTracks[i].UncheckedAt(itrack), values);
 253       fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
 254     }
 255   }
 256
 257   //Fill Pair information, separately for all pair candidate arrays and the legs
 258   TObjArray arrLegs(100);
 259   for (Int_t i=0; i<10; ++i){
 260     className.Form("Pair_%s",fgkPairClassNames[i]);
 261     className2.Form("Track_Legs_%s",fgkPairClassNames[i]);
 262     Bool_t pairClass=fHistos->GetHistogramList()->FindObject(className.Data())!=0x0;
 263     Bool_t legClass=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
 264     if (!pairClass&&!legClass) continue;
 265     Int_t ntracks=PairArray(i)->GetEntriesFast();
 266     for (Int_t ipair=0; ipair<ntracks; ++ipair){
 267       AliDielectronPair *pair=static_cast<AliDielectronPair*>(PairArray(i)->UncheckedAt(ipair));
 268
 269       //fill pair information
 270       if (pairClass){
 271         AliDielectronVarManager::Fill(pair, values);
 272         fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
 273       }
 274
 275       //fill leg information, don't fill the information twice
 276       if (legClass){
 277         AliVParticle *d1=pair->GetFirstDaughter();
 278         AliVParticle *d2=pair->GetSecondDaughter();
 279         if (!arrLegs.FindObject(d1)){
 280           AliDielectronVarManager::Fill(d1, values);
 281           fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
 282           arrLegs.Add(d1);
 283         }
 284         if (!arrLegs.FindObject(d2)){
 285           AliDielectronVarManager::Fill(d2, values);
 286           fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
 287           arrLegs.Add(d2);
 288         }
 289       }
 290     }
 291     if (legClass) arrLegs.Clear();
 292   }
 293
 294 }
 295 //________________________________________________________________
 296 void AliDielectron::FillHistogramsPair(AliDielectronPair *pair)
 297 {
 298   //
 299   // Fill Histogram information for pairs and the track in the pair
 300   // NOTE: in this funtion the leg information may be filled multiple
 301   //       times. This funtion is used in the track rotation pairing
 302   //       and those legs are not saved!
 303   //
 304   TString  className,className2;
 305   Double_t values[AliDielectronVarManager::kNMaxValues];
 306
 307   //Fill Pair information, separately for all pair candidate arrays and the legs
 308   TObjArray arrLegs(100);
 309   const Int_t type=pair->GetType();
 310   className.Form("Pair_%s",fgkPairClassNames[type]);
 311   className2.Form("Track_Legs_%s",fgkPairClassNames[type]);
 312
 313   Bool_t pairClass=fHistos->GetHistogramList()->FindObject(className.Data())!=0x0;
 314   Bool_t legClass=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
 315
 316   //fill pair information
 317   if (pairClass){
 318     AliDielectronVarManager::Fill(pair, values);
 319     fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
 320   }
 321
 322   if (legClass){
 323     AliVParticle *d1=pair->GetFirstDaughter();
 324     AliDielectronVarManager::Fill(d1, values);
 325     fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
 326
 327     AliVParticle *d2=pair->GetSecondDaughter();
 328     AliDielectronVarManager::Fill(d2, values);
 329     fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
 330   }
 331 }
 332
 333 //________________________________________________________________
 334 void AliDielectron::FillTrackArrays(AliVEvent * const ev, Int_t eventNr)
 335 {
 336   //
 337   // select tracks and fill track candidate arrays
 338   // eventNr = 0: First  event, use track arrays 0 and 1
 339   // eventNr = 1: Second event, use track arrays 2 and 3
 340   //
 341
 342   Int_t ntracks=ev->GetNumberOfTracks();
 343   UInt_t selectedMask=(1<<fTrackFilter.GetCuts()->GetEntries())-1;
 344   for (Int_t itrack=0; itrack<ntracks; ++itrack){
 345     //get particle
 346     AliVParticle *particle=ev->GetTrack(itrack);
 347     //TODO: temporary solution, perhaps think about a better implementation
 348     //      This is needed to use AliESDpidCuts, which relies on the ESD event
 349     //      is set as a AliESDtrack attribute... somehow ugly!
 350     if (ev->IsA()==AliESDEvent::Class()){
 351       AliESDtrack *track=static_cast<AliESDtrack*>(particle);
 352       track->SetESDEvent(static_cast<AliESDEvent*>(ev)); //only in trunk...
 353     }
 354
 355     //apply track cuts
 356     if (fTrackFilter.IsSelected(particle)!=selectedMask) continue;
 357
 358     //fill selected particle into the corresponding track arrays
 359     Short_t charge=particle->Charge();
 360     if (charge>0)      fTracks[eventNr*2].Add(particle);
 361     else if (charge<0) fTracks[eventNr*2+1].Add(particle);
 362   }
 363 }
 364
 365 //________________________________________________________________
 366 void AliDielectron::PairPreFilter(Int_t arr1, Int_t arr2, TObjArray &arrTracks1, TObjArray &arrTracks2)
 367 {
 368   //
 369   // Prefilter tracks from pairs
 370   // Needed for datlitz rejections
 371   // remove all tracks from the Single track arrays that pass the cuts in this filter
 372   //
 373   Int_t pairIndex=GetPairIndex(arr1,arr2);
 374
 375   Int_t ntrack1=arrTracks1.GetEntriesFast();
 376   Int_t ntrack2=arrTracks2.GetEntriesFast();
 377
 378   AliDielectronPair candidate;
 379
 380   UInt_t selectedMask=(1<<fPairPreFilter.GetCuts()->GetEntries())-1;
 381
 382   for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1){
 383     Int_t end=ntrack2;
 384     if (arr1==arr2) end=itrack1;
 385     Bool_t accepted=kFALSE;
 386     for (Int_t itrack2=0; itrack2<end; ++itrack2){
 387       AliVTrack *track1=static_cast<AliVTrack*>(arrTracks1.UncheckedAt(itrack1));
 388       AliVTrack *track2=static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack2));
 389       if (!track1 || !track2) continue;
 390       //create the pair
 391       candidate.SetTracks(track1, fPdgLeg1,
 392                           track2, fPdgLeg2);
 393       candidate.SetType(pairIndex);
 394       candidate.SetLabel(AliDielectronMC::Instance()->GetLabelMotherWithPdg(&candidate,fPdgMother));
 395
 396       //pair cuts
 397       UInt_t cutMask=fPairPreFilter.IsSelected(&candidate);
 398
 399       //apply cut
 400       if (cutMask!=selectedMask) continue;
 401       accepted=kTRUE;
 402       //remove the tracks from the Track arrays
 403       arrTracks2.AddAt(0x0,itrack2);
 404       //in case of like sign remove the track from both arrays!
 405       if (arr1==arr2) arrTracks1.AddAt(0x0, itrack2);
 406     }
 407     if ( accepted ) arrTracks1.AddAt(0x0,itrack1);
 408   }
 409   //compress the track arrays
 410   arrTracks1.Compress();
 411   arrTracks2.Compress();
 412
 413   //apply leg cuts after the pre filter
 414   if ( fPairPreFilterLegs.GetCuts()->GetEntries()>0 ) {
 415     selectedMask=(1<<fPairPreFilterLegs.GetCuts()->GetEntries())-1;
 416     //loop over tracks from array 1
 417     for (Int_t itrack=0; itrack<arrTracks1.GetEntriesFast();++itrack){
 418       //test cuts
 419       UInt_t cutMask=fPairPreFilterLegs.IsSelected(arrTracks1.UncheckedAt(itrack));
 420       //apply cut
 421       if (cutMask!=selectedMask) arrTracks1.AddAt(0x0,itrack);;
 422     }
 423     arrTracks1.Compress();
 424
 425     //in case of like sign don't loop over second array
 426     if (arr1==arr2) {
 427       arrTracks2=arrTracks1;
 428     } else {
 429
 430       //loop over tracks from array 2
 431       for (Int_t itrack=0; itrack<arrTracks2.GetEntriesFast();++itrack){
 432       //test cuts
 433         UInt_t cutMask=fPairPreFilterLegs.IsSelected(arrTracks2.UncheckedAt(itrack));
 434       //apply cut
 435         if (cutMask!=selectedMask) arrTracks2.AddAt(0x0,itrack);
 436       }
 437       arrTracks2.Compress();
 438
 439     }
 440   }
 441 }
 442
 443 //________________________________________________________________
 444 void AliDielectron::FillPairArrays(Int_t arr1, Int_t arr2)
 445 {
 446   //
 447   // select pairs and fill pair candidate arrays
 448   //
 449
 450   TObjArray arrTracks1=fTracks[arr1];
 451   TObjArray arrTracks2=fTracks[arr2];
 452
 453   //process pre filter if set
 454   if ( fPairPreFilter.GetCuts()->GetEntries()>0 ) PairPreFilter(arr1, arr2, arrTracks1, arrTracks2);
 455
 456   Int_t pairIndex=GetPairIndex(arr1,arr2);
 457
 458   Int_t ntrack1=arrTracks1.GetEntriesFast();
 459   Int_t ntrack2=arrTracks2.GetEntriesFast();
 460
 461   AliDielectronPair *candidate=new AliDielectronPair;
 462
 463   UInt_t selectedMask=(1<<fPairFilter.GetCuts()->GetEntries())-1;
 464
 465   for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1){
 466     Int_t end=ntrack2;
 467     if (arr1==arr2) end=itrack1;
 468     for (Int_t itrack2=0; itrack2<end; ++itrack2){
 469       //create the pair
 470       candidate->SetTracks(static_cast<AliVTrack*>(arrTracks1.UncheckedAt(itrack1)), fPdgLeg1,
 471                            static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack2)), fPdgLeg2);
 472       candidate->SetType(pairIndex);
 473       candidate->SetLabel(AliDielectronMC::Instance()->GetLabelMotherWithPdg(candidate,fPdgMother));
 474
 475       //pair cuts
 476       UInt_t cutMask=fPairFilter.IsSelected(candidate);
 477
 478       //CF manager for the pair
 479       if (fCfManagerPair) fCfManagerPair->Fill(cutMask,candidate);
 480
 481       //apply cut
 482       if (cutMask!=selectedMask) continue;
 483
 484       //add the candidate to the candidate array
 485       PairArray(pairIndex)->Add(candidate);
 486       //get a new candidate
 487       candidate=new AliDielectronPair;
 488     }
 489   }
 490   //delete the surplus candidate
 491   delete candidate;
 492 }
 493
 494 //________________________________________________________________
 495 void AliDielectron::FillPairArrayTR()
 496 {
 497   //
 498   // select pairs and fill pair candidate arrays
 499   //
 500   UInt_t selectedMask=(1<<fPairFilter.GetCuts()->GetEntries())-1;
 501
 502   while ( fTrackRotator->NextCombination() ){
 503     AliDielectronPair candidate;
 504     candidate.SetTracks(fTrackRotator->GetTrackP(), fPdgLeg1, fTrackRotator->GetTrackN(), fPdgLeg2);
 505     candidate.SetType(10);
 506
 507     //pair cuts
 508     UInt_t cutMask=fPairFilter.IsSelected(&candidate);
 509
 510     //CF manager for the pair
 511     if (fCfManagerPair) fCfManagerPair->Fill(cutMask,&candidate);
 512
 513     //apply cut
 514     if (cutMask==selectedMask&&fHistos) FillHistogramsPair(&candidate);
 515   }
 516 }
 517
 518 //________________________________________________________________
 519 void AliDielectron::FillDebugTree()
 520 {
 521   //
 522   // Fill Histogram information for tracks and pairs
 523   //
 524
 525   //Fill Debug tree
 526   for (Int_t i=0; i<10; ++i){
 527     Int_t ntracks=PairArray(i)->GetEntriesFast();
 528     for (Int_t ipair=0; ipair<ntracks; ++ipair){
 529       fDebugTree->Fill(static_cast<AliDielectronPair*>(PairArray(i)->UncheckedAt(ipair)));
 530     }
 531   }
 532 }
 533
 534 //________________________________________________________________
 535 void AliDielectron::SaveDebugTree()
 536 {
 537   //
 538   // delete the debug tree, this will also write the tree
 539   //
 540   if (fDebugTree) fDebugTree->DeleteStreamer();
 541 }
 542