PWGPP/TRD/AliTRDonlineTrackletQA.cxx

   1 // AliTRDonlineTrackletQA implements the standard QA for the TRD
   2 // on-line tracklets. It reads the tracklets provided by
   3 // AliTRDonlineTrackletFilter and analyses them. The QA can be run
   4 // on tracklets from simulation or from real data. In addition it is
   5 // possible to compare the tracklets from real data to re-simulation
   6 // based on the ADC data.
   7
   8 #include "TFile.h"
   9 #include "TTree.h"
  10 #include "TChain.h"
  11 #include "TH1F.h"
  12 #include "TH2F.h"
  13 #include "TProfile.h"
  14 #include "TCanvas.h"
  15 #include "Math/Factory.h"
  16
  17 #include "AliAnalysisManager.h"
  18 #include "AliESDEvent.h"
  19 #include "AliAODEvent.h"
  20 #include "AliMCEvent.h"
  21 #include "AliESDInputHandler.h"
  22 #include "AliAODHandler.h"
  23 #include "AliMCEventHandler.h"
  24 #include "AliLog.h"
  25 #include "AliESDTrdTrack.h"
  26 #include "AliESDTrdTracklet.h"
  27 #include "AliTRDtrackletMCM.h"
  28 #include "AliTRDtrackletWord.h"
  29 #include "AliTRDtrapConfig.h"
  30 #include "AliTRDmcmSim.h"
  31 #include "AliTRDcalibDB.h"
  32
  33 #include "AliVParticle.h"
  34 #include "AliMCParticle.h"
  35 #include "AliTrackReference.h"
  36 #include "AliMagF.h"
  37 #include "TGeoGlobalMagField.h"
  38
  39 #include "AliTRDonlineTrackletQA.h"
  40
  41 ClassImp(AliTRDonlineTrackletQA)
  42
  43 AliTRDonlineTrackletQA::AliTRDonlineTrackletQA(const char *name) :
  44   AliAnalysisTask(name, ""),
  45   fESD(0x0),
  46   fInputHandler(0x0),
  47   fInputEvent(0x0),
  48   fOutputAOD(0x0),
  49   fMCEvent(0x0),
  50   fTrackletsRaw(0x0),
  51   fTrackletsSim(0x0),
  52   fOutputList(0x0),
  53   fHistYpos(0x0),
  54   fHistYres(0x0),
  55   fHistYresDy(0x0),
  56   fHistdY(0x0),
  57   fHistdYres(0x0),
  58   fHistYresESD(0x0),
  59   fHistdYresESD(0x0),
  60   fHistCanddY(0x0),
  61   fHistFounddY(0x0),
  62   fHistTrklPerRef(0x0),
  63   fHistdYdYref(0x0),
  64   fHistYposRaw(0x0),
  65   fHistdYRaw(0x0),
  66   fHistAlphaRaw(0x0),
  67   fHistYdYRaw(0x0),
  68   fHistZrow(0x0),
  69   fHistZrowRaw(0x0),
  70   fHistPid(0x0),
  71   fHistPidRaw(0x0),
  72   fHistPidDiff(0x0),
  73   fHistYdiff(0x0),
  74   fHistdYdiff(0x0),
  75   fHistdYdYraw(0x0),
  76   fHistFitYres(0x0),
  77   fHistFitDyresEven(0x0),
  78   fHistFitDyresOdd(0x0),
  79   fHistNoMatchSim(0x0),
  80   fHistNoMatchRaw(0x0),
  81   fHistResY(0x0),
  82   fHistResZ(0x0),
  83   fTreeTracklets(0x0),
  84   fY(0.),
  85   fDY(0.),
  86   fYdiff(0.),
  87   fDYdiff(0.),
  88   fQ0(0),
  89   fQ1(0),
  90   fNHits(0),
  91   fMinPt(1.),
  92   fGeo(new AliTRDgeometry),
  93   fNevent(0),
  94   fTrackletTree(0x0)
  95 {
  96   // ctor
  97
  98   DefineInput(0, TChain::Class());
  99   DefineInput(1, TTree::Class());
 100
 101   DefineOutput(0, TTree::Class());
 102   DefineOutput(1, TList::Class());
 103 }
 104
 105 AliTRDonlineTrackletQA::~AliTRDonlineTrackletQA()
 106 {
 107   // dtor
 108
 109   delete fGeo;
 110 }
 111
 112 void AliTRDonlineTrackletQA::ConnectInputData(const Option_t * /* option */)
 113 {
 114   // connect input data
 115
 116   fInputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
 117   if (fInputHandler)
 118     fInputEvent = fInputHandler->GetEvent();
 119
 120   AliMCEventHandler *mcH = (AliMCEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler();
 121   if (mcH)
 122     fMCEvent = mcH->MCEvent();
 123 }
 124
 125 void AliTRDonlineTrackletQA::CreateOutputObjects()
 126 {
 127   // create output objects
 128
 129   OpenFile(1);
 130
 131   fOutputList = new TList();
 132   fOutputList->SetOwner(kTRUE);
 133
 134   fHistYpos        = new TH1F("ypos",
 135                               "Tracklet (sim) y-position;y (cm);count",
 136                               8192/32, -4096*160e-4, 4095*160e-4);
 137   fHistYposRaw     = new TH1F("ypos_raw",
 138                               "Tracklet (raw) y-position;y (cm);count",
 139                               8192/32, -4096*160e-4, 4095*160e-4);
 140   fHistYres        = new TH1F("yres",
 141                               "Tracklet (sim) #Deltay;y_{tracklet}-y_{MC} (cm);count",
 142                               8192/32, -4096/32*160e-4, 4095/32*160e-4);
 143   fHistYresDy      = new TH2F("yresdy",
 144                               "Tracklet (sim) #Deltay;y_{tracklet}-y_{MC} (cm);deflection (bin)",
 145                               8192/32, -4096/32*160e-4, 4095/32*160e-4,
 146                               128, -64.5, 63.5);
 147   fHistYresESD     = new TH1F("yresesd",
 148                               "Tracklet #Deltay;y (cm);count",
 149                               100, -10, 10);
 150   fHistYdiff       = new TH1F("ydiff",
 151                               "Tracklet #Deltay (sim - raw);y_{sim}-y_{raw} (160 #mum);count",
 152                               200, -100, 100);
 153   for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
 154     fHistYlocal[iLayer] = new TH2F(Form("ylocal_%i", iLayer),
 155                                    Form("Tracklet local y, layer %i;y_{MC} (pad width);y_{trkl} (pad width)", iLayer),
 156                                    100, -1, 1, 100, -1, 1);
 157   }
 158
 159   fHistdY          = new TH1F("dy",
 160                               "deflection (sim);dy (140 #mum)",
 161                               128, -64.5, 63.5);
 162   fHistdYRaw       = new TH1F("dy_raw",
 163                               "deflection (raw);dy (140 #mum)",
 164                               128, -64.5, 63.5);
 165   fHistAlphaRaw    = new TH1F("fHistAlphaRaw",
 166                               "angle w.r.t. to straight line",
 167                               256, -128.5, 127.5);
 168   fHistdYres       = new TH1F("dyres",
 169                               "deflection residual;dy (cm)",
 170                               128, -1., 1.);
 171   fHistdYresESD    = new TH1F("dyresesd",
 172                               "deflection residual;dy (cm)",
 173                               128, -1., 1.);
 174   fHistCanddY      = new TH1F("dycand",
 175                               "deflection;dy (140 #mum)",
 176                               128, -64.5, 63.5);
 177   fHistFounddY     = new TH1F("dyfound",
 178                               "deflection;dy (140 #mum)",
 179                               128, -64.5, 63.5);
 180   fHistdYdiff      = new TH1F("dydiff",
 181                               "deflection #Deltady;dy_{sim}-dy_{raw} (140 #mum)",
 182                               100, -2., 2.);
 183   fHistdYdYraw     = new TH2F("dydyraw",
 184                               "deflection from sim. vs raw;dy_{sim} (140 #mum);dy_{raw} (140 #mum)",
 185                               128, -64.5, 63.5, 128, -64.5, 63.5);
 186
 187   fHistTrklPerRef  = new TH1F("trklperref",
 188                               "No. of tracklets per track reference;no. of tracklets",
 189                               10, -0.5, 9.5);
 190
 191   fHistdYdYref     = new TH2F("dydyref",
 192                               "deflection vs. deflection from track reference;dy_{ref} (140 #mum);dy (140 #mum)",
 193                               128, -64.5, 63.5, 128, -64.5, 63.5);
 194
 195   fHistZrow        = new TH1F("zrow",
 196                               "z-position;pad row",
 197                               16, -0.5, 15.5);
 198   fHistZrowRaw     = new TH1F("zrow-raw",
 199                               "z-position;pad row",
 200                               16, -0.5, 15.5);
 201
 202   fHistPid         = new TH1F("pid",
 203                               "pid",
 204                               256, -0.5, 255.5);
 205   fHistPidRaw      = new TH1F("pid-raw",
 206                               "pid",
 207                               256, -0.5, 255.5);
 208   fHistPidDiff     = new TH1F("piddiff",
 209                               "piddiff",
 210                               256, -127.5, 128.5);
 211
 212   fHistYdYRaw      = new TH2F("ydyraw",
 213                               "y vs dy (raw tracklets);y (cm);dy (140 #mum)",
 214                               8192/32, -4096*160e-4, 4095*160e-4,
 215                               128, -64.5, 63.5);
 216
 217   fHistFitYres     = new TH1F("fityres",
 218                               "Tracklet #Deltay;y_{tracklet}-y_{track} (cm);count",
 219                               8192/32, -4096/32*160e-4, 4095/32*160e-4);
 220
 221   fHistFitDyresEven= new TH1F("fitdyreseven",
 222                               "Tracklet #Deltady;dy_{tracklet}-dy_{track} (cm);count",
 223                               100, -2., 2.);
 224   fHistFitDyresOdd = new TH1F("fitdyresodd",
 225                               "Tracklet #Deltady;dy_{tracklet}-dy_{track} (cm);count",
 226                               100, -2., 2.);
 227
 228   fHistNoMatchSim  = new TH2F("nomatchsim",
 229                               "Unmatched tracklets from Simulation",
 230                               8192/32, -4096*160e-4, 4095*160e-4,
 231                               540, -0.5, 539.5);
 232   fHistNoMatchRaw  = new TH2F("nomatchraw",
 233                               "Unmatched tracklets from raw data",
 234                               8192/32, -4096*160e-4, 4095*160e-4,
 235                               540, -0.5, 539.5);
 236
 237   fHistResY        = new TH1F("resy",
 238                               "Residuals to GTU track",
 239                               100, -10., 10.);
 240   fHistResZ        = new TH1F("resz",
 241                               "Residuals to GTU track",
 242                               100, -20., 20.);
 243
 244   fTreeTracklets   = new TTree("trkl", "trkl");
 245   fTreeTracklets->Branch("y", &fY);
 246   fTreeTracklets->Branch("dy", &fDY);
 247   fTreeTracklets->Branch("ydiff", &fYdiff);
 248   fTreeTracklets->Branch("dydiff", &fDYdiff);
 249   fTreeTracklets->Branch("q0", &fQ0);
 250   fTreeTracklets->Branch("q1", &fQ1);
 251   fTreeTracklets->Branch("nhits", &fNHits);
 252
 253   fOutputList->Add(fHistYpos);
 254   fOutputList->Add(fHistdY);
 255   fOutputList->Add(fHistZrow);
 256   fOutputList->Add(fHistPid);
 257
 258   fOutputList->Add(fHistYres);
 259   fOutputList->Add(fHistYresDy);
 260   fOutputList->Add(fHistCanddY);
 261   fOutputList->Add(fHistFounddY);
 262   fOutputList->Add(fHistTrklPerRef);
 263   fOutputList->Add(fHistdYres);
 264   fOutputList->Add(fHistYresESD);
 265   fOutputList->Add(fHistdYresESD);
 266   fOutputList->Add(fHistdYdYref);
 267
 268   for (Int_t iLayer = 0; iLayer < 6; iLayer++)
 269     fOutputList->Add(fHistYlocal[iLayer]);
 270
 271   fOutputList->Add(fHistYposRaw);
 272   fOutputList->Add(fHistdYRaw);
 273   fOutputList->Add(fHistAlphaRaw);
 274   fOutputList->Add(fHistZrowRaw);
 275   fOutputList->Add(fHistPidRaw);
 276   fOutputList->Add(fHistYdYRaw);
 277
 278   fOutputList->Add(fHistYdiff);
 279   fOutputList->Add(fHistdYdiff);
 280   fOutputList->Add(fHistPidDiff);
 281   fOutputList->Add(fHistdYdYraw);
 282
 283   fOutputList->Add(fHistFitYres);
 284   fOutputList->Add(fHistFitDyresEven);
 285   fOutputList->Add(fHistFitDyresOdd);
 286
 287   fOutputList->Add(fHistNoMatchSim);
 288   fOutputList->Add(fHistNoMatchRaw);
 289
 290   fOutputList->Add(fHistResY);
 291   fOutputList->Add(fHistResZ);
 292
 293   fOutputList->Add(fTreeTracklets);
 294
 295   PostData(1, fOutputList);
 296 }
 297
 298 void AliTRDonlineTrackletQA::Exec(const Option_t * /* option */)
 299 {
 300   // execute this for each event
 301
 302   // connect input data
 303   fTrackletTree = (TTree*) GetInputData(1);
 304   if (fTrackletTree) {
 305     fTrackletTree->SetBranchAddress("tracklets_sim", &fTrackletsSim);
 306     fTrackletTree->SetBranchAddress("tracklets_raw", &fTrackletsRaw);
 307     fTrackletTree->GetEntry(fTrackletTree->GetEntriesFast()-1);
 308   }
 309
 310   fESD = dynamic_cast<AliESDEvent*> (fInputEvent);
 311
 312   TList trackletsSim[540];
 313   TList trackletsRaw[540];
 314
 315   // prepare raw tracklets for comparison
 316   Int_t detRaw;
 317   Int_t robRaw;
 318   Int_t mcmRaw;
 319   Int_t yRaw;
 320   Int_t dyRaw;
 321   TTree trklRaw("raw tracklets", "raw tracklets");
 322   trklRaw.Branch("det", &detRaw);
 323   trklRaw.Branch("rob", &robRaw);
 324   trklRaw.Branch("mcm", &mcmRaw);
 325   trklRaw.Branch("y", &yRaw);
 326   trklRaw.Branch("dy", &dyRaw);
 327   trklRaw.SetDirectory(0x0);
 328   // prepare simulated tracklets for comparison
 329   Int_t detSim;
 330   Int_t robSim;
 331   Int_t mcmSim;
 332   Int_t ySim;
 333   Int_t dySim;
 334   TTree trklSim("sim tracklets", "sim tracklets");
 335   trklSim.Branch("det", &detSim);
 336   trklSim.Branch("rob", &robSim);
 337   trklSim.Branch("mcm", &mcmSim);
 338   trklSim.Branch("y", &ySim);
 339   trklSim.Branch("dy", &dySim);
 340   trklSim.SetDirectory(0x0);
 341
 342   // ----- read ESD tracklets -----
 343   if (fESD) {
 344     Int_t nTracklets = fESD->GetNumberOfTrdTracklets();
 345     Int_t nTrackletsSim = 0;
 346     for (Int_t iTracklet = 0; iTracklet < nTracklets; iTracklet++) {
 347       AliESDTrdTracklet *tracklet = fESD->GetTrdTracklet(iTracklet);
 348
 349       if (tracklet) {
 350         // check for tracklet from simulation (label >= -1)
 351         if (tracklet->GetLabel() >= -1) {
 352           fHistYpos->Fill(tracklet->GetLocalY());
 353           fHistdY->Fill(tracklet->GetBinDy());
 354           fHistZrow->Fill(tracklet->GetBinZ());
 355           fHistPid->Fill(tracklet->GetPID());
 356
 357           detSim = tracklet->GetDetector();
 358           robSim = tracklet->GetROB();
 359           mcmSim = tracklet->GetMCM();
 360           ySim   = tracklet->GetBinY();
 361           dySim  = tracklet->GetBinDy();
 362           trklSim.Fill();
 363
 364           AliDebug(1, Form("trkl sim 0x%08x in %4i (ESD)", tracklet->GetTrackletWord(), tracklet->GetHCId()));
 365           ++nTrackletsSim;
 366
 367           if (tracklet->GetDetector() > -1) {
 368             trackletsSim[tracklet->GetDetector()].Add(tracklet);
 369           }
 370
 371           PlotMC(tracklet);
 372           PlotESD(tracklet);
 373         }
 374         // otherwise raw tracklet
 375         else {
 376           if (tracklet->GetDetector() > -1)
 377             trackletsRaw[tracklet->GetDetector()].Add(tracklet);
 378
 379           fHistYposRaw->Fill(tracklet->GetLocalY());
 380           fHistdYRaw->Fill(tracklet->GetBinDy());
 381           Float_t alpha = tracklet->GetBinDy() - 3./140.e-4*tracklet->GetLocalY()/GetX(tracklet);
 382           fHistAlphaRaw->Fill(alpha);
 383           fHistZrowRaw->Fill(tracklet->GetBinZ());
 384           fHistPidRaw->Fill(tracklet->GetPID());
 385           fHistYdYRaw->Fill(tracklet->GetLocalY(), tracklet->GetBinDy());
 386
 387           detRaw = tracklet->GetDetector();
 388           robRaw = tracklet->GetROB();
 389           mcmRaw = tracklet->GetMCM();
 390           yRaw   = tracklet->GetBinY();
 391           dyRaw  = tracklet->GetBinDy();
 392           trklRaw.Fill();
 393         }
 394       }
 395     }
 396     AliDebug(1, Form("no. of simulated tracklets in ESDs: %i", nTrackletsSim));
 397   }
 398
 399   // ----- simulated tracklet from TRD.Tracklets.root -----
 400   if (fTrackletsSim) {
 401     AliDebug(1, Form("no. of simulated tracklets in TRD.Tracklets.root: %i", fTrackletsSim->GetEntries()));
 402     for (Int_t iTracklet = 0; iTracklet < fTrackletsSim->GetEntries(); iTracklet++) {
 403       AliTRDtrackletMCM *trkl = (AliTRDtrackletMCM*) (*fTrackletsSim)[iTracklet];
 404       AliDebug(1, Form("trkl sim 0x%08x in %4i (TRD)", trkl->GetTrackletWord(), trkl->GetHCId()));
 405     }
 406   }
 407
 408   // ----- match simulated and raw tracklets (detector-wise) -----
 409   for (Int_t iDetector = 0; iDetector < 540; iDetector++) {
 410     for (Int_t iTracklet = 0; iTracklet < trackletsRaw[iDetector].GetEntries(); iTracklet++) {
 411       AliESDTrdTracklet *trackletRaw = (AliESDTrdTracklet*) trackletsRaw[iDetector].At(iTracklet);
 412
 413       // search for the best match
 414       AliESDTrdTracklet *trackletMatch = 0x0;
 415       Int_t dyMax = 20;
 416       Int_t yMax = 100;
 417       for (Int_t iTrackletSim = 0; iTrackletSim < trackletsSim[iDetector].GetEntries(); iTrackletSim++) {
 418         AliESDTrdTracklet *trackletSim = (AliESDTrdTracklet*) trackletsSim[iDetector].At(iTrackletSim);
 419         AliDebug(10, Form("comparing 0x%08x with 0x%08x, det: %i %i",
 420                           trackletRaw->GetTrackletWord(), trackletSim->GetTrackletWord(),
 421                           trackletRaw->GetDetector(), trackletSim->GetDetector() ));
 422         // require equal z
 423         if (trackletRaw->GetBinZ() != trackletSim->GetBinZ())
 424           continue;
 425         // match on y
 426         if (TMath::Abs(trackletRaw->GetBinY() - trackletSim->GetBinY()) > yMax)
 427           continue;
 428
 429         trackletMatch = trackletSim;
 430         yMax  = TMath::Abs(trackletRaw->GetBinY()  - trackletSim->GetBinY());
 431         dyMax = TMath::Abs(trackletRaw->GetBinDy() - trackletSim->GetBinDy());
 432       }
 433       if (trackletMatch) {
 434         fHistYdiff->Fill(trackletRaw->GetBinY()  - trackletMatch->GetBinY());
 435         fHistdYdiff->Fill(3.*trackletRaw->GetDyDx() - 3.*trackletMatch->GetDyDx());
 436         fHistdYdYraw->Fill(trackletMatch->GetBinDy(), trackletRaw->GetBinDy());
 437         fHistPidDiff->Fill(trackletRaw->GetPID() - trackletMatch->GetPID());
 438         trackletsSim[iDetector].Remove(trackletMatch);
 439       }
 440       else {
 441         // store unmatched raw tracklets
 442         fHistNoMatchRaw->Fill(trackletRaw->GetLocalY(), trackletRaw->GetDetector());
 443         // AliInfo(Form("unmatched raw tracklet: 0x%08x", trackletRaw->GetTrackletWord()));
 444       }
 445     }
 446     // store the unmatched sim tracklets
 447     for (Int_t iTracklet = 0; iTracklet < trackletsSim[iDetector].GetEntries(); iTracklet++) {
 448       AliESDTrdTracklet *tracklet = (AliESDTrdTracklet*) trackletsSim[iDetector].At(iTracklet);
 449       fHistNoMatchSim->Fill(tracklet->GetLocalY(), tracklet->GetDetector());
 450       // AliInfo(Form("unmatched simulated tracklet: 0x%08x", tracklet->GetTrackletWord()));
 451     }
 452   }
 453
 454   // ----- MC tracks and track references -----
 455   // determine tracklet efficiency
 456   if (fMCEvent) {
 457     Int_t nTracksMC = fMCEvent->GetNumberOfTracks();
 458     // AliInfo(Form("no. of MC tracks: %i", nTracksMC));
 459     for (Int_t iTrack = 0; iTrack < nTracksMC; iTrack++) {
 460       // we want primaries
 461       if (!fMCEvent->IsPhysicalPrimary(iTrack))
 462         continue;
 463
 464       AliMCParticle *mcpart = (AliMCParticle*) fMCEvent->GetTrack(iTrack);
 465
 466       // don't look at tracks with too low pt
 467       if (TMath::Abs(mcpart->Pt()) < fMinPt)
 468         continue;
 469
 470       // look for two track references in a chamber
 471       Int_t nTrackRefs = mcpart->GetNumberOfTrackReferences();
 472       AliTrackReference *tr[2] = { 0x0 };
 473       Int_t nRef = 0;
 474
 475       for (Int_t iTrackRef = 0; iTrackRef < nTrackRefs; iTrackRef++) {
 476         AliTrackReference *trackRef = mcpart->GetTrackReference(iTrackRef);
 477         if (trackRef->DetectorId() != AliTrackReference::kTRD)
 478           continue;
 479         if (trackRef->Pt() < fMinPt)
 480           continue;
 481         Int_t label = trackRef->Label();
 482         if (label < 0)
 483           continue;
 484
 485         // first track reference, remember it
 486         if (nRef == 0) {
 487           tr[nRef] = trackRef;
 488           nRef++;
 489           continue;
 490         }
 491         else {
 492           // next one is too far away, remember it but forget the previous one
 493           if (TMath::Abs(trackRef->LocalX() - tr[0]->LocalX()) > 5.) {
 494             tr[0] = trackRef;
 495             nRef = 1;
 496             continue;
 497           }
 498           // too close to previous track reference
 499           // we don't want it
 500           else if (TMath::Abs(trackRef->LocalX() - tr[0]->LocalX()) < .5) {
 501             continue;
 502           }
 503           // then it must be ok, so we take it
 504           else {
 505             tr[1] = trackRef;
 506             nRef++;
 507           }
 508         }
 509
 510         // calculation deflection from track references
 511         Float_t deflLength = 3. * (tr[1]->LocalY() - tr[0]->LocalY()) / (tr[1]->LocalX() - tr[0]->LocalX());
 512         // if it is too large we reject it
 513         if (deflLength < 1.) {
 514           fHistCanddY->Fill(deflLength/140e-4);
 515         }
 516         else {
 517           nRef = 0;
 518           continue;
 519         }
 520
 521         // now search for tracklets belonging to this track reference
 522         Int_t nTrackletsPerRef(0);
 523         Int_t defl(0);
 524         TIter nextTracklet(trackletsSim);
 525         while (AliESDTrdTracklet *trkl = (AliESDTrdTracklet*) nextTracklet()) {
 526           // they must have the same label
 527           if (trkl->GetLabel() != label)
 528             continue;
 529
 530           // and be close enough in radial position
 531           if (TMath::Abs(trackRef->LocalX() - GetX(trkl)) > 5.)
 532             continue;
 533
 534           // if they are close in position we accept it
 535           if ((TMath::Abs(trackRef->LocalY() - trkl->GetLocalY()) < 5.) &&
 536               (TMath::Abs(trackRef->Z() - GetZ(trkl)) < 5.)) {
 537             defl = trkl->GetDyDx();
 538             nTrackletsPerRef++;
 539           }
 540         }
 541         fHistTrklPerRef->Fill(nTrackletsPerRef);
 542         if (nTrackletsPerRef == 0) {
 543           AliInfo(Form("Track ref without assigned tracklet: x=%4.2f, y=%4.2f, z=%4.2f, pt=%4.2f (%i)",
 544                        trackRef->X(), trackRef->Y(), trackRef->Z(), trackRef->Pt(), trackRef->Label()));
 545         }
 546         if (nTrackletsPerRef == 1) {
 547           fHistdYdYref->Fill(deflLength/140.e-4, defl);
 548           fHistFounddY->Fill(deflLength/140.e-4);
 549         }
 550         nRef = 0;
 551       }
 552     }
 553   }
 554
 555   // ----- ESD tracks -----
 556   if (fESD) {
 557     for (Int_t iTrack = 0; iTrack < fESD->GetNumberOfTracks(); iTrack++) {
 558       AliESDtrack *esdTrack = fESD->GetTrack(iTrack);
 559       AliDebug(1, Form("ESD track pt: %7.2f", esdTrack->Pt()));
 560     }
 561   }
 562
 563   // ----- ESD TRD tracks -----
 564   if (fESD) {
 565     // AliTRDtrapConfig *trapcfg = AliTRDcalibDB::Instance()->GetTrapConfig();
 566     AliDebug(1, Form("no. of TRD tracks: %i", fESD->GetNumberOfTrdTracks()));
 567     for (Int_t iTrack = 0; iTrack < fESD->GetNumberOfTrdTracks(); iTrack++) {
 568       AliESDTrdTrack *trdTrack = fESD->GetTrdTrack(iTrack);
 569       AliInfo(Form("TRD track pt: %7.2f", trdTrack->Pt()));
 570       // AliTRDtrackOnline track;
 571       AliESDTrdTracklet tracklets[6];
 572       for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
 573         AliESDTrdTracklet *trkl = trdTrack->GetTracklet(iLayer);
 574         if (trkl) {
 575           tracklets[iLayer].SetTrackletWord(trkl->GetTrackletWord());
 576           tracklets[iLayer].SetHCId(trkl->GetHCId());
 577           // track.AddTracklet(&(tracklets[iLayer])); //???
 578         }
 579       }
 580       // track.AddParametrization(new AliTRDtrackParametrizationCurved());
 581       // ROOT::Math::Minimizer *minim = ROOT::Math::Factory::CreateMinimizer("Minuit");
 582       // track.Fit(minim);
 583       // const TList &paramList = track.GetParametrizations();
 584       // TIter param(&paramList, kFALSE);
 585       // while (AliTRDtrackParametrization *trackParam = (AliTRDtrackParametrization*) param()) {
 586       //        trackParam->Print();
 587       //        AliTRDtrackParametrizationCurved *curvedParam = dynamic_cast<AliTRDtrackParametrizationCurved*> (trackParam);
 588       //        if (curvedParam) {
 589       //          // calculate residuals
 590       //          for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
 591       //            AliESDTrdTracklet *trkl = trdTrack->GetTracklet(iLayer);
 592       //            if (trkl) {
 593       //              AliTRDtrackPosition pos = curvedParam->ExtrapolateToLayer(iLayer);
 594       //              fHistFitYres->Fill(pos.GetY() - trkl->GetLocalY());
 595       //              Int_t rob = (trkl->GetBinZ()/4) + (trkl->GetHCId() % 0x1 ? 1 : 0);
 596       //              Int_t mcm = trkl->GetBinZ() & 0x3;
 597       //              Float_t deflcorr = ((Int_t) trapcfg->GetDmemUnsigned(AliTRDmcmSim::fgkDmemAddrDeflCorr, trkl->GetDetector(), rob, mcm))
 598       //                * (0.635 + 0.03 * iLayer)/(256.0) / 32.;
 599       //              Float_t scale = 1.;
 600       //              Float_t defl = (3.*trkl->GetDyDx() - deflcorr) * scale + deflcorr;
 601       //              if (iLayer & 0x1) {
 602       //                fHistFitDyresOdd->Fill(pos.GetdY() - defl);
 603       //              }
 604       //              else {
 605       //                fHistFitDyresEven->Fill(pos.GetdY() - defl);
 606       //              }
 607       //            }
 608       //          }
 609       //        }
 610       // }
 611     }
 612   }
 613
 614   PostData(1, fOutputList);
 615 }
 616
 617 void AliTRDonlineTrackletQA::LocalInit()
 618 {
 619
 620 }
 621
 622 void AliTRDonlineTrackletQA::Terminate(const Option_t * /* option */)
 623 {
 624   // upon terminate
 625
 626   fOutputList = dynamic_cast<TList*> (GetOutputData(1));
 627
 628   if (!fOutputList) {
 629     AliError("No output objects found!");
 630     return;
 631   }
 632
 633 //  fHistYpos = dynamic_cast<TH1F*> (fOutputList->FindObject("ypos"));
 634 //  if (fHistYpos) {
 635 //    TCanvas *c = new TCanvas;
 636 //    c->cd();
 637 //    gPad->SetGridx();
 638 //    gPad->SetGridy();
 639 //    fHistYpos->DrawCopy();
 640 //  }
 641 }
 642
 643 void AliTRDonlineTrackletQA::PlotMC(AliESDTrdTracklet *trkl)
 644 {
 645   // compare the given tracklet to the MC information,
 646   // i.e. track references
 647
 648   Int_t label = trkl->GetLabel();
 649
 650   // check for valid label
 651   if (label < 0 ) {
 652     AliDebug(1, "MC tracklet has no label");
 653     return;
 654   }
 655   if (label >= fMCEvent->GetNumberOfTracks()) {
 656     AliError("MC tracklet has invalid label");
 657     return;
 658   }
 659
 660   // some cuts on the tracklet
 661 //  if (!fMCEvent->IsPhysicalPrimary(label))
 662 //    return;
 663 //  if (TMath::Abs(trkl->GetdYdX()) > 0.05)
 664 //    return;
 665 //  if (trkl->GetDetector() % 6 != 5)
 666 //    return;
 667
 668   // get MC particle for this tracklet
 669   AliMCParticle *track = (AliMCParticle*) fMCEvent->GetTrack(label);
 670
 671   // don't look at tracks with too low pt
 672   if (TMath::Abs(track->Pt() < fMinPt))
 673     return;
 674
 675   // serach track references corresponding to the current tracklet
 676   AliTrackReference *tr[2] = { 0x0 };
 677   Int_t nTrackRefs = 0;
 678
 679   for (Int_t iTrackRef = 0; iTrackRef < track->GetNumberOfTrackReferences(); iTrackRef++) {
 680     AliTrackReference *trackRef = track->GetTrackReference(iTrackRef);
 681     if (trackRef->DetectorId() != AliTrackReference::kTRD)
 682       continue;
 683     if (trackRef->Pt() < fMinPt)
 684       continue;
 685
 686     if (TMath::Abs(GetX(trkl) - trackRef->LocalX()) > 5.)
 687       continue;
 688
 689     tr[nTrackRefs++] = trackRef;
 690
 691     if (nTrackRefs == 2)
 692       break;
 693   }
 694
 695   // if there were exactly 2 track references
 696   // (otherwise something is strange and we want to look at clean cases)
 697   // compare tracklet to them
 698   if (nTrackRefs == 2) {
 699     // sanity check to be in the same sector
 700     if ( TMath::Abs((tr[0]->Alpha()*180./TMath::Pi()-10.)/20. - (trkl->GetDetector()/30)) > .1) {
 701       AliError("Track reference in different sector");
 702     }
 703     // require minimal distance in X and maximum deflection in Y
 704     // for the track references
 705     else if ((tr[1]->LocalX() - tr[0]->LocalX()) > 0.1 && TMath::Abs(tr[1]->LocalY() - tr[0]->LocalY()) < 1.) {
 706       // calculate slope from track references
 707       // and check whether it's in the allowed range
 708       Float_t slope = 3. * (tr[1]->LocalY() - tr[0]->LocalY()) / (tr[1]->LocalX() - tr[0]->LocalX());
 709       if (TMath::Abs(slope) < 64*140e-4) {
 710         AliDebug(1,Form("x1: %f, x0: %f, y1: %f, y0:%f",
 711                         tr[1]->LocalX(), tr[0]->LocalX(), tr[1]->LocalY(), tr[0]->LocalY() ));
 712         // calculate y-position scaled to radial position of the tracklet
 713         // and consider the tilting angle of the pads
 714         // since the tracklets are affected by it
 715         Float_t yMC     = (tr[1]->LocalY() + (-0.5+GetX(trkl) - tr[1]->LocalX()) * (tr[1]->LocalY() - tr[0]->LocalY()) / (tr[1]->LocalX() - tr[0]->LocalX()));
 716         Float_t yMCtilt = yMC + (TMath::Tan(TMath::Power(-1, (trkl->GetDetector() % 6))*2.*TMath::Pi()/180.) * (tr[1]->Z() - GetZ(trkl)));
 717         if (TMath::Abs(trkl->GetLocalY() - yMCtilt) > 10.) {
 718           AliError(Form("Deviation too large for tracklet: 0x%08x in det. %i at x = %f, y = %f, z = %f, alpha = %f",
 719                         trkl->GetTrackletWord(), trkl->GetDetector(), GetX(trkl), trkl->GetLocalY(), GetZ(trkl), tr[0]->Alpha()));
 720         }
 721         fHistYres->Fill(trkl->GetLocalY() - yMCtilt);
 722         fHistYresDy->Fill(trkl->GetLocalY() - yMCtilt, trkl->GetDyDx());
 723         // what about tilt correction here ???
 724         fHistdYres->Fill(3. * trkl->GetDyDx() -
 725                          3. * (tr[1]->LocalY() - tr[0]->LocalY()) / (tr[1]->LocalX() - tr[0]->LocalX()));
 726         // plot position deviation in pad-coordinates
 727         // to study the influence of the position LUT
 728         Float_t padWidth = fGeo->GetPadPlane(trkl->GetDetector())->GetWidthIPad();
 729         Float_t yMClocal = yMCtilt/padWidth - floor(yMCtilt/padWidth) - padWidth/2.;
 730         Int_t layer = trkl->GetDetector() % 6;
 731         fHistYlocal[layer]->Fill(yMClocal,
 732                                  trkl->GetLocalY()/padWidth - floor(trkl->GetLocalY()/padWidth) - padWidth/2. - yMClocal);
 733         // and fill everything to the tree
 734         fYdiff = trkl->GetLocalY() - yMCtilt;
 735         fDYdiff = 3. * trkl->GetDyDx() -
 736           3. * (tr[1]->LocalY() - tr[0]->LocalY()) / (tr[1]->LocalX() - tr[0]->LocalX());
 737         fY = trkl->GetLocalY();
 738         fDY = trkl->GetDyDx();
 739         fTreeTracklets->Fill();
 740         // output tracklets with large deviation
 741         if (TMath::Abs(fYdiff) > 0.5) {
 742           AliWarning(Form("tracklet: y=%4.2f, dy=%4.2f, ydiff=%4.2f, dydiff=%4.2f, q0=%5d, q1=%5d, nhits=%2d, label=%i",
 743                           trkl->GetLocalY(), trkl->GetDyDx(), fYdiff, fDYdiff, fQ0, fQ1, fNHits, label));
 744         }
 745       }
 746     }
 747   }
 748 }
 749
 750
 751 void AliTRDonlineTrackletQA::PlotESD(AliESDTrdTracklet *trkl)
 752 {
 753   // plot comparison to ESD
 754
 755   Float_t xTrkl = GetX(trkl);
 756   Float_t yTrkl = trkl->GetLocalY();
 757   Float_t zTrkl = GetZ(trkl);
 758
 759   Float_t alpha = (trkl->GetDetector() / 30) * 20. + 10.;
 760   alpha *= TMath::Pi() / 180.;
 761
 762   AliESDEvent *esdEvent = dynamic_cast<AliESDEvent*> (fInputEvent);
 763   if (!esdEvent)
 764     return;
 765
 766   Float_t mag = 0.5; //((AliMagF*) TGeoGlobalMagField::Instance()->GetField())->SolenoidField();
 767
 768   for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++) {
 769     AliESDtrack *track = esdEvent->GetTrack(iTrack);
 770
 771     if (!track->GetOuterParam())
 772       continue;
 773
 774     AliExternalTrackParam *param = new AliExternalTrackParam(*(track->GetOuterParam()));
 775
 776     AliDebug(10, Form("track %i at x = %f, y = %f",
 777                       iTrack, param->GetX(), param->GetY()));
 778     param->Propagate(alpha, xTrkl, mag);
 779     AliDebug(10, Form("after propagating track %i at x = %f, y = %f",
 780                       iTrack, param->GetX(), param->GetY()));
 781
 782     if ((TMath::Abs(xTrkl - param->GetX()) < 10.) &&
 783         (TMath::Abs(yTrkl - param->GetY()) < 5.) &&
 784         (TMath::Abs(zTrkl - param->GetZ()) < 10.)) {
 785       AliDebug(2, Form("match of tracklet-track: %i <-> %i",
 786                        trkl->GetLabel(), track->GetLabel()));
 787       AliDebug(5, Form("tracklet position: det: %3i  x = %f, y = %f, z = %f, alpha = %f",
 788                        trkl->GetDetector(), GetX(trkl), trkl->GetLocalY(), GetZ(trkl), alpha));
 789       AliDebug(5, Form("after propagating track %i at x = %f, y = %f, z = %f",
 790                       iTrack, param->GetX(), param->GetY(), param->GetZ()));
 791
 792       fHistYresESD->Fill(yTrkl - param->GetY());
 793     }
 794
 795     delete param;
 796   }
 797
 798 }
 799
 800 Int_t AliTRDonlineTrackletQA::GetTrackletsForMC(Int_t /* label */, Int_t /*idx*/ []) const
 801 {
 802   // get tracklets for MC label
 803   // not implemented
 804
 805   return 0;
 806 }