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 <TFile.h>
  26 #include <TSystem.h>
  27 #include <TChain.h>
  28 #include <TNtuple.h>
  29 #include <TList.h>
  30 #include <TH1F.h>
  31 #include <TH2F.h>
  32 #include <TMap.h>
  33 #include <TVector3.h>
  34 #include <TGeoManager.h>
  35
  36 #include "AliLog.h"
  37 #include "AliCDBEntry.h"
  38 #include "AliGeomManager.h"
  39 #include "AliITSReconstructor.h"
  40 #include "AliITSAlignMille2Module.h"
  41 #include "AliITSgeomTGeo.h"
  42 #include "AliTrackPointArray.h"
  43 #include "AliESDInputHandler.h"
  44 #include "AliESDVertex.h"
  45 #include "AliESDtrack.h"
  46 #include "AliESDEvent.h"
  47 #include "AliESDfriend.h"
  48 #include "AliAnalysisTask.h"
  49 #include "AliAnalysisManager.h"
  50 #include "AliAlignmentDataFilterITS.h"
  51
  52
  53 ClassImp(AliAlignmentDataFilterITS)
  54
  55
  56 //________________________________________________________________________
  57 AliAlignmentDataFilterITS::AliAlignmentDataFilterITS(const char *name):
  58 AliAnalysisTask(name,"task"),
  59 fOnlySPDFO(kFALSE),
  60 fGeometryFileName("geometry.root"),
  61 fITSRecoParam(0),
  62 fESD(0),
  63 fESDfriend(0),
  64 fListOfHistos(0),
  65 fspTree(0),
  66 fHistNevents(0),
  67 fHistNpoints(0),
  68 fHistPt(0),
  69 fHistLayer0(0),
  70 fHistLayer1(0),
  71 fHistLayer2(0),
  72 fHistLayer3(0),
  73 fHistLayer4(0),
  74 fHistLayer5(0),
  75 fntExtra(0),
  76 fntCosmicMatching(0)
  77 {
  78   // Constructor
  79
  80   // Define input and output slots here
  81   // Input slot #0 works with a TChain
  82   DefineInput(0, TChain::Class());
  83   // Output slot #0 writes into a TTree
  84   DefineOutput(0,TTree::Class());  //My private output
  85   // Output slot #1 writes into a TList
  86   DefineOutput(1,TList::Class());  //My private output
  87 }
  88
  89 //________________________________________________________________________
  90 AliAlignmentDataFilterITS::~AliAlignmentDataFilterITS()
  91 {
  92   // Destructor
  93   if (fListOfHistos) {
  94     delete fListOfHistos;
  95     fListOfHistos = 0;
  96   }
  97   if (fspTree) {
  98     delete fspTree;
  99     fspTree = 0;
 100   }
 101   if (fHistNevents) {
 102     delete fHistNevents;
 103     fHistNevents = 0;
 104   }
 105   if (fHistNpoints) {
 106     delete fHistNpoints;
 107     fHistNpoints = 0;
 108   }
 109   if (fHistPt) {
 110     delete fHistPt;
 111     fHistPt = 0;
 112   }
 113   if (fHistLayer0) {
 114     delete fHistLayer0;
 115     fHistLayer0 = 0;
 116   }
 117   if (fHistLayer1) {
 118     delete fHistLayer1;
 119     fHistLayer1 = 0;
 120   }
 121   if (fHistLayer2) {
 122     delete fHistLayer2;
 123     fHistLayer2 = 0;
 124   }
 125   if (fHistLayer3) {
 126     delete fHistLayer3;
 127     fHistLayer3 = 0;
 128   }
 129   if (fHistLayer4) {
 130     delete fHistLayer4;
 131     fHistLayer4 = 0;
 132   }
 133   if (fHistLayer5) {
 134     delete fHistLayer5;
 135     fHistLayer5 = 0;
 136   }
 137   if (fntExtra) {
 138     delete fntExtra;
 139     fntExtra = 0;
 140   }
 141   if (fntCosmicMatching) {
 142     delete fntCosmicMatching;
 143     fntCosmicMatching = 0;
 144   }
 145 }
 146
 147 //________________________________________________________________________
 148 void AliAlignmentDataFilterITS::ConnectInputData(Option_t *)
 149 {
 150   // Connect ESD
 151   // Called once
 152
 153   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
 154   if(!tree) {
 155     printf("ERROR: Could not read chain from input slot 0\n");
 156   } else {
 157     // Get the OCDB path and the list of OCDB objects used for reco
 158     TMap *cdbMap = (TMap*)(tree->GetTree()->GetUserInfo())->FindObject("cdbMap");
 159     TList *cdbList = (TList*)(tree->GetTree()->GetUserInfo())->FindObject("cdbList");
 160
 161     //cdbList->Print();
 162     // write the list to the user info of the output tree
 163     if(!fspTree) {
 164       printf("ERROR: fspTree does not exist\n");
 165     } else {
 166       TMap *cdbMapCopy = new TMap(cdbMap->GetEntries());
 167       cdbMapCopy->SetOwner(1);
 168       cdbMapCopy->SetName("cdbMap");
 169       TIter iter1(cdbMap->GetTable());
 170
 171       TPair* pair = 0;
 172       while((pair = dynamic_cast<TPair*> (iter1.Next()))){
 173         TObjString* keyStr = dynamic_cast<TObjString*> (pair->Key());
 174         TObjString* valStr = dynamic_cast<TObjString*> (pair->Value());
 175         cdbMapCopy->Add(new TObjString(keyStr->GetName()), new TObjString(valStr->GetName()));
 176       }
 177
 178       TList *cdbListCopy = new TList();
 179       cdbListCopy->SetOwner(1);
 180       cdbListCopy->SetName("cdbList");
 181
 182       TIter iter2(cdbList);
 183
 184       TObjString* cdbEntry=0;
 185       while((cdbEntry =(TObjString*)(iter2.Next()))) {
 186         cdbListCopy->Add(new TObjString(*cdbEntry));
 187       }
 188       cdbListCopy->Print();
 189
 190
 191       fspTree->GetUserInfo()->Add(cdbMapCopy);
 192       fspTree->GetUserInfo()->Add(cdbListCopy);
 193     }
 194
 195     // Disable all branches and enable only the needed ones
 196     tree->SetBranchStatus("fTriggerMask", 1);
 197     tree->SetBranchStatus("fSPDVertex*", 1);
 198     tree->SetBranchStatus("ESDfriend*", 1);
 199     tree->SetBranchAddress("ESDfriend.",&fESDfriend);
 200     tree->SetBranchStatus("fSPDMult*", 1);
 201
 202     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 203
 204     if(!esdH) {
 205       printf("ERROR: Could not get ESDInputHandler\n");
 206     } else {
 207       fESD = esdH->GetEvent();
 208     }
 209   }
 210
 211   return;
 212 }
 213
 214 //________________________________________________________________________
 215 void AliAlignmentDataFilterITS::Init()
 216 {
 217   // Initialization
 218
 219   return;
 220 }
 221
 222 //________________________________________________________________________
 223 void AliAlignmentDataFilterITS::CreateOutputObjects()
 224 {
 225   // Create the output container
 226   //
 227
 228   // load the geometry
 229   if(!gGeoManager) {
 230     printf("AliAlignmentDataFilterITS::CreateOutputObjects(): loading geometry from %s\n",fGeometryFileName.Data());
 231     AliGeomManager::LoadGeometry(fGeometryFileName.Data());
 232     if(!gGeoManager) {
 233       printf("AliAlignmentDataFilterITS::CreateOutputObjects(): no geometry loaded \n");
 234       return;
 235     }
 236   }
 237
 238   // Several histograms are more conveniently managed in a TList
 239   fListOfHistos = new TList();
 240   fListOfHistos->SetOwner();
 241
 242   fHistNevents = new TH1F("fHistNevents", "Number of processed events; N events; bin",5,-0.5,4.5);
 243   fHistNevents->Sumw2();
 244   fHistNevents->SetMinimum(0);
 245   fListOfHistos->Add(fHistNevents);
 246
 247   fHistNpoints = new TH1F("fHistNpoints", "Number of AliTrackPoints per track; N points; tracks",25,-0.5,24.5);
 248   fHistNpoints->Sumw2();
 249   fHistNpoints->SetMinimum(0);
 250   fListOfHistos->Add(fHistNpoints);
 251
 252   fHistPt = new TH1F("fHistPt", "p_{t} of tracks; p_{t} [GeV/c]; tracks",100,0,50);
 253   fHistPt->Sumw2();
 254   fHistPt->SetMinimum(0);
 255   fListOfHistos->Add(fHistPt);
 256
 257
 258   Float_t zmax=14.;
 259   Int_t nbinsphi=20,nbinsz=4;
 260   fHistLayer0 = new TH2F("fHistLayer0","Points in layer inner SPD; global   #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax);
 261   fListOfHistos->Add(fHistLayer0);
 262   zmax=14.;
 263   nbinsphi=40;nbinsz=4;
 264   fHistLayer1 = new TH2F("fHistLayer1","Points in layer outer SPD; global   #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax);
 265   fListOfHistos->Add(fHistLayer1);
 266   zmax=22.;
 267   nbinsphi=14;nbinsz=6;
 268   fHistLayer2 = new TH2F("fHistLayer2","Points in layer inner SDD; global   #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax);
 269   fListOfHistos->Add(fHistLayer2);
 270   zmax=29.5;
 271   nbinsphi=22;nbinsz=8;
 272   fHistLayer3 = new TH2F("fHistLayer3","Points in layer outer SDD; global   #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax);
 273   fListOfHistos->Add(fHistLayer3);
 274   zmax=45.;
 275   nbinsphi=34;nbinsz=23;
 276   fHistLayer4 = new TH2F("fHistLayer4","Points in layer inner SSD; global   #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax);
 277   fListOfHistos->Add(fHistLayer4);
 278   zmax=51.;
 279   nbinsphi=38;nbinsz=26;
 280   fHistLayer5 = new TH2F("fHistLayer5","Points in layer outer SSD; global   #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax);
 281   fListOfHistos->Add(fHistLayer5);
 282
 283
 284   fntExtra = new TNtuple("fntExtra","extra clusters in ITS","ncls:layer:ladder:volid:phi:x:y:z:xloc:zloc:dxy:dz:d0mu:z0mu");
 285   fListOfHistos->Add(fntExtra);
 286
 287   fntCosmicMatching = new TNtuple("fntCosmicMatching","cosmic tracks matching in ITS","ncls1:ncls2:pt1:pt2:sigmad01:sigmad02:sigmaz01:sigmaz02:dxy:dz:phimu:thetamu:d0mu:z0mu");
 288   fListOfHistos->Add(fntCosmicMatching);
 289
 290   fspTree = new TTree("spTree","Tree with ITS track points");
 291   AliTrackPointArray *array = 0;
 292   Float_t curv=0,curverr=0,runNumber=0;
 293   TObjString *itsaligndata = 0;
 294   TObjString *itscalibrespsdd = 0;
 295   fspTree->Branch("SP","AliTrackPointArray",&array);
 296   fspTree->Branch("curv",&curv);
 297   fspTree->Branch("curverr",&curverr);
 298   fspTree->Branch("run",&runNumber);
 299   fspTree->Branch("ITSAlignData",&itsaligndata);
 300   fspTree->Branch("ITSCalibRespSDD",&itscalibrespsdd);
 301
 302   return;
 303 }
 304
 305 //________________________________________________________________________
 306 void AliAlignmentDataFilterITS::Exec(Option_t */*option*/)
 307 {
 308   // Execute analysis for current event:
 309   // write ITS AliTrackPoints for selected tracks to fspTree
 310
 311   // check the geometry
 312   if(!gGeoManager) {
 313     printf("AliAlignmentDataFilterITS::Exec(): no geometry loaded \n");
 314     return;
 315   }
 316
 317   // check if we have AliITSRecoParam
 318   if(!GetRecoParam()) {
 319     if(!fITSRecoParam) {
 320       printf("AliAlignmentDataFilterITS::Exec(): no AliITSRecoParam\n");
 321       return;
 322     }
 323   }
 324
 325
 326   if(!fESD) {
 327     printf("AliAlignmentDataFilterITS::Exec(): no ESD \n");
 328     return;
 329   }
 330   if(!fESDfriend) {
 331     printf("AliAlignmentDataFilterITS::Exec(): no ESDfriend \n");
 332     return;
 333   }
 334   // attach ESDfriend
 335   fESD->SetESDfriend(fESDfriend);
 336
 337   // Post the data for slot 0
 338   fHistNevents->Fill(0);
 339
 340
 341   // write field value to spTree UserInfo
 342   if(!((fspTree->GetUserInfo())->FindObject("BzkGauss"))) {
 343     Double_t bz=fESD->GetMagneticField();
 344     TString bzString; bzString+=bz;
 345     TObjString *bzObjString = new TObjString(bzString);
 346     TList *bzList = new TList();
 347     bzList->SetOwner(1);
 348     bzList->SetName("BzkGauss");
 349     bzList->Add(bzObjString);
 350     fspTree->GetUserInfo()->Add(bzList);
 351   }
 352
 353
 354   // Process event as Cosmic or Collision
 355   //if(esd->GetEventType()== ???? ) {
 356   printf("AliAlignmentDataFilterITS::Exec(): MOVE ASAP TO esd->GetEventType() !\n");
 357   if(GetRecoParam()->GetAlignFilterCosmics()) {
 358     FilterCosmic(fESD);
 359   } else {
 360     FilterCollision(fESD);
 361   }
 362
 363   PostData(1,fListOfHistos);
 364
 365   return;
 366 }
 367
 368 //________________________________________________________________________
 369 void AliAlignmentDataFilterITS::FilterCosmic(const AliESDEvent *esd)
 370 {
 371   // Extract ITS AliTrackPoints for Cosmics (check angular matching
 372   // of top and bottom track, merge the two tracks, if requested)
 373   //
 374
 375   // Set branch addresses for space points tree
 376   AliTrackPointArray *arrayForTree=0;
 377   Float_t curv,curverr,runNumber;
 378   TObjString *itsaligndata=0;
 379   TObjString *itscalibrespsdd = 0;
 380   fspTree->SetBranchAddress("SP",&arrayForTree);
 381   fspTree->SetBranchAddress("curv",&curv);
 382   fspTree->SetBranchAddress("curverr",&curverr);
 383   fspTree->SetBranchAddress("run",&runNumber);
 384   fspTree->SetBranchAddress("ITSAlignData",&itsaligndata);
 385   fspTree->SetBranchAddress("ITSCalibRespSDD",&itscalibrespsdd);
 386
 387
 388   runNumber = (Float_t)esd->GetRunNumber();
 389
 390   TTree* esdTree = dynamic_cast<TTree*> (GetInputData(0));
 391   // Get the list of OCDB objects used for reco
 392   TList *cdbList = (TList*)(esdTree->GetTree()->GetUserInfo())->FindObject("cdbList");
 393   TIter iter2(cdbList);
 394   TObjString* cdbEntry=0;
 395   TString cdbEntryString;
 396   while((cdbEntry =(TObjString*)(iter2.Next()))) {
 397   cdbEntryString = cdbEntry->GetString();
 398   if(cdbEntryString.Contains("ITS/Align/Data"))
 399     itsaligndata = new TObjString(*cdbEntry);
 400   if(cdbEntryString.Contains("ITS/Calib/RespSDD"))
 401     itscalibrespsdd = new TObjString(*cdbEntry);
 402   }
 403
 404
 405   TString triggeredClass = fESD->GetFiredTriggerClasses();
 406   if(fOnlySPDFO && !triggeredClass.Contains("C0SCO-ABCE-NOPF-CENT")) return;
 407
 408
 409   Int_t ntracks = esd->GetNumberOfTracks();
 410   if(ntracks<2) return;
 411
 412   if(esd->GetPrimaryVertexSPD()->GetNContributors()<0) return;
 413
 414   Double_t vtxpos[3]; esd->GetPrimaryVertexSPD()->GetXYZ(vtxpos);
 415
 416   Int_t *goodtracksArray = new Int_t[ntracks];
 417   Float_t *phiArray = new Float_t[ntracks];
 418   Float_t *thetaArray = new Float_t[ntracks];
 419   Int_t *nclsArray = new Int_t[ntracks];
 420   Int_t ngt=0;
 421   Int_t itrack=0;
 422   for (itrack=0; itrack < ntracks; itrack++) {
 423     AliESDtrack *track = esd->GetTrack(itrack);
 424     if (!track) continue;
 425
 426
 427     if(track->GetNcls(0)<GetRecoParam()->GetAlignFilterMinITSPoints()) continue;
 428
 429     if((GetRecoParam()->GetAlignFilterOnlyITSSATracks()) && track->GetNcls(1)>0) continue;
 430     if((GetRecoParam()->GetAlignFilterOnlyITSTPCTracks()) && track->GetNcls(1)==0) continue;
 431
 432     Float_t phi = track->GetAlpha()+TMath::ASin(track->GetSnp());
 433     Float_t theta = 0.5*TMath::Pi()-TMath::ATan(track->GetTgl());
 434
 435     if(track->Pt()<GetRecoParam()->GetAlignFilterMinPt() ||
 436        track->Pt()>GetRecoParam()->GetAlignFilterMaxPt()) continue;
 437
 438     goodtracksArray[ngt] = itrack;
 439     phiArray[ngt]        = phi;
 440     thetaArray[ngt]      = theta;
 441     nclsArray[ngt]       = track->GetNcls(0);
 442     ngt++;
 443   }
 444
 445   if(ngt<2) {
 446     delete [] goodtracksArray; goodtracksArray=0;
 447     delete [] phiArray; phiArray=0;
 448     delete [] thetaArray; thetaArray=0;
 449     delete [] nclsArray; nclsArray=0;
 450     return;
 451   }
 452
 453   // check matching of the two tracks from the muon
 454   Float_t min = 10000000.;
 455   Int_t maxCls = 0;
 456   Int_t good1 = -1, good2 = -1;
 457   for(Int_t itr1=0; itr1<ngt-1; itr1++) {
 458     for(Int_t itr2=itr1+1; itr2<ngt; itr2++) {
 459       Float_t deltatheta = TMath::Abs(TMath::Pi()-thetaArray[itr1]-thetaArray[itr2]);
 460       if(deltatheta>GetRecoParam()->GetAlignFilterMaxMatchingAngle()) continue;
 461       Float_t deltaphi = TMath::Abs(TMath::Abs(phiArray[itr1]-phiArray[itr2])-TMath::Pi());
 462       if(deltaphi>GetRecoParam()->GetAlignFilterMaxMatchingAngle()) continue;
 463       if(nclsArray[itr1]+nclsArray[itr2] > maxCls) {
 464         maxCls = nclsArray[itr1]+nclsArray[itr2];
 465         min = deltatheta+deltaphi;
 466         good1 = goodtracksArray[itr1];
 467         good2 = goodtracksArray[itr2];
 468       } else if(nclsArray[itr1]+nclsArray[itr2] == maxCls) {
 469         if(deltatheta+deltaphi < min) {
 470           min = deltatheta+deltaphi;
 471           good1 = goodtracksArray[itr1];
 472           good2 = goodtracksArray[itr2];
 473         }
 474       }
 475     }
 476   }
 477
 478   delete [] goodtracksArray; goodtracksArray=0;
 479   delete [] phiArray; phiArray=0;
 480   delete [] thetaArray; thetaArray=0;
 481   delete [] nclsArray; nclsArray=0;
 482
 483   if(good1<0) return;
 484   AliDebug(2,"ok track matching");
 485
 486   // track1 will be the inward track (top)
 487   // track2 the outward (bottom)
 488   AliESDtrack *track1=0;
 489   AliESDtrack *track2=0;
 490   AliESDtrack *track = esd->GetTrack(good1);
 491   if(track->Py()>0) {
 492     track1 = esd->GetTrack(good1);
 493     track2 = esd->GetTrack(good2);
 494   } else {
 495     track1 = esd->GetTrack(good2);
 496     track2 = esd->GetTrack(good1);
 497   }
 498
 499   AliTrackPoint point;
 500   const AliTrackPointArray *array=0;
 501   Int_t ipt,volId,modId,layerId,lay,lad,det;
 502   Int_t jpt=0;
 503   Bool_t layerOK[6][2];
 504   Int_t nclsTrk[2]={0,0};
 505
 506   for(Int_t l1=0;l1<6;l1++) for(Int_t l2=0;l2<2;l2++) layerOK[l1][l2]=kFALSE;
 507
 508   for(itrack=0; itrack<2; itrack++) {
 509     if(itrack==0) {
 510       track = track1;
 511     } else {
 512       track = track2;
 513     }
 514     array = track->GetTrackPointArray();
 515     if(!array) {
 516       AliWarning("No tracks points avaialble");
 517       continue;
 518     }
 519     for(ipt=0; ipt<array->GetNPoints(); ipt++) {
 520       array->GetPoint(point,ipt);
 521       volId = point.GetVolumeID();
 522       if(volId<=0) continue;
 523       layerId = AliGeomManager::VolUIDToLayer(volId,modId);
 524       AliDebug(2,Form("%d %d %d  %f\n",ipt,layerId-1,volId,TMath::Sqrt(point.GetX()*point.GetX()+point.GetY()*point.GetY())));
 525       if(point.IsExtra() &&
 526          (GetRecoParam()->GetAlignFilterSkipExtra())) continue;
 527       if(layerId<1 || layerId>6) continue;
 528       if(!GetRecoParam()->GetAlignFilterUseLayer(layerId-1)) continue;
 529       // check minAngleWrtITSModulePlanes
 530       Double_t p[3]; track->GetDirection(p);
 531       TVector3 pvec(p[0],p[1],p[2]);
 532       Double_t rot[9]; AliGeomManager::GetOrigRotation(volId,rot);
 533       TVector3 normvec(rot[1],rot[4],rot[7]);
 534       Double_t angle = pvec.Angle(normvec);
 535       if(angle>0.5*TMath::Pi()) angle = TMath::Pi()-angle;
 536       angle = 0.5*TMath::Pi()-angle;
 537       if(angle<GetRecoParam()->GetAlignFilterMinAngleWrtModulePlanes()) continue;
 538       layerOK[layerId-1][itrack]=kTRUE;
 539       jpt++;
 540       nclsTrk[itrack]++;
 541     }
 542   }
 543   AliDebug(2,Form("nClsTrk1 %d nClsTrk2 %d\n",nclsTrk[0],nclsTrk[1]));
 544
 545   // read ITS cluster maps
 546   Int_t map1[6],map2[6];
 547   for(Int_t ilay=0;ilay<6;ilay++) {
 548     map1[ilay]=0; map2[ilay]=0;
 549     if(track1->HasPointOnITSLayer(ilay)) map1[ilay]=1;
 550     if(track2->HasPointOnITSLayer(ilay)) map2[ilay]=1;
 551   }
 552   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()));
 553   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()));
 554   Int_t idx1[12],idx2[12];
 555   track1->GetITSclusters(idx1);
 556   track2->GetITSclusters(idx2);
 557   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]));
 558   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]));
 559
 560
 561   if(jpt<GetRecoParam()->GetAlignFilterMinITSPointsMerged()) return;
 562   AliDebug(2,Form(" Total points %d, accepted\n",jpt));
 563   fHistNpoints->Fill(jpt);
 564   fHistPt->Fill(0.5*(track1->Pt()+track2->Pt()));
 565
 566   Float_t d0z0mu[2];
 567   track1->GetDZ(0,0,0,esd->GetMagneticField(),d0z0mu);
 568   //printf("d0mu %f  z0mu %f\n",d0z0mu[0],d0z0mu[1]);
 569
 570   Float_t dzOverlap[2];
 571   Float_t curvArray[2],curverrArray[2];
 572   Double_t globExtra[3],locExtra[3];
 573   if(GetRecoParam()->GetAlignFilterCosmicMergeTracks())
 574     arrayForTree = new AliTrackPointArray(jpt);
 575
 576   jpt=0;
 577   for(itrack=0; itrack<2; itrack++) {
 578     if(itrack==0) {
 579       track = track1;
 580     } else {
 581       track = track2;
 582     }
 583     curvArray[itrack] = track->GetC(esd->GetMagneticField());
 584     curverrArray[itrack] = TMath::Sqrt(track->GetSigma1Pt2())*track->GetC(esd->GetMagneticField())/track->OneOverPt();
 585
 586     if(!(GetRecoParam()->GetAlignFilterCosmicMergeTracks())) {
 587       jpt=0;
 588       arrayForTree = new AliTrackPointArray(nclsTrk[itrack]);
 589     }
 590     array = track->GetTrackPointArray();
 591     for(ipt=0; ipt<array->GetNPoints(); ipt++) {
 592       array->GetPoint(point,ipt);
 593       volId = point.GetVolumeID();
 594       if(volId<=0) continue;
 595       layerId = AliGeomManager::VolUIDToLayer(volId,modId);
 596       if(layerId<1 || layerId>6 || !layerOK[layerId-1][itrack]) continue;
 597       arrayForTree->AddPoint(jpt,&point);
 598       jpt++;
 599       switch(layerId) {
 600       case 1:
 601         fHistLayer0->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 602         break;
 603       case 2:
 604         fHistLayer1->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 605         break;
 606       case 3:
 607         fHistLayer2->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 608         break;
 609       case 4:
 610         fHistLayer3->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 611         break;
 612       case 5:
 613         fHistLayer4->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 614         break;
 615       case 6:
 616         fHistLayer5->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 617         break;
 618       }
 619       // Post the data for slot 0
 620       if(jpt==1) PostData(1,fListOfHistos); // only if this is the first points
 621       if(!point.IsExtra() ||
 622          !(GetRecoParam()->GetAlignFilterFillQANtuples())) continue;
 623       nclsTrk[itrack]--;
 624       for(Int_t ll=1;ll<layerId;ll++) modId+=AliITSgeomTGeo::GetNLadders(ll)*AliITSgeomTGeo::GetNDetectors(ll);
 625       AliITSgeomTGeo::GetModuleId(modId,lay,lad,det);
 626       globExtra[0]=point.GetX();
 627       globExtra[1]=point.GetY();
 628       globExtra[2]=point.GetZ();
 629       AliITSgeomTGeo::GlobalToLocal(lay,lad,det,globExtra,locExtra);
 630       //printf("%d %d %d %d %d  %f %f %f\n",volId,modId,lay,lad,det,locExtra[0],locExtra[1],locExtra[2]);
 631       track->GetDZ(point.GetX(),point.GetY(),point.GetZ(),esd->GetMagneticField(),dzOverlap);
 632       AliTrackPoint pointT;
 633       Float_t radius,radiusT,phiv,phivT,thetav,thetavT;
 634       for(Int_t lll=0;lll<ipt;lll++) {
 635         array->GetPoint(pointT,lll);
 636         if(pointT.GetVolumeID()<=0) continue;
 637         Int_t layerIdT = AliGeomManager::VolUIDToLayer(pointT.GetVolumeID(),modId);
 638         if(layerIdT!=layerId) continue;
 639         radius=TMath::Sqrt((point.GetX()-vtxpos[0])*(point.GetX()-vtxpos[0])+(point.GetY()-vtxpos[1])*(point.GetY()-vtxpos[1]));
 640         radiusT=TMath::Sqrt((pointT.GetX()-vtxpos[0])*(pointT.GetX()-vtxpos[0])+(pointT.GetY()-vtxpos[1])*(pointT.GetY()-vtxpos[1]));
 641         phiv=TMath::ATan2(point.GetY()-vtxpos[1],point.GetX()-vtxpos[0]);
 642         phivT=TMath::ATan2(pointT.GetY()-vtxpos[1],pointT.GetX()-vtxpos[0]);
 643         if(TMath::Abs(point.GetZ()-vtxpos[2])<0.00001 || TMath::Abs(pointT.GetZ()-vtxpos[2])<0.00001) continue;
 644         thetav=TMath::ATan(radius/(point.GetZ()-vtxpos[2]));
 645         thetavT=TMath::ATan(radiusT/(pointT.GetZ()-vtxpos[2]));
 646         dzOverlap[0]=(Float_t)((phivT-phiv)*0.5*(radiusT+radius));
 647         if(TMath::Abs(TMath::Tan(0.5*(thetav+thetavT)))<0.00001) continue;
 648         dzOverlap[1]=(Float_t)((pointT.GetZ()-point.GetZ())-(radiusT-radius)/TMath::Tan(0.5*(thetav+thetavT)));
 649         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]);
 650       }
 651     }
 652
 653     if(!(GetRecoParam()->GetAlignFilterCosmicMergeTracks())) {
 654       curv = curvArray[itrack];
 655       curverr = curverrArray[itrack];
 656       fspTree->Fill();
 657     }
 658   }
 659
 660   if(GetRecoParam()->GetAlignFilterCosmicMergeTracks()) {
 661     curv = 0.5*(curvArray[0]-curvArray[1]);  // the "-" is because the two tracks have opposite curvature!
 662     curverr = 0.5*TMath::Sqrt(curverrArray[0]*curverrArray[0]+curverrArray[1]*curverrArray[1]);
 663     fspTree->Fill();
 664   }
 665   PostData(0,fspTree);
 666
 667   if(!(GetRecoParam()->GetAlignFilterFillQANtuples())) return;
 668   // fill ntuple with track-to-track matching
 669   Float_t phimu,thetamu,phiout,thetaout,dphi,dtheta,rotymu,rotyout,droty;
 670   Float_t d0[2],z0[2];
 671   Float_t sigmad0[2],sigmaz0[2];
 672   phimu = track1->GetAlpha()+TMath::ASin(track1->GetSnp());
 673   thetamu = 0.5*TMath::Pi()-TMath::ATan(track1->GetTgl());
 674   phiout = track2->GetAlpha()+TMath::ASin(track2->GetSnp());
 675   thetaout = 0.5*TMath::Pi()-TMath::ATan(track2->GetTgl());
 676   rotymu = TMath::ATan2(track1->Px(),track1->Pz());
 677   rotyout = TMath::ATan2(track2->Px(),track2->Pz());
 678
 679   dphi = phimu - (phiout+TMath::Pi());
 680   dtheta = thetamu - (TMath::Pi()-thetaout);
 681   if(rotymu>0) {
 682     droty = rotymu - (rotyout+TMath::Pi());
 683   } else {
 684     droty = rotymu - (rotyout-TMath::Pi());
 685   }
 686
 687   Double_t alpha = TMath::ATan2(track1->Py(),track1->Px());
 688
 689   track1->Propagate(alpha,0.,esd->GetMagneticField());
 690   track2->Propagate(alpha,0.,esd->GetMagneticField());
 691   d0[0] = track1->GetY();
 692   z0[0] = track1->GetZ();
 693   d0[1] = track2->GetY();
 694   z0[1] = track2->GetZ();
 695   Float_t dxy = -(d0[0]-d0[1]);
 696   Float_t dz  = z0[0]-z0[1];
 697   sigmad0[0] = TMath::Sqrt(track1->GetSigmaY2());
 698   sigmaz0[0] = TMath::Sqrt(track1->GetSigmaZ2());
 699   sigmad0[1] = TMath::Sqrt(track2->GetSigmaY2());
 700   sigmaz0[1] = TMath::Sqrt(track2->GetSigmaZ2());
 701   /*
 702   Double_t xyz1atxl0[3],xyz1atxl1[3],xyz2atxl0[3],xyz2atxl1[3];
 703   track1->GetXYZAt(0.,esd->GetMagneticField(),xyz1atxl0);
 704   track1->GetXYZAt(1.,esd->GetMagneticField(),xyz1atxl1);
 705   track2->GetXYZAt(0.,esd->GetMagneticField(),xyz2atxl0);
 706   track2->GetXYZAt(1.,esd->GetMagneticField(),xyz2atxl1);
 707   Float_t x1aty0 = (xyz1atxl0[0]*xyz1atxl1[1]-xyz1atxl0[1]*xyz1atxl1[0])/(xyz1atxl1[1]-xyz1atxl0[1]);
 708   Float_t x2aty0 = (xyz2atxl0[0]*xyz2atxl1[1]-xyz2atxl0[1]*xyz2atxl1[0])/(xyz2atxl1[1]-xyz2atxl0[1]);
 709   Float_t dxaty0 = x1aty0-x2aty0;
 710   */
 711   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]);
 712
 713   return;
 714 }
 715
 716 //________________________________________________________________________
 717 void AliAlignmentDataFilterITS::FilterCollision(const AliESDEvent *esd)
 718 {
 719   // Extract ITS AliTrackPoints for Cosmics (check angular matching
 720   // of top and bottom track, merge the two tracks, if requested)
 721   //
 722
 723   // Set branch addresses for space points tree
 724   AliTrackPointArray *arrayForTree=0;
 725   Float_t curv,curverr,runNumber;
 726   TObjString *itsaligndata=0;
 727   TObjString *itscalibrespsdd = 0;
 728   fspTree->SetBranchAddress("SP",&arrayForTree);
 729   fspTree->SetBranchAddress("curv",&curv);
 730   fspTree->SetBranchAddress("curverr",&curverr);
 731   fspTree->SetBranchAddress("run",&runNumber);
 732   fspTree->SetBranchAddress("ITSAlignData",&itsaligndata);
 733   fspTree->SetBranchAddress("ITSCalibRespSDD",&itscalibrespsdd);
 734
 735
 736   runNumber = (Float_t)esd->GetRunNumber();
 737
 738   TTree* esdTree = dynamic_cast<TTree*> (GetInputData(0));
 739   // Get the list of OCDB objects used for reco
 740   TList *cdbList = (TList*)(esdTree->GetTree()->GetUserInfo())->FindObject("cdbList");
 741   TIter iter2(cdbList);
 742   TObjString* cdbEntry=0;
 743   TString cdbEntryString;
 744   while((cdbEntry =(TObjString*)(iter2.Next()))) {
 745   cdbEntryString = cdbEntry->GetString();
 746   if(cdbEntryString.Contains("ITS/Align/Data"))
 747     itsaligndata = new TObjString(*cdbEntry);
 748   if(cdbEntryString.Contains("ITS/Calib/RespSDD"))
 749     itscalibrespsdd = new TObjString(*cdbEntry);
 750   }
 751
 752   Int_t ntracks = esd->GetNumberOfTracks();
 753
 754   if(ntracks==0) return;
 755
 756   if(esd->GetPrimaryVertexSPD()->GetNContributors()<=0) return;
 757   TString vtitle = esd->GetPrimaryVertexSPD()->GetTitle();
 758   if(!vtitle.Contains("3D")) return;
 759
 760   Double_t vtxpos[3]; esd->GetPrimaryVertexSPD()->GetXYZ(vtxpos);
 761
 762   Int_t ncls=0;
 763   Double_t pt=-10000.;
 764   Double_t d0z0[2],covd0z0[3];
 765   const AliTrackPointArray *array = 0;
 766
 767   for (Int_t itrack=0; itrack < ntracks; itrack++) {
 768     AliESDtrack * track = esd->GetTrack(itrack);
 769     if (!track) continue;
 770
 771     if(track->GetNcls(0)<GetRecoParam()->GetAlignFilterMinITSPoints()) continue;
 772
 773     if((GetRecoParam()->GetAlignFilterOnlyITSSATracks()) && track->GetNcls(1)>0) continue;
 774     if((GetRecoParam()->GetAlignFilterOnlyITSTPCTracks()) && track->GetNcls(1)==0) continue;
 775
 776     if(track->Pt()<GetRecoParam()->GetAlignFilterMinPt() ||
 777        track->Pt()>GetRecoParam()->GetAlignFilterMaxPt()) continue;
 778
 779     pt = track->Pt();
 780     ncls = track->GetNcls(0);
 781     Double_t maxd=10000.;
 782     track->PropagateToDCA(esd->GetPrimaryVertex(),esd->GetMagneticField(),maxd,d0z0,covd0z0);
 783
 784     // read ITS cluster map
 785     Int_t map[6];
 786     for(Int_t ilay=0;ilay<6;ilay++) {
 787       map[ilay]=0;
 788       if(track->HasPointOnITSLayer(ilay)) map[ilay]=1;
 789     }
 790     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()));
 791     Int_t idx[12];
 792     track->GetITSclusters(idx);
 793     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]));
 794
 795
 796     AliTrackPoint point;
 797     Int_t ipt,volId,modId,layerId,lay,lad,det;
 798     Int_t jpt=0;
 799     Bool_t layerOK[6]; for(Int_t l1=0;l1<6;l1++) layerOK[l1]=kFALSE;
 800
 801     array = track->GetTrackPointArray();
 802     if(!array) continue;
 803     for(ipt=0; ipt<array->GetNPoints(); ipt++) {
 804       array->GetPoint(point,ipt);
 805       volId = point.GetVolumeID();
 806       if(volId<=0) continue;
 807       layerId = AliGeomManager::VolUIDToLayer(volId,modId);
 808       if(layerId<1 || layerId>6) continue;
 809       if(point.IsExtra() &&
 810          (GetRecoParam()->GetAlignFilterSkipExtra())) continue;
 811       layerOK[layerId-1]=kTRUE;
 812       jpt++;
 813     }
 814
 815     if(jpt < GetRecoParam()->GetAlignFilterMinITSPoints()) continue;
 816
 817     fHistNpoints->Fill(jpt);
 818     fHistPt->Fill(pt);
 819     PostData(1,fListOfHistos);
 820
 821     Float_t dzOverlap[2];
 822     Double_t globExtra[3],locExtra[3];
 823     arrayForTree = new AliTrackPointArray(jpt);
 824     jpt=0;
 825     array = track->GetTrackPointArray();
 826     if(!array) continue;
 827     for(ipt=0; ipt<array->GetNPoints(); ipt++) {
 828       array->GetPoint(point,ipt);
 829       volId = point.GetVolumeID();
 830       layerId = AliGeomManager::VolUIDToLayer(volId,modId);
 831       if(layerId<1 || layerId>6 || !layerOK[layerId-1]) continue;
 832       arrayForTree->AddPoint(jpt,&point);
 833       switch(layerId) {
 834       case 1:
 835         fHistLayer0->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 836         break;
 837       case 2:
 838         fHistLayer1->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 839         break;
 840       case 3:
 841         fHistLayer2->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 842         break;
 843       case 4:
 844         fHistLayer3->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 845         break;
 846       case 5:
 847         fHistLayer4->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 848         break;
 849       case 6:
 850         fHistLayer5->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 851         break;
 852       }
 853       jpt++;
 854       if(!point.IsExtra() ||
 855          !(GetRecoParam()->GetAlignFilterFillQANtuples())) continue;
 856       ncls--;
 857       for(Int_t ll=1;ll<layerId;ll++) modId+=AliITSgeomTGeo::GetNLadders(ll)*AliITSgeomTGeo::GetNDetectors(ll);
 858       AliITSgeomTGeo::GetModuleId(modId,lay,lad,det);
 859       globExtra[0]=point.GetX();
 860       globExtra[1]=point.GetY();
 861       globExtra[2]=point.GetZ();
 862       AliITSgeomTGeo::GlobalToLocal(lay,lad,det,globExtra,locExtra);
 863       track->GetDZ(point.GetX(),point.GetY(),point.GetZ(),esd->GetMagneticField(),dzOverlap);
 864       AliTrackPoint pointT;
 865       Float_t radius,radiusT,phiv,phivT,thetav,thetavT;
 866       for(Int_t lll=0;lll<ipt;lll++) {
 867         array->GetPoint(pointT,lll);
 868         Int_t layerIdT = AliGeomManager::VolUIDToLayer(pointT.GetVolumeID(),modId);
 869         if(layerIdT!=layerId) continue;
 870         radius=TMath::Sqrt((point.GetX()-vtxpos[0])*(point.GetX()-vtxpos[0])+(point.GetY()-vtxpos[1])*(point.GetY()-vtxpos[1]));
 871         radiusT=TMath::Sqrt((pointT.GetX()-vtxpos[0])*(pointT.GetX()-vtxpos[0])+(pointT.GetY()-vtxpos[1])*(pointT.GetY()-vtxpos[1]));
 872         phiv=TMath::ATan2(point.GetY()-vtxpos[1],point.GetX()-vtxpos[0]);
 873         phivT=TMath::ATan2(pointT.GetY()-vtxpos[1],pointT.GetX()-vtxpos[0]);
 874         if(TMath::Abs(point.GetZ()-vtxpos[2])<0.00001 || TMath::Abs(pointT.GetZ()-vtxpos[2])<0.00001) continue;
 875         thetav=TMath::ATan(radius/(point.GetZ()-vtxpos[2]));
 876         thetavT=TMath::ATan(radiusT/(pointT.GetZ()-vtxpos[2]));
 877         dzOverlap[0]=(Float_t)((phivT-phiv)*0.5*(radiusT+radius));
 878         if(TMath::Abs(TMath::Tan(0.5*(thetav+thetavT)))<0.00001) continue;
 879         dzOverlap[1]=(Float_t)((pointT.GetZ()-point.GetZ())-(radiusT-radius)/TMath::Tan(0.5*(thetav+thetavT)));
 880         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]);
 881       }
 882     }
 883
 884     curv = track->GetC(esd->GetMagneticField());
 885     curverr = TMath::Sqrt(track->GetSigma1Pt2())*track->GetC(esd->GetMagneticField())/track->OneOverPt();
 886
 887     fspTree->Fill();
 888
 889   } // end of tracks loop
 890
 891   PostData(0,fspTree);
 892
 893   return;
 894 }
 895
 896 //________________________________________________________________________
 897 void AliAlignmentDataFilterITS::Terminate(Option_t */*option*/)
 898 {
 899   // Terminate analysis
 900   //
 901   AliDebug(2,"AliITSAlignmentDataFiler: Terminate() \n");
 902
 903   fspTree = dynamic_cast<TTree*> (GetOutputData(0));
 904   if (!fspTree) {
 905     printf("ERROR: fspTree not available\n");
 906     return;
 907   }
 908
 909   fListOfHistos = dynamic_cast<TList*> (GetOutputData(1));
 910   if (!fListOfHistos) {
 911     printf("ERROR: fListOfHistos not available\n");
 912     return;
 913   }
 914
 915   fHistNevents = dynamic_cast<TH1F*>(fListOfHistos->FindObject("fHistNevents"));
 916   fHistNpoints = dynamic_cast<TH1F*>(fListOfHistos->FindObject("fHistNpoints"));
 917   fHistPt = dynamic_cast<TH1F*>(fListOfHistos->FindObject("fHistPt"));
 918   fHistLayer0 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer0"));
 919   fHistLayer1 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer1"));
 920   fHistLayer2 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer2"));
 921   fHistLayer3 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer3"));
 922   fHistLayer4 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer4"));
 923   fHistLayer5 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer5"));
 924   fntExtra = dynamic_cast<TNtuple*>(fListOfHistos->FindObject("fntExtra"));
 925   fntCosmicMatching = dynamic_cast<TNtuple*>(fListOfHistos->FindObject("fntCosmicMatching"));
 926
 927
 928
 929   return;
 930 }
 931 //-------------------------------------------------------------------------------
 932 const AliITSRecoParam *AliAlignmentDataFilterITS::GetRecoParam() const
 933 {
 934   //
 935   // Return the ITSRecoParam object
 936   //
 937   if(AliITSReconstructor::GetRecoParam()) {
 938     return AliITSReconstructor::GetRecoParam();
 939   } else if(fITSRecoParam) {
 940     return fITSRecoParam;
 941   } else return NULL;
 942 }
 943 //--------------------------------------------------------------------------------
 944 Int_t AliAlignmentDataFilterITS::WriteTrackPointsInIdealGeom(Char_t *fin,
 945                                                              Char_t *fout,
 946                                                              Char_t *fmis,
 947                                                              Char_t *fgeo,
 948                                                              Bool_t prn)
 949 {
 950   //
 951   // Convert AliTrackPoints in fin, reconstructed with fmis, back
 952   // to ideal geometry
 953   //
 954   // M. Lunardon
 955   //
 956
 957
 958   TGeoHMatrix deltahm;
 959
 960   // Load geometry
 961   if (gSystem->AccessPathName(fgeo)) {
 962     printf("couldn't find geometry file %s - skipping...\n",fmis);
 963     return -1;
 964   }
 965
 966   TFile *geofile=TFile::Open(fgeo);
 967   TGeoManager *fgGeometry=NULL;
 968
 969   fgGeometry=(TGeoManager*)geofile->Get("ALICE");
 970
 971   if (!fgGeometry)
 972     fgGeometry=(TGeoManager*)geofile->Get("Geometry");
 973
 974   if (!fgGeometry) {
 975     AliCDBEntry *entry = (AliCDBEntry*)geofile->Get("AliCDBEntry");
 976     if (entry)
 977       fgGeometry = (TGeoManager*)entry->GetObject();
 978   }
 979
 980   if (!fgGeometry) return -1;
 981   AliGeomManager::SetGeometry(fgGeometry);
 982   if(!AliGeomManager::GetGeometry()) return -1;
 983
 984
 985   // open alignment file
 986   if (gSystem->AccessPathName(fmis)) {
 987     printf("couldn't open alignment file %s - skipping...\n",fmis);
 988     return -2;
 989   }
 990   TFile *pref = TFile::Open(fmis);
 991   if (!pref->IsOpen()) return -2;
 992
 993
 994   /// apply alignment to ideal geometry
 995   TClonesArray *prea=(TClonesArray*)pref->Get("ITSAlignObjs");
 996   if (!prea) {
 997     if (pref->Get("AliCDBEntry"))
 998       prea = (TClonesArray*) ((AliCDBEntry*)pref->Get("AliCDBEntry"))->GetObject();
 999   }
1000   if (!prea) return -3;
1001   Int_t nprea=prea->GetEntriesFast();
1002   printf("Array of input misalignments with %d entries\n",nprea);
1003   AliGeomManager::ApplyAlignObjsToGeom(*prea); // apply all levels of objs
1004
1005   AliTrackPointArray *tpain=NULL;
1006   TFile *tpainfile=NULL;
1007   TTree *treein=NULL;
1008   AliTrackPoint point;
1009   AliITSAlignMille2Module *m2[2200];
1010   for (Int_t i=0; i<2198; i++)
1011     m2[i]=new AliITSAlignMille2Module(AliITSAlignMille2Module::GetVolumeIDFromIndex(i));
1012
1013   // open input file
1014   if (gSystem->AccessPathName(fin)) {
1015     printf("couldn't open file %s - skipping...\n",fin);
1016     return -4;
1017   }
1018   tpainfile = TFile::Open(fin);
1019   if (!tpainfile->IsOpen()) return -4;
1020
1021   treein=(TTree*)tpainfile->Get("spTree");
1022   if (!treein) return -5;
1023   Float_t curv,curverr,runNumber;
1024   TObjString *itsaligndata=0;
1025   TObjString *itscalibrespsdd = 0;
1026   treein->SetBranchAddress("SP", &tpain);
1027   treein->SetBranchAddress("curv", &curv);
1028   treein->SetBranchAddress("curverr", &curverr);
1029   treein->SetBranchAddress("run",&runNumber);
1030   treein->SetBranchAddress("ITSAlignData",&itsaligndata);
1031   treein->SetBranchAddress("ITSCalibRespSDD",&itscalibrespsdd);
1032
1033   int ntrks=treein->GetEntries();
1034   printf("Reading %d tracks from %s\n",ntrks,fin);
1035
1036
1037   // open output file
1038   TFile *pointsFile = TFile::Open(fout,"RECREATE");
1039   if (!pointsFile || !pointsFile->IsOpen()) {
1040     printf("Can't open output file %s !",fout);
1041     return -6;
1042   }
1043   AliTrackPointArray *array = new AliTrackPointArray();
1044
1045   // new!
1046   TTree *treeout=(TTree*)treein->Clone("spTree");
1047   treeout->Reset();
1048   treeout->SetBranchAddress("SP", &array);
1049   treeout->SetBranchAddress("curv", &curv);
1050   treeout->SetBranchAddress("curverr", &curverr);
1051   treeout->SetBranchAddress("run",&runNumber);
1052   treeout->SetBranchAddress("ITSAlignData",&itsaligndata);
1053   treeout->SetBranchAddress("ITSCalibRespSDD",&itscalibrespsdd);
1054
1055   // tracks main loop
1056   for (Int_t it=0; it<ntrks; it++) {
1057     if (!(it%5000) ) printf("...processing track n. %d\n",it);
1058
1059     treein->GetEvent(it);
1060
1061     //////////////////////////////
1062
1063     AliTrackPointArray *atp=tpain;
1064     AliTrackPointArray *atps=NULL;
1065     Int_t npts=atp->GetNPoints();
1066
1067     AliTrackPoint p;
1068     // check points in specific places
1069
1070     // build a new track
1071     atps=new AliTrackPointArray(npts);
1072
1073     Int_t npto=0;
1074     for (int i=0; i<npts; i++) {
1075       atp->GetPoint(p,i);
1076
1077       UShort_t volid=atp->GetVolumeID()[i];
1078       Int_t index=AliITSAlignMille2Module::GetIndexFromVolumeID(volid);
1079
1080
1081       // dealign point
1082       // get MODIFIED matrix
1083       TGeoHMatrix *svMatrix = m2[index]->GetSensitiveVolumeMatrix(p.GetVolumeID());
1084       //TGeoHMatrix *svOrigMatrix = mm->GetSensitiveVolumeOrigGlobalMatrix(p.GetVolumeID());
1085
1086       Double_t pg[3],pl[3];
1087       pg[0]=p.GetX();
1088       pg[1]=p.GetY();
1089       pg[2]=p.GetZ();
1090       if (prn) printf("Global coordinates of measured point : X=%f  Y=%f  Z=%f \n",pg[0],pg[1],pg[2]);
1091       svMatrix->MasterToLocal(pg,pl);
1092
1093       // check that things went OK: local y should be 0.
1094       if(TMath::Abs(pl[1])>1.e-6) {
1095         printf("AliAlignmentDataFilterITS::WriteTrackPointsInIdealGeom: ERROR, local y = %f (should be zero)\n",pl[1]);
1096         return -7;
1097       }
1098
1099       if (prn) printf("Local coordinates of measured point : X=%f  Y=%f  Z=%f \n",pl[0],pl[1],pl[2]);
1100
1101       // update covariance matrix
1102       TGeoHMatrix hcov;
1103       Double_t hcovel[9];
1104       hcovel[0]=(Double_t)(p.GetCov()[0]);
1105       hcovel[1]=(Double_t)(p.GetCov()[1]);
1106       hcovel[2]=(Double_t)(p.GetCov()[2]);
1107       hcovel[3]=(Double_t)(p.GetCov()[1]);
1108       hcovel[4]=(Double_t)(p.GetCov()[3]);
1109       hcovel[5]=(Double_t)(p.GetCov()[4]);
1110       hcovel[6]=(Double_t)(p.GetCov()[2]);
1111       hcovel[7]=(Double_t)(p.GetCov()[4]);
1112       hcovel[8]=(Double_t)(p.GetCov()[5]);
1113       hcov.SetRotation(hcovel);
1114       // now rotate in local system
1115       hcov.Multiply(svMatrix);
1116       hcov.MultiplyLeft(&svMatrix->Inverse());
1117       // now hcov is LOCAL COVARIANCE MATRIX
1118
1119       /// get original matrix of sens. vol.
1120       TGeoHMatrix *svOrigMatrix = m2[index]->GetSensitiveVolumeOrigGlobalMatrix(p.GetVolumeID());
1121       // modify global coordinates according with pre-aligment
1122       svOrigMatrix->LocalToMaster(pl,pg);
1123       // now rotate in local system
1124       hcov.Multiply(&svOrigMatrix->Inverse());
1125       hcov.MultiplyLeft(svOrigMatrix);
1126       // hcov is back in GLOBAL RF
1127       Float_t pcov[6];
1128       pcov[0]=hcov.GetRotationMatrix()[0];
1129       pcov[1]=hcov.GetRotationMatrix()[1];
1130       pcov[2]=hcov.GetRotationMatrix()[2];
1131       pcov[3]=hcov.GetRotationMatrix()[4];
1132       pcov[4]=hcov.GetRotationMatrix()[5];
1133       pcov[5]=hcov.GetRotationMatrix()[8];
1134
1135       p.SetXYZ(pg[0],pg[1],pg[2],pcov);
1136       if (prn) printf("New global coordinates of measured point : X=%f  Y=%f  Z=%f \n",pg[0],pg[1],pg[2]);
1137       atps->AddPoint(npto,&p);
1138       if (prn) printf("Adding point[%d] = ( %f , %f , %f )     volid = %d\n",npto,atps->GetX()[npto],atps->GetY()[npto],atps->GetZ()[npto],atps->GetVolumeID()[npto] );
1139       if (prn) p.Print("");
1140
1141       npto++;
1142     }
1143
1144
1145     ////////////////////////////////////////////////////////////
1146     array = atps;
1147     treeout->Fill();
1148
1149     delete atps;
1150     atps=NULL;
1151
1152   } // end loop on tracks
1153
1154   pointsFile->Write();
1155   pointsFile->Close();
1156   tpainfile->Close();
1157
1158   return 0;
1159 }