PWGLF/SPECTRA/PiKaPr/HMPID/AliHMPIDAnalysisTask.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, 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 // AliHMPIDAnalysysTask - Class representing a basic analysis tool of HMPID data
  18 // A set of histograms is created.
  19 //==============================================================================
  20 //
  21 // By means of AliHMPIDAnalysisTask.C macro it is possible to use this class
  22 // to perform the analysis on local data, on data on alien using local machine
  23 // and on CAF.
  24
  25 #ifndef AliHMPIDAnalysisTASK_CXX
  26 #define AliHMPIDAnalysisTASK_CXX
  27
  28
  29 #include "TH1.h"
  30 #include "TH2.h"
  31 #include "TFile.h"
  32 #include "TCanvas.h"
  33 #include "TGraphErrors.h"
  34 #include "AliPIDResponse.h"
  35 #include "AliPID.h"
  36 #include "AliVEvent.h"
  37 #include "AliVParticle.h"
  38 #include "AliVTrack.h"
  39 #include "AliESDtrackCuts.h"
  40 #include "AliAnalysisFilter.h"
  41 #include "AliAnalysisManager.h"
  42 #include "AliESDInputHandler.h"
  43 #include "AliMCEventHandler.h"
  44 #include "AliMCEvent.h"
  45 #include "AliESDtrack.h"
  46 #include "AliPID.h"
  47 #include "AliLog.h"
  48 #include "AliHMPIDAnalysisTask.h"
  49
  50 ClassImp(AliHMPIDAnalysisTask)
  51
  52 //__________________________________________________________________________
  53 AliHMPIDAnalysisTask::AliHMPIDAnalysisTask() :
  54   fESD(0x0),fMC(0x0),fUseMC(kTRUE),
  55   fHmpHistList(0x0),
  56   fHmpNevents(0x0),
  57   fZvertex(0x0),
  58   fPIDResponse(0x0),
  59   fTrackCuts(0x0),
  60   fTrackFilter(0x0),
  61   fTree(0x0)
  62 {
  63   //
  64   //Default ctor
  65   //
  66   for (Int_t i=0; i<69; i++) fVar[i]=0;
  67 }
  68
  69 //___________________________________________________________________________
  70 AliHMPIDAnalysisTask::AliHMPIDAnalysisTask(const Char_t* name) :
  71   AliAnalysisTaskSE(name),
  72   fESD(0x0), fMC(0x0), fUseMC(kTRUE),
  73   fHmpHistList(0x0),
  74   fHmpNevents(0x0),
  75   fZvertex(0x0),
  76   fPIDResponse(0x0),
  77   fTrackCuts(0x0),
  78   fTrackFilter(0x0),
  79   fTree(0x0)
  80 {
  81   //
  82   // Constructor. Initialization of Inputs and Outputs
  83   //
  84   for (Int_t i=0; i<69; i++) fVar[i]=0;
  85
  86   DefineOutput(1,TList::Class());
  87   DefineOutput(2,TTree::Class());
  88 }
  89
  90 //___________________________________________________________________________
  91 AliHMPIDAnalysisTask& AliHMPIDAnalysisTask::operator=(const AliHMPIDAnalysisTask& c)
  92 {
  93   //
  94   // Assignment operator
  95   //
  96   if (this!=&c) {
  97     AliAnalysisTaskSE::operator=(c);
  98     fESD             = c.fESD;
  99     fMC              = c.fMC;
 100     fUseMC           = c.fUseMC;
 101     fHmpHistList     = c.fHmpHistList;
 102     fHmpNevents      = c.fHmpNevents;
 103     fZvertex         = c.fZvertex;
 104     fPIDResponse     = c.fPIDResponse;
 105     fTrackCuts       = c.fTrackCuts;
 106     fTrackFilter     = c.fTrackFilter;
 107     fTree            = c.fTree;
 108     for (Int_t i=0; i<69; i++) fVar[i]=c.fVar[i];
 109   }
 110   return *this;
 111 }
 112
 113 //___________________________________________________________________________
 114 AliHMPIDAnalysisTask::AliHMPIDAnalysisTask(const AliHMPIDAnalysisTask& c) :
 115   AliAnalysisTaskSE(c),
 116   fESD(c.fESD),fMC(c.fMC),fUseMC(c.fUseMC),
 117   fHmpHistList(c.fHmpHistList),
 118   fHmpNevents(c.fHmpNevents),
 119   fZvertex(c.fZvertex),
 120   fPIDResponse(c.fPIDResponse),
 121   fTrackCuts(c.fTrackCuts),
 122   fTrackFilter(c.fTrackFilter),
 123   fTree(c.fTree)
 124 {
 125   //
 126   // Copy Constructor
 127   //
 128   for (Int_t i=0; i<69; i++) fVar[i]=c.fVar[i];
 129 }
 130
 131 //___________________________________________________________________________
 132 AliHMPIDAnalysisTask::~AliHMPIDAnalysisTask() {
 133   //
 134   //destructor
 135   //
 136   Info("~AliHMPIDAnalysisTask","Calling Destructor");
 137   if (fHmpHistList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fHmpHistList;
 138 }
 139
 140 //___________________________________________________________________________
 141 void AliHMPIDAnalysisTask::ConnectInputData(Option_t *)
 142 {
 143   // Connect ESD here
 144
 145   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 146   if (!esdH) {
 147     AliDebug(2,Form("ERROR: Could not get ESDInputHandler"));
 148   } else
 149     fESD = esdH->GetEvent();
 150     fPIDResponse = esdH->GetPIDResponse();
 151
 152   if (fUseMC){
 153     // Connect MC
 154     AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 155     if (!mcH) {
 156       AliDebug(2,Form("ERROR: Could not get MCEventHandler"));
 157       fUseMC = kFALSE;
 158     } else
 159       fMC = mcH->MCEvent();
 160       if (!fMC) AliDebug(2,Form("ERROR: Could not get MCEvent"));
 161   }
 162
 163   fTrackCuts = new AliESDtrackCuts("fTrackCuts", "Standard");
 164   fTrackCuts->SetAcceptKinkDaughters(kFALSE);
 165   fTrackCuts->SetMinNClustersTPC(80);
 166   fTrackCuts->SetMaxChi2PerClusterTPC(4);
 167   fTrackCuts->SetMaxDCAToVertexXY(3);
 168   fTrackCuts->SetMaxDCAToVertexZ(3);
 169   fTrackCuts->SetRequireTPCRefit(kTRUE);
 170   fTrackCuts->SetRequireITSRefit(kTRUE);
 171   fTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
 172
 173   fTrackFilter = new AliAnalysisFilter("trackFilter");
 174   fTrackFilter->AddCuts(fTrackCuts);
 175
 176 }
 177 //***************************************************************************************************************************************************************************
 178 void AliHMPIDAnalysisTask::UserExec(Option_t *)
 179 {
 180
 181   fHmpNevents->Fill(0);
 182
 183   const AliESDVertex *vertex = fESD->GetPrimaryVertexTracks();
 184
 185   if(!vertex || !vertex->GetStatus() || vertex->GetNContributors()<1) {
 186
 187     // SPD vertex
 188     vertex = fESD->GetPrimaryVertexSPD();
 189     if(!vertex) {
 190       PostData(1,fHmpHistList);
 191       PostData(2,fTree);
 192       return;
 193     }
 194     if(!vertex->GetStatus()){
 195       PostData(1,fHmpHistList);
 196       PostData(2,fTree);
 197       return;
 198     }
 199     if(vertex->GetNContributors()<1){
 200       PostData(1,fHmpHistList);
 201       PostData(2,fTree);
 202       return;
 203     }
 204   }
 205
 206   Double_t vtxPos[3] = {999, 999, 999};
 207   if(vertex) vertex->GetXYZ(vtxPos);
 208
 209   fHmpNevents->Fill(1);
 210
 211   fZvertex->Fill(vtxPos[2]);
 212
 213   if(!fPIDResponse){
 214       PostData(1,fHmpHistList);
 215       PostData(2,fTree);
 216       return;
 217    }
 218
 219   Double_t probsHMP[5], probsTPC[5], probsTOF[5];
 220   Float_t  nSigmasTOF[5], nSigmasTPC[5];
 221
 222   AliESDtrack *track = 0;
 223
 224   Double_t ktol = 0.001;
 225   Double_t r[3], rin[3], rout[3];
 226
 227   //
 228   // Main loop function, executed on Event basis
 229   //
 230   for (Int_t iTrack = 0; iTrack<fESD->GetNumberOfTracks(); iTrack++) {
 231
 232     track = fESD->GetTrack(iTrack);
 233     if(!track) continue;
 234
 235     //if(!fTrackFilter->IsSelected(track)) continue;
 236     //if(!track->GetStatus() || !AliESDtrack::kTOFout || !AliESDtrack::kTIME || track->GetTOFsignal()<12000 || track->GetTOFsignal()>200000 || track->GetIntegratedLength()<365) continue;
 237
 238     track->GetXYZ(r);
 239     track->GetInnerXYZ(rin);
 240     track->GetOuterXYZ(rout);
 241
 242     if(Equal(track->GetHMPIDsignal(),-20.,ktol)) continue;
 243     if(track->GetHMPIDcluIdx() < 0) continue;
 244
 245     Int_t q, nph;
 246     Float_t x, y;
 247     Float_t xpc, ypc, th, ph;
 248     track->GetHMPIDmip(x,y,q,nph);
 249     track->GetHMPIDtrk(xpc,ypc,th,ph);
 250
 251     Int_t ITSrefit = track->IsOn(AliESDtrack::kITSrefit);
 252     Int_t TPCrefit = track->IsOn(AliESDtrack::kTPCrefit);
 253
 254     if(Equal(x,0.,ktol) && Equal(y,0.,ktol) && Equal(xpc,0.,ktol) && Equal(ypc,0.,ktol)) continue;
 255
 256     Double_t pHmp[3] = {0}, pHmp3 = 0;
 257     if (track->GetOuterHmpPxPyPz(pHmp)) pHmp3 = TMath::Sqrt(pHmp[0]*pHmp[0]+pHmp[1]*pHmp[1]+pHmp[2]*pHmp[2]);
 258
 259     Float_t b[2];
 260     Float_t bCov[3];
 261     track->GetImpactParameters(b,bCov);
 262
 263     track->GetHMPIDpid(probsHMP);
 264     fPIDResponse->ComputeTOFProbability(track,AliPID::kSPECIES,probsTOF);
 265     fPIDResponse->ComputeTPCProbability(track,AliPID::kSPECIES,probsTPC);
 266
 267     for(Int_t ispecie=0; ispecie<AliPID::kSPECIES; ispecie++){
 268       nSigmasTOF[ispecie] = fPIDResponse->NumberOfSigmasTOF((AliVTrack*)track,(AliPID::EParticleType)ispecie);
 269       nSigmasTPC[ispecie] = fPIDResponse->NumberOfSigmasTPC((AliVTrack*)track,(AliPID::EParticleType)ispecie);
 270     }
 271
 272     fVar[0] = track->GetHMPIDcluIdx()/1000000;
 273     fVar[1] = pHmp3;
 274     fVar[2] = (Float_t)track->P();
 275     fVar[3] = xpc;
 276     fVar[4] = ypc;
 277     fVar[5] = x;
 278     fVar[6] = y;
 279     fVar[7] = (Float_t)track->GetHMPIDsignal();
 280     fVar[8] = q;
 281     fVar[9] = th;
 282     fVar[10] = ph;
 283     fVar[11] = (Float_t)track->GetSign();
 284     fVar[12] = (Float_t)nph;
 285     fVar[13] = (Float_t)track->GetNcls(1);
 286     fVar[14] = (Float_t)probsHMP[0];
 287     fVar[15] = (Float_t)probsHMP[1];
 288     fVar[16] = (Float_t)probsHMP[2];
 289     fVar[17] = (Float_t)probsHMP[3];
 290     fVar[18] = (Float_t)probsHMP[4];
 291     fVar[19] = (Float_t)track->GetTOFsignal();
 292     fVar[20] = (Float_t)track->GetKinkIndex(0);
 293     fVar[21] = (Float_t)track->Xv();
 294     fVar[22] = (Float_t)track->Yv();
 295     fVar[23] = (Float_t)track->Zv();
 296     fVar[24] = (Float_t)track->GetTPCchi2();
 297     fVar[25] = b[0];
 298     fVar[26] = b[1];
 299     fVar[27] = track->GetHMPIDcluIdx()%1000000/1000;
 300     fVar[28] = vtxPos[0];
 301     fVar[29] = vtxPos[1];
 302     fVar[30] = vtxPos[2];
 303     fVar[31] = (Float_t)ITSrefit;
 304     fVar[32] = (Float_t)TPCrefit;
 305     fVar[33] = (Float_t)track->Eta();
 306     fVar[34] = (Float_t)r[0];
 307     fVar[35] = (Float_t)r[1];
 308     fVar[36] = (Float_t)r[2];
 309     fVar[37] = (Float_t)rout[0];
 310     fVar[38] = (Float_t)rout[1];
 311     fVar[39] = (Float_t)rout[2];
 312     fVar[40] = (Float_t)pHmp[0];
 313     fVar[41] = (Float_t)pHmp[1];
 314     fVar[42] = (Float_t)pHmp[2];
 315     fVar[43] = (Float_t)track->Px();
 316     fVar[44] = (Float_t)track->Py();
 317     fVar[45] = (Float_t)track->Pz();
 318     fVar[46] = (Float_t)track->GetTPCClusterInfo(2,1);
 319     fVar[47] = (Float_t)track->GetTPCNclsF();
 320     fVar[48] = nSigmasTOF[0];
 321     fVar[49] = nSigmasTOF[1];
 322     fVar[50] = nSigmasTOF[2];
 323     fVar[51] = nSigmasTOF[3];
 324     fVar[52] = nSigmasTOF[4];
 325     fVar[53] = nSigmasTPC[0];
 326     fVar[54] = nSigmasTPC[1];
 327     fVar[55] = nSigmasTPC[2];
 328     fVar[56] = nSigmasTPC[3];
 329     fVar[57] = nSigmasTPC[4];
 330     fVar[58] = (Float_t)probsTOF[0];
 331     fVar[59] = (Float_t)probsTOF[1];
 332     fVar[60] = (Float_t)probsTOF[2];
 333     fVar[61] = (Float_t)probsTOF[3];
 334     fVar[62] = (Float_t)probsTOF[4];
 335     fVar[63] = (Float_t)probsTPC[0];
 336     fVar[64] = (Float_t)probsTPC[1];
 337     fVar[65] = (Float_t)probsTPC[2];
 338     fVar[66] = (Float_t)probsTPC[3];
 339     fVar[67] = (Float_t)probsTPC[4];
 340     fVar[68] = (Float_t)track->GetHMPIDchi2();
 341
 342     fTree->Fill();
 343
 344   }//track loop
 345
 346   /* PostData(0) is taken care of by AliAnalysisTaskSE */
 347   PostData(1,fHmpHistList);
 348   PostData(2,fTree);
 349 }
 350 //___________________________________________________________________________
 351 void AliHMPIDAnalysisTask::Terminate(Option_t*)
 352 {
 353   // The Terminate() function is the last function to be called during
 354   // a query. It always runs on the client, it can be used to present
 355   // the results graphically or save the results to file.
 356
 357   Info("Terminate"," ");
 358
 359   if (!fUseMC) return;
 360
 361   fHmpHistList = dynamic_cast<TList*> (GetOutputData(1));
 362
 363   if (!fHmpHistList) {
 364     AliError("Histogram List is not available");
 365     return;
 366   }
 367
 368
 369   AliAnalysisTaskSE::Terminate();
 370
 371 }
 372 //___________________________________________________________________________
 373 void AliHMPIDAnalysisTask::UserCreateOutputObjects() {
 374   //
 375   //HERE ONE CAN CREATE OUTPUT OBJECTS
 376   //TO BE SET BEFORE THE EXECUTION OF THE TASK
 377   //
 378
 379   //slot #1
 380 //   OpenFile(1);
 381    fHmpHistList = new TList();
 382    fHmpHistList->SetOwner();
 383
 384    fHmpNevents = new TH1F("fHmpNevents","Number of events",2,0,2);
 385    fHmpHistList->Add(fHmpNevents);
 386
 387    fZvertex = new TH1F("fZvertex","Z primary vertex distribution",4000,-20,20);
 388    fHmpHistList->Add(fZvertex);
 389
 390 //   OpenFile(2);
 391    fTree = new TTree("Tree","Tree with data");
 392    fTree->Branch("Chamber",&fVar[0]);
 393    fTree->Branch("pHmp3",&fVar[1]);
 394    fTree->Branch("P",&fVar[2]);
 395    fTree->Branch("Xpc",&fVar[3]);
 396    fTree->Branch("Ypc",&fVar[4]);
 397    fTree->Branch("X",&fVar[5]);
 398    fTree->Branch("Y",&fVar[6]);
 399    fTree->Branch("HMPIDsignal",&fVar[7]);
 400    fTree->Branch("Charge",&fVar[8]);
 401    fTree->Branch("Theta",&fVar[9]);
 402    fTree->Branch("Phi",&fVar[10]);
 403    fTree->Branch("Sign",&fVar[11]);
 404    fTree->Branch("NumPhotons",&fVar[12]);
 405    fTree->Branch("NumTPCclust",&fVar[13]);
 406    fTree->Branch("ProbHMP0",&fVar[14]);
 407    fTree->Branch("ProbHMP1",&fVar[15]);
 408    fTree->Branch("ProbHMP2",&fVar[16]);
 409    fTree->Branch("ProbHMP3",&fVar[17]);
 410    fTree->Branch("ProbHMP4",&fVar[18]);
 411    fTree->Branch("TOFsignal",&fVar[19]);
 412    fTree->Branch("KinkIndex",&fVar[20]);
 413    fTree->Branch("Xv",&fVar[21]);
 414    fTree->Branch("Yv",&fVar[22]);
 415    fTree->Branch("Zv",&fVar[23]);
 416    fTree->Branch("TPCchi2",&fVar[24]);
 417    fTree->Branch("b0",&fVar[25]);
 418    fTree->Branch("b1",&fVar[26]);
 419    fTree->Branch("ClustSize",&fVar[27]);
 420    fTree->Branch("PrimVertexX",&fVar[28]);
 421    fTree->Branch("PrimVertexY",&fVar[29]);
 422    fTree->Branch("PrimVertexZ",&fVar[30]);
 423    fTree->Branch("ITSrefit",&fVar[31]);
 424    fTree->Branch("TPCrefit",&fVar[32]);
 425    fTree->Branch("Eta",&fVar[33]);
 426    fTree->Branch("xTrack",&fVar[34]);
 427    fTree->Branch("yTrack",&fVar[35]);
 428    fTree->Branch("zTrack",&fVar[36]);
 429    fTree->Branch("xOuterTrack",&fVar[37]);
 430    fTree->Branch("yOuterTrack",&fVar[38]);
 431    fTree->Branch("zOuterTrack",&fVar[39]);
 432    fTree->Branch("pHmpX",&fVar[40]);
 433    fTree->Branch("pHmpY",&fVar[41]);
 434    fTree->Branch("pHmpZ",&fVar[42]);
 435    fTree->Branch("Px",&fVar[43]);
 436    fTree->Branch("Py",&fVar[44]);
 437    fTree->Branch("Pz",&fVar[45]);
 438    fTree->Branch("nCrossedRowsTPC",&fVar[46]);
 439    fTree->Branch("findableClustTPC",&fVar[47]);
 440    fTree->Branch("NSigmasTOF0",&fVar[48]);
 441    fTree->Branch("NSigmasTOF1",&fVar[49]);
 442    fTree->Branch("NSigmasTOF2",&fVar[50]);
 443    fTree->Branch("NSigmasTOF3",&fVar[51]);
 444    fTree->Branch("NSigmasTOF4",&fVar[52]);
 445    fTree->Branch("NSigmasTPC0",&fVar[53]);
 446    fTree->Branch("NSigmasTPC1",&fVar[54]);
 447    fTree->Branch("NSigmasTPC2",&fVar[55]);
 448    fTree->Branch("NSigmasTPC3",&fVar[56]);
 449    fTree->Branch("NSigmasTPC4",&fVar[57]);
 450    fTree->Branch("ProbTOF0",&fVar[58]);
 451    fTree->Branch("ProbTOF1",&fVar[59]);
 452    fTree->Branch("ProbTOF2",&fVar[60]);
 453    fTree->Branch("ProbTOF3",&fVar[61]);
 454    fTree->Branch("ProbTOF4",&fVar[62]);
 455    fTree->Branch("ProbTPC0",&fVar[63]);
 456    fTree->Branch("ProbTPC1",&fVar[64]);
 457    fTree->Branch("ProbTPC2",&fVar[65]);
 458    fTree->Branch("ProbTPC3",&fVar[66]);
 459    fTree->Branch("ProbTPC4",&fVar[67]);
 460    fTree->Branch("HmpSigma",&fVar[68]);
 461
 462    PostData(1,fHmpHistList);
 463    PostData(2,fTree);
 464 }
 465
 466 //____________________________________________________________________________________________________________________________________
 467 Bool_t AliHMPIDAnalysisTask::Equal(Double_t x, Double_t y, Double_t tolerance)
 468 {
 469  return abs(x - y) <= tolerance ;
 470 }
 471
 472 #endif