PWG2/RESONANCES/AliRsnAnalysisEffSE.cxx

   1 //
   2 // Class AliRsnAnalysisEffSE
   3 //
   4 // TODO
   5 // TODO
   6 //
   7 // authors: Alberto Pulvirenti (alberto.pulvirenti@ct.infn.it)
   8 //          Martin Vala (martin.vala@cern.ch)
   9 //
  10 #include <Riostream.h>
  11 #include "AliStack.h"
  12 #include "AliESDEvent.h"
  13 #include "AliMCEvent.h"
  14
  15 #include "AliCFContainer.h"
  16 #include "AliRsnCutMgr.h"
  17 #include "AliRsnFunctionAxis.h"
  18 #include "AliRsnAnalysisEffSE.h"
  19 #include "AliRsnPairDef.h"
  20 #include "AliRsnCutSet.h"
  21
  22 ClassImp(AliRsnAnalysisEffSE)
  23
  24 //_____________________________________________________________________________
  25 AliRsnAnalysisEffSE::AliRsnAnalysisEffSE(const char *name) :
  26   AliRsnVAnalysisTaskSE(name),
  27   fEventCuts(0x0),
  28   fStepListMC(0),
  29   fStepListESD(0),
  30   fAxisList(0),
  31   fPairDefList(0),
  32   fContainerList(0x0),
  33   fVar(0),
  34   fPair()
  35 {
  36 //
  37 // Default constructor.
  38 //
  39
  40   AliDebug(AliLog::kDebug+2,"<-");
  41
  42   DefineOutput(2, TList::Class());
  43
  44   AliDebug(AliLog::kDebug+2,"->");
  45 }
  46
  47 //_____________________________________________________________________________
  48 AliRsnAnalysisEffSE::AliRsnAnalysisEffSE(const AliRsnAnalysisEffSE& copy) :
  49   AliRsnVAnalysisTaskSE(copy),
  50   fEventCuts(copy.fEventCuts),
  51   fStepListMC(copy.fStepListMC),
  52   fStepListESD(copy.fStepListESD),
  53   fAxisList(copy.fAxisList),
  54   fPairDefList(copy.fPairDefList),
  55   fContainerList(copy.fContainerList),
  56   fVar(0),
  57   fPair()
  58 {
  59 //
  60 // Copy constrtuctor
  61 //
  62 }
  63
  64 //_____________________________________________________________________________
  65 void AliRsnAnalysisEffSE::RsnUserCreateOutputObjects()
  66 {
  67 //
  68 // Creation of output objects.
  69 // These are created through the utility methods in the analysis manager,
  70 // which produces a list of histograms for each specified set of pairs.
  71 // Each of these lists is added to the main list of this task.
  72 //
  73
  74   AliDebug(AliLog::kDebug+2,"<-");
  75
  76   // get number of steps and axes
  77   Int_t iaxis;
  78   Int_t nAxes  = fAxisList.GetEntries();
  79   Int_t nSteps = (Int_t)fStepListMC.GetEntries() + (Int_t)fStepListESD.GetEntries();
  80
  81   if (!nSteps) {
  82     AliError("No steps defined");
  83     return;
  84   }
  85   if (!nAxes) {
  86     AliError("No axes defined");
  87     return;
  88   }
  89
  90   // initialize variable list
  91   fVar.Set(nAxes);
  92
  93   // retrieve number of bins for each axis
  94   Int_t   *nBins     = new Int_t[nAxes];
  95   TArrayD *binLimits = new TArrayD[nAxes];
  96   for (iaxis = 0; iaxis < nAxes; iaxis++) {
  97     AliRsnFunctionAxis *fcnAxis = (AliRsnFunctionAxis*)fAxisList.At(iaxis);
  98     binLimits[iaxis] = fcnAxis->GetArray();
  99     nBins[iaxis] = binLimits[iaxis].GetSize() - 1;
 100   }
 101
 102   // create ouput list of containers
 103   fContainerList = new TList;
 104   fContainerList->SetOwner();
 105   fContainerList->SetName(Form("%s_containers", GetName()));
 106
 107   // initialize output list
 108   OpenFile(2);
 109   fOutList[1] = new TList();
 110   fOutList[1]->SetOwner();
 111
 112   // create the containers
 113   Int_t i, nDef = (Int_t)fPairDefList.GetEntries();
 114   for (i = 0; i < nDef; i++) {
 115     AliRsnPairDef *def = (AliRsnPairDef*)fPairDefList[i];
 116     AliCFContainer *cont = new AliCFContainer(Form("%s", def->GetPairName().Data()), "", nSteps, nAxes, nBins);
 117     // set the bin limits for each axis
 118     for (iaxis = 0; iaxis < nAxes; iaxis++) cont->SetBinLimits(iaxis, binLimits[iaxis].GetArray());
 119     // add the container to output list
 120     fContainerList->Add(cont);
 121   }
 122
 123   fOutList[1]->Add(fContainerList);
 124
 125   AliDebug(AliLog::kDebug+2,"->");
 126 }
 127
 128 //_____________________________________________________________________________
 129 void AliRsnAnalysisEffSE::RsnUserExec(Option_t*)
 130 {
 131 //
 132 // Execution of the analysis task.
 133 // Recovers the input event and processes it with all included pair objects.
 134 // In this case, we NEED to have ESD and MC, otherwise we cannod do anything.
 135 //ULTIMO UNDO
 136 /*
 137   AliRsnDaughter trk;
 138   for (Int_t i = 0; i <= AliPID::kSPECIES; i++)
 139   {
 140     cout << AliPID::ParticleName((AliPID::EParticleType)i) << ": " << endl;
 141     for (Int_t m = 0; m < AliRsnDaughter::kMethods; m++)
 142     {
 143       cout << "-- method: " << AliRsnDaughter::MethodName((AliRsnDaughter::EPIDMethod)m) << endl;
 144       Char_t   sign[2] = {'+', '-'};
 145       for (Int_t s = 0; s < 2; s++)
 146       {
 147         TArrayI *a = fRsnPIDIndex.GetTracksArray((AliRsnDaughter::EPIDMethod)m, sign[s], (AliPID::EParticleType)i);
 148         Int_t n = a->GetSize();
 149         for (Int_t j = 0; j < n; j++)
 150         {
 151           Int_t k = a->At(j);
 152           cout << "-- -- track " << Form("%4d ", k) << ": ";
 153           fRsnEvent.SetDaughter(trk, k);
 154           cout << "charge = " << (trk.IsPos() ? "+ " : (trk.IsNeg() ? "- " : "0 "));
 155           cout << "truePID = " << Form("%10s ", AliPID::ParticleName(trk.PerfectPID()));
 156           cout << "realPID = " << Form("%10s ", AliPID::ParticleName(trk.RealisticPID()));
 157           cout << endl;
 158           cout << "-- -- -- weights (computed): ";
 159           for (Int_t q = 0; q < AliPID::kSPECIES; q++)
 160             cout << Form("%15.12f", trk.ComputedWeights()[q]) << ' ';
 161           cout << endl;
 162           cout << "-- -- -- weights (original): ";
 163           for (Int_t q = 0; q < AliPID::kSPECIES; q++)
 164             cout << Form("%15.12f", trk.GetRef()->PID()[q]) << ' ';
 165           cout << endl;
 166         }
 167       }
 168     }
 169   } return;
 170   */
 171
 172   AliDebug(AliLog::kDebug+2,"<-");
 173   if (!fESDEvent || !fMCEvent) {
 174     AliError("This task can process only ESD + MC events");
 175     return;
 176   }
 177   fRsnEvent.SetRef(fESDEvent, fMCEvent);
 178
 179   // if general event cuts are added to the task (recommended)
 180   // they are checked here on the RSN event interface and,
 181   // if the event does not pass them, it is skipped and ProcessInfo
 182   // is updated accordingly
 183   if (fEventCuts) {
 184     if (!fEventCuts->IsSelected(AliRsnCut::kEvent, &fRsnEvent)) {
 185       fTaskInfo.SetEventUsed(kFALSE);
 186       return;
 187     }
 188   }
 189
 190   // if cuts are passed or not cuts were defined,
 191   // update the task info before processing the event
 192   fTaskInfo.SetEventUsed(kTRUE);
 193
 194   // process first MC steps and then ESD steps
 195   AliRsnPairDef *pairDef = 0;
 196   TObjArrayIter iter(&fPairDefList);
 197   while ( (pairDef = (AliRsnPairDef*)iter.Next()) )
 198   {
 199     ProcessEventMC(pairDef);
 200     ProcessEventESD(pairDef);
 201   }
 202
 203   // Post the data
 204   PostData(2, fOutList[1]);
 205
 206   AliDebug(AliLog::kDebug+2,"->");
 207 }
 208
 209 //_____________________________________________________________________________
 210 void AliRsnAnalysisEffSE::ProcessEventMC(AliRsnPairDef *pairDef)
 211 {
 212 //
 213 // Process current event with the definitions of the specified step in MC list
 214 // and store results in the container slot defined in second argument
 215 //
 216
 217   AliStack      *stack = fMCEvent->Stack();
 218   AliMCParticle *mother, *daughter;
 219
 220   if (!pairDef) return;
 221   AliCFContainer *cont = (AliCFContainer*)fContainerList->FindObject(pairDef->GetPairName().Data());
 222   if (!cont) return;
 223
 224   // other utility variables
 225   Int_t i[2], j, ipart;
 226
 227   // in this case, we first take the resonance from MC
 228   // and then we find its daughters and compute cuts on them
 229   for (ipart = 0; ipart < stack->GetNprimary(); ipart++) {
 230     mother = (AliMCParticle*) fMCEvent->GetTrack(ipart);
 231     if (mother->Particle()->GetNDaughters() != 2) continue;
 232
 233     i[0] = mother->Particle()->GetFirstDaughter();
 234     i[1] = mother->Particle()->GetLastDaughter();
 235
 236     for (j = 0; j < 2; j++) {
 237       daughter = (AliMCParticle*) fMCEvent->GetTrack(i[j]);
 238       fDaughter[j].SetRef(daughter);
 239       fDaughter[j].SetParticle(daughter->Particle());
 240       fDaughter[j].FindMotherPDG(stack);
 241     }
 242
 243     if (fDaughter[0].ChargeC() != pairDef->GetCharge(0)) continue;
 244     if (fDaughter[1].ChargeC() != pairDef->GetCharge(1)) continue;
 245     if (fDaughter[0].PerfectPID() != pairDef->GetType(0)) continue;
 246     if (fDaughter[1].PerfectPID() != pairDef->GetType(1)) continue;
 247
 248     fPair.SetPair(&fDaughter[0], &fDaughter[1]);
 249
 250     // create pair
 251     FillContainer(cont, &fStepListMC, pairDef, 0);
 252   }
 253 }
 254
 255 //_____________________________________________________________________________
 256 void AliRsnAnalysisEffSE::ProcessEventESD(AliRsnPairDef *pairDef)
 257 {
 258 //
 259 // Process current event with the definitions of the specified step in ESD list
 260 // and store results in the container slot defined in second argument
 261 //
 262
 263   Int_t i0, i1, first = (Int_t)fStepListMC.GetEntries();
 264
 265   if (!pairDef) return;
 266   AliCFContainer *cont = (AliCFContainer*)fContainerList->FindObject(pairDef->GetPairName().Data());
 267   if (!cont) return;
 268
 269   // get arrays of all charged tracks
 270   TArrayI *a0 = fRsnPIDIndex.GetTracksArray(AliRsnDaughter::kPerfect, pairDef->GetCharge(0), pairDef->GetType(0));
 271   TArrayI *a1 = fRsnPIDIndex.GetTracksArray(AliRsnDaughter::kPerfect, pairDef->GetCharge(1), pairDef->GetType(1));
 272
 273   // external loop on tracks
 274   for (i0 = 0; i0 < a0->GetSize(); i0++) {
 275     // connect interface
 276     fRsnEvent.SetDaughter(fDaughter[0], a0->At(i0));
 277     if (!fDaughter[0].IsOK()) continue;
 278     fDaughter[0].SetRequiredPID(pairDef->GetType(0));
 279
 280     // internal loop on tracks
 281     for (i1 = 0; i1 < a1->GetSize(); i1++) {
 282       // connect interface
 283       fRsnEvent.SetDaughter(fDaughter[1], a1->At(i1));
 284       if (!fDaughter[1].IsOK()) continue;
 285       fDaughter[1].SetRequiredPID(pairDef->GetType(1));
 286       // build pair
 287       fPair.SetPair(&fDaughter[0], &fDaughter[1]);
 288       if (TMath::Abs(fPair.CommonMother()) != pairDef->GetMotherPDG()) continue;
 289       // fill containers
 290       FillContainer(cont, &fStepListESD, pairDef, first);
 291     }
 292   }
 293 }
 294
 295 //_____________________________________________________________________________
 296 void AliRsnAnalysisEffSE::FillContainer(AliCFContainer *cont, const TObjArray *stepList, AliRsnPairDef *pd, Int_t firstOutStep)
 297 {
 298 //
 299 // Fill the containers
 300 //
 301
 302   Int_t iaxis, nAxes  = fAxisList.GetEntries();
 303   Int_t istep, nSteps = stepList->GetEntries();
 304
 305   // compute values for all axes
 306   for (iaxis = 0; iaxis < nAxes; iaxis++) {
 307     AliRsnFunctionAxis *fcnAxis = (AliRsnFunctionAxis*)fAxisList.At(iaxis);
 308     switch (fcnAxis->GetAxisObject()) {
 309     case AliRsnFunctionAxis::kPair:
 310       fVar[iaxis] = (Double_t)fcnAxis->Eval(&fPair, pd);
 311       break;
 312     case AliRsnFunctionAxis::kEvent:
 313       fVar[iaxis] = (Double_t)fcnAxis->Eval(&fRsnEvent);
 314       break;
 315     default:
 316       fVar[iaxis] = 0.0;
 317     }
 318   }
 319
 320   // fill all steps
 321   for (istep = 0; istep < nSteps; istep++) {
 322     AliRsnCutMgr *cutMgr = (AliRsnCutMgr*)stepList->At(istep);
 323     if (!cutMgr->IsSelected(AliRsnCut::kParticle, &fDaughter[0])) break;
 324     if (!cutMgr->IsSelected(AliRsnCut::kParticle, &fDaughter[1])) break;
 325     if (!cutMgr->IsSelected(AliRsnCut::kPair,     &fPair)) break;
 326     cont->Fill(fVar.GetArray(), istep + firstOutStep);
 327   }
 328 }
 329
 330 //_____________________________________________________________________________
 331 void AliRsnAnalysisEffSE::RsnTerminate(Option_t*)
 332 {
 333 //
 334 // Termination.
 335 // Could be added some monitor histograms here.
 336 //
 337
 338   AliDebug(AliLog::kDebug+2,"<-");
 339   AliDebug(AliLog::kDebug+2,"->");
 340 }
 341
 342 //_____________________________________________________________________________
 343 void AliRsnAnalysisEffSE::AddPairDef(AliRsnPairDef* pairDef)
 344 {
 345 //
 346 //  Adds pair definition
 347 //
 348   fPairDefList.AddLast(pairDef);
 349 }