Updated version of the fast TPC simulation which works with the ESD IO (A.Dainese)
authorhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 23 Mar 2006 14:56:27 +0000 (14:56 +0000)
committerhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 23 Mar 2006 14:56:27 +0000 (14:56 +0000)
TPC/AliBarrelRec_TPCparam.C
TPC/AliTPCtrackerParam.cxx
TPC/AliTPCtrackerParam.h
TPC/AliTPCtrackingParamDB.C

index d9db3e4..b04dd08 100644 (file)
-/****************************************************************************
- * This macro performs track and vertex reconstruction in TPC and ITS.      *
- * The ITS Kalman tracker V2 is feeded "with" parameterized TPC tracks.     * 
- *                                                                          *
- * Reconstruction is performed in the following steps:                      *
- *             1) TPC tracking parameterization                             *
- *             2) ITS clusters: slow or fast                                *
- *             3) Primary vertex reconstruction                             *
- *                - read from event header for Pb-Pb events                 *
- *                - determined using points in pixels for pp/pA events      *
- *             4) ITS track finding V2                                      *
- *                - in pp/pA, redetermine the position of primary vertex    *
- *                  using the reconstructed tracks                          *
- *             5) Create a reference file with simulation info (p,PDG...)   *
- *                                                                          *
- * If mode='A' all 5 steps are executed                                     *
- * If mode='B' only steps 4-5 are executed                                  *
- *                                                                          *  
- *  Origin: A.Dainese, Padova,   andrea.dainese@pd.infn.it                  * 
- *  (from AliTPCtest.C & AliITStestV2.C by I.Belikov)                       *
- ****************************************************************************/
-
-// structure for track references
-typedef struct {
-  Int_t lab;
-  Int_t pdg;
-  Int_t mumlab;
-  Int_t mumpdg;
-  Float_t Vx,Vy,Vz;
-  Float_t Px,Py,Pz;
-} RECTRACK;
-
-//===== Functions definition ================================================= 
-
-void CopyVtx(const Char_t *inName,const Char_t *outName);
-
-void ITSFindClustersV2(Char_t SlowOrFast);
-
-void ITSFindTracksV2(Int_t *skipEvt);
-
-void ITSMakeRefFile(Int_t *skipEvt);
-
-void MarkEvtsToSkip(const Char_t *evtsName,Int_t *skipEvt);
-
-void PrimaryVertex(const Char_t *outName,Char_t vtxMode);
-
-void TPCParamTracks(Int_t coll,Double_t Bfield);
-
-Int_t UpdateEvtsToSkip(const Char_t *logName,const Char_t *evtsName);
-
-void VtxFromHeader(const Char_t *outName,Bool_t smear);
-
-void VtxFromTracks(const Char_t *outName);
-
-void ZvtxFromSPD(const Char_t *outName);
-
-//=============================================================================
-
-// number of events to be processed
-Int_t    gNevents;
-// magnetic field
-Double_t gBfieldValue;
-
-void AliBarrelRec_TPCparam(Int_t n=-1,Char_t mode='A') {
-
-  //---------------------------------------------------------------------
-  //                    CONFIGURATION
+void AliBarrelRec_TPCparam(Int_t firstEvent=0,Int_t lastEvent=0) {
   //
-  // _Magnetic_field_
-  gBfieldValue = 0.4;
+  // Macro to create AliESDs.root using parametrized TPC tracking
+  // and AliITStrackerV2
   //
-  // _Type_of_collision_ (needed for TPC tracking parameterization) 
-  // Available choices:   !!! ONLY B = 0.4 TESLA !!!
-  //    collcode = 0  ->   PbPb6000 (HIJING with b<2fm) 
-  //    collcode = 1  ->   low multiplicity: pp or pA
-  Int_t collcode = 1;  
-  // 
-  // _ITS_clusters_reconstruction_
-  // Available choices:  (from AliITStestV2.C)
-  //    SlowOrFast = 's'    slow points
-  //    SlowOrFast = 'f'    fast points
-  Char_t SlowOrFast = 'f';
-  //
-  // _Primary_vertex_for_ITS_tracking_
-  // Available choices:
-  //    Vtx4Tracking = 'H'   from event Header
-  //    --- for Pb-Pb ---
-  //    Vtx4Tracking = 'S'   from event header + Smearing 
-  //                                           (x=15,y=15,z=10) micron
-  //    --- for pp/pA ---
-  //    Vtx4Tracking = 'P'   z from pixels, x,y in(0,0)
-  Char_t Vtx4Tracking = 'P';
-  // _Primary_vertex_for_analysis_ (AliESDVertex stored in tracks file)
-  // Available choices:
-  //    Vtx4Analysis = 'C'   Copy the same used for tracking
-  //    --- for pp/pA ---
-  //    Vtx4Analysis = 'T'   x,y,z from Tracks
-  Char_t Vtx4Analysis = 'T';
-  //
-  //                  END CONFIGURATION
-  //---------------------------------------------------------------------
-
-  const Char_t *name=" AliBarrelRec_TPCparam";
-  printf("\n %s\n",name);
-  gBenchmark->Start(name);
-
-  if(n==-1) { // read number of events to be processed from file
-    TFile *f = new TFile("galice.root");
-    gAlice = (AliRun*)f->Get("gAlice");
-    n = gAlice->GetEventsPerRun();
-    delete gAlice; 
-    gAlice=0;
-    f->Close(); 
-    delete f;
-    printf(" All %d events in file will be processed\n",n);
-  }
-  gNevents = n;
-
-
-  // conversion constant for kalman tracks 
-  AliKalmanTrack::SetConvConst(100/0.299792458/gBfieldValue);
-
-  // Selection of execution mode
-  switch(mode) {
-  case 'A':
-    // Build TPC tracks with parameterization
-    TPCParamTracks(collcode,gBfieldValue);
-  
-    // ITS clusters
-    ITSFindClustersV2(SlowOrFast);
-
-    // Vertex for ITS tracking
-    PrimaryVertex("Vtx4Tracking.root",Vtx4Tracking);
-
-    break;
-  
-  case 'B':
-    printf("       ---> only tracking in ITS <---\n");
-
-    // Update list of events to be skipped
-    if(!UpdateEvtsToSkip("itstracking.log","evtsToSkip.dat")) return;
-
-    break;
-  }
-
-  // Mark events that have to be skipped (if any)
-  Int_t *skipEvt = new Int_t[gNevents];
-  for(Int_t i=0; i<gNevents; i++) skipEvt[i] = 0;
-  if(!gSystem->AccessPathName("evtsToSkip.dat",kFileExists)) 
-    MarkEvtsToSkip("evtsToSkip.dat",skipEvt);
-    
-  // Tracking in ITS
-  ITSFindTracksV2(skipEvt); 
-
-  // Vertex for analysis 
-  PrimaryVertex("AliITStracksV2.root",Vtx4Analysis);
-
-  // Create ITS tracks reference file
-  ITSMakeRefFile(skipEvt);
-  delete [] skipEvt;
-
-
-
-  gBenchmark->Stop(name);
-  gBenchmark->Show(name);
-
-  return;
-}
-//-----------------------------------------------------------------------------
-void CopyVtx(const Char_t *inName,const Char_t *outName) {
-
-  // Open input and output files
-  TFile *inFile = new TFile(inName);
-  TFile *outFile = new TFile(outName,"update");
-
-  TDirectory *curdir;
-  Char_t vname[20];
-
-
-  for(Int_t ev=0; ev<gNevents; ev++) {
-    sprintf(vname,"Vertex_%d",ev);
-    AliESDVertex *vertex = (AliESDVertex*)inFile->Get(vname);
-    if(!vertex) continue;
-    curdir = gDirectory;
-    outFile->cd();
-    vertex->Write();
-    curdir->cd();
-    vertex = 0;
-  }
-
-  inFile->Close();
-  outFile->Close();
-  delete inFile;
-  delete outFile;
-
-  return;
-}
-//-----------------------------------------------------------------------------
-void ITSFindClustersV2(Char_t SlowOrFast) {
-
-  printf("\n------------------------------------\n");
-
-  const Char_t *name="ITSFindClustersV2";
-  printf("\n %s\n",name);
-  gBenchmark->Start(name);
-
-  //---  taken from AliITStestV2.C--------------------------------------
+  // Input files:
+  // - galice.root
+  // - Kinematics.root
+  // - TrackReferences.root
+  // - ITS.RecPoints.root (use AliRecontruction class)
+  // - ITS.Vertex.root (use $ALICE_ROOT/ITS/AliITSVertexerZTest.C)
   //
-  if (SlowOrFast=='f') {
-    //cerr<<"Fast AliITSRecPoint(s) !\n";
-    //gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSHits2FastRecPoints.C");
-    //AliITSHits2FastRecPoints();
-  } else {
-    gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSHits2SDigits.C");
-    AliITSHits2SDigits();
-    gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSSDigits2Digits.C");
-    AliITSSDigits2Digits();
-    //gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSDigits2RecPoints.C");
-    //AliITSDigits2RecPoints();
-  }
-  gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSFindClustersV2.C");
-  AliITSFindClustersV2(SlowOrFast,gNevents);
+  // A. Dainese - LNL
   //
-  //--------------------------------------------------------------------
 
 
-  gBenchmark->Stop(name);
-  gBenchmark->Show(name);
+   /**** Initialization of the NewIO *******/
+
+   if (gAlice) {
+      delete gAlice->GetRunLoader();
+      delete gAlice; 
+      gAlice=0;
+   }
+
+   AliRunLoader *rl = AliRunLoader::Open("galice.root");
+   if (rl == 0x0) {
+      cerr<<"Can not open session"<<endl;
+      return;
+   }
+   Int_t retval = rl->LoadgAlice();
+   if (retval) {
+      cerr<<"LoadgAlice returned error"<<endl;
+      delete rl;
+      return;
+   }
+   retval = rl->LoadHeader();
+   if (retval) {
+      cerr<<"LoadHeader returned error"<<endl;
+      delete rl;
+      return;
+   }
+   gAlice=rl->GetAliRun();
+       
+
+   TDatabasePDG *DataBase = TDatabasePDG::Instance();
+
+   // Get field from galice.root
+   AliMagF *fiel = (AliMagF*)gAlice->Field();
+   Double_t fieval=TMath::Abs((Double_t)fiel->SolenoidField()/10.);
+   // Set the conversion constant between curvature and Pt
+   AliTracker::SetFieldMap(fiel,kTRUE);
+
+   /**** The TPC corner ********************/
+
+   Int_t collcode = 1; // pp collisions
+   AliTPCtrackerParam tpcTrackerPar(collcode,fieval);
+   tpcTrackerPar.Init();
+
+   //**** Switch on the PID class (mandatory!)
+   AliPID pid;
+
+   /**** The ITS corner ********************/
+
+   AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
+   if (itsl == 0x0) {
+      cerr<<"Cannot get the ITS loader"<<endl;
+      return;
+   }
+   itsl->LoadRecPoints("read");
+
+   AliITS *dITS = (AliITS*)gAlice->GetDetector("ITS");
+   if (!dITS) {
+      cerr<<"Cannot find the ITS detector !"<<endl;
+      return;
+   }
+   AliITSgeom *geom = dITS->GetITSgeom();
+
+   //An instance of the ITS tracker
+   AliITStrackerV2 itsTracker(geom);
+   Int_t ITSclusters[6] = {1,1,1,1,1,1};
+   itsTracker.SetLayersNotToSkip(ITSclusters);
+
+   /***** The TREE is born *****/
+   
+   TTree *esdTree=new TTree("esdTree","Tree with ESD objects");
+   AliESD *event=0;
+   esdTree->Branch("ESD","AliESD",&event);
+   
+   if(firstEvent>rl->GetNumberOfEvents()) firstEvent=rl->GetNumberOfEvents()-1;
+   if(lastEvent>rl->GetNumberOfEvents())  lastEvent=rl->GetNumberOfEvents()-1;
+   cout<<" Number of events: "<<1+lastEvent-firstEvent<<endl;
+   
+   TFile *ppZ = TFile::Open("ITS.Vertex.root"); // z vertices from SPD
+   AliESDVertex *myvertex = new AliESDVertex();
+   Char_t zver[100];
+   Double_t vtx[3];
+   Double_t cvtx[6];
+
+
+   //<----------------------------------The Loop over events begins
+   TStopwatch timer;
+   Int_t trc;
+   for(Int_t i=firstEvent; i<=lastEvent; i++) { 
+     
+     cerr<<" Processing event number : "<<i<<endl;
+     AliESD *event = new AliESD(); 
+     event->SetRunNumber(gAlice->GetRunNumber());
+     event->SetEventNumber(i);
+     event->SetMagneticField(gAlice->Field()->SolenoidField());
+     rl->GetEvent(i);
+
+    //***** Primary vertex 
+     sprintf(zver,"Event%d/Vertex",i);
+     myvertex = (AliESDVertex*)ppZ->Get(zver);
+     if(!myvertex) {
+       esdTree->Fill(); delete event;
+       continue;
+     }      
+     event->SetVertex(myvertex);
+     myvertex->GetXYZ(vtx);
+     myvertex->GetCovMatrix(cvtx);
+
+     //***** TPC tracking
+     if ( (trc=tpcTrackerPar.BuildTPCtracks(event)) ) {
+       printf("exiting TPC tracker with code %d in event %d\n",trc,i);
+       esdTree->Fill(); delete event;
+       continue;
+     }
+
+     //***** ITS tracking
+     itsTracker.SetVertex(vtx,cvtx);
+     TTree *itsTree=itsl->TreeR();
+     if (!itsTree) {
+        cerr<<"Can't get the ITS cluster tree !\n";
+       esdTree->Fill(); delete event;
+        return;
+     }     
+     itsTracker.UnloadClusters();
+     itsTracker.LoadClusters(itsTree);
+     if ( (trc=itsTracker.Clusters2Tracks(event)) ) {
+       printf("exiting ITS tracker with code %d in event %d\n",trc,i);
+       esdTree->Fill(); delete event;
+       continue;
+     }
+
+
+     esdTree->Fill();
+     delete event;
+
+   }//<-----------------------------------The Loop over events ends here
+   timer.Stop(); timer.Print();
+
+   //        The AliESDs.root is born
+   TFile *ef = TFile::Open("AliESDs.root","RECREATE"); 
+   if (!ef || !ef->IsOpen()) {cerr<<"Can't open AliESDs.root !\n"; return;}
+
+   esdTree->Write(); //Write the TREE and close everything
+   delete esdTree;
+   ef->Close();
+
+   delete rl;
 
    return;
 }
-//-----------------------------------------------------------------------------
-Int_t ITSFindTracksV2(Int_t *skipEvt) {
-
-  printf("\n------------------------------------\n");
-  
-  const Char_t *name="ITSFindTracksV2";
-  printf("\n %s\n",name);
-  gBenchmark->Start(name);
-
-  TFile *outFile     = new TFile("AliITStracksV2.root","recreate");
-  TFile *inTPCtrks   = new TFile("AliTPCtracksParam.root");
-  TFile *inVertex    = new TFile("Vtx4Tracking.root");
-  TFile *inClusters  = new TFile("AliITSclustersV2.root");
-  
-  AliITSgeom *geom=(AliITSgeom*)inClusters->Get("AliITSgeom");
-  if(!geom) { printf("can't get ITS geometry !\n"); return;}
-  
-  Double_t vtx[3];
-  Int_t flag1stPass,flag2ndPass;
-  Char_t vname[20];
-
-  // open logfile for done events
-  FILE *logfile = fopen("itstracking.log","w");
-
-  // Instantiate AliITStrackerV2
-  AliITStrackerV2 tracker(geom);
-
-  // loop on events
-  for(Int_t ev=0; ev<gNevents; ev++){
-    // write to logfile of begun events
-    fprintf(logfile,"%d\n",ev);
-
-    if(skipEvt[ev]) continue;
-    printf(" --- Processing event %d ---\n",ev);
-
-    // pass event number to the tracker
-    tracker.SetEventNumber(ev);
-
-    // set position of primary vertex
-    sprintf(vname,"Vertex_%d",ev);
-    AliESDVertex *vertex = (AliESDVertex*)inVertex->Get(vname);
-    if(vertex) {
-      vertex->GetXYZ(vtx);
-      delete vertex;
-    } else {
-      printf(" AliESDVertex not found for event %d\n",ev);
-      printf(" Using (0,0,0) for ITS tracking\n");
-      vtx[0] = vtx[1] = vtx[2] = 0.;
-    }
-
-    flag1stPass=1; // vtx constraint
-    flag2ndPass=0; // no vtx constraint
-
-    // no vtx constraint if vertex not found
-    if(vtx[2]<-999.) {
-      flag1stPass=0;
-      vtx[2]=0.;
-    }
-
-    tracker.SetVertex(vtx);  
-
-    // setup vertex constraint in the two tracking passes
-    Int_t flags[2];
-    flags[0]=flag1stPass;
-    tracker.SetupFirstPass(flags);
-    flags[0]=flag2ndPass;
-    tracker.SetupSecondPass(flags);
-    
-    // find the tracks
-    tracker.Clusters2Tracks(inTPCtrks,outFile);
-
-  } // loop on events
-
-  fprintf(logfile,"%d\n",gNevents); //this means all evts are successfully completed
-  fclose(logfile);
-
-  delete geom;
-
-  inTPCtrks->Close();
-  inClusters->Close();
-  inVertex->Close();
-  outFile->Close();
-
-  gBenchmark->Stop(name);
-  gBenchmark->Show(name);
-  
-  return;
-}
-//-----------------------------------------------------------------------------
-void ITSMakeRefFile(Int_t *skipEvt) {
-
-  printf("\n------------------------------------\n");
-
-  const Char_t *name="ITSMakeRefFile";
-  printf("\n %s\n",name);
-  gBenchmark->Start(name);
-  
-  
-  TFile *out = TFile::Open("ITStracksRefFile.root","recreate");
-  TFile *trk = TFile::Open("AliITStracksV2.root");
-  TFile *kin = TFile::Open("galice.root");
-
-  
-  // Get gAlice object from file
-  gAlice=(AliRun*)kin->Get("gAlice");
-  
-  Int_t label;
-  TParticle *Part;  
-  TParticle *Mum;
-  RECTRACK rectrk;
-  
-
-  for(Int_t ev=0; ev<gNevents; ev++){
-    if(skipEvt[ev]) continue;
-    printf(" --- Processing event %d ---\n",ev);
-
-    gAlice->GetEvent(ev);  
-
-    trk->cd();
-
-    // Tree with ITS tracks
-    char tname[100];
-    sprintf(tname,"TreeT_ITS_%d",ev);
-
-    TTree *tracktree=(TTree*)trk->Get(tname);
-    if(!tracktree) continue;
-    AliITStrackV2 *itstrack=new AliITStrackV2; 
-    tracktree->SetBranchAddress("tracks",&itstrack);
-    Int_t nentr=(Int_t)tracktree->GetEntries();
-
-    // Tree for true track parameters
-    char ttname[100];
-    sprintf(ttname,"Tree_Ref_%d",ev);
-    TTree *reftree = new TTree(ttname,"Tree with true track params");
-    reftree->Branch("rectracks",&rectrk,"lab/I:pdg:mumlab:mumpdg:Vx/F:Vy:Vz:Px:Py:Pz");
-
-    for(Int_t i=0; i<nentr; i++) {
-      tracktree->GetEvent(i);
-      label = TMath::Abs(itstrack->GetLabel());
-
-      Part = (TParticle*)gAlice->Particle(label);
-      rectrk.lab=label;
-      rectrk.pdg=Part->GetPdgCode();
-      rectrk.mumlab = Part->GetFirstMother();
-      if(Part->GetFirstMother()>=0) {
-       Mum = (TParticle*)gAlice->Particle(Part->GetFirstMother());
-       rectrk.mumpdg=Mum->GetPdgCode();
-      } else {
-       rectrk.mumpdg=-1;
-      }
-      rectrk.Vx=Part->Vx();
-      rectrk.Vy=Part->Vy();
-      rectrk.Vz=Part->Vz();
-      rectrk.Px=Part->Px();
-      rectrk.Py=Part->Py();
-      rectrk.Pz=Part->Pz();
-      
-      reftree->Fill();
-    } // loop on tracks   
-
-    out->cd();
-    reftree->Write();
-
-    delete itstrack;
-    delete reftree;
-  } // loop on events
-
-  trk->Close();
-  kin->Close();
-  out->Close();
-  
-  gBenchmark->Stop(name);
-  gBenchmark->Show(name);
-  
-
-  return;
-}
-//-----------------------------------------------------------------------------
-void MarkEvtsToSkip(const Char_t *evtsName,Int_t *skipEvt) {
-
-  printf("\n------------------------------------\n");
-  printf("\nChecking for events to skip...\n");
-
-  Int_t evt,ncol;
-
-  FILE *f = fopen(evtsName,"r");
-  while(1) {
-    ncol = fscanf(f,"%d",&evt);
-    if(ncol<1) break;
-    skipEvt[evt] = 1;
-    printf(" event %d will be skipped\n",evt);
-  }
-  fclose(f);
-
-  return;
-}
-//-----------------------------------------------------------------------------
-void PrimaryVertex(const Char_t *outName,Char_t vtxMode) {
-
-  printf("\n------------------------------------\n");
-
-  const Char_t *name="PrimaryVertex";
-  printf("\n %s\n",name);
-  gBenchmark->Start(name);
-
-  switch(vtxMode) {
-  case 'H':
-    printf(" ... from event header\n");
-    VtxFromHeader(outName,kFALSE);
-    break;
-  case 'S':
-    printf(" ... from event header + smearing\n");
-    VtxFromHeader(outName,kTRUE);
-    break;
-  case 'P':
-    printf(" ... z from pixels for pp/pA\n");
-    ZvtxFromSPD(outName);
-    break;
-  case 'T':
-    printf(" ... from tracks for pp/pA\n");
-    VtxFromTracks(outName);
-    break;
-  case 'C':
-    printf(" ... copied from Vtx4Tracking.root to AliITStracksV2.root\n");
-    CopyVtx("Vtx4Tracking.root",outName);
-    break;
-  }
-
-  gBenchmark->Stop(name);
-  gBenchmark->Show(name);
-
-  return;
-}
-//-----------------------------------------------------------------------------
-void TPCParamTracks(Int_t coll,Double_t Bfield) {
-
-  printf("\n------------------------------------\n");
-
-  const Char_t *name="TPCParamTracks";
-  printf("\n %s\n",name);
-  gBenchmark->Start(name);
-
-  TFile *outFile=TFile::Open("AliTPCtracksParam.root","recreate");
-  TFile *inFile =TFile::Open("galice.root");
-  AliTPCtrackerParam tracker(coll,Bfield,gNevents);
-  tracker.BuildTPCtracks(inFile,outFile);
-
-  delete gAlice; gAlice=0;
-
-  inFile->Close();
-  outFile->Close();
-
-  gBenchmark->Stop(name);
-  gBenchmark->Show(name);
-
-  return;
-}
-//-----------------------------------------------------------------------------
-Int_t UpdateEvtsToSkip(const Char_t *logName,const Char_t *evtsName) {
-
-    if(!gSystem->AccessPathName(logName,kFileExists)) { 
-      FILE *ifile = fopen(logName,"r");
-      Int_t lEvt=0,nCol=1;
-      while(nCol>0) {
-       nCol = fscanf(ifile,"%d",&lEvt);
-      }
-      fclose(ifile);
-      if(lEvt==gNevents) { 
-       printf(" All events already reconstructed\n"); 
-       return 0;  
-      } else {
-       FILE *ofile = fopen("evtsToSkip.dat","a");
-       fprintf(ofile,"%d\n",lEvt);
-       fclose(ofile);
-      }
-    } else { 
-      printf("File itstracking.log not found\n");
-    }
-
-    return 1;
-}
-//-----------------------------------------------------------------------------
-void VtxFromHeader(const Char_t *outName,Bool_t smear) {
-
-  TDatime t;
-  UInt_t seed = t.Get();
-  gRandom->SetSeed(seed);
-
-  TFile *galice  = new TFile("galice.root");  
-  TFile *outFile = new TFile(outName,"update");
-
-  TDirectory *curdir;
-  Double_t pos[3],sigma[3];
-  if(smear) {
-    sigma[0]=15.e-4;
-    sigma[1]=15.e-4;
-    sigma[2]=10.e-4;
-  } else {
-    sigma[0]=0.;
-    sigma[1]=0.;
-    sigma[2]=0.;
-  }
-  Char_t vname[20];
-
-  galice->cd();
-
-  for(Int_t ev=0; ev<gNevents; ev++){
-    printf(" event %d\n",ev);
-    sprintf(vname,"Vertex_%d",ev);
-    TArrayF o = 0;
-    o.Set(3);
-    AliHeader* header = 0;
-    TTree* treeE = (TTree*)gDirectory->Get("TE");
-    treeE->SetBranchAddress("Header",&header);
-    treeE->GetEntry(ev);
-    AliGenEventHeader* genHeader = header->GenEventHeader();
-    if(genHeader) {
-      // get primary vertex position
-      genHeader->PrimaryVertex(o);
-      pos[0] = (Double_t)o[0];
-      pos[1] = (Double_t)o[1];
-      pos[2] = (Double_t)o[2];
-      if(smear) {
-       pos[0] = gRandom->Gaus(pos[0],sigma[0]);
-       pos[1] = gRandom->Gaus(pos[1],sigma[1]);
-       pos[2] = gRandom->Gaus(pos[2],sigma[2]);
-      }
-      // create AliESDVertex
-      AliESDVertex *vertex = new AliESDVertex(pos,sigma,vname);
-    } else {
-      printf(" ! event header not found : setting vertex to (0,0,0) !");
-      pos[0] = 0.;
-      pos[1] = 0.;
-      pos[2] = 0.;
-      // create AliESDVertex
-      AliESDVertex *vertex = new AliESDVertex(pos,sigma,vname);
-    }    
-    delete header;
-    // write AliESDVertex to file
-    curdir = gDirectory;
-    outFile->cd();
-    if(smear) {
-      vertex->SetTitle("vertex from header, smeared");
-    } else {
-      vertex->SetTitle("vertex from header");
-    }
-    vertex->Write();
-    curdir->cd();
-    vertex = 0;
-  }
-
-  outFile->Close();
-  galice->Close();
-
-  delete outFile;
-  delete galice;
-
-  return;
-}
-//-----------------------------------------------------------------------------
-void VtxFromTracks(const Char_t *outName) {
-
-  // Open input and output files
-  TFile *inFile  = new TFile("AliITStracksV2.root");
-  TFile *outFile = new TFile(outName,"update");
-
-  // set AliRun object to 0
-  if(gAlice) gAlice = 0;
-
-  // Create vertexer
-  AliITSVertexerTracks *vertexer = 
-    new AliITSVertexerTracks(inFile,outFile,gBfieldValue);
-  vertexer->SetFirstEvent(0);
-  vertexer->SetLastEvent(gNevents-1);
-  vertexer->SetDebug(0);
-  vertexer->PrintStatus();
-  // Find vertices
-  vertexer->FindVertices();
-
-  delete vertexer;
-
-  inFile->Close();
-  outFile->Close();
-  delete inFile;
-  delete outFile;
-
-  return;
-}
-//-----------------------------------------------------------------------------
-void ZvtxFromSPD(const Char_t *outName) {
-
-  // create fast RecPoints, which are used for vertex finding
-  cerr<<"Fast AliITSRecPoint(s) !\n";
-  gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSHits2FastRecPoints.C");
-  AliITSHits2FastRecPoints(0,gNevents-1);
-
-  // delphi ---> azimuthal range to accept tracklets
-  // window ---> window in Z around the peak of tracklets proj. in mm
-  Float_t delphi=0.05;
-  Float_t window=3.;
-  Float_t initx=0.;
-  Float_t inity=0.;
-
-  TFile *infile = new TFile("galice.root");
-  TFile *outfile = new TFile(outName,"update");
-
-  AliITSVertexerPPZ *vertexer = new AliITSVertexerPPZ(infile,outfile,initx,inity);
-  vertexer->SetFirstEvent(0);
-  vertexer->SetLastEvent(gNevents-1);
-  vertexer->SetDebug(0);
-  vertexer->SetDiffPhiMax(delphi);
-  vertexer->SetWindow(window);
-  vertexer->PrintStatus();
-  vertexer->FindVertices();
-  delete vertexer;
-  vertexer=0;
-
-  outfile->Close();
-  infile->Close();
-  delete infile;
-  delete outfile;
-
-
-  return;
-}
-//-----------------------------------------------------------------------------
-
-
-
-
-
-
-
-
-
index 00acbc6..3067362 100644 (file)
  * - All the code necessary to create a BataBase has been included in     *
  *   class (see the macro AliTPCtrackingParamDB.C for the details).       *
  *                                                                        *
- *  Origin: Andrea Dainese, Padova - e-mail: andrea.dainese@pd.infn.it    *
+ * 2006/03/16: Adapted to ESD input/output                                *
  *                                                                        *
+ *  Origin: Andrea Dainese, Padova - e-mail: andrea.dainese@pd.infn.it    *
+ *       adapted to ESD output by Marcello Lunardon, Padova               *
  **************************************************************************/
 //  *
 // This is a dummy comment
@@ -69,6 +71,7 @@
 #include <TH2F.h>
 #include <TLegend.h>
 #include <TLine.h>
+#include <TMatrixD.h>
 #include <TParticle.h>
 #include <TStyle.h>
 #include <TSystem.h>
@@ -117,25 +120,21 @@ ClassImp(AliTPCtrackerParam)
 
 //-----------------------------------------------------------------------------
 AliTPCtrackerParam::AliTPCtrackerParam(Int_t kcoll, Double_t kBz,
-                                      Int_t kn, const char* evfoldname):
+                                      const char* evfoldname):
   fEvFolderName(evfoldname) {
 //-----------------------------------------------------------------------------
 // This is the class conctructor 
 //-----------------------------------------------------------------------------
 
-  fNevents = kn;         // events to be processed
   fBz = kBz;             // value of the z component of L3 field (Tesla)
   fColl = kcoll;         // collision code (0: PbPb6000; 1: pp)
   fSelAndSmear = kTRUE; // by default selection and smearing are done
 
-  if(fBz!=0.4) {
-    cerr<<"AliTPCtrackerParam::AliTPCtrackerParam:  Invalid field!\n";
-    cerr<<"      Available:  0.4"<<endl;
+  if(fBz!=0.4 && fBz!=0.5) {
+    Fatal("AliTPCtrackerParam","AliTPCtrackerParam::AliTPCtrackerParam:  Invalid field!\n      Available:  0.4 or 0.5");
   }
-  if(fColl!=0 && fColl!=1) { 
-    cerr<<"AliTPCtrackerParam::AliTPCtrackerParam:  Invalid collision!\n";
-    cerr<<"      Available:  0   ->   PbPb6000"<<endl;
-    cerr<<"                  1   ->   pp"<<endl; 
+  if(fColl!=0 && fColl!=1) {
+    Fatal("AliTPCtrackerParam","AliTPCtrackerParam::AliTPCtrackerParam:  Invalid collision!\n      Available:  0   ->   PbPb6000\n                  1   ->   pp"); 
   }
 
   fDBfileName = gSystem->Getenv("ALICE_ROOT");  
@@ -144,6 +143,8 @@ AliTPCtrackerParam::AliTPCtrackerParam(Int_t kcoll, Double_t kBz,
   if(fColl==0) fDBfileName.Append("PbPb6000");
   if(fColl==1) fDBfileName.Append("pp");
   if(fBz==0.4) fDBfileName.Append("_B0.4T.root");
+  // use same DB for 0.4 and 0.5 T; for 0.5 T, correction done in CookTrack()
+  if(fBz==0.5) fDBfileName.Append("_B0.4T.root");
 }
 //-----------------------------------------------------------------------------
 AliTPCtrackerParam::~AliTPCtrackerParam() {}
@@ -175,307 +176,244 @@ AliTPCtrackerParam::AliTPCseedGeant::AliTPCseedGeant(
       fAlpha *= TMath::Pi();
 }
 //-----------------------------------------------------------------------------
-Int_t AliTPCtrackerParam::BuildTPCtracks(const TFile *inp, TFile *out) {
+Int_t AliTPCtrackerParam::Init() {
 //-----------------------------------------------------------------------------
-// This function creates the TPC parameterized tracks
+// This function reads the DB from file
 //-----------------------------------------------------------------------------
 
-  Error("BuildTPCtracks","in and out parameters ignored. new io");
-  
-/********************************************/
-  AliRunLoader* rl = AliRunLoader::GetRunLoader(fEvFolderName);
-  if (rl == 0x0)
-   {
-     Error("BuildTPCtracks","Can not get Run Loader from event folder named %s.",
-           fEvFolderName.Data());
-     return 2;
-   }
-  AliLoader* tpcloader = rl->GetLoader("TPCLoader");
-  if (tpcloader == 0x0)
-   {
-     Error("BuildTPCtracks","Can not get TPC Loader from Run Loader.");
-     return 3;
-   }
-  
-/********************************************/
-
-  TFile *fileDB=0;
-  TTree *covTreePi[50];
-  TTree *covTreeKa[50];
-  TTree *covTreePr[50];
-  TTree *covTreeEl[50];
-  TTree *covTreeMu[50];
-
   if(fSelAndSmear) {
-    cerr<<"+++\n+++ Reading DataBase from:\n+++ "<<
-      fDBfileName.Data()<<"\n+++\n"; 
+    printf("+++\n+++ Reading DataBase from:%s\n+++\n+++\n",fDBfileName.Data()); 
     // Read paramters from DB file
     if(!ReadAllData(fDBfileName.Data())) {
-      cerr<<"AliTPCtrackerParam::BuildTPCtracks: \
-             Could not read data from DB\n\n"; return 1; 
+      printf("AliTPCtrackerParam::BuildTPCtracks: \
+             Could not read data from DB\n\n"); return 1; 
     }
+    
+  } else printf("\n ! Creating ALL TRUE tracks at TPC inner radius !\n\n");
 
-    // Read the trees with regularized cov. matrices from DB
-    TString str;
-    fileDB = TFile::Open(fDBfileName.Data());
-    Int_t nBinsPi = fDBgridPi.GetTotBins();
-    for(Int_t l=0; l<nBinsPi; l++) {
-      str = "/CovMatrices/Pions/CovTreePi_bin";
-      str+=l;
-      covTreePi[l] = (TTree*)fileDB->Get(str.Data());
-    }
-    Int_t nBinsKa = fDBgridKa.GetTotBins();
-    for(Int_t l=0; l<nBinsKa; l++) {
-      str = "/CovMatrices/Kaons/CovTreeKa_bin";
-      str+=l;
-      covTreeKa[l] = (TTree*)fileDB->Get(str.Data());
-    }
-    Int_t nBinsPr = fDBgridPr.GetTotBins();
-    for(Int_t l=0; l<nBinsPr; l++) {
-      if(fColl==0) {      
-       str = "/CovMatrices/Pions/CovTreePi_bin";
-      } else {
-       str = "/CovMatrices/Protons/CovTreePr_bin";
-      }
-      str+=l;
-      covTreePr[l] = (TTree*)fileDB->Get(str.Data());
-    }
-    Int_t nBinsEl = fDBgridEl.GetTotBins();
-    for(Int_t l=0; l<nBinsEl; l++) {
-      str = "/CovMatrices/Electrons/CovTreeEl_bin";
-      str+=l;
-      covTreeEl[l] = (TTree*)fileDB->Get(str.Data());
-    }
-    Int_t nBinsMu = fDBgridMu.GetTotBins();
-    for(Int_t l=0; l<nBinsMu; l++) {
-      if(fColl==0) {      
-       str = "/CovMatrices/Pions/CovTreePi_bin";
-      } else {
-       str = "/CovMatrices/Muons/CovTreeMu_bin";
-      }
-      str+=l;
-      covTreeMu[l] = (TTree*)fileDB->Get(str.Data());
-    }
 
-  } else cerr<<"\n ! Creating ALL TRUE tracks at TPC 1st hit !\n\n";
+  // Check if value in the galice file is equal to selected one (fBz)
+  AliMagF *fiel = (AliMagF*)gAlice->Field();
+  Double_t fieval=TMath::Abs((Double_t)fiel->SolenoidField()/10.);
+  printf("Magnetic field is %6.2f Tesla\n",fieval);
+  if(fBz!=fieval) {
+    printf("AliTPCtrackerParam::BuildTPCtracks:  Invalid field!");
+    printf("Field selected is: %f T\n",fBz);
+    printf("Field found on file is: %f T\n",fieval);
+    return 1;
+  }
+
+  // Set the conversion constant between curvature and Pt
+  AliTracker::SetFieldMap(fiel,kTRUE);
+
+  return 0;
+}
+//-----------------------------------------------------------------------------
+Int_t AliTPCtrackerParam::BuildTPCtracks(AliESD *event) {
+//-----------------------------------------------------------------------------
+// This function creates the TPC parameterized tracks and writes them
+// as AliESDtrack objects in the ESD event
+//-----------------------------------------------------------------------------
 
-  TFile *infile=(TFile*)inp;
-  infile->cd();
 
-  // Get gAlice object from file
+  if(!event) return -1;
+
+  AliRunLoader* rl = AliRunLoader::GetRunLoader(fEvFolderName);
+  if (rl == 0x0) {
+    Error("BuildTPCtracks","Can not get Run Loader from event folder named %s",
+         fEvFolderName.Data());
+    return 2;
+  }
+  AliLoader* tpcloader = rl->GetLoader("TPCLoader");
+  if (tpcloader == 0x0) {
+    Error("BuildTPCtracks","Can not get TPC Loader from Run Loader.");
+    return 3;
+  }
   
-  rl->LoadgAlice();
-  rl->LoadHeader();
+  // Get gAlice object from file  
+  if(!gAlice) rl->LoadgAlice();
+  //rl->LoadHeader();
+  rl->LoadKinematics();
   tpcloader->LoadHits("read");
   
   if(!(gAlice=rl->GetAliRun())) {
-    cerr<<"Can not get gAlice from Run Loader !\n";
+    printf("Cannot get gAlice from Run Loader !\n");
     return 1;
   }
 
-  // Check if value in the galice file is equal to selected one (fBz)
-  AliMagF *fiel = (AliMagF*)gAlice->Field();
-  Double_t fieval=(Double_t)fiel->SolenoidField()/10.;
-  printf("Magnetic field is %6.2f Tesla\n",fieval);
-  if(fBz!=fieval) {
-    cerr<<"AliTPCtrackerParam::BuildTPCtracks:  Invalid field!"<<endl;
-    cerr<<"Field selected is: "<<fBz<<" T\n";
-    cerr<<"Field found on file is: "<<fieval<<" T\n";
-    return 0;
-  }
-
   // Get TPC detector 
   AliTPC *tpc=(AliTPC*)gAlice->GetDetector("TPC");
-  Int_t ver = tpc->IsVersion(); 
-  cerr<<"+++ TPC version "<<ver<<" has been found !\n";
 
   rl->CdGAFile();
-  AliTPCParam *digp=(AliTPCParam*)infile->Get("75x40_100x60");
+  AliTPCParam *digp=(AliTPCParam*)gDirectory->Get("75x40_100x60");
   if(digp){
     delete digp;
     digp = new AliTPCParamSR();
   }
-  else digp=(AliTPCParam*)infile->Get("75x40_100x60_150x60");
+  else digp=(AliTPCParam*)gDirectory->Get("75x40_100x60_150x60");
   
   if(!digp) { cerr<<"TPC parameters have not been found !\n"; return 1; }
   tpc->SetParam(digp);
 
-  // Set the conversion constant between curvature and Pt
-  AliTracker::SetFieldMap(fiel,kTRUE);
-
   TParticle       *part=0;
   AliTPCseedGeant *seed=0;
   AliTPCtrack     *tpctrack=0;
   Double_t     sPt,sEta;
-  Int_t        cl=0,bin,label,pdg,charge;
+  Int_t        bin,label,pdg,charge;
   Int_t        tracks;
   Int_t        nParticles,nSeeds,arrentr;
-  Char_t       tname[100];
   //Int_t nSel=0,nAcc=0;
 
+  Int_t evt=event->GetEventNumber();
+  
+  tracks=0;
 
-  // loop over first n events in file
-  for(Int_t evt=0; evt<fNevents; evt++){
-    cerr<<"+++\n+++ Processing event "<<evt<<"\n+++\n";
-    tracks=0;
-
-    // tree for TPC tracks
-    sprintf(tname,"TreeT_TPC_%d",evt);
-    TTree *tracktree = new TTree(tname,"Tree with TPC tracks");
-    tracktree->Branch("tracks","AliTPCtrack",&tpctrack,20000,0);
-
-    // array for TPC tracks
-    TObjArray tArray(20000);
-
-    // array for TPC seeds with geant info
-    TObjArray sArray(20000);
-
-    // get the particles stack
-    nParticles = (Int_t)gAlice->GetEvent(evt);
+  // array for TPC tracks
+  TObjArray tArray(20000);
+  
+  // array for TPC seeds with geant info
+  TObjArray sArray(20000);
+  
+  // get the particles stack
+  nParticles = (Int_t)gAlice->GetEvent(evt);
     
-    Bool_t   *done     = new Bool_t[nParticles];
-    Int_t    *pdgCodes = new Int_t[nParticles];
-    Double_t *ptkine   = new Double_t[nParticles];
-    Double_t *pzkine   = new Double_t[nParticles];
-
-    // loop on particles and store pdg codes
-    for(Int_t l=0; l<nParticles; l++) {
-      part        = (TParticle*)gAlice->GetMCApp()->Particle(l);
-      pdgCodes[l] = part->GetPdgCode();
-      ptkine[l]   = part->Pt();
-      pzkine[l]   = part->Pz();
-      done[l]     = kFALSE;
-    }
-    cerr<<"+++\n+++ Number of particles in event "<<evt<<":  "<<nParticles<<
-      "\n+++\n";
-
-    cerr<<"\n ********** MAKING SEEDS *****************"<<endl<<endl;
-
-
-    // Create the seeds for the TPC tracks at the inner radius of TPC
-    if(fColl==0) {
-      // Get TreeH with hits
-      TTree *th = tpcloader->TreeH(); 
-      MakeSeedsFromHits(tpc,th,sArray);
-    } else {
-      // Get TreeTR with track references
-      TTree *ttr = rl->TreeTR();
-      MakeSeedsFromRefs(ttr,sArray);
-    }
-
-
-    nSeeds = sArray.GetEntries();
-    cerr<<"\n\n+++\n+++ Number of seeds: "<<nSeeds<<"\n+++\n";
-
-
-    cerr<<"\n ********** BUILDING TRACKS **************"<<endl<<endl;
-
-    // loop over entries in sArray
-    for(Int_t l=0; l<nSeeds; l++) {
-      //if(l%1000==0) cerr<<"  --- Processing seed "
-      //               <<l<<" of "<<nSeeds<<" ---\r";
-
-      seed = (AliTPCseedGeant*)sArray.At(l);
-
-      // this is TEMPORARY: only for reconstruction of pp production for charm
-      if(fColl==1) cl = CheckLabel(seed,nParticles,ptkine,pzkine);
-      if(cl) continue;
+  Bool_t   *done     = new Bool_t[nParticles];
+  Int_t    *pdgCodes = new Int_t[nParticles];
+  
+  // loop on particles and store pdg codes
+  for(Int_t l=0; l<nParticles; l++) {
+    part        = (TParticle*)gAlice->GetMCApp()->Particle(l);
+    pdgCodes[l] = part->GetPdgCode();
+    done[l]     = kFALSE;
+  }
 
-      // Get track label
-      label = seed->GetLabel();
-
-      // check if this track has already been processed
-      if(done[label]) continue;
-      // PDG code & electric charge
-      pdg = pdgCodes[label];
-      if(pdg>200 || pdg==-11 || pdg==-13) { charge=1; }
-      else if(pdg<-200 || pdg==11 || pdg==13) { charge=-1; }
-      else continue;
-      pdg = TMath::Abs(pdg);
-      if(pdg>3000) pdg=211;
-      if(fSelAndSmear) SetParticle(pdg);
+  printf("+++ Number of particles:  %d\n",nParticles);
 
+  // Create the seeds for the TPC tracks at the inner radius of TPC
+  if(fColl==0) {
+    // Get TreeH with hits
+    TTree *th = tpcloader->TreeH(); 
+    MakeSeedsFromHits(tpc,th,sArray);
+  } else {
+    // Get TreeTR with track references
+    rl->LoadTrackRefs();
+    TTree *ttr = rl->TreeTR();
+    if (!ttr) return -3;
+    MakeSeedsFromRefs(ttr,sArray);
+  }
 
-      sPt  = seed->GetPt();
-      sEta = seed->GetEta();
-      
-      // Apply selection according to TPC efficiency
-      //if(TMath::Abs(pdg)==211) nAcc++;
-      if(fSelAndSmear && !SelectedTrack(sPt,sEta)) continue; 
-      //if(TMath::Abs(pdg)==211) nSel++;
-
-      // create AliTPCtrack object
-      BuildTrack(seed,charge);
-       
-      if(fSelAndSmear) {
-       bin = fDBgrid->GetBin(sPt,sEta);
-       switch (pdg) {
-       case 211:
-         fCovTree = covTreePi[bin];
-         break;
-       case 321:
-         fCovTree = covTreeKa[bin];
-         break;
-       case 2212:
-         fCovTree = covTreePr[bin];
-         break;
-       case 11:
-         fCovTree = covTreeEl[bin];
-         break;
-       case 13:
-         fCovTree = covTreeMu[bin];
-         break;
-       }
-       // deal with covariance matrix and smearing of parameters
-       CookTrack(sPt,sEta);
-
-       // assign the track a dE/dx and make a rough PID
-       CookdEdx(sPt,sEta);
+  nSeeds = sArray.GetEntries();
+  printf("+++ Number of seeds: %d\n",nSeeds);
+    
+  // loop over entries in sArray
+  for(Int_t l=0; l<nSeeds; l++) {
+    //if(l%1==0) printf("  --- Processing seed %d of %d ---\n",l,nSeeds);
+    
+    seed = (AliTPCseedGeant*)sArray.At(l);
+    
+    // Get track label
+    label = seed->GetLabel();
+    
+    // check if this track has already been processed
+    if(done[label]) continue;
+
+    // PDG code & electric charge
+    pdg = pdgCodes[label];
+    if(pdg>200 || pdg==-11 || pdg==-13) { charge=1; }
+    else if(pdg<-200 || pdg==11 || pdg==13) { charge=-1; }
+    else continue;
+    pdg = TMath::Abs(pdg);
+    if(pdg>3000) pdg=211;
+    
+    if(fSelAndSmear) SetParticle(pdg);
+    
+    sPt  = seed->GetPt();
+    sEta = seed->GetEta();
+    
+    // Apply selection according to TPC efficiency
+    //if(TMath::Abs(pdg)==211) nAcc++;
+    if(fSelAndSmear && !SelectedTrack(sPt,sEta)) continue; 
+    //if(TMath::Abs(pdg)==211) nSel++;
+
+    // create AliTPCtrack object
+    BuildTrack(seed,charge);
+
+    if(fSelAndSmear) {
+      bin = fDBgrid->GetBin(sPt,sEta);
+      switch (pdg) {
+      case 211:
+       //fCovTree = &(fCovTreePi[bin]);
+       fCovTree = fCovTreePi[bin];
+       break;
+      case 321:
+       //fCovTree = &(fCovTreeKa[bin]);
+       fCovTree = fCovTreeKa[bin];
+       break;
+      case 2212:
+       //fCovTree = &(fCovTreePr[bin]);
+       fCovTree = fCovTreePr[bin];
+       break;
+      case 11:
+       //fCovTree = &(fCovTreeEl[bin]);
+       fCovTree = fCovTreeEl[bin];
+       break;
+      case 13:
+       //fCovTree = &(fCovTreeMu[bin]);
+       fCovTree = fCovTreeMu[bin];
+       break;
       }
+      // deal with covariance matrix and smearing of parameters
+      CookTrack(sPt,sEta);
 
-      // put track in array
-      AliTPCtrack *iotrack = new AliTPCtrack(fTrack);
-      iotrack->SetLabel(label);
-      tArray.AddLast(iotrack);
-      // Mark track as "done" and register the pdg code
-      done[label] = kTRUE; 
-      tracks++;
-    } // loop over entries in sArray
-
-
-    // sort array with TPC tracks (decreasing pT)
-    tArray.Sort();
-
-    arrentr = tArray.GetEntriesFast();
-    for(Int_t l=0; l<arrentr; l++) {
-      tpctrack=(AliTPCtrack*)tArray.UncheckedAt(l);
-      tracktree->Fill();
+      // assign the track a dE/dx and make a rough PID
+      CookdEdx(sPt,sEta);
     }
-
-
-    // write the tree with tracks in the output file
-    out->cd();
-    tracktree->Write(); 
     
-    delete tracktree;
-    delete [] done;
-    delete [] pdgCodes;
-    delete [] ptkine;
-    delete [] pzkine;    
-
-    printf("\n\n+++\n+++ Number of TPC tracks: %d\n+++\n",tracks);
-    //cerr<<"Average Eff: "<<(Float_t)nSel/nAcc<<endl;
-
-    sArray.Delete();
-    tArray.Delete();
+    // put track in array
+    AliTPCtrack *iotrack = new AliTPCtrack(fTrack);
+    iotrack->SetLabel(label);
+    tArray.AddLast(iotrack);
+    // Mark track as "done" and register the pdg code
+    done[label] = kTRUE; 
+    tracks++;
+    
+  } // loop over entries in sArray
+  
+  // sort array with TPC tracks (decreasing pT)
+  tArray.Sort();
+  
+  // convert to AliESDtrack and write to AliESD
+  arrentr = tArray.GetEntriesFast(); 
+  Int_t idx;
+  Double_t wgts[5];
+  for(Int_t l=0; l<arrentr; l++) {
+    tpctrack=(AliTPCtrack*)tArray.UncheckedAt(l);
+    AliESDtrack ioESDtrack;
+    ioESDtrack.UpdateTrackParams(tpctrack,AliESDtrack::kTPCin);
+    // rough PID
+    wgts[0]=0.; wgts[1]=0.; wgts[2]=0.; wgts[3]=0.; wgts[4]=0.;
+    if(TMath::Abs(tpctrack->GetMass()-0.9)<0.1) {
+      idx = 4; // proton
+    } else if(TMath::Abs(tpctrack->GetMass()-0.5)<0.1) {
+      idx = 3; // kaon
+    } else {
+      idx = 2; // pion
+    }
+    wgts[idx] = 1.;
+    ioESDtrack.SetESDpid(wgts);
+    event->AddTrack(&ioESDtrack);
+  }
+  
+  
+  delete [] done;
+  delete [] pdgCodes;
+  
+  printf("+++ Number of TPC tracks: %d\n",tracks);
+  //cerr<<"Average Eff: "<<(Float_t)nSel/nAcc<<endl;
+  
+  sArray.Delete();
+  tArray.Delete();
   
-    infile->cd();
-  } // loop on events
-
-  if(fileDB) fileDB->Close();
-
   return 0;
 }
 //-----------------------------------------------------------------------------
@@ -1029,11 +967,11 @@ void AliTPCtrackerParam::BuildTrack(AliTPCseedGeant *s,Int_t ch) {
   cc[1]=cc[3]=cc[4]=cc[6]=cc[7]=cc[8]=cc[10]=cc[11]=cc[12]=cc[13]=0.;
   
   // Magnetic field
-  TVector3 bfield(0.,0.,fBz);
+  TVector3 bfield(0.,0.,-fBz);
   
   
   // radius [cm] of track projection in (x,y) 
-  Double_t rho = s->GetPt()*100./0.299792458/bfield.Z();
+  Double_t rho = s->GetPt()*100./0.299792458/TMath::Abs(bfield.Z());
   // center of track projection in local reference frame
   TVector3 sMom,sPos;
 
@@ -1060,7 +998,7 @@ void AliTPCtrackerParam::BuildTrack(AliTPCseedGeant *s,Int_t ch) {
   xx[0] = s->GetYL();
   xx[1] = s->GetZL();
   xx[3] = s->GetPz()/s->GetPt();
-  xx[4] = -ch/rho;
+  xx[4] = ch/rho;
   xx[2] = xx[4]*x0;
 
   // create the object AliTPCtrack    
@@ -1070,40 +1008,8 @@ void AliTPCtrackerParam::BuildTrack(AliTPCseedGeant *s,Int_t ch) {
   return;
 }
 //-----------------------------------------------------------------------------
-Int_t AliTPCtrackerParam::CheckLabel(AliTPCseedGeant *s,Int_t nPart,
-                                    Double_t *ptkine,Double_t *pzkine) const {
-//-----------------------------------------------------------------------------
-// This function checks if the label of the seed has been correctly 
-// assigned (to do only for pp charm production with AliRoot v3-08-02)
-//-----------------------------------------------------------------------------
-
-  Int_t sLabel = s->GetLabel();
-  Double_t sPt = s->GetPt();
-  Double_t sPz = s->GetPz();
-  
-  // check if the label is correct (comparing momentum)
-  if(sLabel<nPart && 
-     TMath::Abs(sPt-ptkine[sLabel])*
-     TMath::Abs(sPz-pzkine[sLabel])<0.001) return 0;
-
-  if((sLabel-30)>=nPart) return 1;  
-
-  Double_t diff=0,mindiff=1000.;
-  Int_t bestLabel=0;
-
-  for(Int_t i=sLabel-30; i<sLabel; i++) {
-    if(i<0 || i>=nPart) continue;
-    diff = TMath::Abs(sPt-ptkine[i])*TMath::Abs(sPz-pzkine[i]); 
-    if(diff<mindiff) { mindiff = diff; bestLabel = i; }
-  }
-
-  if(mindiff>0.001) return 1;
-  s->SetLabel(bestLabel);
-
-  return 0;
-}
-//-----------------------------------------------------------------------------
 void AliTPCtrackerParam::CompareTPCtracks(
+                          Int_t nEvents,
                           const Char_t* galiceName,
                           const Char_t* trkGeaName,
                           const Char_t* trkKalName,
@@ -1168,7 +1074,7 @@ void AliTPCtrackerParam::CompareTPCtracks(
   Char_t tname[100];
 
   // loop on events in file
-  for(Int_t evt=0; evt<fNevents; evt++) {  
+  for(Int_t evt=0; evt<nEvents; evt++) {  
     cerr<<"\n  --- Reading tracks for event "<<evt<<" ---\n\n";
     sprintf(tname,"TreeT_TPC_%d",evt);
     
@@ -1426,24 +1332,28 @@ void AliTPCtrackerParam::CookdEdx(Double_t pt,Double_t eta) {
   //Very rough PID
   Double_t p = TMath::Sqrt(1.+t.GetTgl()*t.GetTgl())*pt;
 
+  Double_t massPi = (Double_t)TDatabasePDG::Instance()->GetParticle(211)->Mass();
+  Double_t massKa = (Double_t)TDatabasePDG::Instance()->GetParticle(321)->Mass();
+  Double_t massPr = (Double_t)TDatabasePDG::Instance()->GetParticle(2212)->Mass();
+
   if (p<0.6) {
     if (dEdx < 39.+ 12./(p+0.25)/(p+0.25)) { 
-      t.AssignMass(AliPID::ParticleMass(AliPID::kPion)); new(&fTrack) AliTPCtrack(t); return;
+      t.AssignMass(massPi); new(&fTrack) AliTPCtrack(t); return;
     }
     if (dEdx < 39.+ 12./p/p) { 
-      t.AssignMass(AliPID::ParticleMass(AliPID::kKaon)); new(&fTrack) AliTPCtrack(t); return;
+      t.AssignMass(massKa); new(&fTrack) AliTPCtrack(t); return;
     }
-    t.AssignMass(AliPID::ParticleMass(AliPID::kProton)); new(&fTrack) AliTPCtrack(t); return;
+    t.AssignMass(massPr); new(&fTrack) AliTPCtrack(t); return;
   }
 
   if (p<1.2) {
     if (dEdx < 39.+ 12./(p+0.25)/(p+0.25)) { 
-      t.AssignMass(AliPID::ParticleMass(AliPID::kPion)); new(&fTrack) AliTPCtrack(t); return;
+      t.AssignMass(massPi); new(&fTrack) AliTPCtrack(t); return;
     }
-    t.AssignMass(AliPID::ParticleMass(AliPID::kProton)); new(&fTrack) AliTPCtrack(t); return;
+    t.AssignMass(massPr); new(&fTrack) AliTPCtrack(t); return;
   }
 
-  t.AssignMass(AliPID::ParticleMass(AliPID::kPion)); new(&fTrack) AliTPCtrack(t); return;
+  t.AssignMass(massPi); new(&fTrack) AliTPCtrack(t); return;
 }
 //-----------------------------------------------------------------------------
 void AliTPCtrackerParam::CookTrack(Double_t pt,Double_t eta) {
@@ -1493,10 +1403,10 @@ void AliTPCtrackerParam::CookTrack(Double_t pt,Double_t eta) {
   cc[12]= covmat.c42*RegFunc(trkKine,trkRegPar);
   cc[13]= covmat.c43;
   for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(8,l);
-  cc[14]= covmat.c44*RegFunc(trkKine,trkRegPar);
-   
-  TMatrixD covMatSmear(5,5);
-    
+  cc[14]= covmat.c44*RegFunc(trkKine,trkRegPar);  
+
+
+  TMatrixD covMatSmear(5,5);    
   covMatSmear = GetSmearingMatrix(cc,pt,eta);
 
   // get track original parameters
@@ -1509,6 +1419,7 @@ void AliTPCtrackerParam::CookTrack(Double_t pt,Double_t eta) {
   xx[4]=fTrack.GetC();
     
   // use smearing matrix to smear the original parameters
+  xxsm[0]=xref;
   SmearTrack(xx,xxsm,covMatSmear);
     
   AliTPCtrack track(0,xxsm,cc,xref,alpha);
@@ -2093,8 +2004,8 @@ void AliTPCtrackerParam::MakeSeedsFromHits(AliTPC *tpc,TTree *th,
   Int_t label;
   Int_t nTracks=(Int_t)th->GetEntries();
 
-  cerr<<"+++\n+++ Number of \"primary tracks\"(entries in TreeH): "<<nTracks<<
-         "\n+++\n\n";
+  cout<<"+++ Number of \"primary tracks\"(entries in TreeH): "<<nTracks<<
+         "\n";
    
   AliTPChit *tpcHit=0;
 
@@ -2155,13 +2066,12 @@ void AliTPCtrackerParam::MakeSeedsFromRefs(TTree *ttr,
   if(!b) {cerr<<"TPC branch of TreeTR not found"<<endl; return; }
   b->SetAddress(&tkRefArray);
   Int_t nTkRef = (Int_t)b->GetEntries();
-  cerr<<"+++\n+++ Number of entries in TreeTR(TPC): "<<nTkRef<<
-         "\n+++\n\n";
+  cerr<<"+++ Number of entries in TreeTR(TPC): "<<nTkRef<<"\n";
 
   // loop on track references
   for(Int_t l=0; l<nTkRef; l++){
-    if(l%1000==0) cerr<<"  --- Processing primary track "
-                     <<l<<" of "<<nTkRef<<" ---\r";
+    //if(l%1000==0) cerr<<"  --- Processing primary track "
+    //       <<l<<" of "<<nTkRef<<" ---\r";
     if(!b->GetEvent(l)) continue;
     nnn = tkRefArray->GetEntriesFast();
 
@@ -2185,7 +2095,8 @@ void AliTPCtrackerParam::MakeSeedsFromRefs(TTree *ttr,
       AliTPCseedGeant *ioseed = new AliTPCseedGeant(xg,yg,zg,px,py,pz,label);
 
       // reject if not at the inner part of TPC
-      if(TMath::Abs(ioseed->GetXL()-82.9701) > 0.1) { 
+      if(TMath::Abs(ioseed->GetXL()-83.8) > 0.2) { 
+       //cerr<<"not at TPC inner part: XL = "<<ioseed->GetXL()<<endl;
        delete ioseed; continue; 
       }
 
@@ -2349,6 +2260,75 @@ Int_t AliTPCtrackerParam::ReadAllData(const Char_t* inName) {
   if(!ReaddEdx(inName,11)) return 0;
   if(!ReaddEdx(inName,13)) return 0;
   if(!ReadDBgrid(inName)) return 0;
+  if(!ReadCovTrees(inName)) return 0;
+
+  return 1;
+}
+//-----------------------------------------------------------------------------
+Int_t AliTPCtrackerParam::ReadCovTrees(const Char_t* inName) {
+//-----------------------------------------------------------------------------
+// This function reads the covariance matrix trees from the DB
+//-----------------------------------------------------------------------------
+
+  if(gSystem->AccessPathName(inName,kFileExists)) { 
+    cerr<<"AliTPCtrackerParam::ReaddEdx: "<<inName<<" not found\n"; 
+    return 0; }
+
+  TString str;
+
+  TFile *inFile = TFile::Open(inName);
+
+
+  Int_t nBinsPi = fDBgridPi.GetTotBins();
+  for(Int_t l=0; l<nBinsPi; l++) {
+    str = "/CovMatrices/Pions/CovTreePi_bin";
+    str+=l;
+    fCovTreePi[l] = (TTree*)inFile->Get(str.Data());
+    //fCovTree = (TTree*)inFile->Get(str.Data());
+    //fCovTreePi[l] = new TTree(*fCovTree);
+  }
+  Int_t nBinsKa = fDBgridKa.GetTotBins();
+  for(Int_t l=0; l<nBinsKa; l++) {
+    str = "/CovMatrices/Kaons/CovTreeKa_bin";
+    str+=l;
+    fCovTreeKa[l] = (TTree*)inFile->Get(str.Data());
+    //fCovTree = (TTree*)inFile->Get(str.Data());
+    //fCovTreeKa[l] = new TTree(*fCovTree);
+  }
+  Int_t nBinsPr = fDBgridPr.GetTotBins();
+  for(Int_t l=0; l<nBinsPr; l++) {
+    if(fColl==0) {      
+      str = "/CovMatrices/Pions/CovTreePi_bin";
+    } else {
+      str = "/CovMatrices/Protons/CovTreePr_bin";
+    }
+    str+=l;
+    fCovTreePr[l] = (TTree*)inFile->Get(str.Data());
+    //fCovTree = (TTree*)inFile->Get(str.Data());
+    //fCovTreePr[l] = new TTree(*fCovTree);
+  }
+  Int_t nBinsEl = fDBgridEl.GetTotBins();
+  for(Int_t l=0; l<nBinsEl; l++) {
+    str = "/CovMatrices/Electrons/CovTreeEl_bin";
+    str+=l;
+    fCovTreeEl[l] = (TTree*)inFile->Get(str.Data());
+    //fCovTree = (TTree*)inFile->Get(str.Data());
+    //fCovTreeEl[l] = new TTree(*fCovTree);
+  }
+  Int_t nBinsMu = fDBgridMu.GetTotBins();
+  for(Int_t l=0; l<nBinsMu; l++) {
+    if(fColl==0) {      
+      str = "/CovMatrices/Pions/CovTreePi_bin";
+    } else {
+      str = "/CovMatrices/Muons/CovTreeMu_bin";
+    }
+    str+=l;
+    fCovTreeMu[l] = (TTree*)inFile->Get(str.Data());
+    //fCovTree = (TTree*)inFile->Get(str.Data());
+    //fCovTreeMu[l] = new TTree(*fCovTree);
+  }
+
+  //inFile->Close();
 
   return 1;
 }
@@ -2548,16 +2528,83 @@ Int_t AliTPCtrackerParam::ReadRegParams(const Char_t* inName,Int_t pdg) {
   if(gSystem->AccessPathName(inName,kFileExists)) { 
     cerr<<"AliTPCtrackerParam::ReadRegParams: "<<inName<<" not found\n"; 
     return 0; }
+  // Introduced change to "reverse" the TMatrixD read from file.
+  // Needed because storage mode of TMatrixD changed from column-wise
+  // to rwo-wise in ROOT.
+  //
+  // A.Dainese 03/06/05 
+
+  TMatrixD dummyMatrix(9,3);
 
   TFile *inFile = TFile::Open(inName);
   switch (pdg) {
   case 211:
     fRegPar = (TMatrixD*)inFile->Get("/RegParams/Pions/RegPions");
     new(&fRegParPi) TMatrixD(*fRegPar);
+    //printf("before\n");
+    //for(Int_t jj=0;jj<9;jj++) printf("%13.10f   %13.10f  %13.10f\n",fRegParPi(jj,0),fRegParPi(jj,1),fRegParPi(jj,2));
+    dummyMatrix = fRegParPi;  
+    fRegParPi(0,0) = dummyMatrix(0,0);
+    fRegParPi(0,1) = dummyMatrix(0,1);
+    fRegParPi(0,2) = dummyMatrix(0,2);
+    fRegParPi(1,0) = dummyMatrix(3,0);
+    fRegParPi(1,1) = dummyMatrix(1,1);
+    fRegParPi(1,2) = dummyMatrix(1,2);
+    fRegParPi(2,0) = dummyMatrix(6,0);
+    fRegParPi(2,1) = dummyMatrix(3,2);
+    fRegParPi(2,2) = dummyMatrix(2,2);
+    fRegParPi(3,0) = dummyMatrix(1,0);
+    fRegParPi(3,1) = dummyMatrix(4,0);
+    fRegParPi(3,2) = dummyMatrix(7,0);
+    fRegParPi(4,0) = dummyMatrix(3,1);
+    fRegParPi(4,1) = dummyMatrix(4,1);
+    fRegParPi(4,2) = dummyMatrix(7,1);
+    fRegParPi(5,0) = dummyMatrix(6,1);
+    fRegParPi(5,1) = dummyMatrix(4,2);
+    fRegParPi(5,2) = dummyMatrix(7,2);
+    fRegParPi(6,0) = dummyMatrix(2,0);
+    fRegParPi(6,1) = dummyMatrix(5,0);
+    fRegParPi(6,2) = dummyMatrix(8,0);
+    fRegParPi(7,0) = dummyMatrix(2,1);
+    fRegParPi(7,1) = dummyMatrix(5,1);
+    fRegParPi(7,2) = dummyMatrix(8,1);
+    fRegParPi(8,0) = dummyMatrix(6,2);
+    fRegParPi(8,1) = dummyMatrix(5,2);
+    fRegParPi(8,2) = dummyMatrix(8,2);
+    //printf("after\n");
+    //for(Int_t jj=0;jj<9;jj++) printf("%13.10f   %13.10f  %13.10f\n",fRegParPi(jj,0),fRegParPi(jj,1),fRegParPi(jj,2));
     break;
   case 321:
     fRegPar = (TMatrixD*)inFile->Get("/RegParams/Kaons/RegKaons");
     new(&fRegParKa) TMatrixD(*fRegPar);
+    dummyMatrix = fRegParKa;  
+    fRegParKa(0,0) = dummyMatrix(0,0);
+    fRegParKa(0,1) = dummyMatrix(0,1);
+    fRegParKa(0,2) = dummyMatrix(0,2);
+    fRegParKa(1,0) = dummyMatrix(3,0);
+    fRegParKa(1,1) = dummyMatrix(1,1);
+    fRegParKa(1,2) = dummyMatrix(1,2);
+    fRegParKa(2,0) = dummyMatrix(6,0);
+    fRegParKa(2,1) = dummyMatrix(3,2);
+    fRegParKa(2,2) = dummyMatrix(2,2);
+    fRegParKa(3,0) = dummyMatrix(1,0);
+    fRegParKa(3,1) = dummyMatrix(4,0);
+    fRegParKa(3,2) = dummyMatrix(7,0);
+    fRegParKa(4,0) = dummyMatrix(3,1);
+    fRegParKa(4,1) = dummyMatrix(4,1);
+    fRegParKa(4,2) = dummyMatrix(7,1);
+    fRegParKa(5,0) = dummyMatrix(6,1);
+    fRegParKa(5,1) = dummyMatrix(4,2);
+    fRegParKa(5,2) = dummyMatrix(7,2);
+    fRegParKa(6,0) = dummyMatrix(2,0);
+    fRegParKa(6,1) = dummyMatrix(5,0);
+    fRegParKa(6,2) = dummyMatrix(8,0);
+    fRegParKa(7,0) = dummyMatrix(2,1);
+    fRegParKa(7,1) = dummyMatrix(5,1);
+    fRegParKa(7,2) = dummyMatrix(8,1);
+    fRegParKa(8,0) = dummyMatrix(6,2);
+    fRegParKa(8,1) = dummyMatrix(5,2);
+    fRegParKa(8,2) = dummyMatrix(8,2);
     break;
   case 2212:
     if(fColl==0) {
@@ -2566,10 +2613,66 @@ Int_t AliTPCtrackerParam::ReadRegParams(const Char_t* inName,Int_t pdg) {
       fRegPar = (TMatrixD*)inFile->Get("/RegParams/Protons/RegProtons");
     }
     new(&fRegParPr) TMatrixD(*fRegPar);
+    dummyMatrix = fRegParPr;  
+    fRegParPr(0,0) = dummyMatrix(0,0);
+    fRegParPr(0,1) = dummyMatrix(0,1);
+    fRegParPr(0,2) = dummyMatrix(0,2);
+    fRegParPr(1,0) = dummyMatrix(3,0);
+    fRegParPr(1,1) = dummyMatrix(1,1);
+    fRegParPr(1,2) = dummyMatrix(1,2);
+    fRegParPr(2,0) = dummyMatrix(6,0);
+    fRegParPr(2,1) = dummyMatrix(3,2);
+    fRegParPr(2,2) = dummyMatrix(2,2);
+    fRegParPr(3,0) = dummyMatrix(1,0);
+    fRegParPr(3,1) = dummyMatrix(4,0);
+    fRegParPr(3,2) = dummyMatrix(7,0);
+    fRegParPr(4,0) = dummyMatrix(3,1);
+    fRegParPr(4,1) = dummyMatrix(4,1);
+    fRegParPr(4,2) = dummyMatrix(7,1);
+    fRegParPr(5,0) = dummyMatrix(6,1);
+    fRegParPr(5,1) = dummyMatrix(4,2);
+    fRegParPr(5,2) = dummyMatrix(7,2);
+    fRegParPr(6,0) = dummyMatrix(2,0);
+    fRegParPr(6,1) = dummyMatrix(5,0);
+    fRegParPr(6,2) = dummyMatrix(8,0);
+    fRegParPr(7,0) = dummyMatrix(2,1);
+    fRegParPr(7,1) = dummyMatrix(5,1);
+    fRegParPr(7,2) = dummyMatrix(8,1);
+    fRegParPr(8,0) = dummyMatrix(6,2);
+    fRegParPr(8,1) = dummyMatrix(5,2);
+    fRegParPr(8,2) = dummyMatrix(8,2);
     break;
   case 11:
     fRegPar = (TMatrixD*)inFile->Get("/RegParams/Electrons/RegElectrons");
     new(&fRegParEl) TMatrixD(*fRegPar);
+    dummyMatrix = fRegParEl;  
+    fRegParEl(0,0) = dummyMatrix(0,0);
+    fRegParEl(0,1) = dummyMatrix(0,1);
+    fRegParEl(0,2) = dummyMatrix(0,2);
+    fRegParEl(1,0) = dummyMatrix(3,0);
+    fRegParEl(1,1) = dummyMatrix(1,1);
+    fRegParEl(1,2) = dummyMatrix(1,2);
+    fRegParEl(2,0) = dummyMatrix(6,0);
+    fRegParEl(2,1) = dummyMatrix(3,2);
+    fRegParEl(2,2) = dummyMatrix(2,2);
+    fRegParEl(3,0) = dummyMatrix(1,0);
+    fRegParEl(3,1) = dummyMatrix(4,0);
+    fRegParEl(3,2) = dummyMatrix(7,0);
+    fRegParEl(4,0) = dummyMatrix(3,1);
+    fRegParEl(4,1) = dummyMatrix(4,1);
+    fRegParEl(4,2) = dummyMatrix(7,1);
+    fRegParEl(5,0) = dummyMatrix(6,1);
+    fRegParEl(5,1) = dummyMatrix(4,2);
+    fRegParEl(5,2) = dummyMatrix(7,2);
+    fRegParEl(6,0) = dummyMatrix(2,0);
+    fRegParEl(6,1) = dummyMatrix(5,0);
+    fRegParEl(6,2) = dummyMatrix(8,0);
+    fRegParEl(7,0) = dummyMatrix(2,1);
+    fRegParEl(7,1) = dummyMatrix(5,1);
+    fRegParEl(7,2) = dummyMatrix(8,1);
+    fRegParEl(8,0) = dummyMatrix(6,2);
+    fRegParEl(8,1) = dummyMatrix(5,2);
+    fRegParEl(8,2) = dummyMatrix(8,2);
     break;
   case 13:
     if(fColl==0) {
@@ -2578,6 +2681,34 @@ Int_t AliTPCtrackerParam::ReadRegParams(const Char_t* inName,Int_t pdg) {
       fRegPar = (TMatrixD*)inFile->Get("/RegParams/Muons/RegMuons");
     }
     new(&fRegParMu) TMatrixD(*fRegPar);
+    dummyMatrix = fRegParMu;  
+    fRegParMu(0,0) = dummyMatrix(0,0);
+    fRegParMu(0,1) = dummyMatrix(0,1);
+    fRegParMu(0,2) = dummyMatrix(0,2);
+    fRegParMu(1,0) = dummyMatrix(3,0);
+    fRegParMu(1,1) = dummyMatrix(1,1);
+    fRegParMu(1,2) = dummyMatrix(1,2);
+    fRegParMu(2,0) = dummyMatrix(6,0);
+    fRegParMu(2,1) = dummyMatrix(3,2);
+    fRegParMu(2,2) = dummyMatrix(2,2);
+    fRegParMu(3,0) = dummyMatrix(1,0);
+    fRegParMu(3,1) = dummyMatrix(4,0);
+    fRegParMu(3,2) = dummyMatrix(7,0);
+    fRegParMu(4,0) = dummyMatrix(3,1);
+    fRegParMu(4,1) = dummyMatrix(4,1);
+    fRegParMu(4,2) = dummyMatrix(7,1);
+    fRegParMu(5,0) = dummyMatrix(6,1);
+    fRegParMu(5,1) = dummyMatrix(4,2);
+    fRegParMu(5,2) = dummyMatrix(7,2);
+    fRegParMu(6,0) = dummyMatrix(2,0);
+    fRegParMu(6,1) = dummyMatrix(5,0);
+    fRegParMu(6,2) = dummyMatrix(8,0);
+    fRegParMu(7,0) = dummyMatrix(2,1);
+    fRegParMu(7,1) = dummyMatrix(5,1);
+    fRegParMu(7,2) = dummyMatrix(8,1);
+    fRegParMu(8,0) = dummyMatrix(6,2);
+    fRegParMu(8,1) = dummyMatrix(5,2);
+    fRegParMu(8,2) = dummyMatrix(8,2);
     break;
   }
   inFile->Close();
@@ -3221,15 +3352,25 @@ void AliTPCtrackerParam::SmearTrack(Double_t *xx,Double_t *xxsm,TMatrixD cov)
 //-----------------------------------------------------------------------------
 // This function smears track parameters according to streched cov. matrix
 //-----------------------------------------------------------------------------
+  Double_t xref=xxsm[0]; xxsm[0]=0;
+
   AliGausCorr *corgen = new AliGausCorr(cov,5);
   TArrayD corr(5);
   corgen->GetGaussN(corr);
+  // check on fP2(ESD) = fX*fP4-fP2(TPC)
+  // with fX=xref (not smeared), fP4=xx[4]+corr[4] e fP2=xx[2]+corr[2];
+  // if fP2(ESD)^2 > 1 -> resmear...
+  Double_t fp2esd=xref*(xx[4]+corr[4])-(xx[2]+corr[2]);
+  while ( (fp2esd*fp2esd) > 1.0 ) {
+    Warning("SmearTrack()","Badly smeared track, retrying...");
+    corgen->GetGaussN(corr);
+    fp2esd=xref*(xx[4]+corr[4])-(xx[2]+corr[2]);
+  }
+
   delete corgen;
   corgen = 0;
 
-  for(Int_t l=0;l<5;l++) {
-    xxsm[l] = xx[l]+corr[l];
-  }
+  for(Int_t l=0;l<5;l++) xxsm[l] = xx[l]+corr[l];
 
   return;
 }
index 3e52725..d790347 100644 (file)
@@ -19,7 +19,7 @@
 //                                  
 //                                                                        
 //  Origin: Andrea Dainese, Padova - e-mail: andrea.dainese@pd.infn.it     
-//                                                                        
+//  adapted to ESD output by: Marcello Lunardon, Padova
 /////////////////////////////////////////////////////////////////////////
 
 
@@ -30,6 +30,7 @@ class TTree;
 #include "AliConfig.h"
 #include "AliTPCkineGrid.h"
 #include "AliTPCtrack.h"
+#include "AliESD.h"
 //----------------------------
 
 class AliTPC;
@@ -39,7 +40,7 @@ class AliTPCtrackerParam:
 {
    
  public:
-  AliTPCtrackerParam(Int_t coll=0, Double_t Bz=0.4, Int_t n=1,
+  AliTPCtrackerParam(Int_t coll=0, Double_t Bz=0.4,
                     const char* evfoldname = AliConfig::GetDefaultEventFolderName());
   virtual ~AliTPCtrackerParam();
 
@@ -47,14 +48,16 @@ class AliTPCtrackerParam:
   //
   AliTPCtrackerParam(const AliTPCtrackerParam& p);
   //
-  Int_t BuildTPCtracks(const TFile *inp, TFile *out);
+  Int_t Init();
+  Int_t BuildTPCtracks(AliESD* event);
   // these functions are used to create a DB of cov. matrices,
   // including regularization, efficiencies and dE/dx
   void  AllGeantTracks() { fSelAndSmear=kFALSE; return; }
   void  AnalyzedEdx(const Char_t *outName,Int_t pdg);
   void  AnalyzePulls(const Char_t *outName);
   void  AnalyzeResolutions(Int_t pdg);
-  void  CompareTPCtracks(const Char_t *galiceName="galice.root",
+  void  CompareTPCtracks(Int_t nEvents=1,
+                        const Char_t *galiceName="galice.root",
                         const Char_t *trkGeaName="AliTPCtracksGeant.root",
                         const Char_t *trkKalName="AliTPCtracksSorted.root",
                         const Char_t *covmatName="CovMatrix.root",
@@ -87,6 +90,9 @@ class AliTPCtrackerParam:
                    Int_t lab=0);
     Int_t    GetLabel() const { return fLabel; }
     Double_t GetAlpha() const { return fAlpha; }      
+    Double_t GetXG() const { return fXg; }
+    Double_t GetYG() const { return fYg; }
+    Double_t GetZG() const { return fZg; }
     Double_t GetXL() const 
       { return fXg*TMath::Cos(fAlpha)+fYg*TMath::Sin(fAlpha); }
     Double_t GetYL() const 
@@ -122,7 +128,6 @@ class AliTPCtrackerParam:
 
   TString fEvFolderName;//! name of data folder
 
-  Int_t           fNevents;     // number of events in the file to be processed
   Double_t        fBz;          // value of the z component of L3 field (Tesla)
   Int_t           fColl;        // collision code (0: PbPb6000; 1: pp)
   Bool_t          fSelAndSmear; // if kFALSE returns GEANT tracks 
@@ -133,6 +138,11 @@ class AliTPCtrackerParam:
 
   TTree          *fCovTree;  // tree with regularized cov matrices 
                              // for the current track
+  TTree           *fCovTreePi[30];
+  TTree           *fCovTreeKa[30];
+  TTree           *fCovTreePr[30];
+  TTree           *fCovTreeEl[30];
+  TTree           *fCovTreeMu[30];
 
   AliTPCkineGrid *fDBgrid;   // grid for the cov matrix look-up table  
   AliTPCkineGrid  fDBgridPi; //               "                  for pions  
@@ -178,8 +188,6 @@ class AliTPCtrackerParam:
 
 
   void     BuildTrack(AliTPCseedGeant *s,Int_t ch);
-  Int_t    CheckLabel(AliTPCseedGeant *s,Int_t nPart,
-                     Double_t *ptkine,Double_t *pzkine) const;
   void     CookdEdx(Double_t pt,Double_t eta);
   void     CookTrack(Double_t pt,Double_t eta);
   TMatrixD GetSmearingMatrix(Double_t *cc, Double_t pt,Double_t eta) const;
@@ -193,6 +201,7 @@ class AliTPCtrackerParam:
   Int_t    ReadEffs(const Char_t *inName);
   Int_t    ReadPulls(const Char_t *inName);
   Int_t    ReadRegParams(const Char_t *inName,Int_t pdg);
+  Int_t    ReadCovTrees(const Char_t* inName); 
   Bool_t   SelectedTrack(Double_t pt, Double_t eta) const;
   void     SetParticle(Int_t pdg);  
   void     SmearTrack(Double_t *xx,Double_t *xxsm,TMatrixD cov) const;
index 77a6973..1783593 100644 (file)
@@ -22,6 +22,7 @@
 #ifndef __CINT__
 #include "Riostream.h"
 #include <TFile.h>
+#include <TTree.h>
 #include <TStopwatch.h>
 #include <TObject.h>
 #include "alles.h"
 #include "AliITStrackerV2.h"
 #include "AliKalmanTrack.h"
 #include "AliTPCtrackerParam.h"
+#include "AliTracker.h"
+#include "AliESD.h"
+#include "AliRun.h"
+#include "AliRunLoader.h"
 #endif
 
-void TPCParamTracks(const Char_t *galice,const Char_t *outname,const Int_t coll,const Double_t Bfield,Int_t n);
-void CreateAllGeantTracks(const Char_t *galice, const Char_t *outname,const Int_t coll,const Double_t Bfield,Int_t n);
+Int_t TPCParamTracks(const Int_t coll=1,Int_t firstEvent=0,Int_t lastEvent=0);
+void CreateAllGeantTracks(const Int_t coll=1,Int_t nev=1);
 void TrackCompare(const Int_t coll,const Double_t Bfield,Int_t n);
 void BuildDataBase(const Int_t coll,const Double_t Bfield);
 
-void AliTPCtrackingParamDB(Int_t n=1) {
-
-  const Char_t *name=" AliTPCtrackingParamTest";
-  cerr<<'\n'<<name<<"...\n";
-  gBenchmark->Start(name);
-
-
-  const Char_t *TPCtrkName="AliTPCtracksParam.root";
-  const Char_t *TPCgeatrkName="AliTPCtracksGeant.root";
-  const Char_t *galiceName="galice.root";
-
-  // set here the code for the type of collision (needed for TPC tracking
-  // parameterization). available collisions:
-  //
-  // coll = 0 ->   PbPb6000 (HIJING with b<2fm) 
-  const Int_t    collcode = 0;  
-  // set here the value of the magnetic field
-  const Double_t BfieldValue = 0.4;
-
-  AliKalmanTrack::SetConvConst(100/0.299792458/BfieldValue);
-
-  // ********** Build TPC tracks with parameterization *********** // 
-  TPCParamTracks(galiceName,TPCtrkName,collcode,BfieldValue,n);
-
-  // ********** Build All True TPC tracks *********** //
-  CreateAllGeantTracks(galiceName,TPCgeatrkName,collcode,BfieldValue,n);
-    
-  // ********** Compare Kalman tracks with tracks at 1st hit *********** //
-  TrackCompare(collcode,BfieldValue);
-  
-   
-  // ********** Merge files, compute pulls, and dEdx *********** //
-  BuildDataBase(collcode,BfieldValue);
-  
-
-
-
-  gBenchmark->Stop(name);
-  gBenchmark->Show(name);
-
-  return;
-}
-
 //_____________________________________________________________________________
-void TPCParamTracks(const Char_t *galice, const Char_t *outname,
-                   const Int_t coll,const Double_t Bfield,Int_t n) {
+Int_t TPCParamTracks(const Int_t coll,Int_t firstEvent,Int_t lastEvent) {
 //
 // Ordinary TPC tracking parameterization
 //
-  cerr<<"\n*******************************************************************\n";
 
-  const Char_t *name="TPCParamTracks";
-  cerr<<'\n'<<name<<"...\n";
-  gBenchmark->Start(name);
-
-  TFile *outfile=TFile::Open(outname,"recreate");
-  TFile *infile =TFile::Open(galice);
-
-  AliTPCtrackerParam tracker(coll,Bfield,n);
-  tracker.BuildTPCtracks(infile,outfile);
-
-  delete gAlice; gAlice=0;
-
-  infile->Close();
-  outfile->Close();
-
-  gBenchmark->Stop(name);
-  gBenchmark->Show(name);
-
-  return;
+   /**** Initialization of the NewIO *******/
+
+   if (gAlice) {
+      delete gAlice->GetRunLoader();
+      delete gAlice; 
+      gAlice=0;
+   }
+
+   AliRunLoader *rl = AliRunLoader::Open("galice.root");
+   if (rl == 0x0) {
+      cerr<<"Can not open session"<<endl;
+      return;
+   }
+   Int_t retval = rl->LoadgAlice();
+   if (retval) {
+      cerr<<"LoadgAlice returned error"<<endl;
+      delete rl;
+      return;
+   }
+   retval = rl->LoadHeader();
+   if (retval) {
+      cerr<<"LoadHeader returned error"<<endl;
+      delete rl;
+      return;
+   }
+   gAlice=rl->GetAliRun();
+       
+
+   TDatabasePDG *DataBase = TDatabasePDG::Instance();
+
+   // Get field from galice.root
+   AliMagF *fiel = (AliMagF*)gAlice->Field();
+   Double_t fieval=TMath::Abs((Double_t)fiel->SolenoidField()/10.);
+   // Set the conversion constant between curvature and Pt
+   AliTracker::SetFieldMap(fiel,kTRUE);
+
+   /**** The TPC corner ********************/
+
+   AliTPCtrackerParam tpcTrackerPar(coll,fieval);
+   tpcTrackerPar.Init();
+
+   /***** The TREE is born *****/
+   
+   TTree *esdTree=new TTree("esdTree","Tree with ESD objects");
+   AliESD *event=0;
+   esdTree->Branch("ESD","AliESD",&event);
+   
+   if(firstEvent>rl->GetNumberOfEvents()) firstEvent=rl->GetNumberOfEvents()-1;
+   if(lastEvent>rl->GetNumberOfEvents())  lastEvent=rl->GetNumberOfEvents()-1;
+   cout<<" Number of events: "<<1+lastEvent-firstEvent<<endl;
+   
+   //<----------------------------------The Loop over events begins
+   TStopwatch timer;
+   Int_t trc;
+   for(Int_t i=firstEvent; i<=lastEvent; i++) { 
+     
+     cerr<<" Processing event number : "<<i<<endl;
+     AliESD *event = new AliESD(); 
+     event->SetRunNumber(gAlice->GetRunNumber());
+     event->SetEventNumber(i);
+     event->SetMagneticField(gAlice->Field()->SolenoidField());
+     rl->GetEvent(i);
+
+     if ( (trc=tpcTrackerPar.BuildTPCtracks(event)) ) {
+       printf("exiting tracker with code %d in event %d\n",trc,i);
+       esdTree->Fill(); delete event;
+       continue;
+     }
+
+     esdTree->Fill();
+     delete event;
+
+   }//<-----------------------------------The Loop over events ends here
+   timer.Stop(); timer.Print();
+
+   //        The AliESDs.root is born
+   TFile *ef = TFile::Open("AliESDs.root","RECREATE"); 
+   if (!ef || !ef->IsOpen()) {cerr<<"Can't open AliESDs.root !\n"; return;}
+
+   esdTree->Write(); //Write the TREE and close everything
+   delete esdTree;
+   ef->Close();
+
+   delete rl;
+
+   return;
 }
 //_____________________________________________________________________________
-void CreateAllGeantTracks(const Char_t *galice, const Char_t *outname,
-                         const Int_t coll,const Double_t Bfield,Int_t n) {
+void CreateAllGeantTracks(const Int_t coll,Int_t nev) {
 //
 // Get all tracks at TPC 1st hit w/o selection and smearing
 //