TRD/qaRec/AliTRDpidRefMaker.cxx

   1 #include "TPDGCode.h"
   2 #include "TFile.h"
   3 #include "TTree.h"
   4 #include "TEventList.h"
   5 #include "TObjArray.h"
   6 #include "TH2.h"
   7
   8 #include "AliLog.h"
   9 #include "AliESDtrack.h"
  10 #include "AliTrackReference.h"
  11
  12 #include "AliTRDReconstructor.h"
  13 #include "AliTRDtrackV1.h"
  14 #include "AliTRDseedV1.h"
  15 #include "AliTRDpidRefMaker.h"
  16 #include "info/AliTRDv0Info.h"
  17
  18
  19 // Defines and implements basic functionality for building reference data for TRD PID.
  20 //
  21 // Here is the list of functionality provided by this class
  22 // 1.
  23 // 2.
  24 //
  25 // Authors:
  26 //   Alex Bercuci <A.Bercuci@gsi.de>
  27 //   Alex Wilk    <wilka@uni-muenster.de>
  28 //   Markus Fasel <mfasel@gsi.de>
  29 //   Markus Heide <mheide@uni-muenster.de>
  30 //
  31
  32 ClassImp(AliTRDpidRefMaker)
  33
  34 //________________________________________________________________________
  35 AliTRDpidRefMaker::AliTRDpidRefMaker(const char *name, const char *title)
  36   :AliTRDrecoTask(name, title)
  37   ,fReconstructor(0x0)
  38   ,fV0s(0x0)
  39   ,fData(0x0)
  40   ,fRefPID(kMC)
  41   ,fRefP(kMC)
  42   ,fTrainFreq(1.)
  43   ,fTestFreq(0.)
  44   ,fLy(-1)
  45   ,fP(-1.)
  46 {
  47   //
  48   // Default constructor
  49   //
  50
  51   fReconstructor = new AliTRDReconstructor();
  52   fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
  53   memset(fdEdx, 0, 10*sizeof(Float_t));
  54   memset(fPID, 0, AliPID::kSPECIES*sizeof(Float_t));
  55
  56   DefineInput(1, TObjArray::Class());
  57   DefineOutput(1, TTree::Class());
  58 }
  59
  60
  61 //________________________________________________________________________
  62 AliTRDpidRefMaker::~AliTRDpidRefMaker()
  63 {
  64   if(fReconstructor) delete fReconstructor;
  65 }
  66
  67
  68 //________________________________________________________________________
  69 void AliTRDpidRefMaker::ConnectInputData(Option_t *opt)
  70 {
  71   AliTRDrecoTask::ConnectInputData(opt);
  72   fV0s = dynamic_cast<TObjArray*>(GetInputData(1));
  73 }
  74
  75 //________________________________________________________________________
  76 void AliTRDpidRefMaker::CreateOutputObjects()
  77 {
  78   // Create histograms
  79   // Called once
  80
  81   OpenFile(0, "RECREATE");
  82   fContainer = new TObjArray();
  83   fContainer->SetName(Form("Moni%s", GetName()));
  84
  85   TH2 *h2 = new TH2I("hPDG","Particle abundance", AliPID::kSPECIES, -0.5, 4.5, AliTRDCalPID::kNMom, -0.5, AliTRDCalPID::kNMom-0.5);
  86   TAxis *ax = h2->GetXaxis();
  87   ax->SetNdivisions(505);
  88   ax->SetTitle("Particle species");
  89   for(Int_t is=AliPID::kSPECIES; is--;) ax->SetBinLabel(is+1, AliPID::ParticleShortName(is));
  90   h2->GetYaxis()->SetTitle("P bins");
  91   h2->GetYaxis()->SetNdivisions(511);
  92
  93   fContainer->AddAt(h2, 0);
  94 }
  95
  96 //________________________________________________________________________
  97 void AliTRDpidRefMaker::Exec(Option_t *)
  98 {
  99   // Main loop
 100   // Called for each event
 101
 102
 103
 104
 105
 106 //   DEBUG ?!
 107 //   for(Int_t iv0=0; iv0<fV0s->GetEntriesFast(); iv0++){
 108 //     v0 = dynamic_cast<AliTRDv0Info*>(fV0s->At(iv0));
 109 //     v0->Print();
 110 //   }
 111
 112   Int_t labelsacc[10000]; // MC labels
 113   memset(labelsacc, 0, sizeof(Int_t) * 10000);
 114
 115   AliTRDtrackInfo     *track = 0x0;
 116   //AliTRDv0Info           *v0 = 0x0;
 117   AliTRDtrackV1    *TRDtrack = 0x0;
 118   AliTrackReference     *ref = 0x0;
 119   //AliExternalTrackParam *esd = 0x0;
 120   AliTRDseedV1 *TRDtracklet = 0x0;
 121   for(Int_t itrk=0, nTRD=0; itrk<fTracks->GetEntriesFast(); itrk++){
 122     track = (AliTRDtrackInfo*)fTracks->UncheckedAt(itrk);
 123     if(!track->HasESDtrack()) continue;
 124     ULong_t status = track->GetStatus();
 125     if(!(status&AliESDtrack::kTPCout)) continue;
 126
 127     if(!(TRDtrack = track->GetTrack())) continue;
 128     //&&(track->GetNumberOfClustersRefit()
 129
 130     // TOO STRONG and might introduce a bias if short
 131     // tracks are to be analysed (A.Bercuci 23.09.09)
 132     // use only tracks that hit 6 chambers
 133     //if(!(TRDtrack->GetNumberOfTracklets() == AliTRDgeometry::kNlayer)) continue;
 134
 135
 136     if(HasMCdata()) labelsacc[nTRD++] = track->GetLabel();
 137
 138     // fill the pid information
 139     memset(fPID, 0, AliPID::kSPECIES*sizeof(Float_t));
 140     switch(fRefPID){
 141     case kV0: SetRefPID(kV0, 0x0/*v0*/, fPID); break;
 142     case kMC: SetRefPID(kMC, track, fPID); break;
 143     case kRec: SetRefPID(kRec, TRDtrack, fPID); break;
 144     }
 145
 146     // fill the momentum and dE/dx information
 147     TRDtrack -> SetReconstructor(fReconstructor);
 148     for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++){
 149       if(!(TRDtracklet = TRDtrack -> GetTracklet(ily))) continue;
 150       if(!CookdEdx(TRDtracklet)) continue;
 151       switch(fRefP){
 152       case kMC:
 153         if(!HasMCdata()){
 154           AliError("Could not retrive reference momentum from MC");
 155           return;
 156         }
 157         if(!(ref = track->GetTrackRef(TRDtracklet))) continue;
 158         fP = ref->P();
 159         break;
 160       case kRec:
 161         fP = TRDtracklet->GetMomentum();
 162         break;
 163       default:
 164         AliWarning("Momentum reference source not implemented");
 165         return;
 166       }
 167       fLy = ily;
 168       Fill();
 169     }
 170
 171     for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
 172       AliDebug(4, Form("PDG is %d %f", iPart, fPID[iPart]));
 173     }
 174   }
 175
 176   PostData(0, fContainer);
 177   PostData(1, fData);
 178 }
 179
 180
 181 //________________________________________________________________________
 182 void AliTRDpidRefMaker::Fill()
 183 {
 184   fData->Fill();
 185 }
 186
 187 //________________________________________________________________________
 188 void AliTRDpidRefMaker::Terminate(Option_t *)
 189 {
 190   // Draw result to the screen
 191   // Called once at the end of the query
 192
 193   fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
 194   if (!fContainer) {
 195     AliError("List not available");
 196     return;
 197   }
 198 }
 199
 200
 201 //________________________________________________________________________
 202 void AliTRDpidRefMaker::SetRefPID(ETRDpidRefMakerSource select, void *source, Float_t *pid)
 203 {
 204 // !!!! PREMILMINARY FUNCTION !!!!
 205 //
 206 // this is the place for the V0 procedure
 207 // as long as there is no one implemented,
 208 // just the probabilities
 209 // of the TRDtrack are used!
 210
 211   switch(select){
 212   case kV0:
 213     {
 214       AliTRDv0Info *v0 = static_cast<AliTRDv0Info*>(source);
 215       v0->Print(); // do something
 216     }
 217     break;
 218   case kMC:
 219     if(!HasMCdata()){
 220       AliError("Could not retrive reference PID from MC");
 221       return;
 222     }
 223     {
 224       AliTRDtrackInfo *track = static_cast<AliTRDtrackInfo*>(source);
 225       switch(track->GetPDG()){
 226       case kElectron:
 227       case kPositron:
 228         fPID[AliPID::kElectron] = 1.;
 229         break;
 230       case kMuonPlus:
 231       case kMuonMinus:
 232         fPID[AliPID::kMuon] = 1.;
 233         break;
 234       case kPiPlus:
 235       case kPiMinus:
 236         fPID[AliPID::kPion] = 1.;
 237         break;
 238       case kKPlus:
 239       case kKMinus:
 240         fPID[AliPID::kKaon] = 1.;
 241         break;
 242       case kProton:
 243       case kProtonBar:
 244         fPID[AliPID::kProton] = 1.;
 245         break;
 246       }
 247     }
 248     break;
 249   case kRec:
 250     {
 251       AliTRDtrackV1 *TRDtrack = static_cast<AliTRDtrackV1*>(source);
 252       TRDtrack -> SetReconstructor(fReconstructor);
 253       //fReconstructor -> SetOption("nn");
 254       TRDtrack -> CookPID();
 255       for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
 256         pid[iPart] = TRDtrack -> GetPID(iPart);
 257         AliDebug(4, Form("PDG is (in V0info) %d %f", iPart, pid[iPart]));
 258       }
 259     }
 260     break;
 261   default:
 262     AliWarning("PID reference source not implemented");
 263     return;
 264   }
 265 }
 266
 267 //________________________________________________________________________
 268 void AliTRDpidRefMaker::SetAbundance(Float_t train, Float_t test)
 269 {
 270   if(fTrainFreq<0. || fTrainFreq >1. ||
 271      fTestFreq<0.  || fTestFreq >1.){
 272     AliWarning("The input data should be in the interval [0, 1]");
 273     return;
 274   }
 275   if(TMath::Abs(fTrainFreq+fTestFreq - 1.) > 0.001){
 276     AliWarning("The sum of the 2 abundances should pe one.");
 277     return;
 278   }
 279
 280   fTrainFreq = train;
 281   fTestFreq = test;
 282 }
 283