TPC/AliBarrelRec_TPCparam.C

   1 void AliBarrelRec_TPCparam(Int_t firstEvent=0,Int_t lastEvent=0) {
   2   //
   3   // Macro to create AliESDs.root using parametrized TPC tracking
   4   // and AliITStrackerSA (MI + SA)
   5   //
   6   // Input files:
   7   // - galice.root
   8   // - Kinematics.root
   9   // - TrackReferences.root
  10   // - ITS.RecPoints.root (use AliRecontruction class)
  11   // - ITS.Vertex.root (use $ALICE_ROOT/ITS/AliITSVertexerZTest.C)
  12   //
  13   // A. Dainese - INFN Legnaro
  14   //
  15
  16   Int_t  collcode = 1; // pp collisions
  17   Bool_t useMeanVtx = kTRUE;
  18
  19
  20   if (gAlice) {
  21     delete gAlice->GetRunLoader();
  22     delete gAlice;
  23     gAlice=0;
  24   }
  25   AliRunLoader *rl = AliRunLoader::Open("galice.root");
  26   if (rl == 0x0) {
  27     cerr<<"Can not open session"<<endl;
  28     return;
  29   }
  30   Int_t retval = rl->LoadgAlice();
  31   if (retval) {
  32     cerr<<"LoadgAlice returned error"<<endl;
  33     delete rl;
  34     return;
  35   }
  36   retval = rl->LoadHeader();
  37   if (retval) {
  38     cerr<<"LoadHeader returned error"<<endl;
  39     delete rl;
  40     return;
  41   }
  42   gAlice=rl->GetAliRun();
  43
  44
  45
  46   // Get field from galice.root
  47   AliMagF *fiel = (AliMagF*)gAlice->Field();
  48   Double_t fieval=TMath::Abs((Double_t)fiel->SolenoidField()/10.);
  49   // Set the conversion constant between curvature and Pt
  50   AliTracker::SetFieldMap(fiel,kTRUE);
  51
  52   /**** The TPC corner ********************/
  53
  54   AliTPCtrackerParam tpcTrackerPar(collcode,fieval);
  55   tpcTrackerPar.Init();
  56
  57   //**** Switch on the PID class (mandatory!)
  58   AliPID pid;
  59
  60   /**** The ITS corner ********************/
  61
  62   AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
  63   if (itsl == 0x0) {
  64     cerr<<"Cannot get the ITS loader"<<endl;
  65     return;
  66   }
  67   itsl->LoadRecPoints("read");
  68
  69   AliITS *dITS = (AliITS*)gAlice->GetDetector("ITS");
  70   if (!dITS) {
  71     cerr<<"Cannot find the ITS detector !"<<endl;
  72     return;
  73   }
  74   AliITSgeom *geom = dITS->GetITSgeom();
  75
  76   // Instance of the ITS tracker
  77   AliITStrackerSA itsTracker(geom);
  78   Int_t ITSclusters[6] = {1,1,1,1,1,1};
  79   itsTracker.SetLayersNotToSkip(ITSclusters);
  80
  81
  82   // Primary vertex reconstruction in pp
  83   AliESDVertex *initVertex = 0;
  84   TFile *invtx = new TFile("AliESDVertexMean.root");
  85   if(collcode==1 && useMeanVtx) {
  86     // open the mean vertex
  87     initVertex = (AliESDVertex*)invtx->Get("vtxmean");
  88   } else {
  89     Double_t pos[3]={0.5,0.5,0.};
  90     Double_t err[3]={3.,3.,5.};
  91     initVertex = new AliESDVertex(pos,err);
  92   }
  93   invtx->Close();
  94   delete invtx;
  95   AliVertexerTracks *vertexer = new AliVertexerTracks();
  96   vertexer->SetVtxStart(initVertex);
  97   vertexer->SetDebug(0);
  98   delete initVertex;
  99   initVertex=0;
 100
 101   /***** The TREE for ESD is born *****/
 102   TTree *esdTree=new TTree("esdTree","Tree with ESD objects");
 103   AliESD *event=0;
 104   esdTree->Branch("ESD","AliESD",&event);
 105
 106   if(firstEvent>rl->GetNumberOfEvents()) firstEvent=rl->GetNumberOfEvents()-1;
 107   if(lastEvent>rl->GetNumberOfEvents())  lastEvent=rl->GetNumberOfEvents()-1;
 108   cout<<" Number of events: "<<1+lastEvent-firstEvent<<endl;
 109
 110   TFile *ppZ = TFile::Open("ITS.Vertex.root"); // z vertices from SPD
 111   AliESDVertex *vertexSPD = new AliESDVertex();
 112   Char_t zver[100];
 113   Double_t vtx[3]={0,0,0};
 114   Double_t sigmavtx[3]={0.07,0.07,0.1};
 115
 116
 117   //<---------------------------------- The Loop over events begins
 118   TStopwatch timer;
 119   Int_t trc;
 120   for(Int_t i=firstEvent; i<=lastEvent; i++) {
 121
 122     cerr<<" Processing event number : "<<i<<endl;
 123     AliESD *event = new AliESD();
 124     event->SetRunNumber(gAlice->GetRunNumber());
 125     event->SetEventNumber(i);
 126     event->SetMagneticField(gAlice->Field()->SolenoidField());
 127     rl->GetEvent(i);
 128
 129     //***** Primary vertex from SPD from file
 130     sprintf(zver,"Event%d/Vertex",i);
 131     vertexSPD = (AliESDVertex*)ppZ->Get(zver);
 132     if(!vertexSPD) {
 133       esdTree->Fill(); delete event;
 134       continue;
 135     }
 136     event->SetVertex(vertexSPD);
 137     vertexSPD->GetXYZ(vtx);
 138     vertexSPD->GetSigmaXYZ(sigmavtx);
 139
 140     //***** TPC tracking
 141     if ( (trc=tpcTrackerPar.BuildTPCtracks(event)) ) {
 142       printf("exiting TPC tracker with code %d in event %d\n",trc,i);
 143       esdTree->Fill(); delete event;
 144       continue;
 145     }
 146
 147     //***** ITS tracking
 148     itsTracker.AliTracker::SetVertex(vtx,sigmavtx);
 149     TTree *itsTree=itsl->TreeR();
 150     if (!itsTree) {
 151       cerr<<"Can't get the ITS cluster tree !\n";
 152       esdTree->Fill(); delete event;
 153       return;
 154     }
 155     itsTracker.UnloadClusters();
 156     itsTracker.LoadClusters(itsTree);
 157     if ( (trc=itsTracker.Clusters2Tracks(event)) ) {
 158       printf("exiting ITS tracker with code %d in event %d\n",trc,i);
 159       esdTree->Fill(); delete event;
 160       continue;
 161     }
 162
 163     //***** Propagate kITSin tracks to local x = 0 (through beam pipe)
 164     ToLocalX0(event);
 165
 166     //***** Vertex from ESD tracks
 167     if(collcode==1) { // pp
 168       AliESDVertex *vertexTrks =
 169         (AliESDVertex*)vertexer->FindPrimaryVertex(event);
 170       event->SetPrimaryVertex(vertexTrks);
 171     }
 172
 173     esdTree->Fill();
 174     delete event;
 175
 176   }//<-----------------------------------The Loop over events ends here
 177   timer.Stop(); timer.Print();
 178
 179   // The AliESDs.root is born
 180   TFile *ef = TFile::Open("AliESDs.root","RECREATE");
 181   if (!ef || !ef->IsOpen()) {cerr<<"Can't open AliESDs.root !\n"; return;}
 182
 183   //Write the TREE and close everything
 184   esdTree->Write();
 185   delete esdTree;
 186   ef->Close();
 187
 188   delete vertexer;
 189   delete rl;
 190
 191   return;
 192 }
 193 //--------------------------------------------------------------------------
 194 void ToLocalX0(AliESD *esd) {
 195
 196   Int_t ntracks = esd->GetNumberOfTracks();
 197   AliESDtrack *esdTrack = 0;
 198
 199   // loop on tracks
 200   for(Int_t tr = 0; tr<ntracks; tr++) {
 201     esdTrack = (AliESDtrack*)esd->GetTrack(tr);
 202     // ask for kITSin
 203     if(!(esdTrack->GetStatus()&AliESDtrack::kITSin)) continue;
 204     AliITStrackV2 *itsTrack = 0;
 205     try {
 206       itsTrack = new AliITStrackV2(*esdTrack);
 207       esdTrack = 0;
 208     }
 209     catch (const Char_t *msg) {
 210         Warning("FindPrimaryVertexForCurrentEvent",msg);
 211         continue;
 212     }
 213     // propagate track to beam pipe (0.8 mm of Be)
 214     if(itsTrack->GetX()>3.) itsTrack->PropagateTo(3.,0.0023,65.19);
 215     // propagate track to (0,0)
 216     itsTrack->PropagateTo(0.,0.,0.);
 217     itsTrack->UpdateESDtrack(AliESDtrack::kITSin);
 218     esdTrack = new AliESDtrack(*(itsTrack->GetESDtrack()));
 219     delete itsTrack;
 220   } // end loop on tracks
 221
 222   return;
 223 }