PWGPP/ITS/AliTrackMatchingTPCITSCosmics.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2009, 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 //
  18 // AliAnalysisTask to study the TPC-ITS track matching
  19 //
  20 // Author: A.Dainese, andrea.dainese@pd.infn.it
  21 /////////////////////////////////////////////////////////////
  22
  23 #include <TTree.h>
  24 #include <TChain.h>
  25 #include <TNtuple.h>
  26 #include <TList.h>
  27 #include <TH1F.h>
  28 #include <TH2F.h>
  29 #include <TVector3.h>
  30 #include <TGeoManager.h>
  31
  32 #include "AliLog.h"
  33 #include "AliGeomManager.h"
  34 #include "AliTrackPointArray.h"
  35 #include "AliESDInputHandler.h"
  36 #include "AliESDVertex.h"
  37 #include "AliESDtrack.h"
  38 #include "AliESDEvent.h"
  39 #include "AliESDfriend.h"
  40 #include "AliAnalysisTask.h"
  41 #include "AliAnalysisManager.h"
  42 #include "AliExternalTrackParam.h"
  43 #include "AliTrackMatchingTPCITSCosmics.h"
  44
  45
  46 ClassImp(AliTrackMatchingTPCITSCosmics)
  47
  48
  49 //________________________________________________________________________
  50 AliTrackMatchingTPCITSCosmics::AliTrackMatchingTPCITSCosmics(const char *name):
  51 AliAnalysisTask(name,"task"),
  52 fOnlySPDFO(kFALSE),
  53 fReadHLTESD(kFALSE),
  54 fGeometryFileName("geometry.root"),
  55 fESD(0),
  56 fList(0),
  57 fHistEvCount(0),
  58 fntTrks(0),
  59 fntPairs(0),
  60 fntITSPairs(0)
  61 {
  62   // Constructor
  63
  64   // Define input and output slots here
  65   // Input slot #0 works with a TChain
  66   DefineInput(0, TChain::Class());
  67   // Output slot #0 writes into a TList
  68   DefineOutput(0,TList::Class());  //My private output
  69 }
  70
  71 //________________________________________________________________________
  72 AliTrackMatchingTPCITSCosmics::~AliTrackMatchingTPCITSCosmics()
  73 {
  74   // Destructor
  75   if (fList) {
  76     delete fList;
  77     fList = 0;
  78   }
  79   if (fHistEvCount) {
  80     delete fHistEvCount;
  81     fHistEvCount = 0;
  82   }
  83   if (fntTrks) {
  84     delete fntTrks;
  85     fntTrks = 0;
  86   }
  87   if (fntPairs) {
  88     delete fntPairs;
  89     fntPairs = 0;
  90   }
  91   if (fntITSPairs) {
  92     delete fntITSPairs;
  93     fntITSPairs = 0;
  94   }
  95 }
  96
  97 //________________________________________________________________________
  98 void AliTrackMatchingTPCITSCosmics::ConnectInputData(Option_t *)
  99 {
 100   // Connect ESD
 101   // Called once
 102
 103   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
 104   if(!tree) {
 105     printf("ERROR: Could not read chain from input slot 0\n");
 106   } else {
 107     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 108
 109     if(!esdH) {
 110       printf("ERROR: Could not get ESDInputHandler\n");
 111     } else {
 112       fESD = esdH->GetEvent();
 113     }
 114   }
 115
 116   return;
 117 }
 118
 119 //________________________________________________________________________
 120 void AliTrackMatchingTPCITSCosmics::Init()
 121 {
 122   // Initialization
 123
 124   return;
 125 }
 126
 127 //________________________________________________________________________
 128 void AliTrackMatchingTPCITSCosmics::CreateOutputObjects()
 129 {
 130   // Create the output container
 131   //
 132
 133   // load the geometry
 134   if(!gGeoManager) {
 135     AliGeomManager::LoadGeometry(fGeometryFileName.Data());
 136     if(!gGeoManager) {
 137       printf("AliTrackMatchingTPCITSCosmics::CreateOutputObjects(): no geometry loaded \n");
 138       return;
 139     }
 140   }
 141
 142   // Several histograms are more conveniently managed in a TList
 143   fList = new TList();
 144   fList->SetOwner();
 145
 146   fHistEvCount = new TH1F("fHistEvCount","0: all, 1: SPDFO, 2: vtx, 3: TPCtrks, 4: ITStrks, 5: HLT, 6: HLT&TPCtrks, 7: HLT&ITStrks",11,-0.5,10.5);
 147   fList->Add(fHistEvCount);
 148
 149   fntTrks = new TNtuple("fntTrks","TPC-ITS matching","ptTPC:nClsTPC:nClsITSSA:nClsITSMI:SSD1:SSD2:phi:z:dx:dy:dz:drphi:dphi:dtgl");
 150   fList->Add(fntTrks);
 151   fntPairs = new TNtuple("fntPairs","pairs at vertex","nClsITSSAin:nClsITSSAout:nClsITSMIin:nClsITSMIout:dxyITSSA:dzITSSA:dxyITSMI:dzITSMI:ptITSSAin:ptITSSAout:ptITSMIin:ptITSMIout:sigmad0MIin:sigmad0MIout");
 152   fList->Add(fntPairs);
 153   fntITSPairs = new TNtuple("fntITSPairs","pairs at vertex","nClsITSSAin:nClsITSSAout:dxyITSSA:dzITSSA:ptITSSAin:ptITSSAout:d0mu");
 154   fList->Add(fntITSPairs);
 155
 156   return;
 157 }
 158
 159 //________________________________________________________________________
 160 void AliTrackMatchingTPCITSCosmics::Exec(Option_t */*option*/)
 161 {
 162   // Execute analysis for current event:
 163   // write ITS AliTrackPoints for selected tracks to fspTree
 164
 165   // check the geometry
 166   if(!gGeoManager) {
 167     printf("AliTrackMatchingTPCITSCosmics::Exec(): no geometry loaded \n");
 168     return;
 169   }
 170
 171   if(!fESD) {
 172     printf("AliTrackMatchingTPCITSCosmics::Exec(): no ESD \n");
 173     return;
 174   }
 175
 176
 177   fHistEvCount->Fill(0);
 178   PostData(0,fList);
 179
 180   // check if event is triggered by HLT
 181   Bool_t hltTrigg=kFALSE;
 182   if(fReadHLTESD) {
 183     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 184     if(!esdH) {
 185       printf("AliTrackMatchingTPCITSCosmics::Exec(): no ESD \n");
 186       return;
 187     }
 188     AliESDEvent *hltESD = esdH->GetHLTEvent();
 189     if(!hltESD) {
 190       printf("AliTrackMatchingTPCITSCosmics::Exec(): no HLT ESD \n");
 191       return;
 192     }
 193     if(hltESD->IsHLTTriggerFired()) {fHistEvCount->Fill(5);hltTrigg=kTRUE;}
 194   }
 195
 196
 197   TString triggeredClass = fESD->GetFiredTriggerClasses();
 198   if(triggeredClass.Contains("C0SCO-ABCE-NOPF-CENT")) {
 199     fHistEvCount->Fill(1);
 200   } else {
 201     if(fOnlySPDFO) {  PostData(0,fList); return; }
 202   }
 203
 204   Int_t ntracks = fESD->GetNumberOfTracks();
 205   printf("CONTR: %d\n",fESD->GetVertex()->GetNContributors());
 206   if(fESD->GetVertex()->GetNContributors()>=0) {
 207     fHistEvCount->Fill(2);
 208   }
 209
 210   Int_t nTrksITSSA=0,nTrksTPC=0,nTrksTPCITS=0;
 211   Int_t idxITSSA[10000],idxTPC[10000];
 212   Int_t idxITSSAin=-1,idxITSSAout=-1;
 213   Int_t idxTPCin=-1,idxTPCout=-1;
 214   for(Int_t itr=0; itr<ntracks; itr++) {
 215     AliESDtrack *t = fESD->GetTrack(itr);
 216     Int_t nClsTPC = t->GetNcls(1);
 217     Int_t nClsITS = t->GetNcls(0);
 218     if(nClsITS>=2 && nClsTPC==0) { // ITS SA
 219       idxITSSA[nTrksITSSA]=itr;
 220       nTrksITSSA++;
 221       //printf("%d\n",itr);
 222     }
 223     if(nClsTPC>=50) {              // TPC
 224       idxTPC[nTrksTPC]=itr;
 225       nTrksTPC++;
 226     }
 227     if(nClsTPC>=50 && nClsITS>=2) {// TPC+ITS
 228       nTrksTPCITS++;
 229       /*printf("  --> TPC+ITS: %d + %d (%d,%d,%d,%d,%d,%d) pt = %f\n",nClsTPC,nClsITS,
 230         (Int_t)(t->HasPointOnITSLayer(0)),
 231         (Int_t)(t->HasPointOnITSLayer(1)),
 232         (Int_t)(t->HasPointOnITSLayer(2)),
 233         (Int_t)(t->HasPointOnITSLayer(3)),
 234         (Int_t)(t->HasPointOnITSLayer(4)),
 235         (Int_t)(t->HasPointOnITSLayer(5)),
 236         t->Pt());*/
 237     }
 238   }
 239
 240   //printf("nTrksTPC %d  nTrksITSSA %d\n",nTrksTPC,nTrksITSSA);
 241   if(nTrksITSSA>0) fHistEvCount->Fill(4);
 242   if(hltTrigg && nTrksITSSA>0) fHistEvCount->Fill(7);
 243   if(nTrksTPC>0)   fHistEvCount->Fill(3);
 244   if(hltTrigg && nTrksTPC>0)   fHistEvCount->Fill(6);
 245   if(nTrksTPCITS>0)   fHistEvCount->Fill(8);
 246   if(hltTrigg && nTrksTPCITS>0)   fHistEvCount->Fill(9);
 247
 248   if(nTrksTPC>2 || nTrksTPC==0 || nTrksITSSA>2 || nTrksITSSA==0) { PostData(0,fList); return; }
 249
 250   for(Int_t itr=0; itr<nTrksTPC; itr++) {
 251     AliESDtrack *t = fESD->GetTrack(idxTPC[itr]);
 252     if(t->Py()>0) idxTPCin=idxTPC[itr];
 253     if(t->Py()<0) idxTPCout=idxTPC[itr];
 254   }
 255
 256
 257   for(Int_t itr=0; itr<nTrksITSSA; itr++) {
 258     AliESDtrack *t = fESD->GetTrack(idxITSSA[itr]);
 259     if(t->Py()>0) idxITSSAin=idxITSSA[itr];
 260     if(t->Py()<0) idxITSSAout=idxITSSA[itr];
 261   }
 262
 263   Double_t xyzITS[3],xyzTPC[3];
 264
 265   // analysis for inward track
 266   if(idxITSSAin>-1 && idxTPCin>-1) {
 267     //printf(" %d\n",idxITSSAin);
 268     AliESDtrack *tITS = fESD->GetTrack(idxITSSAin);
 269     const AliExternalTrackParam *outerITS = tITS->GetOuterParam();
 270     xyzITS[0]=-1000.;xyzITS[1]=-1000.;xyzITS[2]=-1000.;
 271     if(outerITS) outerITS->GetXYZAt(50.,fESD->GetMagneticField(),xyzITS);
 272     AliESDtrack *tTPC = fESD->GetTrack(idxTPCin);
 273     AliExternalTrackParam *innerTPC = new AliExternalTrackParam(*(tTPC->GetInnerParam()));
 274     innerTPC->Rotate(tITS->GetAlpha());
 275     innerTPC->GetXYZAt(50.,fESD->GetMagneticField(),xyzTPC);
 276     fntTrks->Fill(innerTPC->Pt(),
 277                  tTPC->GetNcls(1),
 278                  tITS->GetNcls(0),
 279                  tTPC->GetNcls(0),
 280                  tITS->HasPointOnITSLayer(4),
 281                  tITS->HasPointOnITSLayer(5),
 282                  TMath::ATan2(xyzITS[1],xyzITS[0]),
 283                  xyzITS[2],
 284                  xyzITS[0]-xyzTPC[0],
 285                  xyzITS[1]-xyzTPC[1],
 286                  xyzITS[2]-xyzTPC[2],
 287                  0.5*(TMath::Sqrt(xyzITS[0]*xyzITS[0]+xyzITS[1]*xyzITS[1])+TMath::Sqrt(xyzTPC[0]*xyzTPC[0]+xyzTPC[1]*xyzTPC[1]))*(TMath::ATan2(xyzITS[1],xyzITS[0])-TMath::ATan2(xyzTPC[1],xyzTPC[0])),
 288                  TMath::ATan2(tITS->Py(),tITS->Px())-TMath::ATan2(innerTPC->Py(),innerTPC->Px()),
 289                  tITS->Pz()/tITS->Pt()-innerTPC->Pz()/innerTPC->Pt());
 290     delete innerTPC; innerTPC=NULL;
 291   }
 292
 293   // analysis for outward track
 294   if(idxITSSAout>-1 && idxTPCout>-1) {
 295     //printf(" %d\n",idxITSSAout);
 296     AliESDtrack *tITS = fESD->GetTrack(idxITSSAout);
 297     const AliExternalTrackParam *outerITS = tITS->GetOuterParam();
 298     xyzITS[0]=-1000.;xyzITS[1]=-1000.;xyzITS[2]=-1000.;
 299     if(outerITS) outerITS->GetXYZAt(50.,fESD->GetMagneticField(),xyzITS);
 300     AliESDtrack *tTPC = fESD->GetTrack(idxTPCout);
 301     AliExternalTrackParam *innerTPC = new AliExternalTrackParam(*(tTPC->GetInnerParam()));
 302     innerTPC->Rotate(tITS->GetAlpha());
 303     innerTPC->GetXYZAt(50.,fESD->GetMagneticField(),xyzTPC);
 304     fntTrks->Fill(innerTPC->Pt(),
 305                  tTPC->GetNcls(1),
 306                  tITS->GetNcls(0),
 307                  tTPC->GetNcls(0),
 308                  tITS->HasPointOnITSLayer(4),
 309                  tITS->HasPointOnITSLayer(5),
 310                  TMath::ATan2(xyzITS[1],xyzITS[0]),
 311                  xyzITS[2],
 312                  xyzITS[0]-xyzTPC[0],
 313                  xyzITS[1]-xyzTPC[1],
 314                  xyzITS[2]-xyzTPC[2],
 315                  0.5*(TMath::Sqrt(xyzITS[0]*xyzITS[0]+xyzITS[1]*xyzITS[1])+TMath::Sqrt(xyzTPC[0]*xyzTPC[0]+xyzTPC[1]*xyzTPC[1]))*(TMath::ATan2(xyzITS[1],xyzITS[0])-TMath::ATan2(xyzTPC[1],xyzTPC[0])),
 316                  TMath::ATan2(tITS->Py(),tITS->Px())-TMath::ATan2(innerTPC->Py(),innerTPC->Px()),
 317                  tITS->Pz()/tITS->Pt()-innerTPC->Pz()/innerTPC->Pt()
 318                  );
 319     delete innerTPC; innerTPC=NULL;
 320   }
 321
 322   // inward-outward track ITS only
 323   if(idxITSSAin>-1 && idxITSSAout>-1) {
 324     AliESDtrack *tITSin = fESD->GetTrack(idxITSSAin);
 325     AliESDtrack *tITSout = fESD->GetTrack(idxITSSAout);
 326     if(tITSin->HasPointOnITSLayer(0) &&
 327        tITSin->HasPointOnITSLayer(1) &&
 328        tITSout->HasPointOnITSLayer(0) &&
 329        tITSout->HasPointOnITSLayer(1)) {
 330       Double_t alpha = TMath::ATan2(tITSin->Py(),tITSin->Px());
 331       tITSin->Propagate(alpha,0.,fESD->GetMagneticField());
 332       tITSout->Propagate(alpha,0.,fESD->GetMagneticField());
 333       Float_t d0[2],z0[2];
 334       d0[0] = tITSin->GetY();
 335       z0[0] = tITSin->GetZ();
 336       d0[1] = tITSout->GetY();
 337       z0[1] = tITSout->GetZ();
 338       Float_t dxyITS = -(d0[0]-d0[1]);
 339       Float_t dzITS = z0[0]-z0[1];
 340       fntITSPairs->Fill(tITSin->GetNcls(0),tITSout->GetNcls(0),dxyITS,dzITS,tITSin->Pt(),tITSout->Pt(),TMath::Abs(d0[0]));
 341     }
 342   }
 343
 344   // inward-outward track TPC+ITS
 345   if(idxITSSAin>-1 && idxTPCin>-1 && idxITSSAout>-1 && idxTPCout>-1) {
 346     AliESDtrack *tTPCin = fESD->GetTrack(idxTPCin);
 347     AliESDtrack *tITSin = fESD->GetTrack(idxITSSAin);
 348     AliESDtrack *tTPCout = fESD->GetTrack(idxTPCout);
 349     AliESDtrack *tITSout = fESD->GetTrack(idxITSSAout);
 350     if(tTPCin->HasPointOnITSLayer(0) &&
 351        tTPCin->HasPointOnITSLayer(1) &&
 352        tITSin->HasPointOnITSLayer(0) &&
 353        tITSin->HasPointOnITSLayer(1) &&
 354        tTPCout->HasPointOnITSLayer(0) &&
 355        tTPCout->HasPointOnITSLayer(1) &&
 356        tITSout->HasPointOnITSLayer(0) &&
 357        tITSout->HasPointOnITSLayer(1)) {
 358       Double_t alpha = TMath::ATan2(tTPCin->Py(),tTPCin->Px());
 359       tTPCin->Propagate(alpha,0.,fESD->GetMagneticField());
 360       tTPCout->Propagate(alpha,0.,fESD->GetMagneticField());
 361       Float_t d0[2],z0[2];
 362       d0[0] = tTPCin->GetY();
 363       z0[0] = tTPCin->GetZ();
 364       d0[1] = tTPCout->GetY();
 365       z0[1] = tTPCout->GetZ();
 366       Float_t dxyTPC = -(d0[0]-d0[1]);
 367       Float_t dzTPC = z0[0]-z0[1];
 368       alpha = TMath::ATan2(tITSin->Py(),tITSin->Px());
 369       tITSin->Propagate(alpha,0.,fESD->GetMagneticField());
 370       tITSout->Propagate(alpha,0.,fESD->GetMagneticField());
 371       d0[0] = tITSin->GetY();
 372       z0[0] = tITSin->GetZ();
 373       d0[1] = tITSout->GetY();
 374       z0[1] = tITSout->GetZ();
 375       Float_t dxyITS = -(d0[0]-d0[1]);
 376       Float_t dzITS = z0[0]-z0[1];
 377       fntPairs->Fill(tITSin->GetNcls(0),tITSout->GetNcls(0),tTPCin->GetNcls(0),tTPCout->GetNcls(0),dxyITS,dzITS,dxyTPC,dzTPC,tITSin->Pt(),tITSout->Pt(),tTPCin->Pt(),tTPCout->Pt(),tTPCin->GetSigmaY2(),tTPCout->GetSigmaY2());
 378     }
 379   }
 380
 381   PostData(0,fList);
 382
 383   return;
 384 }
 385
 386 //________________________________________________________________________
 387 void AliTrackMatchingTPCITSCosmics::Terminate(Option_t */*option*/)
 388 {
 389   // Terminate analysis
 390   //
 391   AliDebug(2,"AliTrackMatchingTPCITSCosmics: Terminate() \n");
 392
 393   fList = dynamic_cast<TList*> (GetOutputData(0));
 394   if (!fList) {
 395     printf("ERROR: fList not available\n");
 396     return;
 397   }
 398
 399   return;
 400 }
 401