PWG2/RESONANCES/AliRsnReader.cxx

   1 //
   2 // Class AliRsnReader
   3 //
   4 // This is the universal converter from any kind of source event
   5 // (i.e. ESD, standard AOD, MC) into the internal non-standard
   6 // AOD format used by RSN package.
   7 // ---
   8 // This class reads all tracks in an input event and converts them
   9 // into AliRsnDaughters, and computes preliminarily the PID probabilities
  10 // by doing the Bayesian combination of a set of assigned prior probabilities
  11 // with the PID weights defined in each track.
  12 // ---
  13 // When filling the output event (AliRsnEvent), some arrays of indexes
  14 // are created in order to organize tracks according to their PID and charge,
  15 // which will then be used in further analysis steps.
  16 //
  17 // author: A. Pulvirenti (alberto.pulvirenti@ct.infn.it)
  18 //
  19
  20 #include <TString.h>
  21
  22 #include "AliLog.h"
  23
  24 #include "AliVEvent.h"
  25
  26 #include "AliESDEvent.h"
  27 #include "AliESDtrack.h"
  28 #include "AliESDVertex.h"
  29
  30 #include "AliAODEvent.h"
  31 #include "AliAODTrack.h"
  32 #include "AliAODVertex.h"
  33
  34 #include "AliStack.h"
  35 #include "AliMCEvent.h"
  36 #include "AliMCParticle.h"
  37 #include "AliGenEventHeader.h"
  38
  39 #include "AliRsnMCInfo.h"
  40 #include "AliRsnDaughter.h"
  41 #include "AliRsnEvent.h"
  42 #include "AliRsnPIDWeightsMgr.h"
  43
  44 #include "AliRsnReader.h"
  45
  46 ClassImp(AliRsnReader)
  47
  48 //_____________________________________________________________________________
  49 AliRsnReader::AliRsnReader(AliRsnPIDWeightsMgr *mgr) :
  50     TNamed("RsnReader", ""),
  51     fCheckSplit(kFALSE),
  52     fRejectFakes(kFALSE),
  53     fWeightsMgr(mgr),
  54     fCurrentPIDtype(AliRsnDaughter::kEsd),
  55     fPIDDivValue(0.0),
  56     fITSClusters(0),
  57     fTPCClusters(0),
  58     fTRDClusters(0)
  59 {
  60 //
  61 // Constructor.
  62 //
  63 }
  64
  65 //_____________________________________________________________________________
  66 Bool_t AliRsnReader::Fill
  67 (AliRsnEvent *rsn, AliVEvent *event, AliMCEvent *mc, Bool_t useTPCOnly)
  68 {
  69 //
  70 // According to the class type of event and the selected source
  71 // recalls one of the private reading methods to fill the RsnEvent
  72 // passed by reference as first argument.
  73 //
  74
  75     Bool_t success = kFALSE;
  76     TString str(event->ClassName());
  77
  78     if (!str.CompareTo("AliESDEvent")) {
  79         AliDebug(1, "Reading an ESD event");
  80         success = FillFromESD(rsn, (AliESDEvent*)event, mc, useTPCOnly);
  81     }
  82     else if (!str.CompareTo("AliAODEvent")) {
  83         AliDebug(1, "Reading an AOD event");
  84         success = FillFromAOD(rsn, (AliAODEvent*)event, mc);
  85     }
  86     else if (!str.CompareTo("AliMCEvent")) {
  87         AliDebug(1, "Reading an MC event");
  88         success = FillFromMC(rsn, (AliMCEvent*)event);
  89     }
  90     else {
  91         AliError(Form("ClassName '%s' not recognized as possible source data: aborting.", str.Data()));
  92         return kFALSE;
  93     }
  94
  95     // sort tracks w.r. to Pt (from largest to smallest)
  96     rsn->SortTracks();
  97     return success;
  98 }
  99
 100 //_____________________________________________________________________________
 101 Bool_t AliRsnReader::FillFromESD
 102 (AliRsnEvent *rsn, AliESDEvent *esd, AliMCEvent *mc, Bool_t useTPCOnly)
 103 {
 104 //
 105 // Filler from an ESD event.
 106 // Stores all tracks (if a filter is defined, it will store
 107 // only the ones which survive the cuts).
 108 // If a reference MC event is provided, it is used to store
 109 // the MC informations for each track (true PDG code,
 110 // GEANT label of mother, PDG code of mother, if any).
 111 // When this is used, the 'source' flag of the output
 112 // AliRsnEvent object will be set to 'kESD'.
 113 //
 114
 115     // retrieve stack (if possible)
 116     AliStack *stack = 0x0;
 117     if (mc) stack = mc->Stack();
 118
 119     // get number of tracks
 120     Int_t ntracks = esd->GetNumberOfTracks();
 121
 122     // if required with the flag, scans the event
 123     // and searches all split tracks (= 2 tracks with the same label);
 124     // for each pair of split tracks, only the better (best chi2) is kept
 125     // and the other is rejected: this info is stored into a Boolean array
 126     Int_t i1, i2, lab1, lab2;
 127     Bool_t *accept = new Bool_t[ntracks];
 128     for (i1 = 0; i1 < ntracks; i1++) accept[i1] = kTRUE;
 129     if (fCheckSplit)
 130     {
 131         for (i1 = 0; i1 < ntracks; i1++)
 132         {
 133             AliESDtrack *trk1 = esd->GetTrack(i1);
 134             lab1 = TMath::Abs(trk1->GetLabel());
 135             for (i2 = i1+1; i2 < ntracks; i2++)
 136             {
 137                 AliESDtrack *trk2 = esd->GetTrack(i2);
 138                 lab2 = TMath::Abs(trk2->GetLabel());
 139                 // check if labels are equal
 140                 if (lab1 == lab2)
 141                 {
 142                     if (trk1->GetConstrainedChi2() < trk2->GetConstrainedChi2())
 143                     {
 144                         accept[i1] = kTRUE;
 145                         accept[i2] = kFALSE;
 146                     }
 147                     else
 148                     {
 149                         accept[i1] = kFALSE;
 150                         accept[i2] = kTRUE;
 151                     }
 152                 }
 153             }
 154         }
 155     }
 156
 157     // get primary vertex
 158     Double_t vertex[3];
 159     if (!useTPCOnly)
 160     {
 161         vertex[0] = esd->GetVertex()->GetXv();
 162         vertex[1] = esd->GetVertex()->GetYv();
 163         vertex[2] = esd->GetVertex()->GetZv();
 164     }
 165     else
 166     {
 167         vertex[0] = esd->GetPrimaryVertexTPC()->GetXv();
 168         vertex[1] = esd->GetPrimaryVertexTPC()->GetYv();
 169         vertex[2] = esd->GetPrimaryVertexTPC()->GetZv();
 170     }
 171     rsn->SetPrimaryVertex(vertex[0], vertex[1], vertex[2]);
 172
 173     // store tracks from ESD
 174     Int_t    i, index, label, labmum;
 175     Bool_t check;
 176     AliRsnDaughter temp;
 177     for (index = 0; index < ntracks; index++)
 178     {
 179         // skip track recognized as the worse one in a splitted pair
 180         if (!accept[index])
 181         {
 182             AliInfo(Form("Rejecting split track #%d in this event", index));
 183             continue;
 184         }
 185         // get and check ESD track
 186         AliESDtrack *esdTrack = esd->GetTrack(index);
 187         label = esdTrack->GetLabel();
 188         if (fRejectFakes && (label < 0)) continue;
 189         // copy ESD track data into RsnDaughter
 190         // if unsuccessful, this track is skipped
 191 //         temp.SetCurrentESDPID(fCurrentPIDtype);
 192
 193         if (fITSClusters>0)
 194             if (esdTrack->GetITSclusters(0) < fITSClusters) continue;
 195
 196         if (fTPCClusters>0)
 197             if (esdTrack->GetTPCclusters(0) < fTPCClusters) continue;
 198
 199         if (fTRDClusters>0)
 200             if (esdTrack->GetTRDclusters(0) < fTRDClusters) continue;
 201
 202         check = temp.Adopt(esdTrack,fCurrentPIDtype,fPIDDivValue, useTPCOnly);
 203         if (!check) continue;
 204         // if the AliRsnWeightsMgr object is initialized
 205         // this means that ESD PID weights are not used
 206         // and they are computed according to the Weight manager
 207         if (fWeightsMgr)
 208         {
 209             //AliInfo("Using customized weights");
 210             //AliInfo(Form("ESD weights = %f %f %f %f %f", temp.PID()[0], temp.PID()[1], temp.PID()[2], temp.PID()[3], temp.PID()[4], temp.PID()[5]));
 211             esdTrack->GetITSpid(fWeightsMgr->GetWeightArray(AliRsnPIDWeightsMgr::kITS));
 212             esdTrack->GetTPCpid(fWeightsMgr->GetWeightArray(AliRsnPIDWeightsMgr::kTPC));
 213             esdTrack->GetTRDpid(fWeightsMgr->GetWeightArray(AliRsnPIDWeightsMgr::kTRD));
 214             esdTrack->GetTOFpid(fWeightsMgr->GetWeightArray(AliRsnPIDWeightsMgr::kTOF));
 215             esdTrack->GetHMPIDpid(fWeightsMgr->GetWeightArray(AliRsnPIDWeightsMgr::kHMPID));
 216             for (i = 0; i < AliRsnPID::kSpecies; i++)
 217             {
 218                 temp.SetPIDWeight(i, fWeightsMgr->GetWeight((AliRsnPID::EType)i, temp.Pt()));
 219             }
 220             //AliInfo(Form("Used weights = %f %f %f %f %f", temp.PID()[0], temp.PID()[1], temp.PID()[2], temp.PID()[3], temp.PID()[4], temp.PID()[5]));
 221         }
 222         else
 223         {
 224             //AliInfo("Using standard ESD weights");
 225         }
 226
 227         // if stack is present, copy MC info
 228         if (stack)
 229         {
 230             TParticle *part = stack->Particle(TMath::Abs(label));
 231             if (part)
 232             {
 233                 temp.InitMCInfo(part);
 234                 labmum = part->GetFirstMother();
 235                 if (labmum >= 0)
 236                 {
 237                     TParticle *mum = stack->Particle(labmum);
 238                     temp.GetMCInfo()->SetMotherPDG(mum->GetPdgCode());
 239                 }
 240             }
 241         }
 242         // set index and label and add this object to the output container
 243         temp.SetIndex(index);
 244         temp.SetLabel(label);
 245         AliRsnDaughter *ptr = rsn->AddTrack(temp);
 246         // if problems occurred while storing, that pointer is NULL
 247         if (!ptr) AliWarning(Form("Failed storing track#%d", index));
 248     }
 249
 250     // compute total multiplicity
 251     if (rsn->GetMultiplicity() <= 0)
 252     {
 253         AliDebug(1, "Zero Multiplicity in this event");
 254         return kFALSE;
 255     }
 256
 257     return kTRUE;
 258 }
 259
 260 //_____________________________________________________________________________
 261 Bool_t AliRsnReader::FillFromAOD(AliRsnEvent *rsn, AliAODEvent *aod, AliMCEvent *mc)
 262 {
 263 //
 264 // Filler from an AOD event.
 265 // Stores all tracks (if a filter is defined, it will store
 266 // only the ones which survive the cuts).
 267 // If a reference MC event is provided, it is used to store
 268 // the MC informations for each track (true PDG code,
 269 // GEANT label of mother, PDG code of mother, if any).
 270 // When this is used, the 'source' flag of the output
 271 // AliRsnEvent object will be set to 'kAOD'.
 272 //
 273
 274     // retrieve stack (if possible)
 275     AliStack *stack = 0x0;
 276     if (mc) stack = mc->Stack();
 277
 278     // get number of tracks
 279     Int_t ntracks = aod->GetNTracks();
 280     if (!ntracks)
 281     {
 282         AliWarning("No tracks in this event");
 283         return kFALSE;
 284     }
 285
 286     // get primary vertex
 287     Double_t vertex[3];
 288     vertex[0] = aod->GetPrimaryVertex()->GetX();
 289     vertex[1] = aod->GetPrimaryVertex()->GetY();
 290     vertex[2] = aod->GetPrimaryVertex()->GetZ();
 291     rsn->SetPrimaryVertex(vertex[0], vertex[1], vertex[2]);
 292
 293     // store tracks from ESD
 294     Int_t    index, label, labmum;
 295     Bool_t check;
 296     AliAODTrack *aodTrack = 0;
 297     AliRsnDaughter temp;
 298     TObjArrayIter iter(aod->GetTracks());
 299     while ((aodTrack = (AliAODTrack*)iter.Next()))
 300     {
 301         // retrieve index
 302         index = aod->GetTracks()->IndexOf(aodTrack);
 303         label = aodTrack->GetLabel();
 304         if (fRejectFakes && (label < 0)) continue;
 305         // copy ESD track data into RsnDaughter
 306         // if unsuccessful, this track is skipped
 307         check = temp.Adopt(aodTrack);
 308         if (!check) continue;
 309         // if stack is present, copy MC info
 310         if (stack)
 311         {
 312             TParticle *part = stack->Particle(TMath::Abs(label));
 313             if (part)
 314             {
 315                 temp.InitMCInfo(part);
 316                 labmum = part->GetFirstMother();
 317                 if (labmum >= 0)
 318                 {
 319                     TParticle *mum = stack->Particle(labmum);
 320                     temp.GetMCInfo()->SetMotherPDG(mum->GetPdgCode());
 321                 }
 322             }
 323         }
 324         // set index and label and add this object to the output container
 325         temp.SetIndex(index);
 326         temp.SetLabel(label);
 327         AliRsnDaughter *ptr = rsn->AddTrack(temp);
 328         // if problems occurred while storin, that pointer is NULL
 329         if (!ptr) AliWarning(Form("Failed storing track#%d"));
 330     }
 331
 332     // compute total multiplicity
 333     if (rsn->GetMultiplicity() <= 0)
 334     {
 335         AliDebug(1, "Zero multiplicity in this event");
 336         return kFALSE;
 337     }
 338
 339     return kTRUE;
 340 }
 341
 342 //_____________________________________________________________________________
 343 Bool_t AliRsnReader::FillFromMC(AliRsnEvent *rsn, AliMCEvent *mc)
 344 {
 345 //
 346 // Filler from an ESD event.
 347 // Stores all tracks which generate at least one
 348 // TrackReference (a point in a sensitive volume).
 349 // In this case, the MC info is stored by default and
 350 // perfect particle identification is the unique available.
 351 // When this is used, the 'source' flag of the output
 352 // AliRsnEvent object will be set to 'kMC'.
 353 //
 354
 355     // get number of tracks
 356     Int_t ntracks = mc->GetNumberOfTracks();
 357     if (!ntracks)
 358     {
 359         AliWarning("No tracks in this event");
 360         return kFALSE;
 361     }
 362
 363     AliStack *stack = mc->Stack();
 364
 365     // get primary vertex
 366     TArrayF fvertex(3);
 367     Double_t vertex[3];
 368     mc->GenEventHeader()->PrimaryVertex(fvertex);
 369     vertex[0] = (Double_t)fvertex[0];
 370     vertex[1] = (Double_t)fvertex[1];
 371     vertex[2] = (Double_t)fvertex[2];
 372     rsn->SetPrimaryVertex(vertex[0], vertex[1], vertex[2]);
 373
 374     // store tracks from MC
 375     Int_t    index, labmum;
 376     Bool_t check;
 377     AliRsnDaughter temp;
 378     for (index = 0; index < ntracks; index++)
 379     {
 380         // get and check MC track
 381         AliMCParticle *mcTrack = mc->GetTrack(index);
 382         // if particle has no track references, it is rejected
 383         if (mcTrack->GetNumberOfTrackReferences() <= 0) continue;
 384         // try to insert in the RsnDaughter its data
 385         check = temp.Adopt(mcTrack);
 386         if (!check) continue;
 387         labmum = temp.GetMCInfo()->Mother();
 388         if (labmum >= 0)
 389         {
 390             TParticle *mum = stack->Particle(labmum);
 391             temp.GetMCInfo()->SetMotherPDG(mum->GetPdgCode());
 392         }
 393         // if successful, set other data and stores it
 394         temp.SetIndex(index);
 395         temp.SetLabel(mcTrack->Label());
 396         AliRsnDaughter *ptr = rsn->AddTrack(temp);
 397         // if problems occurred while storin, that pointer is NULL
 398         if (!ptr) AliWarning(Form("Failed storing track#%d", index));
 399     }
 400
 401     // compute total multiplicity
 402     if (rsn->GetMultiplicity() <= 0)
 403     {
 404         AliDebug(1, "Zero multiplicity in this event");
 405         return kFALSE;
 406     }
 407
 408     return kTRUE;
 409 }
 410
 411 void AliRsnReader::SetPIDtype(const AliRsnDaughter::EPIDType & theValue, Double_t divValue)
 412 {
 413     fCurrentPIDtype = theValue;
 414     fPIDDivValue = divValue;
 415 }
 416
 417
 418 void AliRsnReader::SetITSTPCTRDSectors(const Int_t & its, const Int_t & tpc, const Int_t & trd)
 419 {
 420     fITSClusters = its;
 421     fTPCClusters = tpc;
 422     fTRDClusters = trd;
 423 }