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   fPairFilter("PairFilter"),
  92   fPdgMother(443),
  93   fHistos(0x0),
  94   fPairCandidates(new TObjArray(10)),
  95   fCfManagerPair(0x0),
  96   fDebugTree(0x0)
  97 {
  98   //
  99   // Default constructor
 100   //
 101
 102 }
 103
 104 //________________________________________________________________
 105 AliDielectron::AliDielectron(const char* name, const char* title) :
 106   TNamed(name,title),
 107   fEventFilter("EventFilter"),
 108   fTrackFilter("TrackFilter"),
 109   fPairFilter("PairFilter"),
 110   fPdgMother(443),
 111   fHistos(0x0),
 112   fPairCandidates(new TObjArray(10)),
 113   fCfManagerPair(0x0),
 114   fDebugTree(0x0)
 115 {
 116   //
 117   // Named constructor
 118   //
 119
 120 }
 121
 122 //________________________________________________________________
 123 AliDielectron::~AliDielectron()
 124 {
 125   //
 126   // Default destructor
 127   //
 128   if (fHistos) delete fHistos;
 129   if (fPairCandidates) delete fPairCandidates;
 130   if (fDebugTree) delete fDebugTree;
 131 }
 132
 133 //________________________________________________________________
 134 void AliDielectron::Init()
 135 {
 136   //
 137   // Initialise objects
 138   //
 139   if (fCfManagerPair) fCfManagerPair->InitialiseContainer(fPairFilter);
 140
 141 }
 142
 143 //________________________________________________________________
 144 void AliDielectron::Process(AliVEvent *ev1, AliVEvent *ev2)
 145 {
 146   //
 147   // Process the events
 148   //
 149
 150   AliDielectronVarManager::SetEvent(ev1);
 151
 152   //in case we have MC load the MC event and process the MC particles
 153   if (AliDielectronMC::Instance()->HasMC()) {
 154     AliDielectronMC::Instance()->ConnectMCEvent();
 155     ProcessMC();
 156   }
 157
 158   //if candidate array doesn't exist, create it
 159   if (!fPairCandidates->UncheckedAt(0)) {
 160     InitPairCandidateArrays();
 161   } else {
 162     ClearArrays();
 163   }
 164
 165   //mask used to require that all cuts are fulfilled
 166   UInt_t selectedMask=(1<<fEventFilter.GetCuts()->GetEntries())-1;
 167
 168   //apply event cuts
 169     if ((ev1&&fEventFilter.IsSelected(ev1)!=selectedMask) ||
 170         (ev2&&fEventFilter.IsSelected(ev2)!=selectedMask)) return;
 171
 172   AliDielectronVarManager::SetEvent(ev1);
 173
 174   //fill track arrays for the first event
 175   if (ev1) FillTrackArrays(ev1);
 176
 177   //fill track arrays for the second event
 178   if (ev2) FillTrackArrays(ev2,1);
 179
 180   // create pairs and fill pair candidate arrays
 181   for (Int_t itrackArr1=0; itrackArr1<4; ++itrackArr1){
 182     for (Int_t itrackArr2=itrackArr1; itrackArr2<4; ++itrackArr2){
 183       FillPairArrays(itrackArr1, itrackArr2);
 184     }
 185   }
 186
 187   //in case there is a histogram manager, fill the QA histograms
 188   if (fHistos) FillHistograms(ev1);
 189
 190   //fill debug tree if a manager is attached
 191   if (fDebugTree) FillDebugTree();
 192 }
 193
 194 //________________________________________________________________
 195 void AliDielectron::ProcessMC()
 196 {
 197   //
 198   // Process the MC data
 199   //
 200
 201   //loop over all MC data and Fill the CF container if it exist
 202   if (!fCfManagerPair) return;
 203   fCfManagerPair->SetPdgMother(fPdgMother);
 204   AliDielectronMC *dieMC=AliDielectronMC::Instance();
 205   for (Int_t ipart=0; ipart<dieMC->GetNMCTracks();++ipart){
 206     //TODO: MC truth cut properly!!!
 207     AliVParticle *mcPart=dieMC->GetMCTrackFromMCEvent(ipart);
 208     if (!dieMC->IsMCMotherToEE(mcPart, fPdgMother)) continue;
 209     fCfManagerPair->FillMC(mcPart);
 210   }
 211 }
 212
 213 //________________________________________________________________
 214 void AliDielectron::FillHistograms(const AliVEvent *ev)
 215 {
 216   //
 217   // Fill Histogram information for tracks and pairs
 218   //
 219
 220   TString  className;
 221   Double_t values[AliDielectronVarManager::kNMaxValues];
 222   //Fill event information
 223   AliDielectronVarManager::Fill(ev, values);
 224   fHistos->FillClass("Event", AliDielectronVarManager::kNMaxValues, values);
 225
 226   //Fill track information, separately for the track array candidates
 227   for (Int_t i=0; i<4; ++i){
 228     className.Form("Track_%s",fgkTrackClassNames[i]);
 229     Int_t ntracks=fTracks[i].GetEntriesFast();
 230     for (Int_t itrack=0; itrack<ntracks; ++itrack){
 231       AliDielectronVarManager::Fill(fTracks[i].UncheckedAt(itrack), values);
 232       fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
 233     }
 234   }
 235
 236   //Fill Pair information, separately for all pair candidate arrays
 237   for (Int_t i=0; i<10; ++i){
 238     className.Form("Pair_%s",fgkPairClassNames[i]);
 239     Int_t ntracks=PairArray(i)->GetEntriesFast();
 240     for (Int_t ipair=0; ipair<ntracks; ++ipair){
 241       AliDielectronVarManager::Fill(PairArray(i)->UncheckedAt(ipair), values);
 242       fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
 243     }
 244   }
 245
 246 }
 247
 248 //________________________________________________________________
 249 void AliDielectron::FillTrackArrays(AliVEvent * const ev, Int_t eventNr)
 250 {
 251   //
 252   // select tracks and fill track candidate arrays
 253   // eventNr = 0: First  event, use track arrays 0 and 1
 254   // eventNr = 1: Second event, use track arrays 2 and 3
 255   //
 256
 257   Int_t ntracks=ev->GetNumberOfTracks();
 258   UInt_t selectedMask=(1<<fTrackFilter.GetCuts()->GetEntries())-1;
 259   for (Int_t itrack=0; itrack<ntracks; ++itrack){
 260     //get particle
 261     AliVParticle *particle=ev->GetTrack(itrack);
 262     //TODO: temporary solution, perhaps think about a better implementation
 263     //      This is needed to use AliESDpidCuts, which relies on the ESD event
 264     //      is set as a AliESDtrack attribute... somehow ugly!
 265     if (ev->IsA()==AliESDEvent::Class()){
 266       AliESDtrack *track=static_cast<AliESDtrack*>(particle);
 267       track->SetESDEvent(static_cast<AliESDEvent*>(ev)); //only in trunk...
 268     }
 269
 270     //apply track cuts
 271     if (fTrackFilter.IsSelected(particle)!=selectedMask) continue;
 272
 273     //fill selected particle into the corresponding track arrays
 274     Short_t charge=particle->Charge();
 275     if (charge>0)      fTracks[eventNr*2].Add(particle);
 276     else if (charge<0) fTracks[eventNr*2+1].Add(particle);
 277   }
 278 }
 279
 280 //________________________________________________________________
 281 void AliDielectron::FillPairArrays(Int_t arr1, Int_t arr2) {
 282   //
 283   // select pairs and fill pair candidate arrays
 284   //
 285   Int_t pairIndex=GetPairIndex(arr1,arr2);
 286
 287   Int_t ntrack1=fTracks[arr1].GetEntriesFast();
 288   Int_t ntrack2=fTracks[arr2].GetEntriesFast();
 289
 290   AliDielectronPair *candidate=new AliDielectronPair;
 291
 292   UInt_t selectedMask=(1<<fPairFilter.GetCuts()->GetEntries())-1;
 293
 294   for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1){
 295     Int_t end=ntrack2;
 296     if (arr1==arr2) end=itrack1;
 297     for (Int_t itrack2=0; itrack2<end; ++itrack2){
 298       //create the pair TODO: change hardcoded PID of tracks
 299       candidate->SetTracks(static_cast<AliVTrack*>(fTracks[arr1].UncheckedAt(itrack1)), 11,
 300                            static_cast<AliVTrack*>(fTracks[arr2].UncheckedAt(itrack2)), 11);
 301       candidate->SetType(pairIndex);
 302       candidate->SetLabel(AliDielectronMC::Instance()->GetLabelMotherWithPdg(candidate,fPdgMother));
 303
 304       //pair cuts
 305       UInt_t cutMask=fPairFilter.IsSelected(candidate);
 306
 307       //CF manager for the pair
 308       if (fCfManagerPair) fCfManagerPair->Fill(cutMask,candidate);
 309
 310       //apply cut
 311       if (cutMask!=selectedMask) continue;
 312
 313       //add the candidate to the candidate array
 314       PairArray(pairIndex)->Add(candidate);
 315       //get a new candidate
 316       candidate=new AliDielectronPair;
 317     }
 318   }
 319   //delete the surplus candidate
 320   delete candidate;
 321 }
 322
 323 //________________________________________________________________
 324 void AliDielectron::FillDebugTree()
 325 {
 326   //
 327   // Fill Histogram information for tracks and pairs
 328   //
 329
 330   //Fill Debug tree
 331   for (Int_t i=0; i<10; ++i){
 332     Int_t ntracks=PairArray(i)->GetEntriesFast();
 333     for (Int_t ipair=0; ipair<ntracks; ++ipair){
 334       fDebugTree->Fill(static_cast<AliDielectronPair*>(PairArray(i)->UncheckedAt(ipair)));
 335     }
 336   }
 337 }
 338
 339 //________________________________________________________________
 340 void AliDielectron::SaveDebugTree()
 341 {
 342   //
 343   // delete the debug tree, this will also write the tree
 344   //
 345   if (fDebugTree) fDebugTree->DeleteStreamer();
 346 }
 347