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 <TGeoManager.h>
  64
  65
  66 ClassImp(AliTaskITSUPerf)
  67
  68
  69 const char* AliTaskITSUPerf::fgkLabelTypes[AliTaskITSUPerf::kNLabelTypes] = {
  70   "ITSokTPCok"
  71   ,"ITSokTPCfk"
  72   ,"ITSfkTPCok"
  73   ,"ITSfkTPCfk"
  74   ,"ITSTPCmismatch"
  75   ,"ITSTPCnoMatch"
  76 };
  77
  78 //________________________________________________________________________
  79 /*//Default constructor
  80 AliTaskITSUPerf::AliTaskITSUPerf(const char *name)
  81   : AliAnalysisTaskSE(name),
  82 */
  83 //________________________________________________________________________
  84 AliTaskITSUPerf::AliTaskITSUPerf(const char *name)
  85   : AliAnalysisTaskSE(name)
  86     //
  87   ,fOutput(0)
  88   ,fHistosCentMCLb(0)
  89   ,fHistosCent(0)
  90   ,fRPTree(0)
  91   ,fStack(0)
  92   ,fMCEvent(0)
  93   ,fVtxSPD(0)
  94   ,fVtxTrc(0)
  95   ,fESDEvent(0)
  96   ,fUseSpecialOutput(kTRUE)
  97   ,fUseMC(kTRUE)
  98   ,fTrackingCond(0)
  99     //
 100   ,fGeom(0)
 101   ,fITS(0)
 102   ,fNSelTracksMC(0)
 103   ,fMCStatus(0)
 104   ,fNPtBins(20)
 105   ,fNResBins(100)
 106   ,fPtMin(0)
 107   ,fPtMax(10.)
 108   ,fEtaMin(-3.)
 109   ,fEtaMax( 3.)
 110   ,fZVertexMin(-15.)
 111   ,fZVertexMax(15.)
 112     //
 113   ,fCurrCentBin(0)
 114   ,fNCentBins(1)
 115   ,fUseCentralityVar(0)
 116 {
 117   // Constructor
 118
 119   DefineOutput(1, TList::Class());
 120   //
 121 }
 122
 123 //________________________________________________________________________
 124 AliTaskITSUPerf::~AliTaskITSUPerf()
 125 {
 126   // Destructor
 127   // histograms are in the output list and deleted when the output
 128   // list is deleted by the TSelector dtor
 129   if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {  //RRR
 130     printf("Deleteing output\n");
 131     delete fOutput;
 132     fOutput = 0;
 133   }
 134   //
 135   fHistosCentMCLb.Clear();
 136   fHistosCent.Clear();
 137   //
 138   delete fITS;
 139   delete fGeom;
 140   //
 141 }
 142
 143 //________________________________________________________________________
 144 void AliTaskITSUPerf::UserCreateOutputObjects()
 145 {
 146   //
 147   if (fUseSpecialOutput) OpenFile(1);
 148   fOutput = new TList();
 149   fOutput->SetOwner();
 150   //
 151   // init geometry
 152   if (!gGeoManager) {
 153     TString geom = "geometry.root";
 154     AliGeomManager::LoadGeometry(geom.Data());
 155     if (!gGeoManager) AliFatal("Failed to load geometry");
 156   }
 157   //
 158   // create ITS interface
 159   fGeom = new AliITSUGeomTGeo(kTRUE,kTRUE);
 160   AliITSUClusterPix::SetGeom(fGeom);
 161   fITS = new AliITSURecoDet(fGeom,"ITSURecoInterface");
 162   fITS->CreateClusterArrays();
 163   //
 164   // Create histograms
 165   for (int ib=0;ib<fNCentBins;ib++) {
 166     BookHistos(ib);
 167   }
 168   //
 169   int nhist = fOutput->GetEntries();
 170   for (int i=0;i<nhist;i++) {
 171     TObject* hst = fOutput->At(i);
 172     if (!hst || !(hst->InheritsFrom(TH1::Class()))) continue;
 173     ((TH1*)hst)->Sumw2();
 174   }
 175   //
 176   PostData(1, fOutput);
 177   //
 178 }
 179
 180 //________________________________________________________________________
 181 void AliTaskITSUPerf::UserExec(Option_t *)
 182 {
 183   // Main loop
 184   static int evID = 0;
 185   AliInfo(Form("Event: %d",evID++));
 186   //
 187   AliAnalysisManager* anMan = AliAnalysisManager::GetAnalysisManager();
 188   fRPTree = 0;
 189   AliESDInputHandlerRP *handRP = (AliESDInputHandlerRP*)anMan->GetInputEventHandler();
 190   if (!handRP) { AliFatal("No RP handler"); return; }
 191   //
 192   fESDEvent  = handRP->GetEvent();
 193   if (!fESDEvent) { AliFatal("No AliESDEvent"); return; }
 194   //
 195   fRPTree = handRP->GetTreeR("ITS");
 196   if (!fRPTree) { AliFatal("Invalid ITS cluster tree"); return; }
 197   //
 198   AliMCEventHandler* eventHandler = (AliMCEventHandler*)anMan->GetMCtruthEventHandler();
 199   if (!eventHandler) { AliFatal("Could not retrieve MC event handler"); return; }
 200   fMCEvent = eventHandler->MCEvent();
 201   if (!fMCEvent) { AliFatal("Could not retrieve MC event"); return; }
 202   fStack = fMCEvent->Stack();
 203   if (!fStack) { AliFatal("Stack is not available"); return; }
 204   //
 205   // MC Generator info
 206   AliGenEventHeader* mcGenH = 0;
 207   mcGenH = fMCEvent->GenEventHeader();
 208   TArrayF vtmc(3);
 209   mcGenH->PrimaryVertex(vtmc);
 210   for (int i=3;i--;) fVtxMC[i] = vtmc[i];
 211   //
 212   printf("MagField: %f\n",fESDEvent->GetMagneticField());
 213   printf("VtxMC : %+e %+e %+e\n",fVtxMC[0],fVtxMC[1],fVtxMC[2]);
 214   //
 215   fVtxSPD = fESDEvent->GetPrimaryVertexSPD();
 216   //
 217   if (fVtxSPD->GetStatus()) printf("VtxSPD: %+e %+e %+e | Nc=%d\n",fVtxSPD->GetX(),fVtxSPD->GetY(),fVtxSPD->GetZ(),fVtxSPD->GetNContributors());
 218   else printf("VtxSPD: N/A\n");
 219   //
 220   fVtxTrc = fESDEvent->GetPrimaryVertexTracks();
 221   if (fVtxTrc->GetStatus()) printf("VtxTRC: %+e %+e %+e | Nc=%d\n",fVtxTrc->GetX(),fVtxTrc->GetY(),fVtxTrc->GetZ(),fVtxTrc->GetNContributors());
 222   else printf("VtxTRC: N/A\n");
 223   //
 224   BuildMCInfo();
 225   CheckTracks();
 226   //
 227 }
 228
 229
 230 //________________________________________________________________________
 231 void AliTaskITSUPerf::Terminate(Option_t *)
 232 {
 233   Printf("Terminating...");
 234 }
 235
 236 //_________________________________________________________________________
 237 void AliTaskITSUPerf::BookStandardHistosCentMCLb(Int_t bin, Int_t mcLb)
 238 {
 239   // book standard set of histos for specific centrality bin and MC status type,
 240   // adding them to output list and management array
 241   //
 242   const double kMaxDPt = 0.3;
 243   const double kMaxDCA = 0.3;
 244
 245   TH2* h2=0;
 246   // pt resolution
 247   h2 = new TH2F(Form("B%d_%s_%s",bin,fgkLabelTypes[mcLb],"dPT2PTvsPT"),
 248                 Form("B%d %s %s",bin,fgkLabelTypes[mcLb],"#Deltap_{T}/p_{T} vs p_{T}"),
 249                 fNPtBins,fPtMin,fPtMax,fNResBins,-kMaxDPt,kMaxDPt);
 250   AddHisto(&fHistosCentMCLb,h2, GetHistoID(kHResPTvsPTMC,mcLb,bin) );
 251   //
 252   // transverse DCA resolution
 253   h2 = new TH2F(Form("B%d_%s_%s",bin,fgkLabelTypes[mcLb],"DCARvsPT"),
 254                 Form("B%d %s %s",bin,fgkLabelTypes[mcLb],"DCA R vs p_{T}"),
 255                 fNPtBins,fPtMin,fPtMax,fNResBins,-kMaxDCA,kMaxDCA);
 256   AddHisto(&fHistosCentMCLb,h2, GetHistoID(kHResDCARvsPTMC,mcLb,bin) );
 257   //
 258   // Z DCA resolution
 259   h2 = new TH2F(Form("B%d_%s_%s",bin,fgkLabelTypes[mcLb],"DCAZvsPT"),
 260                 Form("B%d %s %s",bin,fgkLabelTypes[mcLb],"DCA Z vs p_{T}"),
 261                 fNPtBins,fPtMin,fPtMax,fNResBins,-kMaxDCA,kMaxDCA);
 262   AddHisto(&fHistosCentMCLb,h2, GetHistoID(kHResDCAZvsPTMC,mcLb,bin) );
 263   //
 264 }
 265
 266 //_________________________________________________________________________
 267 void AliTaskITSUPerf::BookStandardHistosCent(Int_t bin)
 268 {
 269   // book standard set of histos for specific centrality bin
 270   //
 271   TH2* h2=0;
 272   //
 273   // MC labels combination vs pt
 274   h2 = new TH2F(Form("B%d_%s",bin,"MCLabvsPT"),
 275                 Form("B%d %s",bin,"MCLabComb vs p_{T}"),
 276                 fNPtBins,fPtMin,fPtMax,kNLabelTypes,-0.5,kNLabelTypes-0.5);
 277   AddHisto(&fHistosCent,h2, GetHistoID(kHMatchStatus,-1,bin) );
 278   for (int i=0;i<kNLabelTypes;i++) h2->GetYaxis()->SetBinLabel(i+1,fgkLabelTypes[i]);
 279   //
 280
 281 }
 282
 283 //_________________________________________________________________________
 284 void AliTaskITSUPerf::BookHistos(Int_t bin)
 285 {
 286   // book standard set of histos, adding them to output list and management array
 287   //
 288   // standard histos
 289   for (int ilb=0;ilb<kNLabelTypes;ilb++) {  // create similar set of histos for all kind of MC labels
 290     BookStandardHistosCentMCLb(bin,ilb);
 291   }
 292   BookStandardHistosCent(bin);              // st.histos regardless of MClabels comb.
 293 }
 294
 295 //_________________________________________________________________________
 296 void AliTaskITSUPerf::AddHisto(TObjArray* array,TObject* h, Int_t at)
 297 {
 298   // add single histo to the set
 299   if (at>=0) array->AddAtAndExpand(h,at);
 300   else       array->Add(h);
 301   fOutput->Add(h);
 302 }
 303
 304 //_________________________________________________________________________
 305 Int_t AliTaskITSUPerf::GetHistoID(Int_t htype, Int_t mcStat, Int_t centBin) const
 306 {
 307   // generate ID for the histo. htype must be <100, mcStat<kNLabelTypes
 308   //
 309   if (mcStat>=kNLabelTypes) AliFatal(Form("MCStatus ID=%d > max.allowed %d",mcStat,kNLabelTypes));
 310   if (centBin>=fNCentBins)  AliFatal(Form("CentBin  ID=%d > max.allowed %d",centBin,fNCentBins));
 311   //
 312   if (centBin>=0) {
 313     if (mcStat>=0) { // MCtruthness dependent histos
 314       if (htype>=kHNStdHistosCentMC) AliFatal(Form("StdCentMC histo ID=%d > max.allowed %d",htype,kHNStdHistosCentMC));
 315       return htype + kHNStdHistosCentMC*(mcStat + kNLabelTypes*centBin);
 316     }
 317     else {
 318       if (htype>=kHNStdHistosCent) AliFatal(Form("StdCent histo ID=%d > max.allowed %d",htype,kHNStdHistosCent));
 319       return htype + kHNStdHistosCent*centBin;
 320     }
 321   }
 322   // custom histo
 323   return htype;
 324   //
 325 }
 326
 327 //_________________________________________________________________________
 328 void AliTaskITSUPerf::CheckTracks()
 329 {
 330   // build mc truth info for tracks
 331   //
 332   int ntr = fESDEvent->GetNumberOfTracks();
 333   int ntrMC = fStack->GetNtrack();
 334   //
 335   double fldBz = fESDEvent->GetMagneticField();
 336   float dcaRZ[2];
 337   //
 338   for (int itr=0;itr<ntr;itr++) {
 339     AliESDtrack* trc = fESDEvent->GetTrack(itr);
 340     //
 341     // at the moment we consider only TPC/ITS tracks
 342     if (!trc->IsOn(AliESDtrack::kTPCin)) continue;
 343     Int_t labMC    = trc->GetLabel();
 344     Int_t labMCabs = TMath::Abs(labMC);
 345     Int_t labMCTPC = trc->GetTPCLabel();
 346     Int_t labMCITS = trc->GetITSLabel();
 347     Int_t nClTPC   = trc->GetNcls(1);
 348     Int_t nClITS   = trc->GetNcls(0);
 349     //
 350     UInt_t& mcStatus = (UInt_t &)fMCStatus[labMCabs];
 351     Bool_t reject = mcStatus & BIT(kTrCondFail);
 352     if (reject) continue;
 353     //
 354     int mcLabType = GetMCLabType(labMCTPC,labMCITS,nClTPC,nClITS);
 355     //
 356     double pt  = trc->Pt();
 357     double eta = trc->Eta();
 358     //
 359     //    printf("#%3d pt:%5.2f eta:%+5.2f | Lb:%+6d LbTPC:%+6d LbITS:%+6d MCTp:%d | Ntpc:%3d Nits:%3d\n",itr,pt,eta,labMC,labMCTPC,labMCITS,mcLabType, nClTPC,nClITS);
 360     //
 361     GetHisto(&fHistosCent,kHMatchStatus)->Fill(pt,mcLabType);  // matching status
 362     //
 363     AliMCParticle *part  = 0;
 364     //
 365     if (labMCabs<ntrMC) part = (AliMCParticle*)fMCEvent->GetTrack(labMCabs);
 366     //
 367     if (part) {
 368       //
 369       double ptMC = part->Pt();
 370       double etaMC = part->Eta();
 371       //
 372       if ( (mcStatus&BIT(kMCPrimBit)) ) {
 373         // compare MC vs reco track params
 374         if (ptMC>0) GetHisto(&fHistosCentMCLb,kHResPTvsPTMC, mcLabType)->Fill(ptMC, (ptMC-pt)/ptMC);
 375         //
 376         // DCA resolution
 377         trc->GetDZ(fVtxMC[0],fVtxMC[1],fVtxMC[2], fldBz, dcaRZ);
 378         GetHisto(&fHistosCentMCLb,kHResDCARvsPTMC, mcLabType)->Fill(ptMC, dcaRZ[0]);
 379         GetHisto(&fHistosCentMCLb,kHResDCAZvsPTMC, mcLabType)->Fill(ptMC, dcaRZ[1]);
 380         //
 381       }
 382       //
 383     }
 384     //
 385
 386
 387   }
 388   //
 389
 390 }
 391
 392 //_________________________________________________________________________
 393 void AliTaskITSUPerf::BuildMCInfo()
 394 {
 395   // build mc truth info for tracks
 396   //
 397   int ntrMC = fStack->GetNtrack();
 398   //
 399   if (fMCStatus.GetSize()<ntrMC) fMCStatus.Set(ntrMC);
 400   fMCStatus.Reset();
 401   //
 402   for (int itr=0;itr<ntrMC;itr++) {
 403     //
 404     AliMCParticle *part  = (AliMCParticle*)fMCEvent->GetTrack(itr);
 405     if (!part->Charge()) continue;
 406     //
 407     UInt_t& mcStatus = (UInt_t &)fMCStatus[itr];
 408     mcStatus = 0;
 409     //
 410     if (fStack->IsPhysicalPrimary(itr)) mcStatus |= BIT(kMCPrimBit);
 411     //
 412     /*
 413     TParticle* tp = part->Particle();
 414     printf("mc%4d IsPri=%d",itr,fStack->IsPhysicalPrimary(itr)); tp->Print();
 415     */
 416     //
 417   }
 418   //
 419   int nTrCond = fTrackingCond ? fTrackingCond->GetEntriesFast() : 0;
 420   // fill MC clusters info (if available)
 421   if (fRPTree) {
 422     // load clusters
 423     fITS->LoadClusters(fRPTree);
 424     //
 425     for (int ilr=fITS->GetNLayersActive();ilr--;) {
 426       AliITSURecoLayer* lr = fITS->GetLayerActive(ilr);
 427       for (int icl=lr->GetNClusters();icl--;) {
 428         AliCluster* cl = lr->GetCluster(icl);
 429         // check labels
 430         for (int ilb=0;ilb<3;ilb++) {
 431           Int_t lab = cl->GetLabel(ilb);
 432           if (lab<0||lab>=ntrMC) continue;
 433           UInt_t& mcStatus = (UInt_t &)fMCStatus[lab];
 434           mcStatus |= 0x1<<(ilr+kITSHitBits);
 435         }
 436       }
 437     }
 438     //
 439     if (nTrCond) {
 440       for (int itr=0;itr<ntrMC;itr++) {
 441         //
 442         UInt_t& mcStatus = (UInt_t &)fMCStatus[itr];
 443         Bool_t fail = kTRUE;
 444         UShort_t patt = (mcStatus&0xffff);
 445         for (int ic=0;ic<nTrCond;ic++) {
 446           AliITSUTrackCond* cnd = (AliITSUTrackCond*) fTrackingCond->At(ic);
 447           if (cnd->CheckPattern(patt)) {fail=kFALSE; break;}
 448         }
 449         if (fail) mcStatus |= BIT(kTrCondFail);
 450       }
 451     }
 452     //
 453   }
 454   //
 455 }
 456
 457
 458 //_________________________________________________________________________
 459 Int_t AliTaskITSUPerf::GetMCLabType(Int_t labMCTPC,Int_t labMCITS, Int_t nClTPC, Int_t nClITS)
 460 {
 461   // 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
 462   //
 463   int tp = 0;
 464   if (nClTPC>0 && nClITS<1) return kITSTPCNoMatch;
 465   if (TMath::Abs(labMCTPC)!=TMath::Abs(labMCITS)) tp = kITSTPCMismatch;
 466   else {
 467     if   (labMCITS<0) tp = labMCTPC<0 ? kITS0TPC0 : kITS0TPC1;
 468     else              tp = labMCTPC<0 ? kITS1TPC0 : kITS1TPC1;
 469   }
 470   return tp;
 471   //
 472 }
 473
 474 //_________________________________________________________________________
 475 Int_t AliTaskITSUPerf::GetCentralityBin() const
 476 {
 477   // calculate centrality bin
 478   return 0;
 479 }
 480