PWG1/TRD/AliTRDinfoGen.cxx

   1 /**************************************************************************\r
   2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
   3 *                                                                        *\r
   4 * Author: The ALICE Off-line Project.                                    *\r
   5 * Contributors are mentioned in the code where appropriate.              *\r
   6 *                                                                        *\r
   7 * Permission to use, copy, modify and distribute this software and its   *\r
   8 * documentation strictly for non-commercial purposes is hereby granted   *\r
   9 * without fee, provided that the above copyright notice appears in all   *\r
  10 * copies and that both the copyright notice and this permission notice   *\r
  11 * appear in the supporting documentation. The authors make no claims     *\r
  12 * about the suitability of this software for any purpose. It is          *\r
  13 * provided "as is" without express or implied warranty.                  *\r
  14 **************************************************************************/\r
  15 \r
  16 /* $Id: AliTRDinfoGen.cxx 27496 2008-07-22 08:35:45Z cblume $ */\r
  17 \r
  18 ////////////////////////////////////////////////////////////////////////////\r
  19 //\r
  20 //  Tender wagon for TRD performance/calibration train\r
  21 //\r
  22 // In this wagon the information from\r
  23 //   - ESD\r
  24 //   - Friends [if available]\r
  25 //   - MC [if available]\r
  26 // are grouped into AliTRDtrackInfo objects and fed to worker tasks\r
  27 //\r
  28 //  Authors:\r
  29 //    Markus Fasel <M.Fasel@gsi.de>\r
  30 //    Alexandru Bercuci <A.Bercuci@gsi.de>\r
  31 //\r
  32 ////////////////////////////////////////////////////////////////////////////\r
  33 \r
  34 #include <TClonesArray.h>\r
  35 #include <TObjArray.h>\r
  36 #include <TObject.h>\r
  37 #include <TString.h>\r
  38 #include <TH1S.h>\r
  39 #include <TPad.h>\r
  40 #include <TFile.h>\r
  41 #include <TTree.h>\r
  42 #include <TROOT.h>\r
  43 #include <TChain.h>\r
  44 #include <TParticle.h>\r
  45 \r
  46 #include "AliLog.h"\r
  47 #include "AliAnalysisManager.h"\r
  48 #include "AliGeomManager.h"\r
  49 #include "AliCDBManager.h"\r
  50 #include "AliCDBEntry.h"\r
  51 #include "AliCDBPath.h"\r
  52 #include "AliESDEvent.h"\r
  53 #include "AliMCEvent.h"\r
  54 #include "AliESDInputHandler.h"\r
  55 #include "AliMCEventHandler.h"\r
  56 \r
  57 #include "AliESDfriend.h"\r
  58 #include "AliESDfriendTrack.h"\r
  59 #include "AliESDHeader.h"\r
  60 #include "AliESDRun.h"\r
  61 #include "AliESDtrack.h"\r
  62 #include "AliESDv0.h"\r
  63 #include "AliESDtrackCuts.h"\r
  64 #include "AliMCParticle.h"\r
  65 #include "AliMultiplicity.h"\r
  66 #include "AliPID.h"\r
  67 #include "AliStack.h"\r
  68 #include "AliTrackReference.h"\r
  69 #include "TTreeStream.h"\r
  70 \r
  71 #include <cstdio>\r
  72 #include <climits>\r
  73 #include <cstring>\r
  74 #include <iostream>\r
  75 #include <memory>\r
  76 \r
  77 #include "AliTRDReconstructor.h"\r
  78 #include "AliTRDrecoParam.h"\r
  79 #include "AliTRDcalibDB.h"\r
  80 #include "AliTRDtrackerV1.h"\r
  81 #include "AliTRDgeometry.h"\r
  82 #include "AliTRDtrackV1.h"\r
  83 #include "AliTRDseedV1.h"\r
  84 #include "AliTRDcluster.h"\r
  85 #include "AliTRDinfoGen.h"\r
  86 #include "AliTRDpwg1Helper.h"\r
  87 #include "info/AliTRDtrackInfo.h"\r
  88 #include "info/AliTRDeventInfo.h"\r
  89 #include "info/AliTRDv0Info.h"\r
  90 #include "info/AliTRDeventCuts.h"\r
  91 \r
  92 ClassImp(AliTRDinfoGen)\r
  93 \r
  94 const Float_t AliTRDinfoGen::fgkITS = 100.; // to be checked\r
  95 const Float_t AliTRDinfoGen::fgkTPC = 290.;\r
  96 const Float_t AliTRDinfoGen::fgkTRD = 365.;\r
  97 \r
  98 const Float_t AliTRDinfoGen::fgkEvVertexZ = 15.;\r
  99 const Int_t   AliTRDinfoGen::fgkEvVertexN = 1;\r
 100 \r
 101 const Float_t AliTRDinfoGen::fgkTrkDCAxy  = 3.;\r
 102 const Float_t AliTRDinfoGen::fgkTrkDCAz   = 10.;\r
 103 const Int_t   AliTRDinfoGen::fgkNclTPC    = 70;\r
 104 const Float_t AliTRDinfoGen::fgkPt        = 0.2;\r
 105 const Float_t AliTRDinfoGen::fgkEta       = 0.9;\r
 106 AliTRDReconstructor* AliTRDinfoGen::fgReconstructor(NULL);\r
 107 AliTRDgeometry* AliTRDinfoGen::fgGeo(NULL);\r
 108 //____________________________________________________________________\r
 109 AliTRDinfoGen::AliTRDinfoGen()\r
 110   :AliAnalysisTaskSE()\r
 111   ,fEvTrigger(NULL)\r
 112   ,fESDev(NULL)\r
 113   ,fMCev(NULL)\r
 114   ,fEventCut(NULL)\r
 115   ,fTrackCut(NULL)\r
 116   ,fV0Cut(NULL)\r
 117   ,fOCDB("local://$ALICE_ROOT/OCDB")\r
 118   ,fTrackInfo(NULL)\r
 119   ,fEventInfo(NULL)\r
 120   ,fV0Info(NULL)\r
 121   ,fTracksBarrel(NULL)\r
 122   ,fTracksSA(NULL)\r
 123   ,fTracksKink(NULL)\r
 124   ,fV0List(NULL)\r
 125   ,fContainer(NULL)\r
 126   ,fDebugStream(NULL)\r
 127 {\r
 128   //\r
 129   // Default constructor\r
 130   //\r
 131   SetNameTitle("TRDinfoGen", "MC-REC TRD-track list generator");\r
 132 }\r
 133 \r
 134 //____________________________________________________________________\r
 135 AliTRDinfoGen::AliTRDinfoGen(char* name)\r
 136   :AliAnalysisTaskSE(name)\r
 137   ,fEvTrigger(NULL)\r
 138   ,fESDev(NULL)\r
 139   ,fMCev(NULL)\r
 140   ,fEventCut(NULL)\r
 141   ,fTrackCut(NULL)\r
 142   ,fV0Cut(NULL)\r
 143   ,fOCDB("local://$ALICE_ROOT/OCDB")\r
 144   ,fTrackInfo(NULL)\r
 145   ,fEventInfo(NULL)\r
 146   ,fV0Info(NULL)\r
 147   ,fTracksBarrel(NULL)\r
 148   ,fTracksSA(NULL)\r
 149   ,fTracksKink(NULL)\r
 150   ,fV0List(NULL)\r
 151   ,fContainer(NULL)\r
 152   ,fDebugStream(NULL)\r
 153 {\r
 154   //\r
 155   // Default constructor\r
 156   //\r
 157   SetTitle("MC-REC TRD-track list generator");\r
 158   DefineOutput(AliTRDpwg1Helper::kTracksBarrel, TObjArray::Class());\r
 159   DefineOutput(AliTRDpwg1Helper::kTracksSA, TObjArray::Class());\r
 160   DefineOutput(AliTRDpwg1Helper::kTracksKink, TObjArray::Class());\r
 161   DefineOutput(AliTRDpwg1Helper::kEventInfo, AliTRDeventInfo::Class());\r
 162   DefineOutput(AliTRDpwg1Helper::kV0List, TObjArray::Class());\r
 163   DefineOutput(AliTRDpwg1Helper::kMonitor, TObjArray::Class()); // histogram list\r
 164 }\r
 165 \r
 166 //____________________________________________________________________\r
 167 AliTRDinfoGen::~AliTRDinfoGen()\r
 168 {\r
 169 // Destructor\r
 170   if(fgGeo) delete fgGeo;\r
 171   if(fgReconstructor) delete fgReconstructor;\r
 172   if(fDebugStream) delete fDebugStream;\r
 173   if(fEvTrigger) delete fEvTrigger;\r
 174   if(fV0Cut) delete fV0Cut;\r
 175   if(fTrackCut) delete fTrackCut;\r
 176   if(fEventCut) delete fEventCut;\r
 177   if(fTrackInfo) delete fTrackInfo; fTrackInfo = NULL;\r
 178   if(fEventInfo) delete fEventInfo; fEventInfo = NULL;\r
 179   if(fV0Info) delete fV0Info; fV0Info = NULL;\r
 180   if(fTracksBarrel){ \r
 181     fTracksBarrel->Delete(); delete fTracksBarrel;\r
 182     fTracksBarrel = NULL;\r
 183   }\r
 184   if(fTracksSA){ \r
 185     fTracksSA->Delete(); delete fTracksSA;\r
 186     fTracksSA = NULL;\r
 187   }\r
 188   if(fTracksKink){ \r
 189     fTracksKink->Delete(); delete fTracksKink;\r
 190     fTracksKink = NULL;\r
 191   }\r
 192   if(fV0List){ \r
 193     fV0List->Delete(); \r
 194     delete fV0List;\r
 195     fV0List = NULL;\r
 196   }\r
 197   if(fContainer){ \r
 198     fContainer->Delete(); \r
 199     delete fContainer;\r
 200     fContainer = NULL;\r
 201   }\r
 202 }\r
 203 \r
 204 //____________________________________________________________________\r
 205 Bool_t AliTRDinfoGen::GetRefFigure(Int_t)\r
 206 {\r
 207 // General graphs for PWG1/TRD train\r
 208   if(!gPad){\r
 209     AliWarning("Please provide a canvas to draw results.");\r
 210     return kFALSE;\r
 211   }\r
 212   fContainer->At(0)->Draw("bar");\r
 213   return kTRUE;\r
 214 }\r
 215 \r
 216 //____________________________________________________________________\r
 217 void AliTRDinfoGen::UserCreateOutputObjects()\r
 218 {       \r
 219   //\r
 220   // Create Output Containers (TObjectArray containing 1D histograms)\r
 221   //\r
 222  \r
 223   fTrackInfo = new AliTRDtrackInfo();\r
 224   fEventInfo = new AliTRDeventInfo();\r
 225   fV0Info    = new AliTRDv0Info();\r
 226   fTracksBarrel = new TObjArray(200); fTracksBarrel->SetOwner(kTRUE);\r
 227   fTracksSA = new TObjArray(20); fTracksSA->SetOwner(kTRUE);\r
 228   fTracksKink = new TObjArray(20); fTracksKink->SetOwner(kTRUE);\r
 229   fV0List = new TObjArray(10); fV0List->SetOwner(kTRUE);\r
 230 \r
 231   // define general monitor\r
 232   fContainer = new TObjArray(1); fContainer->SetOwner(kTRUE);\r
 233   TH1 *h=new TH1I("hStat", "Run statistics;Observable;Entries", Int_t(kNObjects), -0.5, Float_t(kNObjects)-0.5);\r
 234   TAxis *ax(h->GetXaxis());\r
 235   ax->SetBinLabel(Int_t(kTracksESD) + 1, "ESD");\r
 236   ax->SetBinLabel(Int_t(kTracksMC) + 1, "MC");\r
 237   ax->SetBinLabel(Int_t(kV0) + 1, "V0");\r
 238   ax->SetBinLabel(Int_t(kTPC) + 1, "TPC");\r
 239   ax->SetBinLabel(Int_t(kTRDin) + 1, "TRDin");\r
 240   ax->SetBinLabel(Int_t(kTRDout) + 1, "TRDout");\r
 241   ax->SetBinLabel(Int_t(kBarrel) + 1, "Barrel");\r
 242   ax->SetBinLabel(Int_t(kBarrelMC) + 1, "BarrelMC");\r
 243   ax->SetBinLabel(Int_t(kSA) + 1, "SA");\r
 244   ax->SetBinLabel(Int_t(kSAMC) + 1, "SAMC");\r
 245   ax->SetBinLabel(Int_t(kKink) + 1, "Kink");\r
 246   ax->SetBinLabel(Int_t(kKinkMC) + 1, "KinkMC");\r
 247   ax->SetBinLabel(Int_t(kBarrelFriend) + 1, "BFriend");\r
 248   ax->SetBinLabel(Int_t(kSAFriend) + 1, "SFriend");\r
 249   fContainer->AddAt(h, 0);\r
 250   PostData(AliTRDpwg1Helper::kMonitor, fContainer);\r
 251 }\r
 252 \r
 253 //____________________________________________________________________\r
 254 Bool_t AliTRDinfoGen::Load(const Char_t *file, const Char_t *dir, const Char_t *name)\r
 255 {\r
 256 // Load data from performance file\r
 257 \r
 258   if(!TFile::Open(file)){\r
 259     AliWarning(Form("Couldn't open file %s.", file));\r
 260     return kFALSE;\r
 261   }\r
 262   if(dir){\r
 263     if(!gFile->cd(dir)){\r
 264       AliWarning(Form("Couldn't cd to %s in %s.", dir, file));\r
 265       return kFALSE;\r
 266     }\r
 267   }\r
 268   TObjArray *o(NULL);\r
 269   const Char_t *tn=(name ? name : GetName());\r
 270   if(!(o = (TObjArray*)gDirectory->Get(tn))){\r
 271     AliWarning(Form("Missing histogram container %s.", tn));\r
 272     return kFALSE;\r
 273   }\r
 274   fContainer = (TObjArray*)o->Clone(GetName());\r
 275   gFile->Close();\r
 276   return kTRUE;\r
 277 }\r
 278 \r
 279 //____________________________________________________________________\r
 280 void AliTRDinfoGen::UserExec(Option_t *){\r
 281   //\r
 282   // Run the Analysis\r
 283   //\r
 284 \r
 285   fTracksBarrel->Delete();\r
 286   fTracksSA->Delete();\r
 287   fTracksKink->Delete();\r
 288   fV0List->Delete();\r
 289   fEventInfo->Delete("");\r
 290 \r
 291   fESDev = dynamic_cast<AliESDEvent*>(InputEvent());\r
 292   if(!fESDev){\r
 293     AliError("Failed retrieving ESD event");\r
 294     return;\r
 295   }\r
 296   // WARNING\r
 297   // This part may conflict with other detectors !!\r
 298   if(!IsInitOCDB()){ \r
 299     AliInfo("Initializing OCDB ...");\r
 300     // prepare OCDB access\r
 301     AliCDBManager* ocdb = AliCDBManager::Instance();\r
 302     ocdb->SetDefaultStorage(fOCDB.Data());\r
 303     ocdb->SetRun(fESDev->GetRunNumber());\r
 304     // create geo manager\r
 305     AliCDBEntry* obj(NULL);\r
 306     if(!(obj = ocdb->Get(AliCDBPath("GRP", "Geometry", "Data")))){\r
 307       AliError("GEOMETRY failed initialization.");\r
 308     } else {\r
 309       AliGeomManager::SetGeometry((TGeoManager*)obj->GetObject());\r
 310       AliGeomManager::GetNalignable("TRD");\r
 311       AliGeomManager::ApplyAlignObjsFromCDB("TRD");\r
 312     }\r
 313     fgGeo = new AliTRDgeometry;\r
 314     fgGeo->CreateClusterMatrixArray();\r
 315     //init magnetic field\r
 316     if(!TGeoGlobalMagField::Instance()->IsLocked() &&\r
 317        !fESDev->InitMagneticField()){\r
 318       AliError("MAGNETIC FIELD failed initialization.");\r
 319     }\r
 320 \r
 321     // set no of time bins\r
 322     AliTRDtrackerV1::SetNTimeBins(AliTRDcalibDB::Instance()->GetNumberOfTimeBinsDCS());\r
 323     AliInfo(Form("OCDB :  Loc[%s] Run[%d] TB[%d]", fOCDB.Data(), ocdb->GetRun(), AliTRDtrackerV1::GetNTimeBins()));\r
 324 \r
 325     // set reco param valid for this run/event\r
 326     AliInfo(Form("Initializing TRD reco params for EventSpecie[%d]...",\r
 327       fESDev->GetEventSpecie()));\r
 328     fgReconstructor = new AliTRDReconstructor();\r
 329     if(!(obj = ocdb->Get(AliCDBPath("TRD", "Calib", "RecoParam")))){\r
 330       AliError("RECO PARAM failed initialization.");\r
 331       fgReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());\r
 332     } else {\r
 333       obj->PrintMetaData();\r
 334       TObjArray *recos((TObjArray*)obj->GetObject());\r
 335       for(Int_t ireco(0); ireco<recos->GetEntriesFast(); ireco++){\r
 336         AliTRDrecoParam *reco((AliTRDrecoParam*)recos->At(ireco));\r
 337         Int_t es(reco->GetEventSpecie());\r
 338         if(!(es&fESDev->GetEventSpecie())) continue;\r
 339         fgReconstructor->SetRecoParam(reco);\r
 340         TString s;\r
 341         if(es&AliRecoParam::kLowMult) s="LowMult";\r
 342         else if(es&AliRecoParam::kHighMult) s="HighMult";\r
 343         else if(es&AliRecoParam::kCosmic) s="Cosmic";\r
 344         else if(es&AliRecoParam::kCalib) s="Calib";\r
 345         else s="Unknown";\r
 346         AliInfo(Form("Using reco params for %s", s.Data()));\r
 347         break;\r
 348       }\r
 349     }\r
 350     SetInitOCDB();\r
 351   }\r
 352 \r
 353   // link MC if available\r
 354   fMCev = MCEvent();\r
 355 \r
 356   // event selection : trigger cut\r
 357   if(UseLocalEvSelection() && fEvTrigger){ \r
 358     Bool_t kTRIGGERED(kFALSE);\r
 359     std::auto_ptr<TObjArray> trig(fEvTrigger->Tokenize(" "));\r
 360     for(Int_t itrig=trig->GetEntriesFast(); itrig--;){\r
 361       const Char_t *trigClass(((TObjString*)(*trig)[itrig])->GetName());\r
 362       if(fESDev->IsTriggerClassFired(trigClass)) {\r
 363         AliDebug(2, Form("Ev[%4d] Trigger[%s]", fESDev->GetEventNumberInFile(), trigClass));\r
 364         kTRIGGERED = kTRUE;\r
 365         break; \r
 366       }\r
 367     }\r
 368     if(!kTRIGGERED){ \r
 369       AliDebug(2, Form("Reject Ev[%4d] Trigger", fESDev->GetEventNumberInFile()));\r
 370       return;\r
 371     }\r
 372     // select only physical events\r
 373     if(fESDev->GetEventType() != 7){ \r
 374       AliDebug(2, Form("Reject Ev[%4d] EvType[%d]", fESDev->GetEventNumberInFile(), fESDev->GetEventType()));\r
 375       return;\r
 376     }\r
 377   }\r
 378 \r
 379   // if the required trigger is a collision trigger then apply event vertex cut\r
 380   if(UseLocalEvSelection() && IsCollision()){\r
 381     const AliESDVertex *vertex = fESDev->GetPrimaryVertex();\r
 382     if(TMath::Abs(vertex->GetZv())<1.e-10 || \r
 383        TMath::Abs(vertex->GetZv())>fgkEvVertexZ || \r
 384        vertex->GetNContributors()<fgkEvVertexN) {\r
 385       AliDebug(2, Form("Reject Ev[%4d] Vertex Zv[%f] Nv[%d]", fESDev->GetEventNumberInFile(), TMath::Abs(vertex->GetZv()), vertex->GetNContributors()));\r
 386       return;\r
 387     }\r
 388   }\r
 389 \r
 390   if(fEventCut && !fEventCut->IsSelected(fESDev, IsCollision())) return;\r
 391 \r
 392   if(!fESDfriend){\r
 393     AliError("Failed retrieving ESD friend event");\r
 394     return;\r
 395   }\r
 396   if(HasMCdata() && !fMCev){\r
 397     AliError("Failed retrieving MC event");\r
 398     return;\r
 399   }\r
 400 \r
 401   new(fEventInfo)AliTRDeventInfo(fESDev->GetHeader(), const_cast<AliESDRun *>(fESDev->GetESDRun()));\r
 402   // Determine centrality\r
 403   // Author: Ionut Arsene <I.C.Arsene@gsi.de>\r
 404   Int_t centralityBin = -1;\r
 405   AliDebug(2, Form("Beam Type: %s", fESDev->GetESDRun()->GetBeamType()));\r
 406   if(TString(fESDev->GetESDRun()->GetBeamType()).Contains("Pb-Pb")){\r
 407     centralityBin = 4;\r
 408     const AliMultiplicity *mult = fESDev->GetMultiplicity();\r
 409     Double_t zdcNeutronEnergy = fESDev->GetZDCN1Energy()+fESDev->GetZDCN2Energy();\r
 410     Double_t itsNTracklets = mult->GetNumberOfTracklets();\r
 411     Double_t centralitySlopes[6] = {0.0, 4.0, 8.0, 20.0, 50.0, 1000000.};\r
 412     AliDebug(1, Form("zdcNeutronEnergy: %f, itsNTracklets: %f\n", zdcNeutronEnergy, itsNTracklets));\r
 413     for(Int_t iCent=1; iCent<=5; ++iCent) {\r
 414       if(zdcNeutronEnergy>centralitySlopes[iCent-1]*itsNTracklets && zdcNeutronEnergy<centralitySlopes[iCent]*itsNTracklets)\r
 415         centralityBin=iCent - 1;\r
 416     }\r
 417     AliDebug(1, Form("Centrality Class: %d", centralityBin));\r
 418   }\r
 419   fEventInfo->SetCentrality(centralityBin);\r
 420 \r
 421   Bool_t *trackMap(NULL);\r
 422   AliStack * mStack(NULL);\r
 423   Int_t nTracksMC = HasMCdata() ? fMCev->GetNumberOfTracks() : 0, nTracksESD = fESDev->GetNumberOfTracks();\r
 424   if(HasMCdata()){\r
 425     mStack = fMCev->Stack();\r
 426     if(!mStack){\r
 427       AliError("Failed retrieving MC Stack");\r
 428       return;\r
 429     }\r
 430     trackMap = new Bool_t[nTracksMC];\r
 431     memset(trackMap, 0, sizeof(Bool_t) * nTracksMC);\r
 432   }\r
 433   \r
 434   Double32_t dedx[100]; Int_t nSlices(0);\r
 435   Int_t nTRDout(0), nTRDin(0), nTPC(0)\r
 436        ,nclsTrklt\r
 437        ,nBarrel(0), nSA(0), nKink(0)\r
 438        ,nBarrelFriend(0), nSAFriend(0)\r
 439        ,nBarrelMC(0), nSAMC(0), nKinkMC(0);\r
 440   AliESDtrack *esdTrack = NULL;\r
 441   AliESDfriendTrack *esdFriendTrack = NULL;\r
 442   TObject *calObject = NULL;\r
 443   AliTRDtrackV1 *track = NULL;\r
 444   AliTRDseedV1 *tracklet = NULL;\r
 445   AliTRDcluster *cl = NULL;\r
 446 \r
 447 \r
 448   // LOOP 0 - over ESD v0s\r
 449   Float_t bField(fESDev->GetMagneticField());\r
 450   AliESDv0 *v0(NULL);\r
 451   Int_t v0pid[AliPID::kSPECIES];\r
 452   for(Int_t iv0(0); iv0<fESDev->GetNumberOfV0s(); iv0++){\r
 453     if(!(v0 = fESDev->GetV0(iv0))) continue;\r
 454     // register v0\r
 455     if(fV0Cut) new(fV0Info) AliTRDv0Info(*fV0Cut);\r
 456     else new(fV0Info) AliTRDv0Info();\r
 457     fV0Info->SetMagField(bField);\r
 458     fV0Info->SetV0tracks(fESDev->GetTrack(v0->GetPindex()), fESDev->GetTrack(v0->GetNindex()));\r
 459     fV0Info->SetV0Info(v0);\r
 460     fV0List->Add(new AliTRDv0Info(*fV0Info));//  kFOUND=kFALSE;\r
 461   }\r
 462 \r
 463 \r
 464   // LOOP 1 - over ESD tracks\r
 465   AliTRDv0Info *v0info=NULL;\r
 466   for(Int_t itrk = 0; itrk < nTracksESD; itrk++){\r
 467     new(fTrackInfo) AliTRDtrackInfo();\r
 468     esdTrack = fESDev->GetTrack(itrk);\r
 469     AliDebug(3, Form("\n%3d ITS[%d] TPC[%d] TRD[%d]\n", itrk, esdTrack->GetNcls(0), esdTrack->GetNcls(1), esdTrack->GetNcls(2)));\r
 470 \r
 471     if(esdTrack->GetStatus()&AliESDtrack::kTPCout) nTPC++;\r
 472     if(esdTrack->GetStatus()&AliESDtrack::kTRDout) nTRDout++;\r
 473     if(esdTrack->GetStatus()&AliESDtrack::kTRDin) nTRDin++;\r
 474 \r
 475     // look at external track param\r
 476     const AliExternalTrackParam *op = esdTrack->GetOuterParam();\r
 477     Double_t xyz[3];\r
 478     if(op){\r
 479       op->GetXYZ(xyz);\r
 480       op->Global2LocalPosition(xyz, op->GetAlpha());\r
 481       AliDebug(3, Form("op @ X[%7.3f]\n", xyz[0]));\r
 482     }\r
 483 \r
 484     // read MC info\r
 485     Int_t fPdg = -1;\r
 486     Int_t label = -1; UInt_t alab=UINT_MAX;\r
 487     if(HasMCdata()){\r
 488       label = esdTrack->GetLabel(); \r
 489       alab = TMath::Abs(label);\r
 490       // register the track\r
 491       if(alab < UInt_t(nTracksMC)){ \r
 492         trackMap[alab] = kTRUE; \r
 493       } else { \r
 494         AliError(Form("MC label[%d] outside scope for Ev[%d] Trk[%d].", label, (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), itrk));\r
 495         continue; \r
 496       }\r
 497       AliMCParticle *mcParticle = NULL; \r
 498       if(!(mcParticle = (AliMCParticle*) fMCev->GetTrack(alab))){\r
 499         AliError(Form("MC particle label[%d] missing for Ev[%d] Trk[%d].", label, (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), itrk));\r
 500         continue;\r
 501       }\r
 502       fPdg = mcParticle->Particle()->GetPdgCode();\r
 503       Int_t nRefs = mcParticle->GetNumberOfTrackReferences();\r
 504       Int_t iref = 0; AliTrackReference *ref = NULL; \r
 505       while(iref<nRefs){\r
 506         ref = mcParticle->GetTrackReference(iref);\r
 507         if(ref->LocalX() > fgkTPC) break;\r
 508         iref++;\r
 509       }\r
 510 \r
 511       fTrackInfo->SetMC();\r
 512       fTrackInfo->SetPDG(fPdg);\r
 513       fTrackInfo->SetPrimary(mcParticle->Particle()->IsPrimary());\r
 514       fTrackInfo->SetLabel(label);\r
 515       Int_t jref = iref;//, kref = 0;\r
 516       while(jref<nRefs){\r
 517         ref = mcParticle->GetTrackReference(jref);\r
 518         if(ref->LocalX() > fgkTRD) break;\r
 519         AliDebug(4, Form("  trackRef[%2d (%2d)] @ %7.3f OK", jref-iref, jref, ref->LocalX()));\r
 520         fTrackInfo->AddTrackRef(ref);\r
 521         jref++;\r
 522       }\r
 523       AliDebug(3, Form("NtrackRefs[%d(%d)]", fTrackInfo->GetNTrackRefs(), nRefs));\r
 524     }\r
 525 \r
 526     // copy some relevant info to TRD track info\r
 527     fTrackInfo->SetStatus(esdTrack->GetStatus());\r
 528     fTrackInfo->SetTrackId(esdTrack->GetID());\r
 529     Double_t p[AliPID::kSPECIES]; esdTrack->GetTRDpid(p);\r
 530     fTrackInfo->SetESDpid(p);\r
 531     fTrackInfo->SetESDpidQuality(esdTrack->GetTRDntrackletsPID());\r
 532     if(!nSlices) nSlices = esdTrack->GetNumberOfTRDslices();\r
 533     memset(dedx, 0, 100*sizeof(Double32_t));\r
 534     Int_t in(0);\r
 535     for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++)\r
 536       for(Int_t is=0; is<nSlices; is++) \r
 537         dedx[in++]=esdTrack->GetTRDslice(il, is);\r
 538     for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++) dedx[in++]=esdTrack->GetTRDmomentum(il);\r
 539     fTrackInfo->SetSlices(in, dedx);\r
 540     fTrackInfo->SetNumberOfClustersRefit(esdTrack->GetNcls(2));\r
 541     // some other Informations which we may wish to store in order to find problematic cases\r
 542     fTrackInfo->SetKinkIndex(esdTrack->GetKinkIndex(0));\r
 543     fTrackInfo->SetTPCncls(static_cast<UShort_t>(esdTrack->GetNcls(1)));\r
 544     nclsTrklt = 0;\r
 545   \r
 546     // set V0pid info\r
 547     for(Int_t iv(0); iv<fV0List->GetEntriesFast(); iv++){\r
 548       if(!(v0info = (AliTRDv0Info*)fV0List->At(iv))) continue;\r
 549       if(!v0info->GetV0Daughter(1) && !v0info->GetV0Daughter(-1)) continue;\r
 550       if(!v0info->HasTrack(fTrackInfo)) continue;\r
 551       memset(v0pid, 0, AliPID::kSPECIES*sizeof(Int_t));\r
 552       fTrackInfo->SetV0();\r
 553       for(Int_t is=AliPID::kSPECIES; is--;){v0pid[is] = v0info->GetPID(is, fTrackInfo);}\r
 554       fTrackInfo->SetV0pid(v0pid);\r
 555       fTrackInfo->SetV0();\r
 556       //const AliTRDtrackInfo::AliESDinfo *ei = fTrackInfo->GetESDinfo();\r
 557       break;\r
 558     }\r
 559 \r
 560     // read REC info\r
 561     esdFriendTrack = (fESDfriend->GetNumberOfTracks() > itrk) ? fESDfriend->GetTrack(itrk): NULL;\r
 562 \r
 563     if(esdFriendTrack){\r
 564       Int_t icalib = 0;\r
 565       while((calObject = esdFriendTrack->GetCalibObject(icalib++))){\r
 566         if(strcmp(calObject->IsA()->GetName(),"AliTRDtrackV1") != 0) continue; // Look for the TRDtrack\r
 567         if(!(track = dynamic_cast<AliTRDtrackV1*>(calObject))) break;\r
 568         AliDebug(4, Form("TRD track OK"));\r
 569         // Set the clusters to unused\r
 570         for(Int_t ipl = 0; ipl < AliTRDgeometry::kNlayer; ipl++){\r
 571           if(!(tracklet = track->GetTracklet(ipl))) continue;\r
 572           tracklet->ResetClusterIter();\r
 573           while((cl = tracklet->NextCluster())) cl->Use(0);\r
 574         }\r
 575         fTrackInfo->SetTrack(track);\r
 576         break;\r
 577       }\r
 578       AliDebug(3, Form("Ntracklets[%d]\n", fTrackInfo->GetNTracklets()));\r
 579     } else AliDebug(3, "No ESD friends");\r
 580     if(op) fTrackInfo->SetOuterParam(op);\r
 581 \r
 582     if(DebugLevel() >= 1){\r
 583       AliTRDtrackInfo info(*fTrackInfo);\r
 584       (*DebugStream()) << "trackInfo"\r
 585       << "TrackInfo.=" << &info\r
 586       << "\n";\r
 587       info.Delete("");\r
 588     }\r
 589 \r
 590     ULong_t status(esdTrack->GetStatus());\r
 591     if((status&AliESDtrack::kTPCout)){\r
 592       if(!esdTrack->GetKinkIndex(0)){ // Barrel  Track Selection\r
 593         Bool_t selected(kTRUE);\r
 594         if(UseLocalTrkSelection()){\r
 595           if(esdTrack->Pt() < fgkPt){ \r
 596             AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] Pt[%5.2f]", fESDev->GetEventNumberInFile(), itrk, esdTrack->Pt()));\r
 597             selected = kFALSE;\r
 598           }\r
 599           if(selected && TMath::Abs(esdTrack->Eta()) > fgkEta){\r
 600             AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] Eta[%5.2f]", fESDev->GetEventNumberInFile(), itrk, TMath::Abs(esdTrack->Eta())));\r
 601             selected = kFALSE;\r
 602           }\r
 603           if(selected && esdTrack->GetTPCNcls() < fgkNclTPC){ \r
 604             AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] NclTPC[%d]", fESDev->GetEventNumberInFile(), itrk, esdTrack->GetTPCNcls()));\r
 605             selected = kFALSE;\r
 606           }\r
 607           Float_t par[2], cov[3];\r
 608           esdTrack->GetImpactParameters(par, cov);\r
 609           if(IsCollision()){ // cuts on DCA\r
 610             if(selected && TMath::Abs(par[0]) > fgkTrkDCAxy){ \r
 611               AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] DCAxy[%f]", fESDev->GetEventNumberInFile(), itrk, TMath::Abs(par[0])));\r
 612               selected = kFALSE;\r
 613             }\r
 614             if(selected && TMath::Abs(par[1]) > fgkTrkDCAz){ \r
 615               AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] DCAz[%f]", fESDev->GetEventNumberInFile(), itrk, TMath::Abs(par[1])));\r
 616               selected = kFALSE;\r
 617             }\r
 618           } else if(selected && fMCev && !fMCev->IsPhysicalPrimary(alab)){;\r
 619             AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] Primary", fESDev->GetEventNumberInFile(), itrk));\r
 620             selected = kFALSE;\r
 621           }\r
 622         }\r
 623         if(fTrackCut && !fTrackCut->IsSelected(esdTrack)) selected = kFALSE;\r
 624         if(selected){ \r
 625           fTracksBarrel->Add(new AliTRDtrackInfo(*fTrackInfo));\r
 626           nBarrel++;\r
 627           if(fTrackInfo->GetTrack()) \r
 628             nBarrelFriend++;\r
 629         }\r
 630       } else {\r
 631         fTracksKink->Add(new AliTRDtrackInfo(*fTrackInfo));\r
 632         nKink++;\r
 633       }\r
 634     } else if((status&AliESDtrack::kTRDout) && !(status&AliESDtrack::kTRDin)){ \r
 635       fTracksSA->Add(new AliTRDtrackInfo(*fTrackInfo));\r
 636       nSA++;\r
 637       if(fTrackInfo->GetTrack()) \r
 638         nSAFriend++;\r
 639     }\r
 640     fTrackInfo->Delete("");\r
 641   }\r
 642 \r
 643   // LOOP 2 - over MC tracks which are passing TRD where the track is not reconstructed\r
 644   if(HasMCdata()){\r
 645     AliDebug(10, "Output of the MC track map:");\r
 646     for(Int_t itk = 0; itk < nTracksMC;  itk++) AliDebug(10, Form("trackMap[%d] = %s", itk, trackMap[itk] == kTRUE ? "TRUE" : "kFALSE"));\r
 647   \r
 648     for(Int_t itk = 0; itk < nTracksMC; itk++){\r
 649       if(trackMap[itk]) continue;\r
 650       AliMCParticle *mcParticle =  (AliMCParticle*) fMCev->GetTrack(TMath::Abs(itk));\r
 651       Int_t fPdg = mcParticle->Particle()->GetPdgCode();\r
 652       Int_t nRefs = mcParticle->GetNumberOfTrackReferences();\r
 653       Int_t iref = 0; AliTrackReference *ref = NULL; \r
 654       Int_t nRefsTRD = 0;\r
 655       new(fTrackInfo) AliTRDtrackInfo();\r
 656       fTrackInfo->SetMC();\r
 657       fTrackInfo->SetPDG(fPdg);\r
 658       while(iref<nRefs){ // count TRD TR\r
 659         Bool_t kIN(kFALSE);\r
 660         ref = mcParticle->GetTrackReference(iref);\r
 661         if(ref->LocalX() > fgkTPC && ref->LocalX() < fgkTRD){\r
 662           fTrackInfo->AddTrackRef(ref);\r
 663           nRefsTRD++;kIN=kTRUE;\r
 664         }\r
 665         AliDebug(4, Form("  trackRef[%2d] @ x[%7.3f] %s", iref, ref->LocalX(), kIN?"IN":"OUT"));\r
 666         iref++;\r
 667       }\r
 668       if(!nRefsTRD){\r
 669         // In this stage we at least require 1 hit inside TRD. What will be done with this tracks is a task for the \r
 670         // analysis job\r
 671         fTrackInfo->Delete("");\r
 672         continue;\r
 673       }\r
 674       fTrackInfo->SetPrimary(mcParticle->Particle()->IsPrimary());\r
 675       fTrackInfo->SetLabel(itk);\r
 676       if(DebugLevel() >= 1){\r
 677         AliTRDtrackInfo info(*fTrackInfo);\r
 678         (*DebugStream()) << "trackInfo"\r
 679         << "TrackInfo.=" << &info\r
 680         << "\n";\r
 681         info.Delete("");\r
 682       }\r
 683       AliDebug(3, Form("Add MC track @ label[%d] nTRDrefs[%d].", itk, nRefsTRD));\r
 684       // check where the track starts\r
 685       ref = mcParticle->GetTrackReference(0);\r
 686       if(ref->LocalX() < fgkITS){ \r
 687         fTracksBarrel->Add(new AliTRDtrackInfo(*fTrackInfo));\r
 688         nBarrelMC++;\r
 689       } else if(ref->LocalX() < fgkTPC) {\r
 690         fTracksKink->Add(new AliTRDtrackInfo(*fTrackInfo));\r
 691         nKinkMC++;\r
 692       } else if(nRefsTRD>6){\r
 693         fTracksSA->Add(new AliTRDtrackInfo(*fTrackInfo));\r
 694         nSAMC++;\r
 695       }\r
 696       fTrackInfo->Delete("");\r
 697     }\r
 698     delete[] trackMap;\r
 699   }\r
 700   AliDebug(1, Form(\r
 701     "\nEv[%3d] Tracks: ESD[%d] MC[%d] V0[%d]\n"\r
 702     "        TPCout[%d] TRDin[%d] TRDout[%d]\n"\r
 703     "        Barrel[%3d+%3d=%3d] SA[%2d+%2d=%2d] Kink[%2d+%2d=%2d]"\r
 704     ,(Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), nTracksESD, nTracksMC, fV0List->GetEntries()\r
 705     , nTPC, nTRDin, nTRDout\r
 706     ,nBarrel, nBarrelMC, fTracksBarrel->GetEntries()\r
 707     ,nSA, nSAMC, fTracksSA->GetEntries()\r
 708     ,nKink, nKinkMC, fTracksKink->GetEntries()\r
 709   ));\r
 710   // save track statistics\r
 711   TH1 *h((TH1S*)fContainer->At(0));\r
 712   h->Fill(Float_t(kTracksESD), nTracksESD);\r
 713   h->Fill(Float_t(kTracksMC), nTracksMC);\r
 714   h->Fill(Float_t(kV0), fV0List->GetEntries());\r
 715   h->Fill(Float_t(kTPC), nTPC);\r
 716   h->Fill(Float_t(kTRDin), nTRDin);\r
 717   h->Fill(Float_t(kTRDout), nTRDout);\r
 718   h->Fill(Float_t(kBarrel), nBarrel);\r
 719   h->Fill(Float_t(kBarrelMC), nBarrelMC);\r
 720   h->Fill(Float_t(kSA), nSA);\r
 721   h->Fill(Float_t(kSAMC), nSAMC);\r
 722   h->Fill(Float_t(kKink), nKink);\r
 723   h->Fill(Float_t(kKinkMC), nKinkMC);\r
 724   h->Fill(Float_t(kBarrelFriend), nBarrelFriend);\r
 725   h->Fill(Float_t(kSAFriend), nSAFriend);\r
 726 \r
 727   PostData(AliTRDpwg1Helper::kTracksBarrel, fTracksBarrel);\r
 728   PostData(AliTRDpwg1Helper::kTracksSA, fTracksSA);\r
 729   PostData(AliTRDpwg1Helper::kTracksKink, fTracksKink);\r
 730   PostData(AliTRDpwg1Helper::kEventInfo, fEventInfo);\r
 731   PostData(AliTRDpwg1Helper::kV0List, fV0List);\r
 732 }\r
 733 \r
 734 //____________________________________________________________________\r
 735 void AliTRDinfoGen::SetLocalV0Selection(const AliTRDv0Info *v0)\r
 736 {\r
 737 // Set V0 cuts from outside\r
 738 \r
 739   if(!fV0Cut) fV0Cut = new AliTRDv0Info(*v0);\r
 740   else new(fV0Cut) AliTRDv0Info(*v0);\r
 741 }\r
 742 \r
 743 //____________________________________________________________________\r
 744 void AliTRDinfoGen::SetTrigger(const Char_t *trigger)\r
 745 {\r
 746   if(!fEvTrigger) fEvTrigger = new TString(trigger);\r
 747   else (*fEvTrigger) = trigger;\r
 748 }\r
 749 \r
 750 //____________________________________________________________________\r
 751 TTreeSRedirector* AliTRDinfoGen::DebugStream()\r
 752 {\r
 753 // Manage debug stream for task\r
 754   if(!fDebugStream){\r
 755     TDirectory *savedir = gDirectory;\r
 756     fDebugStream = new TTreeSRedirector("TRD.DebugInfoGen.root");\r
 757     savedir->cd();\r
 758   }\r
 759   return fDebugStream;\r
 760 }\r
 761 \r
 762 \r
 763 \r