| 8d63376d |
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 extract from ESD tracks the information |
| 19 | // on ITS tracking efficiency and resolutions. |
| 20 | // |
| 21 | // Author: A.Dainese, andrea.dainese@pd.infn.it |
| 22 | ///////////////////////////////////////////////////////////// |
| 23 | |
| 24 | #include <TChain.h> |
| 25 | #include <TTree.h> |
| 26 | #include <TNtuple.h> |
| 27 | #include <TBranch.h> |
| 28 | #include <TClonesArray.h> |
| 29 | #include <TObjArray.h> |
| 30 | #include <TH1F.h> |
| 31 | #include <TH2F.h> |
| 32 | #include <TCanvas.h> |
| 33 | #include <TParticle.h> |
| 34 | |
| 35 | #include "AliAnalysisTask.h" |
| 36 | #include "AliAnalysisManager.h" |
| 37 | |
| 38 | #include "AliMultiplicity.h" |
| 39 | #include "AliVertexerTracks.h" |
| 40 | #include "AliESDtrack.h" |
| 41 | #include "AliExternalTrackParam.h" |
| 42 | #include "AliESDVertex.h" |
| 43 | #include "AliESDEvent.h" |
| 44 | #include "AliESDfriend.h" |
| 45 | #include "AliESDInputHandler.h" |
| 46 | #include "AliESDInputHandlerRP.h" |
| 47 | #include "AliTrackPointArray.h" |
| 48 | #include "AliITSRecPoint.h" |
| 49 | |
| 50 | #include "AliMCEventHandler.h" |
| 51 | #include "AliMCEvent.h" |
| 52 | #include "AliStack.h" |
| 53 | #include "AliLog.h" |
| 54 | |
| 55 | #include "AliGenEventHeader.h" |
| 56 | #include "AliAnalysisTaskITSTrackingCheck.h" |
| 57 | |
| 58 | |
| 59 | ClassImp(AliAnalysisTaskITSTrackingCheck) |
| 60 | |
| 61 | //________________________________________________________________________ |
| 62 | AliAnalysisTaskITSTrackingCheck::AliAnalysisTaskITSTrackingCheck(const char *name) : |
| 63 | AliAnalysisTask(name, "ITSTrackingCheckTask"), |
| 64 | fReadMC(kTRUE), |
| 65 | fReadRPLabels(kFALSE), |
| 66 | fESD(0), |
| 67 | fESDfriend(0), |
| 68 | fOutput(0), |
| 69 | fHistNtracks(0), |
| 70 | fHistNclsITSMI(0), |
| 71 | fHistNclsITSSA(0), |
| 72 | fHistClusterMapITSMI(0), |
| 73 | fHistClusterMapITSMIok(0), |
| 74 | fHistClusterMapITSMIbad(0), |
| 75 | fHistClusterMapITSMIskipped(0), |
| 76 | fHistClusterMapITSMIoutinz(0), |
| 77 | fHistClusterMapITSMInorefit(0), |
| 78 | fHistClusterMapITSMInocls(0), |
| 79 | fHistClusterMapITSSA(0), |
| 80 | fHistClusterMapITSSAok(0), |
| 81 | fHistClusterMapITSSAbad(0), |
| 82 | fHistClusterMapITSSAskipped(0), |
| 83 | fHistClusterMapITSSAoutinz(0), |
| 84 | fHistClusterMapITSSAnorefit(0), |
| 85 | fHistClusterMapITSSAnocls(0), |
| 86 | fHistPhiTPC(0), |
| 87 | fHistPtTPC(0), |
| 88 | fHistPtITSMI2(0), |
| 89 | fHistPtITSMI3(0), |
| 90 | fHistPtITSMI4(0), |
| 91 | fHistPtITSMI5(0), |
| 92 | fHistPtITSMI6(0), |
| 93 | fHistPtITSMISPD(0), |
| 94 | fNtupleESDTracks(0), |
| 95 | fNtupleITSAlignExtra(0), |
| 96 | fNtupleITSAlignSPDTracklets(0) |
| 97 | { |
| 98 | // Constructor |
| 99 | |
| 100 | for(Int_t i=0; i<10; i++) fCountsPerPtBin[i]=0; |
| 101 | |
| 102 | // Define input and output slots here |
| 103 | // Input slot #0 works with a TChain |
| 104 | DefineInput(0, TChain::Class()); |
| 105 | // Output slot #0 writes into a TList container |
| 106 | DefineOutput(0, TList::Class()); //My private output |
| 107 | } |
| 108 | //________________________________________________________________________ |
| 109 | AliAnalysisTaskITSTrackingCheck::~AliAnalysisTaskITSTrackingCheck() |
| 110 | { |
| 111 | // Destructor |
| 112 | |
| 113 | // histograms are in the output list and deleted when the output |
| 114 | // list is deleted by the TSelector dtor |
| 115 | |
| 116 | if (fOutput) { |
| 117 | delete fOutput; |
| 118 | fOutput = 0; |
| 119 | } |
| 120 | } |
| 121 | |
| 122 | //________________________________________________________________________ |
| 123 | void AliAnalysisTaskITSTrackingCheck::ConnectInputData(Option_t *) |
| 124 | { |
| 125 | // Connect ESD or AOD here |
| 126 | // Called once |
| 127 | |
| 128 | TTree* tree = dynamic_cast<TTree*> (GetInputData(0)); |
| 129 | if (!tree) { |
| 130 | Printf("ERROR: Could not read chain from input slot 0"); |
| 131 | } else { |
| 132 | // Disable all branches and enable only the needed ones |
| 133 | // The next two lines are different when data produced as AliESDEvent is read |
| 134 | |
| 135 | tree->SetBranchStatus("ESDfriend*", 1); |
| 136 | tree->SetBranchAddress("ESDfriend.",&fESDfriend); |
| 137 | |
| 138 | AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()); |
| 139 | |
| 140 | |
| 141 | if (!esdH) { |
| 142 | Printf("ERROR: Could not get ESDInputHandler"); |
| 143 | } else { |
| 144 | fESD = esdH->GetEvent(); |
| 145 | |
| 146 | } |
| 147 | } |
| 148 | |
| 149 | return; |
| 150 | } |
| 151 | |
| 152 | //________________________________________________________________________ |
| 153 | void AliAnalysisTaskITSTrackingCheck::CreateOutputObjects() |
| 154 | { |
| 155 | // Create histograms |
| 156 | // Called once |
| 157 | |
| 158 | for(Int_t i=0; i<10; i++) fCountsPerPtBin[i]=0; |
| 159 | |
| 160 | // Several histograms are more conveniently managed in a TList |
| 161 | fOutput = new TList; |
| 162 | fOutput->SetOwner(); |
| 163 | |
| 164 | fHistNtracks = new TH1F("fHistNtracks", "N ESD tracks; N tracks; Events",5000, -0.5, 4999.5); |
| 165 | fHistNtracks->Sumw2(); |
| 166 | fHistNtracks->SetMinimum(0); |
| 167 | fOutput->Add(fHistNtracks); |
| 168 | |
| 169 | fHistNclsITSMI = new TH1F("fHistNclsITSMI", "N ITS clusters per track (MI); N clusters; Counts",7, -0.5, 6.5); |
| 170 | fHistNclsITSMI->Sumw2(); |
| 171 | fHistNclsITSMI->SetMinimum(0); |
| 172 | fOutput->Add(fHistNclsITSMI); |
| 173 | |
| 174 | fHistNclsITSSA = new TH1F("fHistNclsITSSA", "N ITS clusters per track (SA); N clusters; Counts",7, -0.5, 6.5); |
| 175 | fHistNclsITSSA->Sumw2(); |
| 176 | fHistNclsITSSA->SetMinimum(0); |
| 177 | fOutput->Add(fHistNclsITSSA); |
| 178 | |
| 179 | fHistClusterMapITSMI = new TH1F("fHistClusterMapITSMI", "N tracks with point on Layer (MI); Layer; N tracks",6, -0.5, 5.5); |
| 180 | fHistClusterMapITSMI->Sumw2(); |
| 181 | fHistClusterMapITSMI->SetMinimum(0); |
| 182 | fOutput->Add(fHistClusterMapITSMI); |
| 183 | |
| 184 | fHistClusterMapITSSA = new TH1F("fHistClusterMapITSSA", "N tracks with point on Layer (SA); Layer; N tracks",6, -0.5, 5.5); |
| 185 | fHistClusterMapITSSA->Sumw2(); |
| 186 | fHistClusterMapITSSA->SetMinimum(0); |
| 187 | fOutput->Add(fHistClusterMapITSSA); |
| 188 | |
| 189 | fHistClusterMapITSMIok = new TH1F("fHistClusterMapITSMIok", "N tracks with ok on Layer (MI); Layer; N tracks",6, -0.5, 5.5); |
| 190 | fHistClusterMapITSMIok->Sumw2(); |
| 191 | fHistClusterMapITSMIok->SetMinimum(0); |
| 192 | fOutput->Add(fHistClusterMapITSMIok); |
| 193 | |
| 194 | fHistClusterMapITSSAok = new TH1F("fHistClusterMapITSSAok", "N tracks with ok on Layer (SA); Layer; N tracks",6, -0.5, 5.5); |
| 195 | fHistClusterMapITSSAok->Sumw2(); |
| 196 | fHistClusterMapITSSAok->SetMinimum(0); |
| 197 | fOutput->Add(fHistClusterMapITSSAok); |
| 198 | |
| 199 | fHistClusterMapITSMIbad = new TH1F("fHistClusterMapITSMIbad", "N tracks with bad on Layer (MI); Layer; N tracks",6, -0.5, 5.5); |
| 200 | fHistClusterMapITSMIbad->Sumw2(); |
| 201 | fHistClusterMapITSMIbad->SetMinimum(0); |
| 202 | fOutput->Add(fHistClusterMapITSMIbad); |
| 203 | |
| 204 | fHistClusterMapITSSAbad = new TH1F("fHistClusterMapITSSAbad", "N tracks with bad on Layer (SA); Layer; N tracks",6, -0.5, 5.5); |
| 205 | fHistClusterMapITSSAbad->Sumw2(); |
| 206 | fHistClusterMapITSSAbad->SetMinimum(0); |
| 207 | fOutput->Add(fHistClusterMapITSSAbad); |
| 208 | |
| 209 | fHistClusterMapITSMIskipped = new TH1F("fHistClusterMapITSMIskipped", "N tracks with skip on Layer (MI); Layer; N tracks",6, -0.5, 5.5); |
| 210 | fHistClusterMapITSMIskipped->Sumw2(); |
| 211 | fHistClusterMapITSMIskipped->SetMinimum(0); |
| 212 | fOutput->Add(fHistClusterMapITSMIskipped); |
| 213 | |
| 214 | fHistClusterMapITSSAskipped = new TH1F("fHistClusterMapITSSAskipped", "N tracks with skip on Layer (SA); Layer; N tracks",6, -0.5, 5.5); |
| 215 | fHistClusterMapITSSAskipped->Sumw2(); |
| 216 | fHistClusterMapITSSAskipped->SetMinimum(0); |
| 217 | fOutput->Add(fHistClusterMapITSSAskipped); |
| 218 | |
| 219 | fHistClusterMapITSMIoutinz = new TH1F("fHistClusterMapITSMIoutinz", "N tracks out in z on Layer (MI); Layer; N tracks",6, -0.5, 5.5); |
| 220 | fHistClusterMapITSMIoutinz->Sumw2(); |
| 221 | fHistClusterMapITSMIoutinz->SetMinimum(0); |
| 222 | fOutput->Add(fHistClusterMapITSMIoutinz); |
| 223 | |
| 224 | fHistClusterMapITSSAoutinz = new TH1F("fHistClusterMapITSSAoutinz", "N tracks with out in z on Layer (SA); Layer; N tracks",6, -0.5, 5.5); |
| 225 | fHistClusterMapITSSAoutinz->Sumw2(); |
| 226 | fHistClusterMapITSSAoutinz->SetMinimum(0); |
| 227 | fOutput->Add(fHistClusterMapITSSAoutinz); |
| 228 | |
| 229 | fHistClusterMapITSMInorefit = new TH1F("fHistClusterMapITSMInorefit", "N tracks with norefit on Layer (MI); Layer; N tracks",6, -0.5, 5.5); |
| 230 | fHistClusterMapITSMInorefit->Sumw2(); |
| 231 | fHistClusterMapITSMInorefit->SetMinimum(0); |
| 232 | fOutput->Add(fHistClusterMapITSMInorefit); |
| 233 | |
| 234 | fHistClusterMapITSSAnorefit = new TH1F("fHistClusterMapITSSAnorefit", "N tracks with norefit on Layer (SA); Layer; N tracks",6, -0.5, 5.5); |
| 235 | fHistClusterMapITSSAnorefit->Sumw2(); |
| 236 | fHistClusterMapITSSAnorefit->SetMinimum(0); |
| 237 | fOutput->Add(fHistClusterMapITSSAnorefit); |
| 238 | |
| 239 | fHistClusterMapITSMInocls = new TH1F("fHistClusterMapITSMInocls", "N tracks with nocls on Layer (MI); Layer; N tracks",6, -0.5, 5.5); |
| 240 | fHistClusterMapITSMInocls->Sumw2(); |
| 241 | fHistClusterMapITSMInocls->SetMinimum(0); |
| 242 | fOutput->Add(fHistClusterMapITSMInocls); |
| 243 | |
| 244 | fHistClusterMapITSSAnocls = new TH1F("fHistClusterMapITSSAnocls", "N tracks with nocls on Layer (SA); Layer; N tracks",6, -0.5, 5.5); |
| 245 | fHistClusterMapITSSAnocls->Sumw2(); |
| 246 | fHistClusterMapITSSAnocls->SetMinimum(0); |
| 247 | fOutput->Add(fHistClusterMapITSSAnocls); |
| 248 | |
| 249 | fHistPhiTPC = new TH1F("fHistPhiTPC","Azimuthal distribution of TPC tracks; #phi; N tracks",100, -3.1415, 3.1415); |
| 250 | fHistPhiTPC->Sumw2(); |
| 251 | fHistPhiTPC->SetMinimum(0); |
| 252 | fOutput->Add(fHistPhiTPC); |
| 253 | |
| 254 | fHistPtTPC = new TH1F("fHistPtTPC","pt distribution of TPC tracks; p_{t} [GeV/c]; N tracks",80, 0, 40); |
| 255 | fHistPtTPC->Sumw2(); |
| 256 | fHistPtTPC->SetMinimum(0); |
| 257 | fOutput->Add(fHistPtTPC); |
| 258 | |
| 259 | fHistPtITSMI6 = new TH1F("fHistPtITSMI6","pt distribution of ITSMI6 tracks; p_{t} [GeV/c]; N tracks",80, 0, 40); |
| 260 | fHistPtITSMI6->Sumw2(); |
| 261 | fHistPtITSMI6->SetMinimum(0); |
| 262 | fOutput->Add(fHistPtITSMI6); |
| 263 | |
| 264 | fHistPtITSMI5 = new TH1F("fHistPtITSMI5","pt distribution of ITSMI5 tracks; p_{t} [GeV/c]; N tracks",80, 0, 40); |
| 265 | fHistPtITSMI5->Sumw2(); |
| 266 | fHistPtITSMI5->SetMinimum(0); |
| 267 | fOutput->Add(fHistPtITSMI5); |
| 268 | |
| 269 | fHistPtITSMI4 = new TH1F("fHistPtITSMI4","pt distribution of ITSMI4 tracks; p_{t} [GeV/c]; N tracks",80, 0, 40); |
| 270 | fHistPtITSMI4->Sumw2(); |
| 271 | fHistPtITSMI4->SetMinimum(0); |
| 272 | fOutput->Add(fHistPtITSMI4); |
| 273 | |
| 274 | fHistPtITSMI3 = new TH1F("fHistPtITSMI3","pt distribution of ITSMI3 tracks; p_{t} [GeV/c]; N tracks",80, 0, 40); |
| 275 | fHistPtITSMI3->Sumw2(); |
| 276 | fHistPtITSMI3->SetMinimum(0); |
| 277 | fOutput->Add(fHistPtITSMI3); |
| 278 | |
| 279 | fHistPtITSMI2 = new TH1F("fHistPtITSMI2","pt distribution of ITSMI2 tracks; p_{t} [GeV/c]; N tracks",80, 0, 40); |
| 280 | fHistPtITSMI2->Sumw2(); |
| 281 | fHistPtITSMI2->SetMinimum(0); |
| 282 | fOutput->Add(fHistPtITSMI2); |
| 283 | |
| 284 | fHistPtITSMISPD = new TH1F("fHistPtITSMISPD","pt distribution of ITSMISPD tracks; p_{t} [GeV/c]; N tracks",80, 0, 40); |
| 285 | fHistPtITSMISPD->Sumw2(); |
| 286 | fHistPtITSMISPD->SetMinimum(0); |
| 287 | fOutput->Add(fHistPtITSMISPD); |
| 288 | |
| 289 | |
| 290 | |
| 291 | // ntuples |
| 292 | // |
| 293 | fNtupleESDTracks = new TNtuple("fNtupleESDTracks","tracks","pt:eta:phi:d0:z0:sigmad0:sigmaz0:ptMC:pdgMC:d0MC:d0MCv:z0MCv:sigmad0MCv:sigmaz0MCv:ITSflag"); |
| 294 | fOutput->Add(fNtupleESDTracks); |
| 295 | |
| 296 | fNtupleITSAlignExtra = new TNtuple("fNtupleITSAlignExtra","ITS alignment checks: extra clusters","layer:x:y:z:dxy:dz:xloc:zloc:npoints"); |
| 297 | fOutput->Add(fNtupleITSAlignExtra); |
| 298 | |
| 299 | fNtupleITSAlignSPDTracklets = new TNtuple("fNtupleITSAlignSPDTracklets","ITS alignment checks: SPD tracklets wrt SPD vertex","phi:theta:z:dxy:dz"); |
| 300 | fOutput->Add(fNtupleITSAlignSPDTracklets); |
| 301 | |
| 302 | return; |
| 303 | } |
| 304 | |
| 305 | //________________________________________________________________________ |
| 306 | void AliAnalysisTaskITSTrackingCheck::Exec(Option_t *) |
| 307 | { |
| 308 | // Main loop |
| 309 | // Called for each event |
| 310 | |
| 311 | if (!fESD) { |
| 312 | Printf("ERROR: fESD not available"); |
| 313 | return; |
| 314 | } |
| 315 | |
| 316 | |
| 317 | const AliESDVertex *vertexESD = fESD->GetPrimaryVertexTracks(); |
| 318 | if(!vertexESD) return; |
| 319 | if(!(vertexESD->GetStatus())) return; |
| 320 | |
| 321 | // *********** MC info *************** |
| 322 | TArrayF mcVertex(3); |
| 323 | mcVertex[0]=9999.; mcVertex[1]=9999.; mcVertex[2]=9999.; |
| 324 | Float_t dNchdy=-999.; |
| 325 | |
| 326 | TParticle *part=0; |
| 327 | AliESDVertex *vertexMC=0; |
| 328 | AliStack *stack=0; |
| 329 | if (fReadMC) { |
| 330 | AliMCEventHandler *eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()); |
| 331 | if (!eventHandler) { |
| 332 | Printf("ERROR: Could not retrieve MC event handler"); |
| 333 | return; |
| 334 | } |
| 335 | |
| 336 | AliMCEvent* mcEvent = eventHandler->MCEvent(); |
| 337 | if (!mcEvent) { |
| 338 | Printf("ERROR: Could not retrieve MC event"); |
| 339 | return; |
| 340 | } |
| 341 | |
| 342 | stack = mcEvent->Stack(); |
| 343 | if (!stack) { |
| 344 | AliDebug(AliLog::kError, "Stack not available"); |
| 345 | return; |
| 346 | } |
| 347 | |
| 348 | AliHeader* header = mcEvent->Header(); |
| 349 | if (!header) { |
| 350 | AliDebug(AliLog::kError, "Header not available"); |
| 351 | return; |
| 352 | } |
| 353 | AliGenEventHeader* genHeader = header->GenEventHeader(); |
| 354 | genHeader->PrimaryVertex(mcVertex); |
| 355 | |
| 356 | |
| 357 | Int_t ngenpart = (Int_t)stack->GetNtrack(); |
| 358 | //printf("# generated particles = %d\n",ngenpart); |
| 359 | dNchdy=0; |
| 360 | for(Int_t ip=0; ip<ngenpart; ip++) { |
| 361 | part = (TParticle*)stack->Particle(ip); |
| 362 | // keep only electrons, muons, pions, kaons and protons |
| 363 | Int_t apdg = TMath::Abs(part->GetPdgCode()); |
| 364 | if(apdg!=11 && apdg!=13 && apdg!=211 && apdg!=321 && apdg!=2212) continue; |
| 365 | // reject secondaries |
| 366 | if(TMath::Sqrt((part->Vx()-mcVertex[0])*(part->Vx()-mcVertex[0])+(part->Vy()-mcVertex[1])*(part->Vy()-mcVertex[1]))>0.0010) continue; |
| 367 | // reject incoming protons |
| 368 | Double_t energy = part->Energy(); |
| 369 | if(energy>900.) continue; |
| 370 | Double_t pz = part->Pz(); |
| 371 | Double_t y = 0.5*TMath::Log((energy+pz+1.e-13)/(energy-pz+1.e-13)); |
| 372 | if(TMath::Abs(y)<1.0) dNchdy += 0.5; // count 1/2 of particles in |y|<1 |
| 373 | } |
| 374 | //printf("# primary particles = %7.1f\n",dNchdy); |
| 375 | } |
| 376 | // *********** MC info *************** |
| 377 | Double_t mcVtxPos[3]={mcVertex[0],mcVertex[1],mcVertex[2]},mcVtxSigma[3]={0,0,0}; |
| 378 | vertexMC = new AliESDVertex(mcVtxPos,mcVtxSigma); |
| 379 | |
| 380 | // *********** ESD friends *********** |
| 381 | fESD->SetESDfriend(fESDfriend); //Attach the friend to the ESD |
| 382 | // *********** ESD friends *********** |
| 383 | |
| 384 | if(!fESDfriend) printf("no ESD friend\n"); |
| 385 | |
| 386 | // |
| 387 | |
| 388 | // ********** Trigger ***************** |
| 389 | ULong64_t triggerMask; |
| 390 | ULong64_t spdFO = (1 << 14); |
| 391 | ULong64_t v0left = (1 << 11); |
| 392 | ULong64_t v0right = (1 << 12); |
| 393 | |
| 394 | triggerMask=fESD->GetTriggerMask(); |
| 395 | // MB1: SPDFO || V0L || V0R |
| 396 | Bool_t eventTriggered = (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right))); |
| 397 | //MB2: GFO && V0R |
| 398 | //triggerMask & spdFO && ((triggerMask&v0left) || (triggerMask&v0right)) |
| 399 | // ************ Trigger ****************** |
| 400 | if(!eventTriggered) return; |
| 401 | |
| 402 | // SPD vertex |
| 403 | const AliESDVertex *spdv=fESD->GetPrimaryVertexSPD(); |
| 404 | |
| 405 | |
| 406 | Int_t ntracks = fESD->GetNumberOfTracks(); |
| 407 | //printf("Tracks # = %d\n",fESD->GetNumberOfTracks()); |
| 408 | |
| 409 | fHistNtracks->Fill(ntracks); |
| 410 | // Post the data already here |
| 411 | PostData(0, fOutput); |
| 412 | |
| 413 | Int_t idet,status; Float_t xloc,zloc; |
| 414 | |
| 415 | // loop on tracks |
| 416 | for(Int_t itr=0; itr<ntracks; itr++) { |
| 417 | AliESDtrack *track = fESD->GetTrack(itr); |
| 418 | |
| 419 | if(track->GetKinkIndex(0)>0) continue; |
| 420 | |
| 421 | Int_t trkLabel = TMath::Abs(track->GetLabel()); |
| 422 | Int_t nclsITS = track->GetNcls(0); |
| 423 | |
| 424 | |
| 425 | Bool_t itsrefit=kFALSE,tpcin=kFALSE,itsfindableAcc=kFALSE; |
| 426 | if ((track->GetStatus() & AliESDtrack::kITSrefit)) itsrefit=kTRUE; |
| 427 | if ((track->GetStatus() & AliESDtrack::kTPCin)) tpcin=kTRUE; |
| 428 | |
| 429 | if(tpcin) { |
| 430 | fHistNclsITSMI->Fill(nclsITS); |
| 431 | } else { |
| 432 | fHistNclsITSSA->Fill(nclsITS); |
| 433 | } |
| 434 | |
| 435 | // TPC track in ITS acceptance |
| 436 | if(tpcin && TMath::Abs(track->Eta())<0.9 && track->GetNcls(1)>50) { |
| 437 | itsfindableAcc=kTRUE; |
| 438 | fHistPtTPC->Fill(track->Pt()); |
| 439 | fHistPhiTPC->Fill(track->Phi()); |
| 440 | } |
| 441 | |
| 442 | |
| 443 | for(Int_t layer=0; layer<6; layer++) { |
| 444 | if(TESTBIT(track->GetITSClusterMap(),layer)) { |
| 445 | if(tpcin) { |
| 446 | fHistClusterMapITSMI->Fill(layer); |
| 447 | } else { |
| 448 | fHistClusterMapITSSA->Fill(layer); |
| 449 | } |
| 450 | } |
| 451 | track->GetITSModuleIndexInfo(layer,idet,status,xloc,zloc); |
| 452 | if(tpcin) { |
| 453 | if(status==1) fHistClusterMapITSMIok->Fill(layer); |
| 454 | if(status==2) fHistClusterMapITSMIbad->Fill(layer); |
| 455 | if(status==3) fHistClusterMapITSMIskipped->Fill(layer); |
| 456 | if(status==4) fHistClusterMapITSMIoutinz->Fill(layer); |
| 457 | if(status==5) fHistClusterMapITSMInocls->Fill(layer); |
| 458 | if(status==6) fHistClusterMapITSMInorefit->Fill(layer); |
| 459 | } else { |
| 460 | if(status==1) fHistClusterMapITSSAok->Fill(layer); |
| 461 | if(status==2) fHistClusterMapITSSAbad->Fill(layer); |
| 462 | if(status==3) fHistClusterMapITSSAskipped->Fill(layer); |
| 463 | if(status==4) fHistClusterMapITSSAoutinz->Fill(layer); |
| 464 | if(status==5) fHistClusterMapITSSAnocls->Fill(layer); |
| 465 | if(status==6) fHistClusterMapITSSAnorefit->Fill(layer); |
| 466 | } |
| 467 | } |
| 468 | |
| 469 | if(itsfindableAcc) { |
| 470 | if(nclsITS==6) fHistPtITSMI6->Fill(track->Pt()); |
| 471 | if(nclsITS==5) fHistPtITSMI5->Fill(track->Pt()); |
| 472 | if(nclsITS==4) fHistPtITSMI4->Fill(track->Pt()); |
| 473 | if(nclsITS==3) fHistPtITSMI3->Fill(track->Pt()); |
| 474 | if(nclsITS==2) fHistPtITSMI2->Fill(track->Pt()); |
| 475 | if(track->HasPointOnITSLayer(0) && track->HasPointOnITSLayer(1)) |
| 476 | fHistPtITSMISPD->Fill(track->Pt()); |
| 477 | } |
| 478 | |
| 479 | |
| 480 | Int_t iITSflag=0; //ITSflag takes the value 0 if the track has no cluster assigned in the SPDs, 1 (2) if one cluster is assigned in SPD1(2), 3 if two clusters are present. Then the same adding 10,20 or 30 for SDD and 100,200 or 300 for SSD |
| 481 | |
| 482 | if(track->HasPointOnITSLayer(0)) iITSflag+=1; |
| 483 | if(track->HasPointOnITSLayer(1)) iITSflag+=2; |
| 484 | if(track->HasPointOnITSLayer(2)) iITSflag+=10; |
| 485 | if(track->HasPointOnITSLayer(3)) iITSflag+=20; |
| 486 | if(track->HasPointOnITSLayer(4)) iITSflag+=100; |
| 487 | if(track->HasPointOnITSLayer(5)) iITSflag+=200; |
| 488 | |
| 489 | if(iITSflag==333 && track->GetNcls(0)<6) |
| 490 | printf(" ERROR %d %d\n",track->GetNcls(0),track->GetLabel()); |
| 491 | |
| 492 | // number of associated ITS clusters |
| 493 | iITSflag += 1000*track->GetNcls(0); |
| 494 | |
| 495 | // if MC info and is available |
| 496 | // write the number of ITS clusters produced by this track |
| 497 | Int_t nITSclsMC=0; |
| 498 | if(fReadMC && fReadRPLabels) { |
| 499 | nITSclsMC = NumberOfITSClustersMC(trkLabel); |
| 500 | if(nITSclsMC>=0) iITSflag += 10000*nITSclsMC; |
| 501 | } |
| 502 | |
| 503 | // impact parameter to VertexTracks |
| 504 | Double_t d0z0[2],covd0z0[3]; |
| 505 | track->PropagateToDCA(vertexESD,fESD->GetMagneticField(),100.,d0z0,covd0z0); |
| 506 | if(covd0z0[0]<0. || covd0z0[2]<0.) continue; |
| 507 | |
| 508 | // if MC info is available: get particle properties |
| 509 | Float_t ptMC=-999.,pdgMC=-999.,d0MC=-999.; |
| 510 | Double_t d0z0MCv[2]={-999.,-999.},covd0z0MCv[3]={1.,1.,1.}; |
| 511 | if(fReadMC) { |
| 512 | part = (TParticle*)stack->Particle(trkLabel); |
| 513 | ptMC=part->Pt(); |
| 514 | pdgMC=part->GetPdgCode(); |
| 515 | d0MC=ParticleImpParMC(part,vertexMC,0.1*fESD->GetMagneticField()); |
| 516 | track->PropagateToDCA(vertexMC,fESD->GetMagneticField(),100.,d0z0MCv,covd0z0MCv); |
| 517 | if(covd0z0MCv[0]<0. || covd0z0MCv[2]<0.) continue; |
| 518 | |
| 519 | } |
| 520 | |
| 521 | // fill ntuple with track properties |
| 522 | if(SelectPt(track->Pt())) { |
| 523 | Float_t fillArray[19]={track->Pt(),track->Eta(),track->Phi(),d0z0[0],d0z0[1],TMath::Sqrt(covd0z0[0]),TMath::Sqrt(covd0z0[2]),ptMC,pdgMC,d0MC,d0z0MCv[0],d0z0MCv[1],TMath::Sqrt(covd0z0MCv[0]),TMath::Sqrt(covd0z0MCv[2]),(Float_t)iITSflag}; |
| 524 | fNtupleESDTracks->Fill(fillArray); |
| 525 | } |
| 526 | |
| 527 | //--------------------------------------------- |
| 528 | // AliTrackPoints: alignment checks |
| 529 | // |
| 530 | const AliTrackPointArray *array = track->GetTrackPointArray(); |
| 531 | AliTrackPoint point; |
| 532 | Int_t pointOnLayer[6]={0,0,0,0,0,0}; |
| 533 | Int_t indexAssociated[6]={-1,-1,-1,-1,-1,-1},indexExtra=-1; |
| 534 | Bool_t extra=kFALSE; |
| 535 | Int_t layerId,layerExtra=-1; |
| 536 | for(Int_t ipt=0; ipt<array->GetNPoints(); ipt++) { |
| 537 | array->GetPoint(point,ipt); |
| 538 | Float_t r = TMath::Sqrt(point.GetX()*point.GetX()+point.GetY()*point.GetY()); |
| 539 | |
| 540 | if(r>3 && r<6) { |
| 541 | layerId = 0; |
| 542 | } else if(r>6 && r<8) { |
| 543 | layerId = 1; |
| 544 | } else if(r>8 && r<18) { |
| 545 | layerId = 2; |
| 546 | } else if(r>18 && r<25) { |
| 547 | layerId = 3; |
| 548 | } else if(r>25 && r<38) { |
| 549 | layerId = 4; |
| 550 | } else if(r>38 && r<50) { |
| 551 | layerId = 5; |
| 552 | } else { |
| 553 | layerId=100; |
| 554 | } |
| 555 | |
| 556 | // only ITS points |
| 557 | if(layerId>5) continue; |
| 558 | |
| 559 | if(!point.IsExtra()) { |
| 560 | pointOnLayer[layerId]++; |
| 561 | indexAssociated[layerId]=ipt; |
| 562 | } else { |
| 563 | // this is an extra cluster |
| 564 | extra=kTRUE; |
| 565 | layerExtra=layerId; |
| 566 | indexExtra=ipt; |
| 567 | } |
| 568 | } // end loop on AliTrackPoints |
| 569 | |
| 570 | TString vtitle = spdv->GetTitle(); |
| 571 | if(!vtitle.Contains("3D")) continue; |
| 572 | |
| 573 | // SPD tracklet |
| 574 | if(indexAssociated[0]>=0 && indexAssociated[1]>=0) { |
| 575 | AliTrackPoint pointSPD1,pointSPD2; |
| 576 | array->GetPoint(pointSPD1,indexAssociated[0]); |
| 577 | array->GetPoint(pointSPD2,indexAssociated[1]); |
| 578 | Float_t phi=TMath::ATan2(pointSPD2.GetY()-pointSPD1.GetY(),pointSPD2.GetX()-pointSPD1.GetX()); |
| 579 | Float_t lambda=TMath::ATan((pointSPD2.GetZ()-pointSPD1.GetZ())/TMath::Sqrt((pointSPD2.GetX()-pointSPD1.GetX())*(pointSPD2.GetX()-pointSPD1.GetX())+(pointSPD2.GetY()-pointSPD1.GetY())*(pointSPD2.GetY()-pointSPD1.GetY()))); |
| 580 | Float_t theta=0.5*TMath::Pi()-lambda; |
| 581 | TParticle particle(211,0,0,0,0,0,TMath::Cos(phi),TMath::Sin(phi),TMath::Tan(lambda),10.,pointSPD1.GetX(),pointSPD1.GetY(),pointSPD1.GetZ(),0); |
| 582 | AliESDtrack tracklet(&particle); |
| 583 | Float_t dz[2]; |
| 584 | // distance to primary SPD (only if 3D and high multiplicity) |
| 585 | if(spdv->GetNContributors()>10) { |
| 586 | tracklet.GetDZ(spdv->GetXv(),spdv->GetYv(),spdv->GetZv(),0,dz); |
| 587 | fNtupleITSAlignSPDTracklets->Fill(phi,theta,0.5*(pointSPD1.GetZ()+pointSPD2.GetZ()),dz[0],dz[1]); |
| 588 | } |
| 589 | } |
| 590 | |
| 591 | // distance to extra |
| 592 | if(extra && spdv->GetNContributors()>4) { |
| 593 | AliTrackPoint pointExtra,pointAssociated; |
| 594 | array->GetPoint(pointAssociated,indexAssociated[layerExtra]); |
| 595 | array->GetPoint(pointExtra,indexExtra); |
| 596 | Float_t phiExtra = TMath::ATan2(pointExtra.GetY()-spdv->GetYv(),pointExtra.GetX()-spdv->GetXv()); |
| 597 | Float_t phiAssociated = TMath::ATan2(pointAssociated.GetY()-spdv->GetYv(),pointAssociated.GetX()-spdv->GetXv()); |
| 598 | Float_t rExtra = TMath::Sqrt((pointExtra.GetX()-spdv->GetXv())*(pointExtra.GetX()-spdv->GetXv())+(pointExtra.GetY()-spdv->GetYv())*(pointExtra.GetY()-spdv->GetYv())); |
| 599 | Float_t rAssociated = TMath::Sqrt((pointAssociated.GetX()-spdv->GetXv())*(pointAssociated.GetX()-spdv->GetXv())+(pointAssociated.GetY()-spdv->GetYv())*(pointAssociated.GetY()-spdv->GetYv())); |
| 600 | Float_t dzExtra[2]; |
| 601 | dzExtra[0] = (phiExtra-phiAssociated)*0.5*(rExtra+rAssociated); |
| 602 | dzExtra[1] = pointExtra.GetZ()-pointAssociated.GetZ()-(rExtra-rAssociated)*(pointAssociated.GetZ()-spdv->GetZv())/rAssociated; |
| 603 | Float_t xlocExtra=-100.,zlocExtra=-100.; |
| 604 | fNtupleITSAlignExtra->Fill(layerExtra,pointExtra.GetX(),pointExtra.GetY(),pointExtra.GetZ(),xlocExtra,zlocExtra,dzExtra[0],dzExtra[1],nclsITS); |
| 605 | } |
| 606 | |
| 607 | |
| 608 | } // end loop on tracks |
| 609 | |
| 610 | if(vertexMC) { delete vertexMC; vertexMC=0; } |
| 611 | |
| 612 | return; |
| 613 | } |
| 614 | |
| 615 | //________________________________________________________________________ |
| 616 | void AliAnalysisTaskITSTrackingCheck::Terminate(Option_t *) |
| 617 | { |
| 618 | // Draw result to the screen |
| 619 | // Called once at the end of the query |
| 620 | fOutput = dynamic_cast<TList*> (GetOutputData(0)); |
| 621 | if (!fOutput) { |
| 622 | Printf("ERROR: fOutput not available"); |
| 623 | return; |
| 624 | } |
| 625 | |
| 626 | fHistNtracks = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNtracks")); |
| 627 | fHistNclsITSMI = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNclsITSMI")); |
| 628 | fHistNclsITSSA = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNclsITSSA")); |
| 629 | fHistClusterMapITSMI = dynamic_cast<TH1F*>(fOutput->FindObject("fHistClusterMapITSMI")); |
| 630 | fHistClusterMapITSMIbad = dynamic_cast<TH1F*>(fOutput->FindObject("fHistClusterMapITSMIbad")); |
| 631 | fHistClusterMapITSMIskipped = dynamic_cast<TH1F*>(fOutput->FindObject("fHistClusterMapITSMIskipped")); |
| 632 | fHistClusterMapITSMIoutinz = dynamic_cast<TH1F*>(fOutput->FindObject("fHistClusterMapITSMIoutinz")); |
| 633 | fHistClusterMapITSMInorefit = dynamic_cast<TH1F*>(fOutput->FindObject("fHistClusterMapITSMInorefit")); |
| 634 | fHistClusterMapITSMInocls = dynamic_cast<TH1F*>(fOutput->FindObject("fHistClusterMapITSMInocls")); |
| 635 | fHistClusterMapITSSA = dynamic_cast<TH1F*>(fOutput->FindObject("fHistClusterMapITSSA")); |
| 636 | fHistClusterMapITSSAbad = dynamic_cast<TH1F*>(fOutput->FindObject("fHistClusterMapITSSAbad")); |
| 637 | fHistClusterMapITSSAskipped = dynamic_cast<TH1F*>(fOutput->FindObject("fHistClusterMapITSSAskipped")); |
| 638 | fHistClusterMapITSSAoutinz = dynamic_cast<TH1F*>(fOutput->FindObject("fHistClusterMapITSSAoutinz")); |
| 639 | fHistClusterMapITSSAnorefit = dynamic_cast<TH1F*>(fOutput->FindObject("fHistClusterMapITSSAnorefit")); |
| 640 | fHistClusterMapITSSAnocls = dynamic_cast<TH1F*>(fOutput->FindObject("fHistClusterMapITSSAnocls")); |
| 641 | fNtupleESDTracks = dynamic_cast<TNtuple*>(fOutput->FindObject("fNtupleESDTracks")); |
| 642 | fNtupleITSAlignExtra = dynamic_cast<TNtuple*>(fOutput->FindObject("fNtupleITSAlignExtra")); |
| 643 | fNtupleITSAlignSPDTracklets = dynamic_cast<TNtuple*>(fOutput->FindObject("fNtupleITSAlignSPDTracklets")); |
| 644 | |
| 645 | return; |
| 646 | } |
| 647 | //--------------------------------------------------------------------------- |
| 648 | Int_t AliAnalysisTaskITSTrackingCheck::NumberOfITSClustersMC(Int_t label) const |
| 649 | { |
| 650 | // |
| 651 | // Return number of ITS clusters produced by MC particle with given label |
| 652 | // |
| 653 | |
| 654 | AliESDInputHandlerRP *esdHRP = dynamic_cast<AliESDInputHandlerRP*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()); |
| 655 | if(!esdHRP) return -1; |
| 656 | TTree *cTree = (TTree*)esdHRP->GetTreeR("ITS"); |
| 657 | if(!cTree) return -1; |
| 658 | TClonesArray *clusters=0; // new TClonesArray("AliITSRecPoint",10000); |
| 659 | cTree->SetBranchAddress("ITSRecPoints",&clusters); |
| 660 | if(!clusters) return -1; |
| 661 | |
| 662 | AliITSRecPoint *c=0; |
| 663 | Int_t i,n,icl,lay,ilab; |
| 664 | Int_t ncls[6]={0,0,0,0,0,0}; |
| 665 | Int_t nclstot=0; |
| 666 | |
| 667 | for(i=0; i<2198; i++) { |
| 668 | cTree->GetEvent(i); |
| 669 | n=clusters->GetEntriesFast(); |
| 670 | for (icl=0; icl<n; icl++) { |
| 671 | c=(AliITSRecPoint*)clusters->UncheckedAt(icl); |
| 672 | lay=c->GetLayer(); |
| 673 | for(ilab=0;ilab<3;ilab++) { |
| 674 | if(c->GetLabel(ilab)==label) ncls[lay]++; |
| 675 | } |
| 676 | } |
| 677 | } |
| 678 | for(i=0;i<6;i++) { if(ncls[i]) nclstot++; } |
| 679 | |
| 680 | return nclstot; |
| 681 | } |
| 682 | //--------------------------------------------------------------------------- |
| 683 | Double_t AliAnalysisTaskITSTrackingCheck::ParticleImpParMC(TParticle *part, |
| 684 | AliESDVertex *vert, |
| 685 | Double_t bzT) const |
| 686 | { |
| 687 | // |
| 688 | // Return the MC value of the impact parameter |
| 689 | // |
| 690 | |
| 691 | Double_t vx=part->Vx()-vert->GetX(); |
| 692 | Double_t vy=part->Vy()-vert->GetY(); |
| 693 | |
| 694 | Double_t pt=part->Pt(); |
| 695 | Double_t px=part->Px(); |
| 696 | Double_t py=part->Py(); |
| 697 | Double_t charge = (part->GetPdgCode()>0. ? 1. : -1.); |
| 698 | if(TMath::Abs(part->GetPdgCode())<100) charge*=-1.; |
| 699 | |
| 700 | if(px<0.000001) px=0.000001; |
| 701 | Double_t rAnd=((10./2.99792458)*pt/bzT)*100.; |
| 702 | Double_t center[3],d0; |
| 703 | center[0]=vx-(1./charge)*rAnd*(py/pt); |
| 704 | center[1]=vy+(1./charge)*rAnd*(px/pt); |
| 705 | center[2]=TMath::Sqrt(center[0]*center[0]+center[1]*center[1]); |
| 706 | d0 = -center[2]+rAnd; |
| 707 | |
| 708 | return d0; |
| 709 | } |
| 710 | //--------------------------------------------------------------------------- |
| 711 | Bool_t AliAnalysisTaskITSTrackingCheck::SelectPt(Double_t pt) |
| 712 | { |
| 713 | // |
| 714 | // Keep only tracks in given pt bins |
| 715 | // |
| 716 | Double_t ptlower[10]={0.29,0.49,0.75,0.9,1.9,3.5,6.5,9.,19.,27.}; |
| 717 | Double_t ptupper[10]={0.31,0.51,0.85,1.1,2.1,4.5,7.5,11.,21.,33.}; |
| 718 | |
| 719 | for(Int_t i=0; i<10; i++) { |
| 720 | if(pt<ptlower[i] && pt<ptupper[i]) { |
| 721 | fCountsPerPtBin[i]++; |
| 722 | return kTRUE; |
| 723 | } |
| 724 | } |
| 725 | return kFALSE; |
| 726 | } |
| 727 | |
| 728 | |
| 729 | |
| 730 | |
git://git.uio.no / u / mrichter / AliRoot.git / blame