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