PWG1/ITS/AliAlignmentDataFilterITS.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 extract from ESD tracks the AliTrackPointArrays
  19 // with ITS points for selected tracks. This are the input data for alignment
  20 //
  21 // Author: A.Dainese, andrea.dainese@pd.infn.it
  22 /////////////////////////////////////////////////////////////
  23
  24 #include <TTree.h>
  25 #include <TChain.h>
  26 #include <TNtuple.h>
  27 #include <TList.h>
  28 #include <TH1F.h>
  29 #include <TH2F.h>
  30 #include <TMap.h>
  31 #include <TVector3.h>
  32 #include <TGeoManager.h>
  33
  34 #include "AliLog.h"
  35 #include "AliGeomManager.h"
  36 #include "AliITSReconstructor.h"
  37 #include "AliITSgeomTGeo.h"
  38 #include "AliTrackPointArray.h"
  39 #include "AliESDInputHandler.h"
  40 #include "AliESDVertex.h"
  41 #include "AliESDtrack.h"
  42 #include "AliESDEvent.h"
  43 #include "AliESDfriend.h"
  44 #include "AliAnalysisTask.h"
  45 #include "AliAnalysisManager.h"
  46 #include "AliAlignmentDataFilterITS.h"
  47
  48
  49 ClassImp(AliAlignmentDataFilterITS)
  50
  51
  52 //________________________________________________________________________
  53 AliAlignmentDataFilterITS::AliAlignmentDataFilterITS(const char *name):
  54 AliAnalysisTask(name,"task"),
  55 fOnlySPDFO(kFALSE),
  56 fGeometryFileName("geometry.root"),
  57 fITSRecoParam(0),
  58 fESD(0),
  59 fESDfriend(0),
  60 fListOfHistos(0),
  61 fspTree(0),
  62 fHistNevents(0),
  63 fHistNpoints(0),
  64 fHistPt(0),
  65 fHistLayer0(0),
  66 fHistLayer1(0),
  67 fHistLayer2(0),
  68 fHistLayer3(0),
  69 fHistLayer4(0),
  70 fHistLayer5(0),
  71 fntExtra(0),
  72 fntCosmicMatching(0)
  73 {
  74   // Constructor
  75
  76   // Define input and output slots here
  77   // Input slot #0 works with a TChain
  78   DefineInput(0, TChain::Class());
  79   // Output slot #0 writes into a TTree
  80   DefineOutput(0,TTree::Class());  //My private output
  81   // Output slot #1 writes into a TList
  82   DefineOutput(1,TList::Class());  //My private output
  83 }
  84
  85 //________________________________________________________________________
  86 AliAlignmentDataFilterITS::~AliAlignmentDataFilterITS()
  87 {
  88   // Destructor
  89   if (fListOfHistos) {
  90     delete fListOfHistos;
  91     fListOfHistos = 0;
  92   }
  93   if (fspTree) {
  94     delete fspTree;
  95     fspTree = 0;
  96   }
  97   if (fHistNevents) {
  98     delete fHistNevents;
  99     fHistNevents = 0;
 100   }
 101   if (fHistNpoints) {
 102     delete fHistNpoints;
 103     fHistNpoints = 0;
 104   }
 105   if (fHistPt) {
 106     delete fHistPt;
 107     fHistPt = 0;
 108   }
 109   if (fHistLayer0) {
 110     delete fHistLayer0;
 111     fHistLayer0 = 0;
 112   }
 113   if (fHistLayer1) {
 114     delete fHistLayer1;
 115     fHistLayer1 = 0;
 116   }
 117   if (fHistLayer2) {
 118     delete fHistLayer2;
 119     fHistLayer2 = 0;
 120   }
 121   if (fHistLayer3) {
 122     delete fHistLayer3;
 123     fHistLayer3 = 0;
 124   }
 125   if (fHistLayer4) {
 126     delete fHistLayer4;
 127     fHistLayer4 = 0;
 128   }
 129   if (fHistLayer5) {
 130     delete fHistLayer5;
 131     fHistLayer5 = 0;
 132   }
 133   if (fntExtra) {
 134     delete fntExtra;
 135     fntExtra = 0;
 136   }
 137   if (fntCosmicMatching) {
 138     delete fntCosmicMatching;
 139     fntCosmicMatching = 0;
 140   }
 141 }
 142
 143 //________________________________________________________________________
 144 void AliAlignmentDataFilterITS::ConnectInputData(Option_t *)
 145 {
 146   // Connect ESD
 147   // Called once
 148
 149   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
 150   if(!tree) {
 151     printf("ERROR: Could not read chain from input slot 0\n");
 152   } else {
 153     // Get the OCDB path and the list of OCDB objects used for reco
 154     TMap *cdbMap = (TMap*)(tree->GetTree()->GetUserInfo())->FindObject("cdbMap");
 155     TList *cdbList = (TList*)(tree->GetTree()->GetUserInfo())->FindObject("cdbList");
 156
 157     //cdbList->Print();
 158     // write the list to the user info of the output tree
 159     if(!fspTree) {
 160       printf("ERROR: fspTree does not exist\n");
 161     } else {
 162       TMap *cdbMapCopy = new TMap(cdbMap->GetEntries());
 163       cdbMapCopy->SetOwner(1);
 164       cdbMapCopy->SetName("cdbMap");
 165       TIter iter1(cdbMap->GetTable());
 166
 167       TPair* pair = 0;
 168       while((pair = dynamic_cast<TPair*> (iter1.Next()))){
 169         TObjString* keyStr = dynamic_cast<TObjString*> (pair->Key());
 170         TObjString* valStr = dynamic_cast<TObjString*> (pair->Value());
 171         cdbMapCopy->Add(new TObjString(keyStr->GetName()), new TObjString(valStr->GetName()));
 172       }
 173
 174       TList *cdbListCopy = new TList();
 175       cdbListCopy->SetOwner(1);
 176       cdbListCopy->SetName("cdbList");
 177
 178       TIter iter2(cdbList);
 179
 180       TObjString* cdbEntry=0;
 181       while((cdbEntry =(TObjString*)(iter2.Next()))) {
 182         cdbListCopy->Add(new TObjString(*cdbEntry));
 183       }
 184       cdbListCopy->Print();
 185
 186
 187       fspTree->GetUserInfo()->Add(cdbMapCopy);
 188       fspTree->GetUserInfo()->Add(cdbListCopy);
 189     }
 190
 191     // Disable all branches and enable only the needed ones
 192     tree->SetBranchStatus("fTriggerMask", 1);
 193     tree->SetBranchStatus("fSPDVertex*", 1);
 194     tree->SetBranchStatus("ESDfriend*", 1);
 195     tree->SetBranchAddress("ESDfriend.",&fESDfriend);
 196     tree->SetBranchStatus("fSPDMult*", 1);
 197
 198     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 199
 200     if(!esdH) {
 201       printf("ERROR: Could not get ESDInputHandler\n");
 202     } else {
 203       fESD = esdH->GetEvent();
 204     }
 205   }
 206
 207   return;
 208 }
 209
 210 //________________________________________________________________________
 211 void AliAlignmentDataFilterITS::Init()
 212 {
 213   // Initialization
 214
 215   return;
 216 }
 217
 218 //________________________________________________________________________
 219 void AliAlignmentDataFilterITS::CreateOutputObjects()
 220 {
 221   // Create the output container
 222   //
 223
 224   // load the geometry
 225   if(!gGeoManager) {
 226     printf("AliAlignmentDataFilterITS::CreateOutputObjects(): loading geometry from %s\n",fGeometryFileName.Data());
 227     AliGeomManager::LoadGeometry(fGeometryFileName.Data());
 228     if(!gGeoManager) {
 229       printf("AliAlignmentDataFilterITS::CreateOutputObjects(): no geometry loaded \n");
 230       return;
 231     }
 232   }
 233
 234   // Several histograms are more conveniently managed in a TList
 235   fListOfHistos = new TList();
 236   fListOfHistos->SetOwner();
 237
 238   fHistNevents = new TH1F("fHistNevents", "Number of processed events; N events; bin",5,-0.5,4.5);
 239   fHistNevents->Sumw2();
 240   fHistNevents->SetMinimum(0);
 241   fListOfHistos->Add(fHistNevents);
 242
 243   fHistNpoints = new TH1F("fHistNpoints", "Number of AliTrackPoints per track; N points; tracks",25,-0.5,24.5);
 244   fHistNpoints->Sumw2();
 245   fHistNpoints->SetMinimum(0);
 246   fListOfHistos->Add(fHistNpoints);
 247
 248   fHistPt = new TH1F("fHistPt", "p_{t} of tracks; p_{t} [GeV/c]; tracks",100,0,50);
 249   fHistPt->Sumw2();
 250   fHistPt->SetMinimum(0);
 251   fListOfHistos->Add(fHistPt);
 252
 253
 254   Float_t zmax=14.;
 255   Int_t nbinsphi=20,nbinsz=4;
 256   fHistLayer0 = new TH2F("fHistLayer0","Points in layer inner SPD; global   #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax);
 257   fListOfHistos->Add(fHistLayer0);
 258   zmax=14.;
 259   nbinsphi=40;nbinsz=4;
 260   fHistLayer1 = new TH2F("fHistLayer1","Points in layer outer SPD; global   #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax);
 261   fListOfHistos->Add(fHistLayer1);
 262   zmax=22.;
 263   nbinsphi=14;nbinsz=6;
 264   fHistLayer2 = new TH2F("fHistLayer2","Points in layer inner SDD; global   #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax);
 265   fListOfHistos->Add(fHistLayer2);
 266   zmax=29.5;
 267   nbinsphi=22;nbinsz=8;
 268   fHistLayer3 = new TH2F("fHistLayer3","Points in layer outer SDD; global   #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax);
 269   fListOfHistos->Add(fHistLayer3);
 270   zmax=45.;
 271   nbinsphi=34;nbinsz=23;
 272   fHistLayer4 = new TH2F("fHistLayer4","Points in layer inner SSD; global   #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax);
 273   fListOfHistos->Add(fHistLayer4);
 274   zmax=51.;
 275   nbinsphi=38;nbinsz=26;
 276   fHistLayer5 = new TH2F("fHistLayer5","Points in layer outer SSD; global   #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax);
 277   fListOfHistos->Add(fHistLayer5);
 278
 279
 280   fntExtra = new TNtuple("fntExtra","extra clusters in ITS","ncls:layer:ladder:volid:phi:x:y:z:xloc:zloc:dxy:dz:d0mu:z0mu");
 281   fListOfHistos->Add(fntExtra);
 282
 283   fntCosmicMatching = new TNtuple("fntCosmicMatching","cosmic tracks matching in ITS","ncls1:ncls2:pt1:pt2:sigmad01:sigmad02:sigmaz01:sigmaz02:dxy:dz:phimu:thetamu:d0mu:z0mu");
 284   fListOfHistos->Add(fntCosmicMatching);
 285
 286   fspTree = new TTree("spTree","Tree with ITS track points");
 287   const AliTrackPointArray *array = 0;
 288   Float_t curv,curverr,runNumber;
 289   const TObjString *itsaligndata = 0;
 290   fspTree->Branch("SP","AliTrackPointArray",&array);
 291   fspTree->Branch("curv",&curv);
 292   fspTree->Branch("curverr",&curverr);
 293   fspTree->Branch("run",&runNumber);
 294   fspTree->Branch("ITSAlignData",&itsaligndata);
 295
 296
 297   return;
 298 }
 299
 300 //________________________________________________________________________
 301 void AliAlignmentDataFilterITS::Exec(Option_t */*option*/)
 302 {
 303   // Execute analysis for current event:
 304   // write ITS AliTrackPoints for selected tracks to fspTree
 305
 306   // check the geometry
 307   if(!gGeoManager) {
 308     printf("AliAlignmentDataFilterITS::Exec(): no geometry loaded \n");
 309     return;
 310   }
 311
 312   // check if we have AliITSRecoParam
 313   if(!GetRecoParam()) {
 314     if(!fITSRecoParam) {
 315       printf("AliAlignmentDataFilterITS::Exec(): no AliITSRecoParam\n");
 316       return;
 317     }
 318   }
 319
 320
 321   if(!fESD) {
 322     printf("AliAlignmentDataFilterITS::Exec(): no ESD \n");
 323     return;
 324   }
 325   if(!fESDfriend) {
 326     printf("AliAlignmentDataFilterITS::Exec(): no ESDfriend \n");
 327     return;
 328   }
 329   // attach ESDfriend
 330   fESD->SetESDfriend(fESDfriend);
 331
 332   // Post the data for slot 0
 333   fHistNevents->Fill(0);
 334
 335
 336   // write field value to spTree UserInfo
 337   if(!((fspTree->GetUserInfo())->FindObject("BzkGauss"))) {
 338     Double_t bz=fESD->GetMagneticField();
 339     TString bzString; bzString+=bz;
 340     TObjString *bzObjString = new TObjString(bzString);
 341     TList *bzList = new TList();
 342     bzList->SetOwner(1);
 343     bzList->SetName("BzkGauss");
 344     bzList->Add(bzObjString);
 345     fspTree->GetUserInfo()->Add(bzList);
 346   }
 347
 348
 349   // Process event as Cosmic or Collision
 350   //if(esd->GetEventType()== ???? ) {
 351   printf("AliAlignmentDataFilterITS::Exec(): MOVE ASAP TO esd->GetEventType() !\n");
 352   if(GetRecoParam()->GetAlignFilterCosmics()) {
 353     FilterCosmic(fESD);
 354   } else {
 355     FilterCollision(fESD);
 356   }
 357
 358   PostData(1,fListOfHistos);
 359
 360   return;
 361 }
 362
 363 //________________________________________________________________________
 364 void AliAlignmentDataFilterITS::FilterCosmic(const AliESDEvent *esd)
 365 {
 366   // Extract ITS AliTrackPoints for Cosmics (check angular matching
 367   // of top and bottom track, merge the two tracks, if requested)
 368   //
 369
 370   // Set branch addresses for space points tree
 371   AliTrackPointArray *arrayForTree=0;
 372   Float_t curv,curverr,runNumber;
 373   TObjString *itsaligndata=0;
 374   fspTree->SetBranchAddress("SP",&arrayForTree);
 375   fspTree->SetBranchAddress("curv",&curv);
 376   fspTree->SetBranchAddress("curverr",&curverr);
 377   fspTree->SetBranchAddress("run",&runNumber);
 378   fspTree->SetBranchAddress("ITSAlignData",&itsaligndata);
 379
 380
 381   runNumber = (Float_t)esd->GetRunNumber();
 382
 383   TTree* esdTree = dynamic_cast<TTree*> (GetInputData(0));
 384   // Get the list of OCDB objects used for reco
 385   TList *cdbList = (TList*)(esdTree->GetTree()->GetUserInfo())->FindObject("cdbList");
 386   TIter iter2(cdbList);
 387   TObjString* cdbEntry=0;
 388   TString cdbEntryString;
 389   while((cdbEntry =(TObjString*)(iter2.Next()))) {
 390   cdbEntryString = cdbEntry->GetString();
 391   if(cdbEntryString.Contains("ITS/Align/Data"))
 392     itsaligndata = new TObjString(*cdbEntry);
 393   }
 394
 395
 396   TString triggeredClass = fESD->GetFiredTriggerClasses();
 397   if(fOnlySPDFO && !triggeredClass.Contains("C0SCO-ABCE-NOPF-CENT")) return;
 398
 399
 400   Int_t ntracks = esd->GetNumberOfTracks();
 401   if(ntracks<2) return;
 402
 403   if(esd->GetPrimaryVertexSPD()->GetNContributors()<0) return;
 404
 405   Double_t vtxpos[3]; esd->GetPrimaryVertexSPD()->GetXYZ(vtxpos);
 406
 407   Int_t *goodtracksArray = new Int_t[ntracks];
 408   Float_t *phiArray = new Float_t[ntracks];
 409   Float_t *thetaArray = new Float_t[ntracks];
 410   Int_t *nclsArray = new Int_t[ntracks];
 411   Int_t ngt=0;
 412   Int_t itrack=0;
 413   for (itrack=0; itrack < ntracks; itrack++) {
 414     AliESDtrack *track = esd->GetTrack(itrack);
 415     if (!track) continue;
 416
 417
 418     if(track->GetNcls(0)<GetRecoParam()->GetAlignFilterMinITSPoints()) continue;
 419
 420     if((GetRecoParam()->GetAlignFilterOnlyITSSATracks()) && track->GetNcls(1)>0) continue;
 421     if((GetRecoParam()->GetAlignFilterOnlyITSTPCTracks()) && track->GetNcls(1)==0) continue;
 422
 423     Float_t phi = track->GetAlpha()+TMath::ASin(track->GetSnp());
 424     Float_t theta = 0.5*TMath::Pi()-TMath::ATan(track->GetTgl());
 425
 426     if(track->Pt()<GetRecoParam()->GetAlignFilterMinPt() ||
 427        track->Pt()>GetRecoParam()->GetAlignFilterMaxPt()) continue;
 428
 429     goodtracksArray[ngt] = itrack;
 430     phiArray[ngt]        = phi;
 431     thetaArray[ngt]      = theta;
 432     nclsArray[ngt]       = track->GetNcls(0);
 433     ngt++;
 434   }
 435
 436   if(ngt<2) {
 437     delete [] goodtracksArray; goodtracksArray=0;
 438     delete [] phiArray; phiArray=0;
 439     delete [] thetaArray; thetaArray=0;
 440     delete [] nclsArray; nclsArray=0;
 441     return;
 442   }
 443
 444   // check matching of the two tracks from the muon
 445   Float_t min = 10000000.;
 446   Int_t maxCls = 0;
 447   Int_t good1 = -1, good2 = -1;
 448   for(Int_t itr1=0; itr1<ngt-1; itr1++) {
 449     for(Int_t itr2=itr1+1; itr2<ngt; itr2++) {
 450       Float_t deltatheta = TMath::Abs(TMath::Pi()-thetaArray[itr1]-thetaArray[itr2]);
 451       if(deltatheta>GetRecoParam()->GetAlignFilterMaxMatchingAngle()) continue;
 452       Float_t deltaphi = TMath::Abs(TMath::Abs(phiArray[itr1]-phiArray[itr2])-TMath::Pi());
 453       if(deltaphi>GetRecoParam()->GetAlignFilterMaxMatchingAngle()) continue;
 454       if(nclsArray[itr1]+nclsArray[itr2] > maxCls) {
 455         maxCls = nclsArray[itr1]+nclsArray[itr2];
 456         min = deltatheta+deltaphi;
 457         good1 = goodtracksArray[itr1];
 458         good2 = goodtracksArray[itr2];
 459       } else if(nclsArray[itr1]+nclsArray[itr2] == maxCls) {
 460         if(deltatheta+deltaphi < min) {
 461           min = deltatheta+deltaphi;
 462           good1 = goodtracksArray[itr1];
 463           good2 = goodtracksArray[itr2];
 464         }
 465       }
 466     }
 467   }
 468
 469   delete [] goodtracksArray; goodtracksArray=0;
 470   delete [] phiArray; phiArray=0;
 471   delete [] thetaArray; thetaArray=0;
 472   delete [] nclsArray; nclsArray=0;
 473
 474   if(good1<0) return;
 475   AliDebug(2,"ok track matching");
 476
 477   // track1 will be the inward track (top)
 478   // track2 the outward (bottom)
 479   AliESDtrack *track1=0;
 480   AliESDtrack *track2=0;
 481   AliESDtrack *track = esd->GetTrack(good1);
 482   if(track->Py()>0) {
 483     track1 = esd->GetTrack(good1);
 484     track2 = esd->GetTrack(good2);
 485   } else {
 486     track1 = esd->GetTrack(good2);
 487     track2 = esd->GetTrack(good1);
 488   }
 489
 490   AliTrackPoint point;
 491   const AliTrackPointArray *array=0;
 492   Int_t ipt,volId,modId,layerId,lay,lad,det;
 493   Int_t jpt=0;
 494   Bool_t layerOK[6][2];
 495   Int_t nclsTrk[2]={0,0};
 496
 497   for(Int_t l1=0;l1<6;l1++) for(Int_t l2=0;l2<2;l2++) layerOK[l1][l2]=kFALSE;
 498
 499   for(itrack=0; itrack<2; itrack++) {
 500     if(itrack==0) {
 501       track = track1;
 502     } else {
 503       track = track2;
 504     }
 505     array = track->GetTrackPointArray();
 506     if(!array) {
 507       AliWarning("No tracks points avaialble");
 508       continue;
 509     }
 510     for(ipt=0; ipt<array->GetNPoints(); ipt++) {
 511       array->GetPoint(point,ipt);
 512       volId = point.GetVolumeID();
 513       layerId = AliGeomManager::VolUIDToLayer(volId,modId);
 514       AliDebug(2,Form("%d %d\n",ipt,layerId-1));
 515       if(point.IsExtra() &&
 516          (GetRecoParam()->GetAlignFilterSkipExtra())) continue;
 517       if(layerId>6) continue;
 518       if(!GetRecoParam()->GetAlignFilterUseLayer(layerId-1)) continue;
 519       // check minAngleWrtITSModulePlanes
 520       Double_t p[3]; track->GetDirection(p);
 521       TVector3 pvec(p[0],p[1],p[2]);
 522       Double_t rot[9]; AliGeomManager::GetOrigRotation(volId,rot);
 523       TVector3 normvec(rot[1],rot[4],rot[7]);
 524       Double_t angle = pvec.Angle(normvec);
 525       if(angle>0.5*TMath::Pi()) angle = TMath::Pi()-angle;
 526       angle = 0.5*TMath::Pi()-angle;
 527       if(angle<GetRecoParam()->GetAlignFilterMinAngleWrtModulePlanes()) continue;
 528       layerOK[layerId-1][itrack]=kTRUE;
 529       jpt++;
 530       nclsTrk[itrack]++;
 531     }
 532   }
 533   AliDebug(2,Form("nClsTrk1 %d nClsTrk2 %d\n",nclsTrk[0],nclsTrk[1]));
 534
 535   // read ITS cluster maps
 536   Int_t map1[6],map2[6];
 537   for(Int_t ilay=0;ilay<6;ilay++) {
 538     map1[ilay]=0; map2[ilay]=0;
 539     if(track1->HasPointOnITSLayer(ilay)) map1[ilay]=1;
 540     if(track2->HasPointOnITSLayer(ilay)) map2[ilay]=1;
 541   }
 542   AliDebug(2,Form("ITS map 1: %d %d %d %d %d %d pt %f\n",map1[0],map1[1],map1[2],map1[3],map1[4],map1[5],track1->Pt()));
 543   AliDebug(2,Form("ITS map 2: %d %d %d %d %d %d pt %f\n",map2[0],map2[1],map2[2],map2[3],map2[4],map2[5],track2->Pt()));
 544   Int_t idx1[12],idx2[12];
 545   track1->GetITSclusters(idx1);
 546   track2->GetITSclusters(idx2);
 547   AliDebug(2,Form("cls idx 1 %d %d %d %d %d %d %d %d %d %d %d %d\n",idx1[0],idx1[1],idx1[2],idx1[3],idx1[4],idx1[5],idx1[6],idx1[7],idx1[8],idx1[9],idx1[10],idx1[11]));
 548   AliDebug(2,Form("cls idx 2 %d %d %d %d %d %d %d %d %d %d %d %d\n",idx2[0],idx2[1],idx2[2],idx2[3],idx2[4],idx2[5],idx2[6],idx2[7],idx2[8],idx2[9],idx2[10],idx2[11]));
 549
 550
 551   if(jpt<GetRecoParam()->GetAlignFilterMinITSPointsMerged()) return;
 552   AliDebug(2,Form(" Total points %d, accepted\n",jpt));
 553   fHistNpoints->Fill(jpt);
 554   fHistPt->Fill(0.5*(track1->Pt()+track2->Pt()));
 555
 556   Float_t d0z0mu[2];
 557   track1->GetDZ(0,0,0,esd->GetMagneticField(),d0z0mu);
 558   //printf("d0mu %f  z0mu %f\n",d0z0mu[0],d0z0mu[1]);
 559
 560   Float_t dzOverlap[2];
 561   Float_t curvArray[2],curverrArray[2];
 562   Double_t globExtra[3],locExtra[3];
 563   if(GetRecoParam()->GetAlignFilterCosmicMergeTracks())
 564     arrayForTree = new AliTrackPointArray(jpt);
 565
 566   jpt=0;
 567   for(itrack=0; itrack<2; itrack++) {
 568     if(itrack==0) {
 569       track = track1;
 570     } else {
 571       track = track2;
 572     }
 573     curvArray[itrack] = track->GetC(esd->GetMagneticField());
 574     curverrArray[itrack] = TMath::Sqrt(track->GetSigma1Pt2())*track->GetC(esd->GetMagneticField())/track->OneOverPt();
 575
 576     if(!(GetRecoParam()->GetAlignFilterCosmicMergeTracks())) {
 577       jpt=0;
 578       arrayForTree = new AliTrackPointArray(nclsTrk[itrack]);
 579     }
 580     array = track->GetTrackPointArray();
 581     for(ipt=0; ipt<array->GetNPoints(); ipt++) {
 582       array->GetPoint(point,ipt);
 583       volId = point.GetVolumeID();
 584       layerId = AliGeomManager::VolUIDToLayer(volId,modId);
 585       if(layerId>6 || !layerOK[layerId-1][itrack]) continue;
 586       arrayForTree->AddPoint(jpt,&point);
 587       jpt++;
 588       switch(layerId) {
 589       case 1:
 590         fHistLayer0->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 591         break;
 592       case 2:
 593         fHistLayer1->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 594         break;
 595       case 3:
 596         fHistLayer2->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 597         break;
 598       case 4:
 599         fHistLayer3->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 600         break;
 601       case 5:
 602         fHistLayer4->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 603         break;
 604       case 6:
 605         fHistLayer5->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 606         break;
 607       }
 608       // Post the data for slot 0
 609       if(jpt==1) PostData(1,fListOfHistos); // only if this is the first points
 610       if(!point.IsExtra() ||
 611          !(GetRecoParam()->GetAlignFilterFillQANtuples())) continue;
 612       nclsTrk[itrack]--;
 613       for(Int_t ll=1;ll<layerId;ll++) modId+=AliITSgeomTGeo::GetNLadders(ll)*AliITSgeomTGeo::GetNDetectors(ll);
 614       AliITSgeomTGeo::GetModuleId(modId,lay,lad,det);
 615       globExtra[0]=point.GetX();
 616       globExtra[1]=point.GetY();
 617       globExtra[2]=point.GetZ();
 618       AliITSgeomTGeo::GlobalToLocal(lay,lad,det,globExtra,locExtra);
 619       //printf("%d %d %d %d %d  %f %f %f\n",volId,modId,lay,lad,det,locExtra[0],locExtra[1],locExtra[2]);
 620       track->GetDZ(point.GetX(),point.GetY(),point.GetZ(),esd->GetMagneticField(),dzOverlap);
 621       AliTrackPoint pointT;
 622       Float_t radius,radiusT,phiv,phivT,thetav,thetavT;
 623       for(Int_t lll=0;lll<ipt;lll++) {
 624         array->GetPoint(pointT,lll);
 625         Int_t layerIdT = AliGeomManager::VolUIDToLayer(pointT.GetVolumeID(),modId);
 626         if(layerIdT!=layerId) continue;
 627         radius=TMath::Sqrt((point.GetX()-vtxpos[0])*(point.GetX()-vtxpos[0])+(point.GetY()-vtxpos[1])*(point.GetY()-vtxpos[1]));
 628         radiusT=TMath::Sqrt((pointT.GetX()-vtxpos[0])*(pointT.GetX()-vtxpos[0])+(pointT.GetY()-vtxpos[1])*(pointT.GetY()-vtxpos[1]));
 629         phiv=TMath::ATan2(point.GetY()-vtxpos[1],point.GetX()-vtxpos[0]);
 630         phivT=TMath::ATan2(pointT.GetY()-vtxpos[1],pointT.GetX()-vtxpos[0]);
 631         if(TMath::Abs(point.GetZ()-vtxpos[2])<0.00001 || TMath::Abs(pointT.GetZ()-vtxpos[2])<0.00001) continue;
 632         thetav=TMath::ATan(radius/(point.GetZ()-vtxpos[2]));
 633         thetavT=TMath::ATan(radiusT/(pointT.GetZ()-vtxpos[2]));
 634         dzOverlap[0]=(Float_t)((phivT-phiv)*0.5*(radiusT+radius));
 635         if(TMath::Abs(TMath::Tan(0.5*(thetav+thetavT)))<0.00001) continue;
 636         dzOverlap[1]=(Float_t)((pointT.GetZ()-point.GetZ())-(radiusT-radius)/TMath::Tan(0.5*(thetav+thetavT)));
 637         fntExtra->Fill((Float_t)nclsTrk[itrack],(Float_t)(layerId-1),lad,volId,TMath::ATan2(point.GetY(),point.GetX()),point.GetX(),point.GetY(),point.GetZ(),locExtra[0],locExtra[2],dzOverlap[0],dzOverlap[1],d0z0mu[0],d0z0mu[1]);
 638       }
 639     }
 640
 641     if(!(GetRecoParam()->GetAlignFilterCosmicMergeTracks())) {
 642       curv = curvArray[itrack];
 643       curverr = curverrArray[itrack];
 644       fspTree->Fill();
 645     }
 646   }
 647
 648   if(GetRecoParam()->GetAlignFilterCosmicMergeTracks()) {
 649     curv = 0.5*(curvArray[0]+curvArray[1]);
 650     curverr = 0.5*TMath::Sqrt(curverrArray[0]*curverrArray[0]+curverrArray[1]*curverrArray[1]);
 651     fspTree->Fill();
 652   }
 653   PostData(0,fspTree);
 654
 655   if(!(GetRecoParam()->GetAlignFilterFillQANtuples())) return;
 656   // fill ntuple with track-to-track matching
 657   Float_t phimu,thetamu,phiout,thetaout,dphi,dtheta,rotymu,rotyout,droty;
 658   Float_t d0[2],z0[2];
 659   Float_t sigmad0[2],sigmaz0[2];
 660   phimu = track1->GetAlpha()+TMath::ASin(track1->GetSnp());
 661   thetamu = 0.5*TMath::Pi()-TMath::ATan(track1->GetTgl());
 662   phiout = track2->GetAlpha()+TMath::ASin(track2->GetSnp());
 663   thetaout = 0.5*TMath::Pi()-TMath::ATan(track2->GetTgl());
 664   rotymu = TMath::ATan2(track1->Px(),track1->Pz());
 665   rotyout = TMath::ATan2(track2->Px(),track2->Pz());
 666
 667   dphi = phimu - (phiout+TMath::Pi());
 668   dtheta = thetamu - (TMath::Pi()-thetaout);
 669   if(rotymu>0) {
 670     droty = rotymu - (rotyout+TMath::Pi());
 671   } else {
 672     droty = rotymu - (rotyout-TMath::Pi());
 673   }
 674
 675   Double_t alpha = TMath::ATan2(track1->Py(),track1->Px());
 676
 677   track1->Propagate(alpha,0.,esd->GetMagneticField());
 678   track2->Propagate(alpha,0.,esd->GetMagneticField());
 679   d0[0] = track1->GetY();
 680   z0[0] = track1->GetZ();
 681   d0[1] = track2->GetY();
 682   z0[1] = track2->GetZ();
 683   Float_t dxy = -(d0[0]-d0[1]);
 684   Float_t dz  = z0[0]-z0[1];
 685   sigmad0[0] = TMath::Sqrt(track1->GetSigmaY2());
 686   sigmaz0[0] = TMath::Sqrt(track1->GetSigmaZ2());
 687   sigmad0[1] = TMath::Sqrt(track2->GetSigmaY2());
 688   sigmaz0[1] = TMath::Sqrt(track2->GetSigmaZ2());
 689   /*
 690   Double_t xyz1atxl0[3],xyz1atxl1[3],xyz2atxl0[3],xyz2atxl1[3];
 691   track1->GetXYZAt(0.,esd->GetMagneticField(),xyz1atxl0);
 692   track1->GetXYZAt(1.,esd->GetMagneticField(),xyz1atxl1);
 693   track2->GetXYZAt(0.,esd->GetMagneticField(),xyz2atxl0);
 694   track2->GetXYZAt(1.,esd->GetMagneticField(),xyz2atxl1);
 695   Float_t x1aty0 = (xyz1atxl0[0]*xyz1atxl1[1]-xyz1atxl0[1]*xyz1atxl1[0])/(xyz1atxl1[1]-xyz1atxl0[1]);
 696   Float_t x2aty0 = (xyz2atxl0[0]*xyz2atxl1[1]-xyz2atxl0[1]*xyz2atxl1[0])/(xyz2atxl1[1]-xyz2atxl0[1]);
 697   Float_t dxaty0 = x1aty0-x2aty0;
 698   */
 699   fntCosmicMatching->Fill((Float_t)nclsTrk[0],(Float_t)nclsTrk[1],track1->Pt(),track2->Pt(),sigmad0[0],sigmad0[1],sigmaz0[0],sigmaz0[1],dxy,dz,phimu,thetamu,TMath::Abs(d0z0mu[0]),d0z0mu[1]);
 700
 701   return;
 702 }
 703
 704 //________________________________________________________________________
 705 void AliAlignmentDataFilterITS::FilterCollision(const AliESDEvent *esd)
 706 {
 707   // Extract ITS AliTrackPoints for Cosmics (check angular matching
 708   // of top and bottom track, merge the two tracks, if requested)
 709   //
 710
 711   // Set branch addresses for space points tree
 712   AliTrackPointArray *arrayForTree=0;
 713   Float_t curv,curverr,runNumber;
 714   TObjString *itsaligndata=0;
 715   fspTree->SetBranchAddress("SP",&arrayForTree);
 716   fspTree->SetBranchAddress("curv",&curv);
 717   fspTree->SetBranchAddress("curverr",&curverr);
 718   fspTree->SetBranchAddress("run",&runNumber);
 719   fspTree->SetBranchAddress("ITSAlignData",&itsaligndata);
 720
 721
 722   runNumber = (Float_t)esd->GetRunNumber();
 723
 724   TTree* esdTree = dynamic_cast<TTree*> (GetInputData(0));
 725   // Get the list of OCDB objects used for reco
 726   TList *cdbList = (TList*)(esdTree->GetTree()->GetUserInfo())->FindObject("cdbList");
 727   TIter iter2(cdbList);
 728   TObjString* cdbEntry=0;
 729   TString cdbEntryString;
 730   while((cdbEntry =(TObjString*)(iter2.Next()))) {
 731   cdbEntryString = cdbEntry->GetString();
 732   if(cdbEntryString.Contains("ITS/Align/Data"))
 733     itsaligndata = new TObjString(*cdbEntry);
 734   }
 735
 736   Int_t ntracks = esd->GetNumberOfTracks();
 737
 738   if(ntracks==0) return;
 739
 740   if(esd->GetPrimaryVertexTracks()->GetNContributors()<=0) return;
 741
 742   Double_t vtxpos[3]; esd->GetPrimaryVertexTracks()->GetXYZ(vtxpos);
 743
 744   Int_t ncls=0;
 745   Double_t pt=-10000.;
 746   Double_t d0z0[2],covd0z0[3];
 747   const AliTrackPointArray *array = 0;
 748
 749   for (Int_t itrack=0; itrack < ntracks; itrack++) {
 750     AliESDtrack * track = esd->GetTrack(itrack);
 751     if (!track) continue;
 752
 753     if(track->GetNcls(0)<GetRecoParam()->GetAlignFilterMinITSPoints()) continue;
 754
 755     if((GetRecoParam()->GetAlignFilterOnlyITSSATracks()) && track->GetNcls(1)>0) continue;
 756     if((GetRecoParam()->GetAlignFilterOnlyITSTPCTracks()) && track->GetNcls(1)==0) continue;
 757
 758     if(track->Pt()<GetRecoParam()->GetAlignFilterMinPt() ||
 759        track->Pt()>GetRecoParam()->GetAlignFilterMaxPt()) continue;
 760
 761     pt = track->Pt();
 762     ncls = track->GetNcls(0);
 763     Double_t maxd=10000.;
 764     track->PropagateToDCA(esd->GetPrimaryVertex(),esd->GetMagneticField(),maxd,d0z0,covd0z0);
 765
 766     // read ITS cluster map
 767     Int_t map[6];
 768     for(Int_t ilay=0;ilay<6;ilay++) {
 769       map[ilay]=0;
 770       if(track->HasPointOnITSLayer(ilay)) map[ilay]=1;
 771     }
 772     AliDebug(2,Form("ITS map : %d %d %d %d %d %d pt %f\n",map[0],map[1],map[2],map[3],map[4],map[5],track->Pt()));
 773     Int_t idx[12];
 774     track->GetITSclusters(idx);
 775     AliDebug(2,Form("cls idx %d %d %d %d %d %d %d %d %d %d %d %d\n",idx[0],idx[1],idx[2],idx[3],idx[4],idx[5],idx[6],idx[7],idx[8],idx[9],idx[10],idx[11]));
 776
 777
 778     AliTrackPoint point;
 779     Int_t ipt,volId,modId,layerId,lay,lad,det;
 780     Int_t jpt=0;
 781     Bool_t layerOK[6]; for(Int_t l1=0;l1<6;l1++) layerOK[l1]=kFALSE;
 782
 783     array = track->GetTrackPointArray();
 784     if(!array) continue;
 785     for(ipt=0; ipt<array->GetNPoints(); ipt++) {
 786       array->GetPoint(point,ipt);
 787       volId = point.GetVolumeID();
 788       layerId = AliGeomManager::VolUIDToLayer(volId,modId);
 789       if(layerId<1 || layerId>6) continue;
 790       if(point.IsExtra() &&
 791          (GetRecoParam()->GetAlignFilterSkipExtra())) continue;
 792       layerOK[layerId-1]=kTRUE;
 793       jpt++;
 794     }
 795
 796     if(jpt < GetRecoParam()->GetAlignFilterMinITSPoints()) continue;
 797
 798     fHistNpoints->Fill(jpt);
 799     fHistPt->Fill(pt);
 800     PostData(1,fListOfHistos);
 801
 802     Float_t dzOverlap[2];
 803     Double_t globExtra[3],locExtra[3];
 804     arrayForTree = new AliTrackPointArray(jpt);
 805     jpt=0;
 806     array = track->GetTrackPointArray();
 807     if(!array) continue;
 808     for(ipt=0; ipt<array->GetNPoints(); ipt++) {
 809       array->GetPoint(point,ipt);
 810       volId = point.GetVolumeID();
 811       layerId = AliGeomManager::VolUIDToLayer(volId,modId);
 812       if(layerId<1 || layerId>6 || !layerOK[layerId-1]) continue;
 813       if(!point.IsExtra() ||
 814          !(GetRecoParam()->GetAlignFilterFillQANtuples())) continue;
 815       ncls--;
 816       for(Int_t ll=1;ll<layerId;ll++) modId+=AliITSgeomTGeo::GetNLadders(ll)*AliITSgeomTGeo::GetNDetectors(ll);
 817       AliITSgeomTGeo::GetModuleId(modId,lay,lad,det);
 818       globExtra[0]=point.GetX();
 819       globExtra[1]=point.GetY();
 820       globExtra[2]=point.GetZ();
 821       AliITSgeomTGeo::GlobalToLocal(lay,lad,det,globExtra,locExtra);
 822       track->GetDZ(point.GetX(),point.GetY(),point.GetZ(),esd->GetMagneticField(),dzOverlap);
 823       AliTrackPoint pointT;
 824       Float_t radius,radiusT,phiv,phivT,thetav,thetavT;
 825       for(Int_t lll=0;lll<ipt;lll++) {
 826         array->GetPoint(pointT,lll);
 827         Int_t layerIdT = AliGeomManager::VolUIDToLayer(pointT.GetVolumeID(),modId);
 828         if(layerIdT!=layerId) continue;
 829         radius=TMath::Sqrt((point.GetX()-vtxpos[0])*(point.GetX()-vtxpos[0])+(point.GetY()-vtxpos[1])*(point.GetY()-vtxpos[1]));
 830         radiusT=TMath::Sqrt((pointT.GetX()-vtxpos[0])*(pointT.GetX()-vtxpos[0])+(pointT.GetY()-vtxpos[1])*(pointT.GetY()-vtxpos[1]));
 831         phiv=TMath::ATan2(point.GetY()-vtxpos[1],point.GetX()-vtxpos[0]);
 832         phivT=TMath::ATan2(pointT.GetY()-vtxpos[1],pointT.GetX()-vtxpos[0]);
 833         if(TMath::Abs(point.GetZ()-vtxpos[2])<0.00001 || TMath::Abs(pointT.GetZ()-vtxpos[2])<0.00001) continue;
 834         thetav=TMath::ATan(radius/(point.GetZ()-vtxpos[2]));
 835         thetavT=TMath::ATan(radiusT/(pointT.GetZ()-vtxpos[2]));
 836         dzOverlap[0]=(Float_t)((phivT-phiv)*0.5*(radiusT+radius));
 837         if(TMath::Abs(TMath::Tan(0.5*(thetav+thetavT)))<0.00001) continue;
 838         dzOverlap[1]=(Float_t)((pointT.GetZ()-point.GetZ())-(radiusT-radius)/TMath::Tan(0.5*(thetav+thetavT)));
 839         fntExtra->Fill((Float_t)ncls,(Float_t)(layerId-1),lad,volId,TMath::ATan2(point.GetY(),point.GetX()),point.GetX(),point.GetY(),point.GetZ(),locExtra[0],locExtra[2],dzOverlap[0],dzOverlap[1],d0z0[0],d0z0[1]);
 840       }
 841       arrayForTree->AddPoint(jpt,&point);
 842       switch(layerId) {
 843       case 1:
 844         fHistLayer0->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 845         break;
 846       case 2:
 847         fHistLayer1->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 848         break;
 849       case 3:
 850         fHistLayer2->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 851         break;
 852       case 4:
 853         fHistLayer3->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 854         break;
 855       case 5:
 856         fHistLayer4->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 857         break;
 858       case 6:
 859         fHistLayer5->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 860         break;
 861       }
 862       jpt++;
 863     }
 864
 865     curv = track->GetC(esd->GetMagneticField());
 866     curverr = TMath::Sqrt(track->GetSigma1Pt2())*track->GetC(esd->GetMagneticField())/track->OneOverPt();
 867
 868     fspTree->Fill();
 869
 870   } // end of tracks loop
 871
 872   PostData(0,fspTree);
 873
 874   return;
 875 }
 876
 877 //________________________________________________________________________
 878 void AliAlignmentDataFilterITS::Terminate(Option_t */*option*/)
 879 {
 880   // Terminate analysis
 881   //
 882   AliDebug(2,"AliITSAlignmentDataFiler: Terminate() \n");
 883
 884   fspTree = dynamic_cast<TTree*> (GetOutputData(0));
 885   if (!fspTree) {
 886     printf("ERROR: fspTree not available\n");
 887     return;
 888   }
 889
 890   fListOfHistos = dynamic_cast<TList*> (GetOutputData(1));
 891   if (!fListOfHistos) {
 892     printf("ERROR: fListOfHistos not available\n");
 893     return;
 894   }
 895
 896   fHistNevents = dynamic_cast<TH1F*>(fListOfHistos->FindObject("fHistNevents"));
 897   fHistNpoints = dynamic_cast<TH1F*>(fListOfHistos->FindObject("fHistNpoints"));
 898   fHistPt = dynamic_cast<TH1F*>(fListOfHistos->FindObject("fHistPt"));
 899   fHistLayer0 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer0"));
 900   fHistLayer1 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer1"));
 901   fHistLayer2 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer2"));
 902   fHistLayer3 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer3"));
 903   fHistLayer4 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer4"));
 904   fHistLayer5 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer5"));
 905   fntExtra = dynamic_cast<TNtuple*>(fListOfHistos->FindObject("fntExtra"));
 906   fntCosmicMatching = dynamic_cast<TNtuple*>(fListOfHistos->FindObject("fntCosmicMatching"));
 907
 908
 909
 910   return;
 911 }
 912