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