ITS/UPGRADE/testITSU/AliTaskITSUPerf.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 ///////////////////////////////////////////////////////////////////////////
  17 // Class AliTaskITSUPerf                                                 //
  18 // Analysis task to produce data and MC histos needed for tracklets      //
  19 // dNdEta extraction in multiple bins in one go                          //
  20 // Author:  ruben.shahoyan@cern.ch                                       //
  21 ///////////////////////////////////////////////////////////////////////////
  22
  23 #include "TFile.h"
  24 #include "TChain.h"
  25 #include "TTree.h"
  26 #include "TRandom.h"
  27 #include "TH1F.h"
  28 #include "TH2F.h"
  29 #include "TList.h"
  30 #include "TObjArray.h"
  31 #include "TGeoGlobalMagField.h"
  32
  33 #include "AliAnalysisManager.h"
  34
  35 #include "AliESDEvent.h"
  36 #include "AliESDInputHandler.h"
  37 #include "AliESDInputHandlerRP.h"
  38 #include "AliCDBPath.h"
  39 #include "AliCDBManager.h"
  40 #include "AliCDBEntry.h"
  41 #include "AliCDBStorage.h"
  42 #include "AliGeomManager.h"
  43 #include "AliMagF.h"
  44 #include "AliMCEventHandler.h"
  45 #include "AliMCEvent.h"
  46 #include "AliMCParticle.h"
  47 #include "AliStack.h"
  48 #include "AliGenEventHeader.h"
  49 #include "AliCentrality.h"
  50 #include "AliLog.h"
  51 #include "AliPhysicsSelection.h"
  52 #include "AliGenEventHeader.h"
  53 #include "AliGenHijingEventHeader.h"
  54 #include "AliGenDPMjetEventHeader.h"
  55 #include "AliGenCocktailEventHeader.h"
  56 #include "AliESDtrackCuts.h"
  57
  58 #include "AliTaskITSUPerf.h"
  59 #include "../ITS/UPGRADE/AliITSURecoDet.h"
  60 #include "../ITS/UPGRADE/AliITSUGeomTGeo.h"
  61 #include "../ITS/UPGRADE/AliITSUClusterPix.h"
  62 #include "../ITS/UPGRADE/AliITSUTrackCond.h"
  63 #include "../ITS/UPGRADE/AliITSUAux.h"
  64 #include <TGeoManager.h>
  65 #include <TTree.h>
  66
  67
  68 ClassImp(AliTaskITSUPerf)
  69
  70
  71 const char* AliTaskITSUPerf::fgkLabelTypes[AliTaskITSUPerf::kNLabelTypes] = {
  72   "ITSokTPCok"
  73   ,"ITSokTPCfk"
  74   ,"ITSfkTPCok"
  75   ,"ITSfkTPCfk"
  76   ,"ITSTPCmismatch"
  77   ,"ITSTPCnoMatch"
  78 };
  79
  80
  81 typedef struct {
  82   Bool_t  prim;
  83   Bool_t  rcbl;
  84   Char_t  nClITS;
  85   Char_t  nClTPC;
  86   Char_t  nClITSMC;
  87   Char_t  mcCl[7];
  88   Char_t  rcCl[7];
  89   Char_t  fcCl[7];
  90   Char_t  charge;
  91   Float_t zV;
  92   Float_t ptMC;
  93   Float_t etaMC;
  94   Float_t pt;
  95   Float_t eta;
  96   Float_t ptTPC;
  97   Float_t etaTPC;
  98   Float_t dcaR;
  99   Float_t dcaZ;
 100   Float_t dcaRE;
 101   Float_t dcaZE;
 102   Float_t ptE;
 103   Float_t ptTPCE;
 104   Float_t phi;
 105   Int_t   pdg;
 106   Int_t   lbl;
 107   Int_t   lblTPC;
 108   Int_t   spl;
 109   Int_t   evID;
 110 } trInfo_t;
 111
 112 trInfo_t trackInfo;
 113
 114 //________________________________________________________________________
 115 /*//Default constructor
 116 AliTaskITSUPerf::AliTaskITSUPerf(const char *name)
 117   : AliAnalysisTaskSE(name),
 118 */
 119 //________________________________________________________________________
 120 AliTaskITSUPerf::AliTaskITSUPerf(const char *name)
 121   : AliAnalysisTaskSE(name)
 122     //
 123   ,fOutput(0)
 124   ,fHistosCentMCLb(0)
 125   ,fHistosCent(0)
 126   ,fRPTree(0)
 127   ,fStack(0)
 128   ,fMCEvent(0)
 129   ,fVtxSPD(0)
 130   ,fVtxTrc(0)
 131   ,fESDEvent(0)
 132   ,fUseSpecialOutput(kTRUE)
 133   ,fUseMC(kTRUE)
 134   ,fTrackingCond(0)
 135     //
 136   ,fGeom(0)
 137   ,fITS(0)
 138   ,fNSelTracksMC(0)
 139   ,fMCStatus(0)
 140     //
 141 #ifdef _CLUS_LIST_
 142   ,fClLists(0)
 143 #endif
 144     //
 145   ,fNPtBins(20)
 146   ,fNResBins(100)
 147   ,fMinTPCclusters(60)
 148   ,fPtMin(0)
 149   ,fPtMax(10.)
 150   ,fEtaMin(-3.)
 151   ,fEtaMax( 3.)
 152   ,fZVertexMin(-15.)
 153   ,fZVertexMax(15.)
 154     //
 155   ,fCurrCentBin(0)
 156   ,fNCentBins(1)
 157   ,fUseCentralityVar(0)
 158   ,fTPCCut(0)
 159   ,fTree(0)
 160 {
 161   // Constructor
 162
 163   DefineOutput(1, TList::Class());
 164   //
 165 }
 166
 167 //________________________________________________________________________
 168 AliTaskITSUPerf::~AliTaskITSUPerf()
 169 {
 170   // Destructor
 171   // histograms are in the output list and deleted when the output
 172   // list is deleted by the TSelector dtor
 173   if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {  //RRR
 174     printf("Deleteing output\n");
 175     delete fOutput;
 176     fOutput = 0;
 177   }
 178   //
 179   fHistosCentMCLb.Clear();
 180   fHistosCent.Clear();
 181   //
 182   delete fITS;
 183   delete fGeom;
 184   delete fTPCCut;
 185   //
 186 }
 187
 188 //________________________________________________________________________
 189 void AliTaskITSUPerf::UserCreateOutputObjects()
 190 {
 191   //
 192   if (fUseSpecialOutput) OpenFile(1);
 193   fOutput = new TList();
 194   fOutput->SetOwner();
 195   //
 196   // init geometry
 197   if (!gGeoManager) {
 198     TString geom = "geometry.root";
 199     AliGeomManager::LoadGeometry(geom.Data());
 200     if (!gGeoManager) AliFatal("Failed to load geometry");
 201   }
 202   //
 203   // create ITS interface
 204   fGeom = new AliITSUGeomTGeo(kTRUE,kTRUE);
 205   AliITSUClusterPix::SetGeom(fGeom);
 206   fITS = new AliITSURecoDet(fGeom,"ITSURecoInterface");
 207   fITS->CreateClusterArrays();
 208   //
 209   // Create histograms
 210   for (int ib=0;ib<fNCentBins;ib++) {
 211     BookHistos(ib);
 212   }
 213   //
 214   int nhist = fOutput->GetEntries();
 215   for (int i=0;i<nhist;i++) {
 216     TObject* hst = fOutput->At(i);
 217     if (!hst || !(hst->InheritsFrom(TH1::Class()))) continue;
 218     ((TH1*)hst)->Sumw2();
 219   }
 220   //
 221   fTree = new TTree("trinfo","trinfo");
 222   fTree->Branch("prim",&trackInfo.prim,"prim/O");
 223   fTree->Branch("rcbl",&trackInfo.rcbl,"rcbl/O");
 224   fTree->Branch("nClITS",&trackInfo.nClITS,"nClITS/b");
 225   fTree->Branch("nClTPC",&trackInfo.nClTPC,"nClTPC/b");
 226   fTree->Branch("nClITSMC",&trackInfo.nClITSMC,"nClITSMC/b");
 227   fTree->Branch("mcCl",&trackInfo.mcCl,"mcCl[7]/b");
 228   fTree->Branch("rcCl",&trackInfo.rcCl,"rcCl[7]/b");
 229   fTree->Branch("fcCl",&trackInfo.fcCl,"fcCl[7]/b");
 230   fTree->Branch("charge",&trackInfo.charge,"charge/B");
 231   fTree->Branch("zV",    &trackInfo.zV,"zV/F");
 232   fTree->Branch("ptMC",  &trackInfo.ptMC,"ptMC/F");
 233   fTree->Branch("etaMC", &trackInfo.etaMC,"etaMC/F");
 234   fTree->Branch("pt",  &trackInfo.pt,"pt/F");
 235   fTree->Branch("eta", &trackInfo.eta,"eta/F");
 236   fTree->Branch("ptTPC",  &trackInfo.ptTPC,"ptTPC/F");
 237   fTree->Branch("etaTPC", &trackInfo.etaTPC,"etaTPC/F");
 238   fTree->Branch("dcaR", &trackInfo.dcaR,"dcaR/F");
 239   fTree->Branch("dcaZ", &trackInfo.dcaZ,"dcaZ/F");
 240   fTree->Branch("dcaRE", &trackInfo.dcaRE,"dcaRE/F");
 241   fTree->Branch("dcaZE", &trackInfo.dcaZE,"dcaZE/F");
 242   fTree->Branch("ptE",  &trackInfo.ptE,"ptE/F");
 243   fTree->Branch("ptTPCE",  &trackInfo.ptTPCE,"ptTPCE/F");
 244   fTree->Branch("phi",  &trackInfo.phi,"phi/F");
 245   fTree->Branch("pdg",  &trackInfo.pdg,"pdg/I");
 246   fTree->Branch("lbl",  &trackInfo.lbl,"lbl/I");
 247   fTree->Branch("lblTPC",  &trackInfo.lblTPC,"lblTPC/I");
 248   fTree->Branch("spl",  &trackInfo.spl,"spl/I");
 249   fTree->Branch("evID",  &trackInfo.evID,"evID/I");
 250   fOutput->Add(fTree);
 251   //
 252   fTPCCut = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 253   fTPCCut->SetEtaRange(fEtaMin,fEtaMax);
 254   fTPCCut->SetPtRange(fPtMin,fPtMax);
 255   fTPCCut->SetPtRange(fPtMin,fPtMax);
 256   fTPCCut->SetMinNClustersTPC(fMinTPCclusters);
 257   printf("Cuts: Eta: %f %f  | Pt: %f %f\n",fEtaMin,fEtaMax,fPtMin,fPtMax);
 258   //
 259 #ifdef _CLUS_LIST_
 260   fClLists.SetOwner();
 261 #endif
 262   PostData(1, fOutput);
 263   //
 264 }
 265
 266 //________________________________________________________________________
 267 void AliTaskITSUPerf::UserExec(Option_t *)
 268 {
 269   // Main loop
 270   static int evID = 0;
 271   AliInfo(Form("Event: %d",evID++));
 272   //
 273   AliAnalysisManager* anMan = AliAnalysisManager::GetAnalysisManager();
 274   fRPTree = 0;
 275   AliESDInputHandlerRP *handRP = (AliESDInputHandlerRP*)anMan->GetInputEventHandler();
 276   if (!handRP) { AliFatal("No RP handler"); return; }
 277   //
 278   fESDEvent  = handRP->GetEvent();
 279   if (!fESDEvent) { AliFatal("No AliESDEvent"); return; }
 280   //
 281   fRPTree = handRP->GetTreeR("ITS");
 282   if (!fRPTree) { AliWarning("Invalid ITS cluster tree"); }
 283   //
 284   AliMCEventHandler* eventHandler = (AliMCEventHandler*)anMan->GetMCtruthEventHandler();
 285   if (!eventHandler) { AliFatal("Could not retrieve MC event handler"); return; }
 286   fMCEvent = eventHandler->MCEvent();
 287   if (!fMCEvent) { AliFatal("Could not retrieve MC event"); return; }
 288   fStack = fMCEvent->Stack();
 289   if (!fStack) { AliFatal("Stack is not available"); return; }
 290   //
 291   // MC Generator info
 292   AliGenEventHeader* mcGenH = 0;
 293   mcGenH = fMCEvent->GenEventHeader();
 294   TArrayF vtmc(3);
 295   mcGenH->PrimaryVertex(vtmc);
 296   for (int i=3;i--;) fVtxMC[i] = vtmc[i];
 297   //
 298   printf("MagField: %f\n",fESDEvent->GetMagneticField());
 299   printf("VtxMC : %+e %+e %+e\n",fVtxMC[0],fVtxMC[1],fVtxMC[2]);
 300   //
 301   fVtxSPD = fESDEvent->GetPrimaryVertexSPD();
 302   //
 303   if (fVtxSPD->GetStatus()) printf("VtxSPD: %+e %+e %+e | Nc=%d\n",fVtxSPD->GetX(),fVtxSPD->GetY(),fVtxSPD->GetZ(),fVtxSPD->GetNContributors());
 304   else printf("VtxSPD: N/A\n");
 305   //
 306   fVtxTrc = fESDEvent->GetPrimaryVertexTracks();
 307   if (fVtxTrc->GetStatus()) printf("VtxTRC: %+e %+e %+e | Nc=%d\n",fVtxTrc->GetX(),fVtxTrc->GetY(),fVtxTrc->GetZ(),fVtxTrc->GetNContributors());
 308   else printf("VtxTRC: N/A\n");
 309   //
 310   BuildMCInfo();
 311   CheckTracks();
 312   //
 313 }
 314
 315
 316 //________________________________________________________________________
 317 void AliTaskITSUPerf::Terminate(Option_t *)
 318 {
 319 #ifdef _CLUS_LIST_
 320   fClLists.Delete();
 321 #endif
 322   Printf("Terminating...");
 323 }
 324
 325 //_________________________________________________________________________
 326 void AliTaskITSUPerf::BookStandardHistosCentMCLb(Int_t bin, Int_t mcLb)
 327 {
 328   // book standard set of histos for specific centrality bin and MC status type,
 329   // adding them to output list and management array
 330   //
 331   const double kMaxDPt = 0.2;
 332   const double kMaxDCA = 0.1;
 333
 334   TH2* h2=0;
 335   // pt resolution
 336   h2 = new TH2F(Form("B%d_%s_%s",bin,fgkLabelTypes[mcLb],"dPT2PTvsPT"),
 337                 Form("B%d %s %s",bin,fgkLabelTypes[mcLb],"#Deltap_{T}/p_{T} vs p_{T}"),
 338                 fNPtBins,fPtMin,fPtMax,fNResBins,-kMaxDPt,kMaxDPt);
 339   AddHisto(&fHistosCentMCLb,h2, GetHistoID(kHResPTvsPTMC,mcLb,bin) );
 340   //
 341   // transverse DCA resolution
 342   h2 = new TH2F(Form("B%d_%s_%s",bin,fgkLabelTypes[mcLb],"DCARvsPT"),
 343                 Form("B%d %s %s",bin,fgkLabelTypes[mcLb],"DCA R vs p_{T}"),
 344                 fNPtBins,fPtMin,fPtMax,fNResBins,-kMaxDCA,kMaxDCA);
 345   AddHisto(&fHistosCentMCLb,h2, GetHistoID(kHResDCARvsPTMC,mcLb,bin) );
 346   //
 347   // Z DCA resolution
 348   h2 = new TH2F(Form("B%d_%s_%s",bin,fgkLabelTypes[mcLb],"DCAZvsPT"),
 349                 Form("B%d %s %s",bin,fgkLabelTypes[mcLb],"DCA Z vs p_{T}"),
 350                 fNPtBins,fPtMin,fPtMax,fNResBins,-kMaxDCA,kMaxDCA);
 351   AddHisto(&fHistosCentMCLb,h2, GetHistoID(kHResDCAZvsPTMC,mcLb,bin) );
 352   //
 353 }
 354
 355 //_________________________________________________________________________
 356 void AliTaskITSUPerf::BookStandardHistosCent(Int_t bin)
 357 {
 358   // book standard set of histos for specific centrality bin
 359   //
 360   TH2* h2=0;
 361   //
 362   // MC labels combination vs pt
 363   h2 = new TH2F(Form("B%d_%s_RCBL",bin,"MCLabvsPT"),
 364                 Form("B%d %s Reconstructable",bin,"MCLabComb vs p_{T}"),
 365                 fNPtBins,fPtMin,fPtMax,kNLabelTypes,-0.5,kNLabelTypes-0.5);
 366   AddHisto(&fHistosCent,h2, GetHistoID(kHMatchStatusRcbl,-1,bin) );
 367   for (int i=0;i<kNLabelTypes;i++) h2->GetYaxis()->SetBinLabel(i+1,fgkLabelTypes[i]);
 368   //
 369   h2 = new TH2F(Form("B%d_%s_NotRCBL_prim",bin,"MCLabvsPT"),
 370                 Form("B%d %s Not Reconstructable, Prim.",bin,"MCLabComb vs p_{T}"),
 371                 fNPtBins,fPtMin,fPtMax,kNLabelTypes,-0.5,kNLabelTypes-0.5);
 372   AddHisto(&fHistosCent,h2, GetHistoID(kHMatchStatusNRcblPrim,-1,bin) );
 373   for (int i=0;i<kNLabelTypes;i++) h2->GetYaxis()->SetBinLabel(i+1,fgkLabelTypes[i]);
 374   //
 375   h2 = new TH2F(Form("B%d_%s_NotRCBL_sec",bin,"MCLabvsPT"),
 376                 Form("B%d %s Not Reconstructable, Sec.",bin,"MCLabComb vs p_{T}"),
 377                 fNPtBins,fPtMin,fPtMax,kNLabelTypes,-0.5,kNLabelTypes-0.5);
 378   AddHisto(&fHistosCent,h2, GetHistoID(kHMatchStatusNRcblSec,-1,bin) );
 379   for (int i=0;i<kNLabelTypes;i++) h2->GetYaxis()->SetBinLabel(i+1,fgkLabelTypes[i]);
 380   //
 381   h2 = new TH2F(Form("B%d_MCLr_Prim",bin),
 382                 Form("B%d MCLr Prim vs pt",bin),
 383                 fNPtBins,fPtMin,fPtMax, 7,-0.5,6.5);
 384   AddHisto(&fHistosCent,h2, GetHistoID(kHMCLrPresPrim,-1,bin) );
 385   //
 386   h2 = new TH2F(Form("B%d_MCLr_Sec",bin),
 387                 Form("B%d MCLr Sec vs pt",bin),
 388                 fNPtBins,fPtMin,fPtMax, 7,-0.5,6.5);
 389   AddHisto(&fHistosCent,h2, GetHistoID(kHMCLrPresSec,-1,bin) );
 390   //
 391 }
 392
 393 //_________________________________________________________________________
 394 void AliTaskITSUPerf::BookHistos(Int_t bin)
 395 {
 396   // book standard set of histos, adding them to output list and management array
 397   //
 398   // standard histos
 399   for (int ilb=0;ilb<kNLabelTypes;ilb++) {  // create similar set of histos for all kind of MC labels
 400     BookStandardHistosCentMCLb(bin,ilb);
 401   }
 402   BookStandardHistosCent(bin);              // st.histos regardless of MClabels comb.
 403 }
 404
 405 //_________________________________________________________________________
 406 void AliTaskITSUPerf::AddHisto(TObjArray* array,TObject* h, Int_t at)
 407 {
 408   // add single histo to the set
 409   if (at>=0) array->AddAtAndExpand(h,at);
 410   else       array->Add(h);
 411   fOutput->Add(h);
 412 }
 413
 414 //_________________________________________________________________________
 415 Int_t AliTaskITSUPerf::GetHistoID(Int_t htype, Int_t mcStat, Int_t centBin) const
 416 {
 417   // generate ID for the histo. htype must be <100, mcStat<kNLabelTypes
 418   //
 419   //  if (mcStat>=kNLabelTypes) AliFatal(Form("MCStatus ID=%d > max.allowed %d",mcStat,kNLabelTypes));
 420   if (centBin>=fNCentBins)  AliFatal(Form("CentBin  ID=%d > max.allowed %d",centBin,fNCentBins));
 421   //
 422   if (centBin>=0) {
 423     if (mcStat>=0) { // MCtruthness dependent histos
 424       if (htype>=kHNStdHistosCentMC) AliFatal(Form("StdCentMC histo ID=%d > max.allowed %d",htype,kHNStdHistosCentMC));
 425       return htype + kHNStdHistosCentMC*(mcStat + kNLabelTypes*centBin);
 426     }
 427     else {
 428       if (htype>=kHNStdHistosCent) AliFatal(Form("StdCent histo ID=%d > max.allowed %d",htype,kHNStdHistosCent));
 429       return htype + kHNStdHistosCent*centBin;
 430     }
 431   }
 432   // custom histo
 433   return htype;
 434   //
 435 }
 436
 437 //_________________________________________________________________________
 438 void AliTaskITSUPerf::CheckTracks()
 439 {
 440   // build mc truth info for tracks
 441   Bool_t dump = kFALSE;
 442   Bool_t dump1 = kFALSE;
 443   static int evid = 0;
 444   trackInfo.evID = evid++;
 445   trackInfo.zV   = fESDEvent->GetPrimaryVertex()->GetZ();
 446   //
 447   int ntr = fESDEvent->GetNumberOfTracks();
 448   int ntrMC = fStack->GetNtrack();
 449   //
 450   double fldBz = fESDEvent->GetMagneticField();
 451   float dcaRZ[2];
 452   //
 453   for (int itr=0;itr<ntr;itr++) {
 454     AliESDtrack* trc = fESDEvent->GetTrack(itr);
 455     //
 456     // at the moment we consider only TPC/ITS tracks
 457     //if (!trc->IsOn(AliESDtrack::kTPCin)) continue;
 458     if (!fTPCCut->IsSelected(trc)) continue;
 459     //
 460     trc->GetDZ(fVtxMC[0],fVtxMC[1],fVtxMC[2], fldBz, dcaRZ);
 461     //
 462     Int_t labMC    = trc->GetLabel();
 463     Int_t labMCabs = TMath::Abs(labMC);
 464     Int_t labMCTPC = trc->GetTPCLabel();
 465     Int_t labMCITS = trc->GetITSLabel();
 466     Int_t nClTPC   = trc->GetNcls(1);
 467     Int_t nClITS   = trc->GetNcls(0);
 468     //
 469     UInt_t& mcStatus = (UInt_t &)fMCStatus[labMCabs];
 470     double pt  = trc->Pt();
 471     double eta = trc->Eta();
 472     //
 473     TH2* hlr = (TH2*)( (mcStatus&BIT(kMCPrimBit)) ? GetHisto(&fHistosCent,kHMCLrPresPrim) : GetHisto(&fHistosCent,kHMCLrPresSec));
 474     if (hlr) {
 475       for (int il=0;il<7;il++) {
 476         if (mcStatus&(0x1<<(il+kITSHitBits))) hlr->Fill(pt,il);
 477       }
 478       hlr->Fill(pt,-99); // count tracks
 479     }
 480     //
 481
 482     Bool_t reject = mcStatus & BIT(kTrCondFail);
 483     //
 484     int mcLabType = GetMCLabType(labMCTPC,labMCITS,nClTPC,nClITS);
 485     //
 486     //    if (eta>fEtaMax || eta<fEtaMin) continue;
 487     //    if (nClTPC<fMinTPCclusters) continue;
 488     //    printf("#%3d pt:%5.2f eta:%+5.2f | Lb:%+6d LbTPC:%+6d LbITS:%+6d MCTp:%d | Ntpc:%3d Nits:%3d | Rej:%d\n",
 489     //     itr,pt,eta,labMC,labMCTPC,labMCITS,mcLabType, nClTPC,nClITS,reject);
 490     AliMCParticle *part  = 0;
 491     if (labMCabs>=ntrMC) continue;
 492     part = (AliMCParticle*)fMCEvent->GetTrack(labMCabs);
 493     double ptMC = part->Pt();
 494     double etaMC = part->Eta();
 495     //
 496     if (fTree) {
 497       trackInfo.lbl = labMC;
 498       trackInfo.lblTPC = labMCTPC;
 499       trackInfo.prim = (mcStatus&BIT(kMCPrimBit)) != 0;
 500       trackInfo.rcbl = !reject;
 501       trackInfo.nClITS = trc->GetNcls(0);
 502       trackInfo.nClTPC = trc->GetNcls(1);
 503       trackInfo.ptMC = ptMC;
 504       trackInfo.etaMC = etaMC;
 505       trackInfo.charge = trc->Charge();
 506       trackInfo.pt = trc->Pt();
 507       trackInfo.eta = trc->Eta();
 508       //
 509
 510       //
 511       trackInfo.spl = trc->GetITSModuleIndex(10);
 512       trackInfo.phi = part->Phi();
 513       trackInfo.pdg = part->PdgCode();
 514       trackInfo.dcaR = dcaRZ[0];
 515       trackInfo.dcaZ = dcaRZ[1];
 516       trackInfo.dcaRE = TMath::Sqrt(trc->GetSigmaY2());
 517       trackInfo.dcaZE = TMath::Sqrt(trc->GetSigmaZ2());
 518       trackInfo.ptE   = TMath::Sqrt(trc->GetSigma1Pt2());
 519
 520       const AliExternalTrackParam* ttpc = trc->GetTPCInnerParam();
 521       if (ttpc) {
 522         trackInfo.ptTPC = ttpc->Pt();
 523         trackInfo.etaTPC = ttpc->Eta();
 524         trackInfo.ptTPCE = TMath::Sqrt(ttpc->GetSigma1Pt2());//*ttpc->Pt()*ttpc->Pt();
 525       }
 526       else {
 527         trackInfo.ptTPC = -1;
 528         trackInfo.etaTPC = 0;
 529         trackInfo.ptTPCE = -1;
 530       }
 531       //
 532       trackInfo.nClITSMC = 0;
 533       for (int il=0;il<7;il++) {
 534         trackInfo.mcCl[il] = (mcStatus & (0x1<<(il+kITSHitBits))) != 0;
 535         if (trackInfo.mcCl[il]) trackInfo.nClITSMC++;
 536         trackInfo.rcCl[il] = trc->HasPointOnITSLayer(il);
 537         trackInfo.fcCl[il] = trc->HasSharedPointOnITSLayer(il);
 538       }
 539       fTree->Fill();
 540     }
 541     if (reject) {
 542       GetHisto(&fHistosCent,(mcStatus & BIT(kMCPrimBit)) ? kHMatchStatusNRcblPrim:kHMatchStatusNRcblSec)->Fill(pt,mcLabType);  // matching status
 543       if (dump) {
 544         printf("--#%4d Pt:%5.2f Eta:%5.2f |",itr,pt,eta);
 545         for (int k=0;k<32;k++) printf("%d", (mcStatus&(0x1<<k)) ? 1:0);
 546         printf("| %+5d %+5d %+5d |%3d %d -> %d\n",labMC,labMCTPC,labMCITS,nClTPC,nClITS,mcLabType);
 547       }
 548       continue;
 549     }
 550     GetHisto(&fHistosCent,kHMatchStatusRcbl)->Fill(pt,mcLabType);  // matching status
 551     //
 552     TObjArray* clList = 0;
 553     //
 554     if (part) {
 555       /*
 556       int pdg = part->PdgCode();
 557       if (TMath::Abs(pdg)!=211) continue;
 558       if (nClITS<7) continue;
 559       */
 560 #ifdef _CLUS_LIST_
 561       //
 562       clList = (TObjArray*)fClLists.At(labMCabs);
 563       if (clList) {
 564         printf("Clusters for track %d MCLabel:%d\n",itr,labMC);
 565         for (int ic=0;ic<clList->GetEntriesFast();ic++) {
 566           AliITSUClusterPix* cl = (AliITSUClusterPix*) clList->At(ic);
 567           cl->Print();
 568         }
 569       }
 570 #endif
 571     }
 572
 573     if (part) {
 574       //
 575       //
 576       if (dump) {
 577         printf("  #%4d Pt:%5.2f Eta:%5.2f |",itr,ptMC,etaMC);
 578         for (int k=0;k<32;k++) printf("%d", (mcStatus&(0x1<<k)) ? 1:0);
 579         printf("| %+5d %+5d %+5d |%3d %d -> %d\n",labMC,labMCTPC,labMCITS,nClTPC,nClITS,mcLabType);
 580       }
 581       //
 582       if ( (mcStatus&BIT(kMCPrimBit)) ) {
 583         if (dump1 && TMath::Abs(ptMC-0.5)<0.1) {
 584           printf("#%4d Pt:%5.2f Eta:%5.2f | ",itr,ptMC,etaMC);
 585           for (int k=0;k<7;k++) printf("(%d/%d)", (mcStatus&(0x1<<k)) ? 1:0, trc->HasPointOnITSLayer(k));
 586           printf("| %+5d %+5d %+5d |%3d %d -> %d\n",labMC,labMCTPC,labMCITS,nClTPC,nClITS,mcLabType);
 587         }
 588         // compare MC vs reco track params
 589         if (ptMC>0) GetHisto(&fHistosCentMCLb,kHResPTvsPTMC, mcLabType)->Fill(ptMC, (ptMC-pt)/ptMC);
 590         //
 591         // DCA resolution
 592         GetHisto(&fHistosCentMCLb,kHResDCARvsPTMC, mcLabType)->Fill(ptMC, dcaRZ[0]);
 593         GetHisto(&fHistosCentMCLb,kHResDCAZvsPTMC, mcLabType)->Fill(ptMC, dcaRZ[1]);
 594         //
 595       }
 596       //
 597     }
 598     //
 599
 600
 601   }
 602   //
 603
 604 }
 605
 606 //_________________________________________________________________________
 607 void AliTaskITSUPerf::BuildMCInfo()
 608 {
 609   // build mc truth info for tracks
 610   //
 611   int ntrMC = fStack->GetNtrack();
 612   //
 613   if (fMCStatus.GetSize()<ntrMC) fMCStatus.Set(ntrMC);
 614   fMCStatus.Reset();
 615 #ifdef _CLUS_LIST_
 616   fClLists.Delete();
 617   if (fClLists.GetSize()<ntrMC) fClLists.Expand(ntrMC);
 618   //
 619   // create empty lists for reconstructed tracks
 620   int ntr = fESDEvent->GetNumberOfTracks();
 621   for (int itr=0;itr<ntr;itr++) {
 622     AliESDtrack* trc = fESDEvent->GetTrack(itr);
 623     int lb = TMath::Abs(trc->GetLabel());
 624     if (lb<ntrMC && !fClLists.At(itr)) {
 625       fClLists.AddAt(new TObjArray(7),lb);
 626     }
 627   }
 628 #endif
 629   //
 630   for (int itr=0;itr<ntrMC;itr++) {
 631     //
 632     AliMCParticle *part  = (AliMCParticle*)fMCEvent->GetTrack(itr);
 633     if (!part->Charge()) continue;
 634     //
 635     UInt_t& mcStatus = (UInt_t &)fMCStatus[itr];
 636     mcStatus = 0;
 637     //
 638     if (fStack->IsPhysicalPrimary(itr)) mcStatus |= BIT(kMCPrimBit);
 639     //
 640     /*
 641     TParticle* tp = part->Particle();
 642     printf("mc%4d IsPri=%d",itr,fStack->IsPhysicalPrimary(itr)); tp->Print();
 643     */
 644     //
 645   }
 646   //
 647   int nTrCond = fTrackingCond ? fTrackingCond->GetEntriesFast() : 0;
 648   // fill MC clusters info (if available)
 649   if (fRPTree) {
 650     // load clusters
 651     fITS->LoadClusters(fRPTree);
 652     fITS->ProcessClusters(); // order in the same way as in reconstruction
 653     //
 654     for (int ilr=0;ilr<fITS->GetNLayersActive();ilr++) {
 655       AliITSURecoLayer* lr = fITS->GetLayerActive(ilr);
 656       for (int icl=lr->GetNClusters();icl--;) {
 657         AliCluster* cl = lr->GetCluster(icl);
 658         // check labels
 659         for (int ilb=0;ilb<3;ilb++) {
 660           Int_t lab = cl->GetLabel(ilb);
 661           if (lab<0||lab>=ntrMC) continue;
 662           UInt_t& mcStatus = (UInt_t &)fMCStatus[lab];
 663           mcStatus |= 0x1<<(ilr+kITSHitBits);
 664 #ifdef _CLUS_LIST_
 665           TObjArray *clList = (TObjArray*)fClLists.At(lab);
 666           if (clList) clList->AddLast(cl);  // create cl.list only for reconstructed tracks
 667 #endif
 668         }
 669       }
 670     }
 671     //
 672     if (nTrCond) {
 673       for (int itr=0;itr<ntrMC;itr++) {
 674         //
 675         UInt_t& mcStatus = (UInt_t &)fMCStatus[itr];
 676         Bool_t fail = kTRUE;
 677         UShort_t patt = (mcStatus&0xffff);
 678         for (int ic=0;ic<nTrCond;ic++) {
 679           AliITSUTrackCond* cnd = (AliITSUTrackCond*) fTrackingCond->At(ic);
 680           if (cnd->CheckPattern(patt)) {fail=kFALSE; break;}
 681         }
 682         if (fail) mcStatus |= BIT(kTrCondFail);
 683         //
 684         /*
 685         int nlr = AliITSUAux::NumberOfBitsSet(patt);
 686         AliMCParticle *part  = (AliMCParticle*)fMCEvent->GetTrack(itr);
 687         double etaMC = part->Eta();
 688         if (etaMC>=fEtaMin && etaMC<=fEtaMax && nlr) {
 689           TH2* hlr = (TH2*)( (mcStatus&BIT(kMCPrimBit)) ? GetHisto(&fHistosCent,kHMCLrPresPrim) : GetHisto(&fHistosCent,kHMCLrPresSec));
 690           if (hlr) {
 691             double ptMC = part->Pt();
 692             for (int il=0;il<7;il++) {
 693               if (mcStatus&(0x1<<(il+kITSHitBits))) hlr->Fill(ptMC,il);
 694             }
 695             hlr->Fill(ptMC,-99); // count tracks
 696           }
 697         }
 698         */
 699       }
 700     }
 701     //
 702   }
 703   //
 704 }
 705
 706
 707 //_________________________________________________________________________
 708 Int_t AliTaskITSUPerf::GetMCLabType(Int_t labMCTPC,Int_t labMCITS, Int_t nClTPC, Int_t nClITS)
 709 {
 710   // build type of track (0: both ITS and TPC correct, 1: ITS corr. TPC fake, 2: ITS fake TPC corr., 3: ITS fake TPC fake, 4: mismatch
 711   //
 712   int tp = 0;
 713   if (nClTPC>0 && nClITS<1) return kITSTPCNoMatch;
 714   if (TMath::Abs(labMCTPC)!=TMath::Abs(labMCITS)) tp = kITSTPCMismatch;
 715   else {
 716     if   (labMCITS<0) tp = labMCTPC<0 ? kITS0TPC0 : kITS0TPC1;
 717     else              tp = labMCTPC<0 ? kITS1TPC0 : kITS1TPC1;
 718   }
 719   return tp;
 720   //
 721 }
 722
 723 //_________________________________________________________________________
 724 Int_t AliTaskITSUPerf::GetCentralityBin() const
 725 {
 726   // calculate centrality bin
 727   return 0;
 728 }
 729